text:
# Magnetism in systems with various dimensionality: A comparison between Fe and Co

## Abstract

A systematic ab initio study is performed for the spin and orbital moments and for the validity of the sum rules for x-ray magnetic circular dichroism for Fe systems with various dimensionality (bulk, Pt-supported monolayers and monatomic wires, free-standing monolayers and monatomic wires). Qualitatively, the results are similar to those for the respective Co systems, with the main difference that for the monatomic Fe wires the $`T_z`$ term in the spin sum rule is much larger than for the Co wires. The spin and orbital moments induced in the Pt substrate are also discussed.

One of the outstanding features of magnetism in low-dimensional spin systems is the importance of the spin-orbit-coupling effects. Whereas the orbital moments in bulk materials of high symmetry are strongly reduced due to orbital quenching, they may attain large values in low-dimensional systems. For instance, in monatomic Co wires at the steps of a vicinal Pt surface an orbital moment of about 0.68$`\mu _\text{B}`$ per Co atom was found Gambardella et al. (2002), which is a factor of about 5 larger than in bulk hcp Co and which represents the largest orbital moment ever reported for a 3$`d`$ itinerant electron system. For the magnetic-anisotropy energy of this Co wire a value was deduced Gambardella et al. (2002) which is about 50 times larger than the one of hcp Co (which is already large). The large magnetic anisotropy is also relevant for the nature of magnetic excitations. It will generate an excitation gap for the linear excitations, i.e., the spin waves, and for one-dimensional systems it will lead to nonlinear excitation modes which possibly have the character of solitons Mikeska and Steiner (1991).

A suitable experimental method to resolve spin and orbital moments is the technique of x-ray magnetic circular dichroism (XMCD). This technique Schütz et al. (1987) is based on the fact that for magnetic systems the absorption coefficient $`\mu (ϵ)`$ as function of the x-ray energy $`ϵ`$ is different for x-rays with left-circular polarization, $`\mu _+(ϵ)`$, and right-circular polarization, $`\mu _{}(ϵ)`$. For example, in a $`3d`$ transition metal system one has to measure these absorption coefficients at the $`L_2`$ and $`L_3`$ edge corresponding to the $`3d`$ transition metal atoms in the investigated system, and one then can obtain the orbital moment $`m_l=\mu _\text{B}l_z`$ and the spin moment $`m_s=\mu _\text{B}\sigma _z`$ for a transition metal atom via the XMCD sum rules Thole et al. (1992); Carra et al. (1993) according to:

$$l_z=\frac{2I_\text{m}N_\text{h}}{I_\text{t}},$$
(1)
$$\sigma _z=\frac{3I_\text{s}N_\text{h}}{I_\text{t}}7T_z,$$
(2)
$$I_\text{m}=\underset{E_\text{F}}{\overset{E_\text{c}}{}}\left[(\mu _\text{c})_{L_3}+(\mu _\text{c})_{L_2}\right]𝑑ϵ,$$
(3)
$$I_\text{s}=\underset{E_\text{F}}{\overset{E_\text{c}}{}}\left[(\mu _\text{c})_{L_3}2(\mu _\text{c})_{L_2}\right]𝑑ϵ,$$
(4)
$$I_\text{t}=\underset{E_\text{F}}{\overset{E_\text{c}}{}}\left[(\mu _\text{t})_{L_3}+(\mu _\text{t})_{L_2}\right]𝑑ϵ,$$
(5)
with the XMCD signal $`\mu _\text{c}=\mu ^+\mu ^{}`$ and with $`\mu _\text{t}=\mu ^++\mu ^{}+\mu ^0`$. Here we have assumed that the $`z`$-axis is parallel to the propagation direction of the x-rays. Then $`\mu _0(ϵ)`$ is the absorption coefficient for linear polarization along the $`z`$-axis. $`N_\text{h}`$ is the number of holes in the $`d`$ part of the valence band and $`T_z`$ is the expectation value of the magnetic dipolar operator

$$\widehat{T}_z=\frac{1}{2}[𝝈3\widehat{𝐫}(\widehat{𝐫}𝝈)]_z,$$
(6)
where $`𝝈`$ denotes the vector of the Pauli matrices. The quantities $`E_\text{F}`$ and $`E_\text{c}`$ denote the Fermi energy and a cutoff energy. For details of such calculations see Ref. Ederer et al., 2002.

For a practical application of these sum rules there are several problems, especially for low-dimensional systems, and we just want to mention two of them. First, the $`T_z`$ term in general is not known from the experiment, and it is therefore often neglected. However, it is a measure of the anisotropy of the spin density in the material, and it is therefore expected that it becomes more and more important when going to systems with more and more reduced dimensionality and hence more and more reduced symmetry. Stöhr and König Stöhr and König (1995) introduced an angle-dependent XMCD technique which allows to eliminate the $`T_z`$ term from the analysis with the spin sum rule. The basic idea is to measure $`\mu _\text{c}(ϵ)`$ for an orientation of the magnetization in $`x`$, $`y`$, and $`z`$ direction. For systems with weak spin-orbit coupling we have $`T_x+T_y+T_z0`$, and if this precondition is fulfilled then the contribution of the magnetic dipolar term drops out when we average the spin sum rule over the three directions of the magnetization. However, for low-dimensional systems the effects of spin-orbit coupling are larger, and it might be that then the precondition for the application of the angle-dependent XMCD is no longer fulfilled. We have shown this explicitely for the extreme situation of monatomic Fe, Co, and Ni wires Ederer et al. (2003). A second problem is that a couple of assumptions and approximations had to be made to derive the XMCD sum rules (see, e.g., Ref. Ederer et al., 2002 and references therein) the validity of which is not guaranteed in all systems. The validity of the XMCD sum rules has been confirmed, e.g., for bulk Fe and Co Chen et al. (1995), but for atoms at free surfaces Wu et al. (1994) and at interfaces Ederer et al. (2002), i.e., for sites with low symmetry, the application of the sum rules may be critical. It is therefore important to figure out the validity of the sum rules for low-dimensional systems.

In Ref. Komelj et al., 2002 we performed a systematic ab initio study for the spin and orbital moments and for the validity of the XMCD sum rules in Co systems with various dimensionality, i.e., for hcp Co, a Co monolayer on Pt (111) and a free standing Co (111) monolayer, a Co monatomic wire at the steps of a Pt (111) surface as well as for a free-standing monatomic Co wire. In the present Brief Report we discuss the respective data for the case of Fe. Corresponding XMCD measurements are on the way for bulk Fe, an Fe monolayer on Pt and a monatomic Fe wire on Pt Enders and Kern .

The validity of the sum rules was tested by calculating the orbital and spin moments on the one hand directly by the ab initio density functional electron theory, on the other hand the absorption spectra were determined by the same ab initio theory and then a second set of orbital and spin moments was obtained via the sum rules. Deviations between the two sets of data show that the sum rules are violated. We performed ab-initio calculations by using the local-spin-density approximation Perdew and Wang (1992) (LSDA) in a combination with spin-orbit coupling, without and with an orbital-polarization (OP) term Eriksson et al. (1990). The orbital polarization term corrects at least in part explicitely for those electronic correlations which are responsible for the orbital polarization effects (in free atoms described by Hund’s second rule) and which are not appropriately described by the LSDA. We used the tight-binding linear-muffin-tin-orbital (LMTO) method in the atomic sphere approximation Andersen and Jepsen (1984), and the WIEN97 code Blaha et al. (1990) which adopts the full-potential linearized-augmented-plane-wave (FLAPW) method Wimmer et al. (1981) in which the spin-orbit coupling and the tools for the calculation of the XMCD spectra Kuneš et al. (2001a, b) and the orbital polarization term Rodriguez et al. (2001) have been implemented. A supercell geometry Komelj et al. (2002) was used for all calculations with perpendicular magnetization for the case of the monolayers and wires. For the monolayer the supercell consists of two Pt and one Fe (111) layers in the fcc stacking and a vacuum sheet on top of the Fe layer corresponding to two empty layers. The vicinal Pt surface with the Fe wires at the steps is modelled by the supercell shown in Fig. 3. To focus on the pure effect of the dimensionality we fixed the nearest-neighbor distances of the atoms for all considered systems to the one of fcc Pt (2.77Å), as it was done in the LSDA study of finite Co chains on a plane Pt (111) surface by Lazarovits et al. Lazarovits et al. (2003). (In contrast, Spišak and Hafner Spišak and Hafner (2002) have taken into account the structural relaxation effects in their LSDA study of monatomic Fe wires on vicinal Cu surfaces.)

Fig. 1 summarizes the results for the Co systems. For the spin moments nearly the same results were obtained when taking into account the OP term or when neglecting it. We therefore only show the data obtained from calculations with the OP term. Six types of calculations were performed:

1. A direct calculation of the spin moments (filled circles connected by full lines).
2. A calculation of the spin moments based on the spin sum rule, thereby including the $`T_z`$ term (filled triangles connected by full lines).
3. A calculation of the spin moments based on the spin sum rule, thereby neglecting the $`T_z`$ term (filled diamonds connected by full lines).
4. A direct calculation of the orbital moments in LSDA (open squares connected by dashed lines).
5. A direct calculation of the orbital moments in LSDA+OP (open circles connected by dashed lines).
6. A calculation of the orbital moments in LSDA+OP based on the orbital sum rule (open triangles connected by dashed lines).

The comparison of the results for the LMTO and the FLAPW calculations given in Figs. 1,2 shows that the general trends are essentially the same for both types of calculations. Quantitatively, there are some differences which result mainly from the fact that the LMTO method adopts a spherical approximation for the effective potential in each atomic sphere (ASA potential) whereas in the FLAPW calculation the full asphericity of the effective potential is taken into account. In the LMTO-ASA method the $`T_z`$ term is determined from the nonspherical charge and spin density obtained after the last iteration step. Because this charge and spin density is calculated from an ASA potential, the influence of the nonspherical parts of the effective potential within the spheres on the asphericity of the charge and spin density is neglected. Therefore, for a quantitative discussion the FLAPW results are more reliable.

The most important results for the Co systems are:

1. The spin moments increase only slightly with decreasing dimensionality.
2. The orbital moments increase strongly with decreasing dimensionality when we take into account the OP term.
3. The $`T_z`$ term appearing in the spin sum rule is only relevant for the monolayers, otherwise it is rather small (see the more realistic FLAPW data), even for the monatomic wires. It has been outlined in Ref. Komelj et al., 2002 that this results from the special band filling for the case of Co wires.
4. Both the spin sum rule and especially the orbital sum rule are rather well fulfilled for all Co systems, irrespective of the dimensionality.

Fig. 2 shows the corresponding results for the case of Fe. Again we want to focus on the pure influence of the dimensionality and we therefore fixed the nearest-neighbor distances of the atoms for all the considered systems to the one of fcc Pt. For bulk Fe we considered the hypothetical ferromagnetic fcc phase. The four most important results are:

1. The spin moments increase only slightly with decreasing dimensionality, as in the case of Co.
2. The orbital moments increase strongly with decreasing dimensionality when we take into account the OP term, as in the case of Co.
3. The contribution of the $`T_z`$ term to the spin sum rule is considerable for the monatomic layers and also for the monatomic wires. This will represent a problem for the determination of the spin moment from experimental XMCD spectra for Fe monolayers and Fe monatomic wires: The $`T_z`$ term cannot be neglected in the spin sum rule, and — at least for the monatomic wires — it cannot be determined safely from angle-dependent XMCD measurements because the above discussed precondition for the application of this technique is not fulfilled Ederer et al. (2003).
4. There is a tendency that the spin and orbital sum rules are a bit less well fulfilled than for the Co systems, but the violations are still moderate.

For an interpretation of the experiments it is also interesting to know the spin and orbital moments induced by the Co or Fe atoms in the Pt substrate, because the Pt atoms contribute to any magnetic property of the system, for instance, to the total magnetization and especially to the magneto-optical properties in the visible regime of light. Experimentally, they can be separated from the spin and orbital moments of the Co and Fe atoms by the XMCD technique Schütz et al. (1987).

Table 1 shows the results obtained from the LSDA calculations for the Pt supported monolayers including the OP term (very similar results are obtained for the pure LSDA calculation). For the Fe monolayer the induced Pt moments are a bit smaller than for the Co monolayer, especially for the second Pt layer.

It should be recalled that no relaxation effects are taken into account and that because of the supercell geometry there is a vacuum sheet below the second Pt layer. Therefore, the quantitative results should not be taken too literally when comparing with experimental data, but the calculations certainly yield the correct order of magnitude for the polarization effect. A similar discussion holds also for the case of the wires on Pt, for which the results for the monatomic Fe wire obtained by the LMTO method are given in Fig. 3. It is interesting to note that for this geometry we get nearly exactly the same results as for the case of a monatomic Co wire.

Altogether, the ab-initio electron theory does not predict a tremendous difference in the physics of Co and Fe monatomic wires although an eventual application of the XMCD spectroscopy on the Fe wires may be less reliable, particularly due to the discussed importance of the magnetic dipolar term.


format: text

text:
# The quest for hot gas in the halo of NGC 1511

## 1 Introduction

The most recent generation of X-ray satellites, XMM-Newton and Chandra, has a much higher sensitivity than previous X-ray missions. This enables us to conduct observations of fainter targets than before and to obtain a much more detailed picture of those known already. Increased sensitivity leads to significant progress in investigations of low surface-brightness emitters, such as for example hot gaseous halos around actively star-forming spiral galaxies. As part of the XMM-Newton Guaranteed Time program we have observed the edge-on galaxy NGC 1511. Here we report on the results from these observations.

The XMM-Newton observations of nearby starburst galaxies form part of a multi-wavelength study of all phases of extraplanar gas in external galaxies, which is conducted in order to assess the importance of halos as repositories of a metal-enriched medium and their significance in terms of galactic chemical evolution and possible metal enrichment of the intergalactic medium (IGM; e.g. Heckman et al. 1990; Dahlem et al. 1995, 1998; Weaver et al. 2000; Heckman 2001). Soft X-ray surveys of nearby spiral galaxies, including starburst systems, based on ROSAT data were published e.g. by Read et al. (1997) and Read & Ponman (2001).

NGC 1511 is an SAa pec: H II type galaxy with warm dust, as traced by its high far-infrared flux ratio. Its salient properties are collected in Table 1 (from Lehnert & Heckman 1995; scaled to $`D=17.5`$ Mpc). These values indicate that NGC 1511 is an intermediate-mass galaxy (more massive than M 82, but less massive than NGC 253) with a high star formation rate (SFR). Optical spectroscopy was published by Lehnert & Heckman (1995) and, recently, by Kewley et al. (2001), establishing its starburst nature. Two extraplanar gas phases were previously detected in NGC 1511: H$`\alpha `$ emission from warm ionized gas (Lehnert & Heckman 1995) and nonthermal radio continuum emission from cosmic-rays and magnetic fields (Dahlem et al. 2001). These detections, together with the high surface brightness of tracers of star formation (SF) in its disk, indicate that it is an actively star-forming galaxy. It is amongst the first starburst galaxies observed with XMM-Newton, together with NGC 253 (Pietsch et al. 2001), NGC 3628 and M 82 (unpublished). Chandra observations of NGC 253, M 82 and NGC 3628 were published by Strickland et al. (2000, 2002), Kaaret et al. (2001) and Matsumoto et al. (2001).

For an overview of papers related to investigations of galactic gaseous halos, refer to Dettmar (1992) and Dahlem (1997). A description of the XMM-Newton observatory has been given by Jansen et al. (2001).

## 2 Observations and data reduction

Our observations of NGC 1511 were carried out in 2000, July 7–8. The total on-source integration time of the European Photon Imaging Camera (EPIC) pn camera (Strüder et al. 2001) is 40.0 ks. The pn camera performed the observation in extended full frame mode and the two MOS cameras (Turner et al. 2001) were put in full frame mode. For all cameras the thin filter position was chosen, offering the highest possible transmission for very soft X-ray emission. The data were calibrated using the XMM-Newton Science Analysis System (SAS), version 5.3.3, and the Current Calibration File (CCF) version of 2002-10-10. Good time intervals were defined when the total count rate, including background, was $`<1`$ cts s-1. After removal of data from time intervals with high background count rates a net integration time of 32.1 ks remains (for the pn camera).

Corrections were applied to the data for the effects of gain variations of the preamplifier channels and charge transfer inefficiency (CTI), bad pixels, chip edges, split events etc. and the counts were converted to the pulse-height invariant (PI) energy scale. Also, in order to fill the gaps in the EPIC pn camera images, images in the corresponding energy bands were produced from the data of the two MOS cameras and all images aligned and superimposed onto each other. The final images presented here are a combination of data from all three EPIC cameras. Photons with energies below 0.2 keV were not used in the data analysis because of the calibration uncertainties in this very soft energy regime. Derived data products (images and spectra) were produced from the final event lists, after application of these corrections.

## 3 X-ray Results and discussion

In this paper we present only X-ray data from the European Photon Imaging Camera (EPIC; from the pn and both MOS cameras). No Reflection Grating Spectrometer (RGS, den Herder et al. 2001) results are used here, because the target is too faint and too extended for high-resolution X-ray spectroscopy.

### 3.1 Soft X-ray imaging

Our observations are the first dedicated X-ray observations of NGC 1511. Using the ROSAT All-Sky Survey (RASS) Browser we verified that it was not detected in the RASS. In the data reduction we followed roughly the procedure adopted by the XMM-Newton Survey Science Centre (SSC, Watson et al. 2001) to split up the total passband into several energy bands. The only difference is that we chose to split the 0.5–2.0 keV band into two bands, see Table 2. In one case we combined two SSC bands to show an image with energies ranging from 2.0 keV to 7.5 keV.

In Fig. 1 we display X-ray images of NGC 1511. The photon distributions which were received in the individual energy bands as defined above were smoothed with an adaptive filter so as to enhance low surface brightness features, while not smearing out the flux of strong point sources. Such images are not suitable for quantitative analyses. The left panel is an overlay of the adaptively-smoothed XMM-Newton 1.3–2.0 keV image on a second generation Digitized Sky Survey (DSS-2) red image. The right panel shows an overlay of the adaptively-smoothed XMM-Newton 2.0–7.5 keV image on the same DSS-2 image of NGC 1511. Contours start at the 2.5-$`\sigma `$ confidence level above the background and are spaced by factors of $`\sqrt{2}`$. The 2.5-$`\sigma `$ levels of the images displayed here are: $`\mathrm{1.14\hspace{0.17em}10}^6`$ cts s-1 arcmin-2 (0.5–1.3 keV), $`\mathrm{0.62\hspace{0.17em}10}^6`$ cts s-1 arcmin-2 (1.3–2.0 keV) and $`\mathrm{1.20\hspace{0.17em}10}^6`$ cts s-1 arcmin-2 (2.0–7.5 keV).

The X-ray emission from NGC 1511 comprises different components, namely hard point sources and also soft extended emission. Extended soft emission is visible across the disk of the galaxy, with a maximum at the position of a prominent SF region on its eastern edge. A weak secondary maximum is visible near the centre of the galaxy. In addition, the soft X-ray emission distribution is extended perpendicular to the disk plane. The strong X-ray point source about $`30^{\prime \prime }`$ north of the centre of NGC 1511, which we will call–according to the IAU convention–XMMU J035936.6–673736, does not have an optical counterpart in either our XMM-Newton Optical Monitor (OM, Mason et al. 2001) images or the DSS. If associated with NGC 1511, it might be an ultra-luminous ($`L_\mathrm{X}=\mathrm{1.18\hspace{0.17em}10}^{40}`$ erg s-1) X-ray source (ULX, see e.g. Strickland et al. 2001, Foschini et al. 2002 and references therein). If not related to NGC 1511, it could be a background source, e.g. an AGN (see Sect. 3.2.4).

The emission distribution will be discussed in more detail below, in comparison with other wavebands.

### 3.2 Integral wide-band X-ray spectroscopy

#### 3.2.1 Spectral extraction

Spectroscopy was performed based on the EPIC pn event list only. Extracting valid (filter #XMMEA\_EP) single and double photon events from NGC 1511, an integral source spectrum was obtained. A separate spectrum was extracted for XMMU J035936.6–673736. The source photon extraction region for NGC 1511 is an ellipse with a major axis radius of $`67^{\prime \prime }`$, a minor axis radius of $`50^{\prime \prime }`$ and a position angle of $`45^{}`$, centred on the nucleus of NGC 1511. From this, a circular region of radius $`15^{\prime \prime }`$ around $`\alpha (2000)`$=03:59:36.6, $`\delta (2000)`$=–67:37:36.8 was excised for separate spectral analysis of the emission from XMMU J035936.6–673736. A second, elliptical region was defined near NGC 1511, on a neighbouring pn chip, to determine the local background.111The background photon extraction region is an ellipse equal in size to that for the source photons, but centred at $`\alpha (2000)`$=03:59:05.40, $`\delta (2000)`$=–67:41:09.9. Extracted spectra were grouped based both on a minimum number of counts (25) for each channel and on the maximum number of spectral channels sampling the pn energy resolution (3). An analysis of the background-subtracted integral spectrum of NGC 1511 over the energy range 0.2–12.0 keV was performed using Xspec v. 11.1.0.

#### 3.2.2 The need for a multi-component fit

Our first attempts, with single emission components plus foreground absorption, for all different kinds of spectral models, provided statistically unacceptable fit results. Therefore, we subsequently tried a combination of a Mewe-Kaastra plasma with Fe-L lines by Liedahl (“Mekal”) model (Mewe et al. 1985, Mewe et al. 1986, Kaastra 1992) plus a power law, also allowing the two components to undergo different absorption. A thermal plasma was chosen to represent the diffuse emission component, while the power law approximates the expected spectrum of the point source(s), which, if associated with NGC 1511, are most likely to be binaries, otherwise background AGN. This improved the quality of the fit (see Table 3), but still some residual emission remained around 0.6, 1.0 and 2.0 keV. Adding to the spectral model a second thermal component leads to a statistically marginally better fit, but improved visibly the fit in the areas that showed residuals before. Thus, the spectrum of NGC 1511 clearly requires a multi-component fit. However, the one presented here is not necessarily the only possible good fit to the spectrum (Table 3). At this time, we can only state that a single thermal component plus power law spectral model cannot explain the entire soft X-ray spectrum of NGC 1511, demonstrating the complexity of the emission spectrum. The nature of the additional third component is as yet uncertain.

#### 3.2.3 Best-fitting spectral model for NGC 1511

The three-component fit with two thermal components described in Table 3 is displayed in Fig. 2. Two thermal components, with their metallicities fixed at 1.0 $`Z_{}`$, were fitted. Fitting the metallicities of the thermal components as a free fit parameter did not statistically improve the fit quality and the $`Z`$ values are not constrained sensibly. Therefore, this fit is not shown here. In the fit with fixed metallicities, which we consider the simplest best-fitting model, the relative contribution of the two thermal components and the power law to the total flux are 12% (0.19 keV), 11% (0.59 keV) and 77%, respectively. These relative flux contributions conceal the fact that diffuse thermal emission dominates the total emission in the soft regime, while at hard energies the power-law component is dominant. We have verified by means of spatially resolved spectroscopy of subregions that the thermal emission corresponds with the extended soft component visible in our images. The best-fitting absorbing column densities are higher than the Galactic foreground ($`N_\mathrm{H}=\mathrm{0.5\hspace{0.17em}10}^{21}`$ cm-2, Dickey & Lockman 1990), suggesting additional internal absorption in NGC 1511. Based on this fit, we determine a total 0.2–12 keV flux of NGC 1511 of $`f_\mathrm{X}=\mathrm{3.04\hspace{0.17em}10}^{13}`$ erg s-1 cm-2. This X-ray flux translates to a 0.2–12 keV X-ray luminosity of NGC 1511 of $`L_\mathrm{X}=\mathrm{1.11\hspace{0.17em}10}^{40}`$ erg s-1. ¿From this value a far-infrared–to–X-ray luminosity ratio of log$`(L_{FIR}/L_\mathrm{X})=3.62`$ follows, which is typical for starburst galaxies (Heckman et al. 1990; Read & Ponman 2001). The total diffuse soft X-ray luminosity of NGC 1511 is $`L_\mathrm{X}=\mathrm{2.55\hspace{0.17em}10}^{39}`$ erg s-1. This leads to a ratio of log$`(L_{\mathrm{FIR}}/L_{\mathrm{X},\mathrm{diff}})=4.25`$. The amount of diffuse soft X-ray emission in NGC 1511, compared to other galaxies, is thus relatively low (Read & Ponman 2001; their Fig. 5, central panel). The temperatures of the thermal plasma components are in line with earlier measurements of other starburst galaxies, based on ROSAT and ASCA data (e.g. Dahlem et al. 1998; Weaver et al. 2000 and references therein). Components at about 0.2–0.3 keV and 0.6 keV were found by us in several starburst galaxies, with the 0.2–0.3 keV component dominant in the halos, while the hotter component is normally observed from gas in the disks. Similar results were also obtained based on previously published data from Chandra and XMM-Newton (see references in Sect. 1).

#### 3.2.4 Model fit to the spectrum of XMMU J035936.6–673736

The best-fitting model spectrum of XMMU J035936.6–673736 (not displayed) is a power-law with $`\mathrm{\Gamma }=1.92_{0.11}^{+0.14}`$ and an absorbing column density of $`N_\mathrm{H}=2.5_{0.3}^{+1.0}\mathrm{\hspace{0.17em}10}^{21}`$ cm-2. The 0.2–12 keV X-ray flux of this source is $`f_\mathrm{X}=\mathrm{3.22\hspace{0.17em}10}^{13}`$ erg s-1 cm-2. The fitted moderately low absorbing column density makes it unlikely that this source is a background AGN, which one would expect to have much higher intrinsic absorption.

### 3.3 UV imaging

The OM onboard XMM-Newton, with its three UV filters in addition to a range of optical filters, opens a second window on the targets observed that is not accessible from the ground, in the UV. As an example, we present in the left panel of Fig. 3 a 230 nm image of NGC 1511 obtained with the UVM2 filter.

This image exhibits the UV emission from massive stars, which are the primary sources of energy heating the warm ionized gas in NGC 1511. This is the reason for the global correspondence between the emission maxima in the UV and H$`\alpha `$ images. However, it is clear from a comparison of the two frames (Fig. 3 \[left panel\] and Fig. 4 \[right panel\]) that there is more diffuse H$`\alpha `$ outside the H II-regions. The presence of such massive stars is yet another indicator of NGC 1511’s activity, as already suggested by its $`L_\mathrm{X}/L_{\mathrm{FIR}}`$ ratio and also by the presence of diffuse halo emission in various wavebands.

In the overlay with the EPIC 0.5–1.3 keV soft X-ray image in the right panel of Fig. 3 one can see that there is basically no spatial correspondence between UV and X-ray emitters in NGC 1511.

The UV point source south-east of the centre of NGC 1511 is probably an unrelated foreground star. It also shows up on an R-band image by Lehnert & Heckman (1995). Its position does not coincide with that of SN1935C or Nova Hydri 1935 (van den Bergh & Hazen 1988).

An OM image obtained with the UVW2 filter (central wavelength about 210 nm), exhibits the same general emission distribution as Fig. 3, however with a lower signal-to-noise ratio. Therefore, it is not displayed here.

## 4 A multi-wavelength picture of NGC 1511

The observations described here are part of a large project aiming at understanding the nature of gaseous halos around actively star-formating galaxies.

In the case of NGC 1511 it is interesting to compare the XMM-Newton images with others obtained in different wavebands. Especially a comparison of the OM image with H$`\alpha `$ and 2 Micron All-Sky Survey (2MASS; Skrutskie 2001) near-infrared frames is useful. Both, the H$`\alpha `$ frame by Lehnert & Heckman (1995) and the 2MASS J-H-Ks-band composite image are displayed in Fig. 4.

### 4.1 Extraplanar emission

The X-ray disk emission follows roughly that of H$`\alpha `$ and radio continuum (Lehnert & Heckman 1995; Dahlem et al. 2001). In particular the correspondence between soft X-ray and radio emission is obvious. The axial ratio of both emission distributions is much smaller than that of the optically visible galaxy and there are two emission maxima in the disk, one near the centre and one in the south-eastern part of the disk. NGC 1511’s disk is disturbed and the disk plane not well-defined, as visible e.g. in the H$`\alpha `$ by Lehnert & Heckman (1995). Also the NIR 2MASS images suggest disturbances of the stellar distribution. This makes it difficult to assess whether emission arises from inside or beyond the disk. Still, the soft X-ray emission to the north of the central region comes from an area beyond the optically visible galaxy, thus possibly from the halo. This corroborates our earlier finding (e.g. Dahlem 1997) that, once gaseous halos are detected around actively star-forming spirals, different phases of the ISM are present. H$`\alpha `$, radio continuum and diffuse soft X-ray emission (i.e., warm ionized gas, cosmic rays, magnetic fields and hot ionized gas) are found to be associated with each other, although not necessarily directly correlated.

### 4.2 Disturbances of NGC 1511

The H$`\alpha `$ frame shows quite impressively how strong the starburst activity in NGC 1511 is and also how disturbed the galaxy is. This disturbance is most likely caused by its companions, NGC 1511a and NGC 1511b, which are observed at projected distances from NGC 1511 of 56.7 kpc and 38.7 kpc, respectively, with very small relative recession velocities of only 38(25) km s-1 with respect to NGC 1511. Using a high-contrast display of the corresponding red DSS-2 plate one can see a tidal tail connected to NGC 1511 (not visible in Fig. 1, because with such a grey-scale no contours could be seen), from the south-eastern end of the galaxy pointing roughly northwards, suggesting an ongoing tidal interaction.

The starburst nature of NGC 1511 is also corroborated by an updated radio spectral index calculation that we performed, using our earlier measurements (Dahlem et al. 2001), plus a 4.85 GHz measurement by Gregory et al. (1994) and an 843 MHz flux density measurement from the Sydney University Molonglo Sky Survey (SUMSS; Bock et al. 1999). $`S_{4.85GHz}=76`$ mJy and $`S_{0.843GHz}=237`$ mJy lead to a new, improved value of $`\overline{\alpha }=0.66\pm 0.08`$ ($`S\nu ^{\overline{\alpha }}`$), a value typical of actively star-forming galaxies (Condon 1992). This new value is, within the ample error bar of our earlier calculations, consistent with these. The Molonglo Telescope 843 MHz map (not displayed) also corroborates the emission distribution of our 1.43 GHz and 2.45 GHz maps (Dahlem et al. 2001) with a low major-to-minor axial ratio, indicating the presence of a synchrotron radio halo.

The 2MASS image suggests that NGC 1511’s appearance as a high-inclination ring system, as seen in the UV image in Fig. 3, might be caused not by such an emission distribution, but by absorption by a massive dust lane that is optically thick even in the NIR, creating a central emission depression.

This multi-wavelength approach demonstrates that NGC 1511 is a heavily disturbed starburst galaxy, which implies that it is not a good candidate for investigations of the dependence of its gaseous halos’ properties on the level and distribution of star formation in the underlying disk, for which one would prefer to study relatively undisturbed systems. Processes other than energy input from stellar sources (stellar winds and supernovae; see Leitherer & Heckman 1995) might influence the gas motions vertical to the disk, such as e.g. a gravitational potential modified by a tidal interaction with a companion galaxy. Still, NGC 1511 is one of a small number of starburst galaxies in which various components of the ISM are observed, within the same volume, above their disk planes. A significant contribution to the heating of its ISM from a hidden AGN can be excluded based on our X-ray data due to the absence of a hard central source, from optical spectra because of the absence of typical AGN emission lines and from our radio data, which also do not show any signatures of AGN activity at the centre.

### 4.3 Related observations

NGC 1511, amongst other objects, is also being observed by us in H I line emission with the Australia Telescope Compact Array to study the properties of its neutral atomic gas. This will lead to new information about the intrinsic absorbing gas distribution and its kinematics. It will, for example, also help to elucidate the issue whether the observed disturbances of NGC 1511’s disk can be explained by interactions with its companions, NGC 1511a and NGC 1511b.

XMM-Newton observations of several more nearby starburst galaxies, including NGC 1808, NGC 3628 and NGC 4666, have been conducted by us. These, together with data from the XMM-Newton archive, will be used to study systematically and in a self-consistent manner the soft X-ray emission properties of galaxies with kpc- to galaxy-scale outflows from their disks.

Acknowledgements: We thank Dr. M. Lehnert for making available to us his H$`\alpha `$ and R-band images of NGC 1511 in digital form. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions.


text:
# An Effective Membrane Model of the Immunological Synapse

## Abstract

The immunological synapse is a patterned collection of different types of receptors and ligands that forms in the intercellular junction between T Cells and antigen presenting cells (APCs) during recognition. The synapse is implicated in information transfer between cells, and is characterized by different spatial patterns of receptors at different stages in the life cycle of T cells. We obtain a minimalist model that captures this experimentally observed phenomenology. A functional RG analysis provides further insights.

T lymphocytes (T cells) are the orchestrators of the adaptive immune response in complex organisms. A key event during activation of the immune response is T cell recognition of cells that display peptides derived from foreign antigens on their surface Abbas et al. (2000). Recent experiments Monks et al. (1998); Grakoui et al. (1999); Krummel et al. (2000); Lee et al. (2002); Stoll et al. (2001) have vividly demonstrated that during this recognition process a highly organized pattern of different types of receptors and ligands forms in the intercellular junction between T cells and antigen presenting cells. This recognition motif is several microns in diameter, and since it is implicated in information transfer between the cells, it is called the immunological synapse. Formation of a synapse is also characteristic of an earlier stage in the life cycle of T cells. Immature T cells (thymocytes) are selected in the thymus so that they are not activated by peptides derived from the organism itself Abbas et al. (2000). In the thymus, thymocytes interact with cells that display self peptides on their surface. Thymocytes that bind strongly are deleted by apoptosis. The synapses formed during thymocyte selection are distinctly different in character Hailman et al. (2002); Ritchie et al. (2002) from those observed during mature T cell activation. Understanding the mechanisms via which synapses form under different circumstances and the biological purpose of creating different spatial patterns of cell surface receptors are active areas of research. In addition to the biological significance, an understanding of these issues may also inspire the creation of synthetic mimics that could carry out biomimetic recognition tasks which could be useful in applications such as targeted drug delivery.

In this letter, starting from a model proposed by Qi et al. Qi et al. (2001), we develop a minimalist model that captures some of the essential physics of synapse formation when apposing membranes contain complementary pairs of receptors and ligands. The model allows us to calculate a phase diagram which delineates the conditions that lead to a transition from synaptic patterns characteristic of mature T cells to those observed during thymocyte selection. This phase diagram may serve as a guide for the design of synthetic analogs.

Consider two membranes containing complementary pairs of receptors and ligands. The intramembrane motion of receptors and ligands is determined either by diffusion or a directed velocity toward the center of the junction Wulfing and Davis (1998). Complementary receptors and ligands can bind to each other if apposed. Different receptor-ligand complexes have different topographical size Dustin and Shaw (1999), and hence the rate of association depends upon local intermembrane separation. The receptor-ligand complexes also dissociate at a prescribed rate. For the case where there are only two types of receptors-ligand pairs (TCR-MHCp and LFA1-ICAM1 in the case of the T cell immunological synapse), the pertinent equations are:

$`{\displaystyle \frac{C_T}{t}}=D_T^2C_Tk_{on}(z)C_TC_M+k_{off}C_{TM}\stackrel{}{}(\stackrel{}{V}C_T),`$
$`{\displaystyle \frac{C_M}{t}}=D_M^2C_Mk_{on}(z)C_TC_M+k_{off}C_{TM},`$
$`{\displaystyle \frac{C_{TM}}{t}}=D_{TM}\left[C_{TM}+{\displaystyle \frac{1}{k_BT}}C_{TM}{\displaystyle \frac{\delta F}{\delta C_{TM}}}\right]+k_{on}(z)C_TC_Mk_{off}C_{TM},`$
$`{\displaystyle \frac{C_L}{t}}=D_L^2C_Lk_{on}^{}(z)C_LC_I+k_{off}^{}C_{LI},`$
$`{\displaystyle \frac{C_I}{t}}=D_I^2C_Ik_{on}^{}(z)C_LC_I+k_{off}^{}C_{LI},`$
$`{\displaystyle \frac{C_{LI}}{t}}=D_{LI}\left[C_{LI}+{\displaystyle \frac{1}{k_BT}}C_{LI}{\displaystyle \frac{\delta F}{\delta C_{LI}}}\right]+k_{on}^{}(z)C_LC_Ik_{off}^{}C_{LI}.`$
Here, $`C`$ represents the concentration of a given molecule, and $`z`$ is the separation of the T cell membrane from the APC membrane. The abbreviations used for different protein molecules are: TCR $``$ T, MHCp $``$ M, TCR-MHCp $``$ TM, LFA1 $``$ L, ICAM1 $``$ I, and LFA1-ICAM1 $``$ LI. The different species diffuse with the corresponding diffusion coefficients $`D`$. $`\stackrel{}{V}`$ is a directed velocity of TCR due to cytoskeletal motion Wulfing and Davis (1998). Since small values of $`\stackrel{}{V}`$ do not change the qualitative physics Hori et al. (2002), for simplicity we do not consider this further in this letter. The binding and dissociation chemical reactions are controlled by the rate constants $`k_{on}(z)`$ and $`k_{off}`$. The local intermembrane separation, $`z`$, evolves according to a time-dependent Ginzburg-Landau equation,

$$M^1\frac{z}{t}=\frac{\delta F}{\delta z}+\zeta ,$$
(2)
where $`M`$ is a phenomenological parameter for the rate at which membrane shapes respond to changes in the free energy, $`F=d^2x\left[\frac{\gamma }{2}(z)^2+\frac{\kappa }{2}(^2z)^2+C_{TM}B_{TM}(z)+C_{LI}B_{LI}(z)\right]`$. The first two terms are due to costs associated with deforming the membranes: If one of the membranes is supported and undeformed, $`\gamma `$ and $`\kappa `$ are the surface tension and bending rigidity of the free membrane; for similar membranes they are half the corresponding parameters for each membrane. $`B_{TM}(z)`$ and $`B_{LI}(z)`$ denote respectively the bond energy gained from forming TM and LI complexes for a specific intermembrane separation $`z`$.

The time scales associated with large-scale membrane shape changes can be much slower than those associated with binding kinetics and protein mobility. In such cases, for a given membrane shape, the local species concentrations are related by

$`{\displaystyle \frac{C_{TM}(z)}{C_TC_M}}={\displaystyle \frac{k_{on}(z)}{k_{off}}}={\displaystyle \frac{K_{on}}{K_{off}}}\mathrm{exp}\left[{\displaystyle \frac{B_{TM}(z)}{k_BT}}\right],`$
$`{\displaystyle \frac{C_{LI}(z)}{C_LC_I}}={\displaystyle \frac{k_{on}^{}(z)}{k_{off}^{}}}={\displaystyle \frac{K_{on}^{}}{K_{off}^{}}}\mathrm{exp}\left[{\displaystyle \frac{B_{LI}(z)}{k_BT}}\right].`$ (3)
Substituting Eq. (An Effective Membrane Model of the Immunological Synapse) in Eq. (2) obtains the following effective dynamical equation for $`z`$

$`M^1{\displaystyle \frac{z}{t}}=`$ $`\gamma ^2z\kappa ^4z`$
$``$ $`{\displaystyle \frac{K_{on}}{K_{off}}}C_TC_M\mathrm{exp}\left({\displaystyle \frac{B_{TM}(z)}{k_BT}}\right)B_{TM}^{}(z)`$
$``$ $`{\displaystyle \frac{K_{on}^{}}{K_{off}^{}}}C_LC_I\mathrm{exp}\left({\displaystyle \frac{B_{LI}(z)}{k_BT}}\right)B_{LI}^{}(z)+\zeta .`$ (4)
This corresponds to model A dynamics (Eq. (2)) for the order parameter $`z`$ with an effective free energy functional

$$\stackrel{~}{F}=d^2x\left[\frac{\gamma }{2}(z)^2+\frac{\kappa }{2}(^2z)^2k_BT\frac{K_{on}}{K_{off}}C_TC_Me^{\frac{B_{TM}(z)}{k_BT}}k_BT\frac{K_{on}^{}}{K_{off}^{}}C_LC_Ie^{\frac{B_{LI}(z)}{k_BT}}\right].$$
(5)
The reduction of the full set of equations to one determining the separation $`z`$ is reminiscent of the elimination of counter-ions in constructing the Poisson–Boltzmann equation, where the concentration of counter-ions is obtained from the Boltzmann weight of a potential, itself determined self-consistently from the counterion charge profile. Here $`z`$ plays the role of the potential, while the protein complexes are analogous to the counterions.

Bound receptor-ligand complexes do not favor membrane shape changes that deform the bond from its natural length. We thus use a harmonic approximation for the functional form of the bond elasticity energy, i.e. $`B_i(z)=B_i^0+\lambda _i(zz_i)^2/2`$, with $`i=1,2`$ for the two possible complexesQi et al. (2001). The maximal bond energies $`B_i^0`$ can be absorbed into $`K_{on}/K_{off}`$, and shall be ignored henceforth. The natural lengths of receptor-ligand proteins are taken to be those appropriate for the T cells ($`z_1=z_{TM}=15nm`$ and $`z_2=z_{LI}=42nm`$) Dustin and Shaw (1999); Garboczi et al. (1996), with corresponding ‘spring constants’ of $`\lambda _i`$ given below. With this approximation, the dynamical Eq. (4) takes the form (measuring time in units such that $`M=1`$)

$$\frac{z}{t}=\frac{\delta }{\delta z}\left\{d^2x\left[\frac{\gamma }{2}(z)^2+\frac{\kappa }{2}(^2z)^2+V(z)\right]\right\}+\zeta ,$$
(6)
with an effective potential

$$V(z)=k_BT\{\begin{array}{cc}C_1^{}\mathrm{exp}\left\{\frac{(zz_1)^2}{2\sigma _1^2}\right\}C_2^{}\mathrm{exp}\left\{\frac{(zz_2)^2}{2\sigma _2^2}\right\},\text{for }z>0\hfill & \\ c,\text{for }z<0.\hfill & \end{array}$$
(7)
In writing Eq. (7) we have assumed that $`\sigma ^2=k_BT/\lambda `$, which is appropriate if the curvature of the receptor-ligand interaction potential is the same in the bound state and at the activation barrier. Other assumptions do not alter the qualitative physics. The parameters $`C_1^{}=(K_{on}/K_{off})C_TC_M`$ and $`C_2^{}=(K_{on}^{}/K_{off}^{})C_LC_I`$ represent the number of complexes (TM and LI, respectively) formed if the intermembrane separation equals the corresponding optimal bond separation. As depicted in Fig. 1, the effective potential has two minima corresponding to the two possible bonds. Negative separations are prevented by a large energy cost $`k_BTc1`$ (hard wall in simulations).

Eq. (6) was solved numerically by a finite difference method. Initially, the upper membrane was assigned a parabolic shape in the contact region, and the intermembrane distance was kept constant ($`50nm`$) outside the contact region. For the membrane parameters we use experimentally determined values Simson et al. (1998) of $`\gamma =700k_BT/\mu m^2`$ and $`\kappa =400k_BT`$, and for the bond elasticity we take $`\sigma _1=13nm`$ and $`\sigma _2=5nm`$ Qi et al. (2001). Simulations were carried out for various ratios of $`C_1^{}/C_2^{}`$, and local concentrations of the TCR-MHCp and LFA1-ICAM1 were obtained from the intermembrane separation using Eqs. (An Effective Membrane Model of the Immunological Synapse).

Zero Noise: In the absence of noise ($`\zeta =0`$), Eq. (6) is a relaxational dynamics toward the minimum of the potential $`V(z)`$, subject to the imposed boundary conditions. For $`C_1^{}>C_2^{}`$ the minimum at $`z_{TM}`$ is deeper, and indeed in the numerical simulations a central cluster of TCR-MHCp complex forms quickly in the intercellular contact area. As the boundary of the contact region is kept at a constant height close to $`z_{LI}`$, the longer LFA1-ICAM1 moves toward the boundary giving rise to a structure very similar to the final synaptic pattern obtained in the experimental and numerical studies Grakoui et al. (1999); Qi et al. (2001); Chakraborty (2002); Weikl et al. (2002).

When the concentration $`C_1^{}`$ is below $`C_2^{}`$, the global minimum of $`V(z)`$ switches to $`z_{LI}`$. In this regime, the receptor-ligand complex LFA1-ICAM1 concentration dominates the entire intercellular junction, and no synaptic pattern is formed. In the noise-less (mean-field) limit the location of the transition between the two patterns is always at $`C_1^{}/C_2^{}=1`$. The patterns obtained and the conditions for transitions between them are quite different when $`\zeta `$ is finite.

Finite Noise: In this case, numerical simulations of Eq. (6) are performed with a random $`\zeta `$ with zero mean, and variance equal to $`k_BT`$ to mimic thermal noise. Starting with sufficiently large values of $`C_1^{}C_2^{}`$, we again find that TCR-MHCp aggregates at the center of the intercellular junction at long times. In Fig. (2a) a cross section of the upper membrane and the associated protein concentrations are shown. Very similar synaptic patterns have been observed Monks et al. (1998); Grakoui et al. (1999); Krummel et al. (2000); Lee et al. (2002); Stoll et al. (2001) for in vivo and in vitro experiments with mature T cells.

For double positive (DP) thymocytes, the TCR concentration is much lower compared to mature T-cells Hailman et al. (2002). This experimental fact corresponds to lower values of $`C_1^{}`$. When we numerically solve Eq. (6) with $`C_1^{}/C_2^{}<f_c`$, the mature synapse is never formed. Rather, we observe (see Fig. 2b) fluctuating patterns with fleeting clusters of TCR-MHCp forming at various locations in a sea of LFA1-ICAM1. Note that this pattern is different in character from the one observed for $`C_1^{}/C_2^{}<1`$ in the noiseless case. Such dynamic patterns have been observed in recent experiments using DP thymocytes Hailman et al. (2002); Ritchie et al. (2002), and can be understood using the simple picture of relaxational motion in a potential with two minima. As we lower the TCR concentration, the minimum at $`z=z_{TM}`$ becomes shallower, and eventually the minimum near $`z=z_{LI}`$ corresponds to the global equilibrium. However, thermal fluctuations allow sampling of the region near the minimum at $`z=z_{TM}`$. This will favor sporadic aggregation of TCR-MHCp at random places. As long as the barrier between the two minima is not very large, such transient fluctuations will be frequent.

The effects of thermal fluctuations can be estimated by expanding $`V(z)`$ to quadratic order around each of its minima. The total energy cost of a deformation of wavenumber $`p`$ is given by $`E_i(p)=\left(k_BTC_i^{}\sigma _i^2+\gamma p^2+\kappa p^4\right)|z(p)|^2`$. Integrating out these (Gaussian) fluctuations yields a free energy per unit area of

$$\frac{F_i^{}}{k_BT}=C_i^{}\frac{1}{2}\frac{dpp}{2\pi }\mathrm{ln}\left(\frac{k_BTC_i^{}\sigma _i^2+\gamma p^2+\kappa p^4}{2\pi k_BT}\right).$$
(8)
The location of the transition between the two classes of patterns can now be estimated from $`F_1^{}=F_2^{}`$. Thermal noise enhances the minimum with larger $`\sigma `$ (width), because of its larger entropy. Hence the more flexible bond is formed more easily when thermal noise is present Chen (2002). Fig. 3 shows that the larger the temperature, the smaller the value of $`f_c`$. This “phase diagram” could serve as a guide for the design of synthetic systems as it shows how to manipulate conditions such that synaptic patterns characteristic of mature T cells or thymocytes are realized.

Self consistency of the above quadratic expansion requires that the extent of the fluctuations should be less than the corresponding variance of the bond. The mean square width of the interface trapped around $`z=z_i`$ is given by

$$W_i^2=\frac{dpp}{2\pi }\frac{k_BT}{k_BTC_i^{}\sigma _i^2+\gamma p^2+\kappa p^4}.$$
(9)
Typical values of $`C_1^{}`$ range from $`O(10)(\mu m)^2`$ for thymocytes to $`O(100)(\mu m)^2`$ in mature T cells. For this range of parameters, the constant term in the denominator of Eq. (9) dominates, until it is cut-off by the bending rigidity at distances of the order of $`1\mu m`$. (The overall size of the synapse is of the order of $`10\mu m`$ Grakoui et al. (1999).) It is then easy to estimate the value of the integral, and we find that the ratios $`W_1/\sigma _1`$ and $`W_2/\sigma _2`$ are of the order of 0.1. While these estimated fluctuations are smaller than the width of the attractive potentials, they are significant enough to make nonlinear corrections to the above quadratic expansion important. Thus the phase diagram in Fig. 3, using Eq. (8) should be regarded as an approximation.

When anharmonic corrections are significant, an alternative approach is to employ a functional renormalization group (RG) scheme used Fisher and Huse (1985); Lipowsky and Fisher (1986, 1987) in the context of wetting transitions. It consists of the following steps: (1) The fluctuating field $`z(x)`$ is divided into two parts corresponding to small wavenumbers ($`0<|p|<\mathrm{\Lambda }/b`$) and large wavenumbers ($`\mathrm{\Lambda }/b<|p|<\mathrm{\Lambda }`$); (2) The large wavenumber fluctuations $`\mathrm{\Lambda }/b<|p|<\mathrm{\Lambda }`$ are integrated out to yield a coarse grained Hamiltonian; and (3) The system is rescaled to new coordinates $`x^{}=x/b`$ and $`z^{}(x^{})=z(x)`$. (The field $`z`$ does not need to be scaled in two dimensions.) The second step can only be done approximately, and a first order expansion of the potential $`V(z)`$ leads to the linear approximationFisher and Huse (1985)

$$V_b(z)[V(z)]=b^2_{\mathrm{}}^{\mathrm{}}\frac{dz^{}}{\sqrt{2\pi }\delta (b)}\mathrm{exp}\left[\frac{(zz^{})^2}{2\delta ^2(b)}\right]V(z^{}).$$
(10)
The width of the convolution, $`\delta (b)`$, is

$$\delta ^2(b)=_{\mathrm{\Lambda }/b}^\mathrm{\Lambda }\frac{dpp}{2\pi }\frac{k_BT}{\gamma p^2+\kappa p^4}=\frac{k_BT}{4\pi \gamma }\left[2\mathrm{ln}b\mathrm{ln}\left(\frac{1+\kappa \mathrm{\Lambda }^2/\gamma }{1+\kappa \mathrm{\Lambda }^2/(b^2\gamma )}\right)\right].$$
(11)
Note that $`\kappa `$ is an irrelevant operator in the RG sense (vanishing as $`b^2`$ under scaling). Indeed the second term in the above equation can be ignored asymptotically for $`b1`$, but may be important for quantitative comparisons when the bending rigidity is large.

Applying Eq. (10) to the potential in Eq. (7) yields the renormalized potential

$$V_b(z)b^2k_BT\left[\frac{\sigma _1C_1^{}}{\sqrt{\sigma _1^2+\delta ^2}}e^{\frac{(zz_1)^2}{2(\sigma _1^2+\delta ^2)}}\frac{\sigma _2C_2^{}}{\sqrt{\sigma _2^2+\delta ^2}}e^{\frac{(zz_2)^2}{2(\sigma _2^2+\delta ^2)}}+c\mathrm{erf}(z,\delta )\right].$$
(12)
The soft wall renormalizes into an error function of width $`\delta `$, while the variance of the Gaussian potential is increased by $`\delta ^2`$. Thus, the bigger the width of the potential, the deeper it gets due to the fluctuations, as it allows more flexibility (entropy) in bond formation. For large values of the width $`\delta `$ (highly fluctuating surfaces), the location of the transition point is shifted to $`C_1^{}/C_2^{}=f_c=\sigma _2/\sigma _10.38`$. Due to the factor of $`b^2`$ in Eq. (12) the net potential gets stronger under RG, and the linear approximation eventually breaks down.

The RG transformation should be extended to a rescaling factor $`b`$ at which the curvature at the minimum of the renormalized potential is such that the resulting fluctuations in width are comparatively small. Beyond this scale, a quadratic expansion around the minimum of the effective potential should be valid. As demonstrated previously, typical parameters for the synapse yield a bare potential at which the above criterion is (barely) satisfied. Thus the linear RG approximation is of limited practical value in this context. It may, however, yield interesting information about the global form of the phase diagram as the concentrations $`C_1^{}`$ and $`C_2^{}`$ are further reduced. From the double well shape of the potential, we anticipate that the overall phase diagram consists of a line of discontinuous transitions separating preferences in the synapse for TM or LI bonds. This line is located at $`C_1^{}/C_2^{}1`$ for large concentrations, and moves to $`C_1^{}/C_2^{}\sigma _1/\sigma _2`$ at lower concentrations. It presumably terminates at an Ising critical point when the fluctuations in width become of the order of the separation $`z_{TM}z_{LI}`$ between the two minima. Sufficiently small values of $`C_1^{}`$ and $`C_2^{}`$ should also result in an unbinding transition which is located at $`C_1^{}\sigma _1+C_2^{}\sigma _2\sqrt{k_BT/\gamma }`$. It is thus useful to make a more detailed study of the phase diagram of the model synapse as well as the dynamic trajectories as various phase boundaries are approached.


text:
# Untitled Document

This paper is withdrawn. A revised version, already accepted for publication in The Astrophysical Journal, will be submitted tomorrow.


text:
# Radiative Transfer on Perturbations in Protoplanetary Disks

## 1 Introduction

Planetary systems are formed from rotating protoplanetary disks, which are often modeled in steady-state, in vertical hydrostatic equilibrium, with gas and dust fully mixed and thermally coupled (Kenyon & Hartmann, 1987). Such disks describe very well the observed properties of T Tauri disks of age $``$1 Myr and typical mass accretion rates of $`10^8M_{}\mathrm{yr}^1`$ (Hartmann et al., 1998; D’Alessio et al., 1998). They are passive disks in the sense modeled by e.g. Chiang & Goldreich (1997), where the main source of photospheric heating is irradiation from the parent star, although viscous heating is still important at small stellar distances and near the midplane. The temperature in such disks is computed under the assumption that the upper surface of the disk is perfectly concave and smooth at all radii, which is a very good description of an unperturbed disk because thermal and gravitational instabilities are damped very efficiently (D’Alessio et al., 1999b).

We are interested in planet formation inside such disks, which means that we might be forced to abandon the above assumption of a perfectly smooth surface. The newly formed planet core will distort it, affect the heating and cooling of the disk locally, and could have significant consequences for the further growth and migration of protoplanets. As discussed in Sasselov & Lecar (2000), the distortion need only be large enough compared to the grazing angle at which the starlight strikes the disk. This small angle has a minimum at $`0.4\mathrm{AU}`$ and increases significantly only at very large distances. The depth of the depression due to the additional mass of the planet, $`m_p`$, will be proportional to $`(r_\mathrm{H}/h)^3`$, where $`r_\mathrm{H}`$ is the Hill radius, and $`h`$ is the local scale height, with a shade area dependent on the grazing angle (Sasselov & Lecar, 2000).

Generally planet formation has been treated numerically in 2-dimensional disks; to date, there is only one recent study of the 3-dimensional effect due to planets embedded in a protoplanetary disk (Bate et al., 2003). Bate et al. model the 3-dimensional response of a gaseous viscous disk with a simple locally isothermal equation of state, without taking into account irradiation from the central star.

In this paper we study in detail the development of a compression in a standard passive disk and the radiative transfer in it. We consider this to be the first step in developing a global 3-dimensional simulation, similar to the one by Bate et al. (2003), but with a realistic equation of state, especially near and inside the protoplanet’s Roche lobe. We are explicitly interested in the planetary growth process before a gap in the disk is formed, hence our approach is limited to $`m_p30\mathrm{M}_{}`$, which results in a small perturbation to the irradiated surfaces of a passive disk.

## 2 Heating Sources

A schematic of our adopted disk model is shown in Figure 1. The disk is flared, so that disk height increases faster than distance to the star, consistent with observational evidence and theoretical models (Kenyon & Hartmann, 1987; Chiang & Goldreich, 1997; D’Alessio et al., 1998). We assume that the dust is well mixed with the gas and is the primary source of opacity at wavelengths characteristic of stellar emission. The dust absorbs incident stellar radiation at the surface and re-emits at longer wavelength, where the opacity is much lower. We are interested in the region of the disk that is optically thick, at distances of $`10`$ AU (D’Alessio et al., 1999a). This gives rise to three different layers with smooth transitions in the disk: the uppermost is optically thin to both stellar and disk radiation, the interior is optically thick to both stellar and disk radiation, and the middle layer is a transition region that is optically thick to stellar radiation, but still optically thin to disk radiation. We shall refer to this optically thick/thin region as the photosphere. The reprocessing of stellar radiation in this region is important for determining the interior temperature.

Chiang & Goldreich (1997) calculate the temperature of only the thin/thin and thick/thick regions using energy conservation considerations, without detailing the temperature structure of the photosphere. D’Alessio et al. (1998,1999a) calculate the detailed temperature structure of the photosphere in disks around young stars in a self-consistent way, including radiative transport and convection. They find that radiative transfer is the dominant mechanism for energy transport in the disk photosphere, and the primary heating sources are stellar irradiation and viscous heating, as shown by comparison to Calvet et al. (1991), who use the Milne-Strittematter treatment of the superposition of solutions (Milne, 1930; Strittmatter, 1974). We shall consider stellar radiation and viscous heating to be the primary sources of heating, where the total temperature is arrived at by calculating the sum of the viscous and radiative fluxes, as $`T^4=T_v^4+T_r^4`$, where $`T_v`$ is the viscous temperature and $`T_r`$ is stellar radiation heating temperature.

### 2.1 Viscous Heating

The temperature due to viscous heating for a constant flux, gray atmosphere is

$$T_v^4=\frac{3F_v}{4\sigma _B}(\tau _d+2/3)$$
(1)
where $`\sigma _B`$ is the Stefan-Boltzmann constant, $`\tau _d`$ is the optical depth of the disk’s own radiation, and the viscous flux $`F_v`$ at a distance $`a`$, for a star of effective temperature $`T_{}`$, mass $`M_{}`$, and radius $`R_{}`$, accreting at a rate $`\dot{M}_a`$ is

$$F_v=\frac{3GM_{}\dot{M}_a}{4\pi a^3}\left[1\left(\frac{R_{}}{a}\right)^{1/2}\right]$$
(2)
(Pringle, 1981).

### 2.2 Stellar Irradiation

We shall use a modified form of the Milne-Strittmatter treatment to calculate the effect of heating by stellar irradiation on a three-dimensional perturbation on a passive accretion disk. We assume that the dust is well-mixed with the gas in the disk and that they are thermodynamically coupled so that the dust temperature and gas temperature are the same. In general, there is no guarantee that these two temperatures are the same, or that the dust is either well mixed or uniform in composition.

Following the Milne-Strittmatter treatment, the incident radiation from the star is considered to be at some characteristic wavelength, typically in the visible, and the radiation emitted by the disk will be at a longer characteristic wavelength, typically in the infrared. The incident, short wavelength radiation will be indicated by the subscript $`s`$ (stellar or scattered), and the re-emitted, long wavelength radiation will be indicated by the subscript $`d`$ (disk or diffuse). Some fraction of the energy of the incident stellar radiation, $`\sigma `$, is scattered at the same frequency. The remaining fraction, $`\alpha =1\sigma `$, is absorbed by dust in the disk and re-emitted at a longer characteristic wavelength, typically in the infrared. This gives the equation of radiative transfer for the stellar short wavelength radiation,

$$\frac{1}{\kappa _s\rho }\widehat{𝐤}I_s=I_s\frac{\sigma E_0e^{\tau _s}}{4\pi }\frac{\sigma I_s𝑑\mathrm{\Omega }}{4\pi }$$
(3)
where $`\rho `$ is the density, $`\kappa _s`$ is the opacity at short wavelengths, and $`E_0`$ is the incident energy flux given by

$$E_0=\sigma _BT_{}^4\left(\frac{R_{}}{a}\right)^2.$$
(4)
The zeroth moment of this equation is

$$\frac{1}{\kappa _s\rho }𝐅_s=\alpha 4\pi J_s\sigma E_0e^{\tau _s},$$
(5)
and the first moment, making use of the Eddington approximation, is

$$\frac{1}{\kappa _s\rho }J_s=\frac{3𝐅_s}{4\pi }.$$
(6)
The radiative transfer equation for the radiation absorbed and re-emitted within the disk is

$$\frac{1}{\kappa _d\rho }\widehat{𝐤}I_d=I_dB$$
(7)
where $`\kappa _d`$ is the optical depth at disk radiation wavelengths. The zeroth and first moments (using the Eddington approximation) are

$$\frac{1}{\kappa _d\rho }𝐅_d=4\pi (J_dB),$$
(8)
$$\frac{1}{\kappa _d\rho }J_d=\frac{3𝐅_d}{4\pi }.$$
(9)
#### 2.2.1 Plane-parallel Disk

In a plane parallel disk atmosphere, this problem is reduced to one dimension, and we can express everything as a function of optical depth normal to the surface, $`\tau _d`$. We can substitute

$$\frac{1}{\kappa _d\rho }\frac{d}{d\tau _d}\text{and}𝐅_{s,d}F_{s,d}.$$

If the angle of incidence of the stellar radiation is $`\mathrm{cos}^1\mu _0`$, then

$$\tau _d=\tau _s\mu _0/q$$
(10)
where $`q`$ is the ratio of opacities

$$q\frac{\kappa _s}{\kappa _d}.$$
(11)
The condition of zero net flux relates the short and long wavelength fluxes as

$$F_s+F_d=E_0\mu _0\mathrm{exp}(q\tau _d/\mu _0)$$
(12)
and the boundary conditions at the surface are that of an isotropic radiation field,

$$2\pi J_s(0)=F_s(0)\text{and}2\pi J_d(0)=F_d(0).$$
(13)
Now, equations (5), (6), (8), (9), and (12) are a closed system of linear differential equations that we can solve directly for the temperature, $`T_r=(\pi B/\sigma _B)^{1/4}`$. The solution is given in Calvet et al. (1991) as

$$B=\frac{\alpha E_0\mu _0}{4\pi }[C_1^{}+C_2^{}\mathrm{exp}(q\tau _d/\mu _0)+C_3^{}\mathrm{exp}(\beta q\tau _d)]$$
(14)
where $`\beta \sqrt{3\alpha }`$,

$`C_1^{}`$ $`=`$ $`(1+C_1)\left(2+{\displaystyle \frac{3\mu _0}{q}}\right)+C_2\left(2+{\displaystyle \frac{3}{\beta q}}\right),`$ (15)
$`C_2^{}`$ $`=`$ $`{\displaystyle \frac{(1+C_1)}{\mu _0}}\left(q{\displaystyle \frac{3\mu _0^2}{q}}\right),`$ (16)
$`C_3^{}`$ $`=`$ $`C_2\beta \left(q{\displaystyle \frac{3}{q\beta ^2}}\right),`$ (17)
and

$`C_1`$ $`=`$ $`{\displaystyle \frac{3\sigma \mu _0^2}{1\beta ^2\mu _0^2}},`$ (18)
$`C_2`$ $`=`$ $`{\displaystyle \frac{\sigma (2+3\mu _0)}{\beta (1+2\beta /3)(1\beta ^2\mu _0^2)}}.`$ (19)
Figure 2 shows how temperature varies with $`\mu _0`$ at the surface ($`\tau _d=0`$) and in the interior ($`\tau _d\mathrm{}`$). The temperatures are normalized to an effective radiant temperature $`T_0=(E_0/\sigma _B)^{1/4}`$, which is the blackbody temperature corresponding to the incident energy flux. Note that at the surface the temperature remains relatively constant with $`\mu _0`$ and is greater than $`T_0`$ due to the absorptive properties of the dust, in accordance with the “superheated” surface layer proposed in Chiang & Goldreich (1997). In contrast, the interior temperature is sensitive to the incident angle, especially at grazing incidence. Thus, perturbations in the surface that change the angle of incidence can significantly affect the interior temperature.

By definition, $`\tau _d=2/3`$ is where the disk becomes optically thick to its own radiation. For $`q1`$, the temperature at $`\tau _d=2/3`$ is can be approximated by the interior temperature, since as $`q\tau _d`$ becomes large, the exponentials in equation (14) vanish and $`B`$ becomes independent of $`\tau _d`$. This indicates that in the absence of viscous heating, the disk interior is vertically isothermal for $`\tau _d2/3`$.

#### 2.2.2 Disk with Perturbation

Now we consider radiative transfer on a perturbed disk. A point $`P(x,y,z)`$ within the disk “sees” radiative flux coming from the surface of the disk. Each area element $`\delta A`$ contributes a solid angle $`\delta \mathrm{\Omega }`$ of flux to $`P`$. Figure 3 shows a schematic of such an area element. The angle of incidence of radiation on this surface element is $`\mathrm{cos}^1\mu `$. If $`𝐤`$ is the vector from $`\delta A`$ to $`P`$ and $`\widehat{𝐧}`$ is the unit normal to $`\delta A`$, then $`\nu `$ is the cos of the angle between $`\widehat{𝐧}`$ and $`𝐤`$.

We approximate the flux contribution from $`\delta A`$ as if it were a surface emitting isotropically with intensity

$$I=\frac{F_d(\tau _d,\mu )}{\pi }.$$

and the corresponding contribution to the total flux is $`d𝐅=I\widehat{𝐤}d\mathrm{\Omega }.`$ Then the total flux at $`P`$ over all surface contributions is

$$𝐅_{\mathrm{tot}}=\frac{1}{\pi }F_d\widehat{𝐤}\delta \mathrm{\Omega }.$$
(20)
Since $``$ is a linear operator,

$$\frac{1}{\kappa _d\rho }𝐅_{\mathrm{tot}}=\frac{1}{\pi }\frac{dF_d}{d\tau _d}\widehat{𝐧}\widehat{𝐤}\delta \mathrm{\Omega }=4\pi \left(\frac{1}{\pi }J_d\nu \delta \mathrm{\Omega }\frac{1}{\pi }B\nu \delta \mathrm{\Omega }\right).$$
(21)
Defining

$`J_{\mathrm{tot}}`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle J_d(\tau _d,\mu )\nu \delta \mathrm{\Omega }}\text{and}`$ (22)
$`B_{\mathrm{tot}}`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle B(\tau _d,\mu )\nu \delta \mathrm{\Omega }}`$ (23)
we see that $`𝐅_{\mathrm{tot}}`$, $`J_{\mathrm{tot}}`$, and $`B_{\mathrm{tot}}`$ satisfy equations (8) and (9), and we can calculate the perturbed temperature as

$$T_r=\left(\frac{\pi B_{\mathrm{tot}}}{\sigma _B}\right)^{1/4}.$$
(24)
For a plane-parallel surface, $`\tau _d`$ and $`\mu `$ are constant, and $`\nu =\mathrm{cos}\theta `$, where $`\theta `$ is the angle with respect to the surface normal. We then recover the solution for the plane-parallel case.

This method of calculating radiative transfer can be applied to any three-dimensional disk configuration, including perturbations induced by protoplanets. Sufficiently massive protoplanets will open a gap in the disk, in which case the effects of shadowing and illumination across the gap can be calculated for all azimuthal angles. In the following section, we will calculate radiative transfer on small perturbations induced by planets too small to open a gap, but large enough to act on the disk locally.

## 3 Hydrostatic Equilibrium

Now we consider a perturbation induced by the gravitational potential of a protoplanet on a circumstellar disk. We shall consider only the gravitational effects of the protoplanet on hydrostatic equilibrium and assume that any resonant effects are washed out by gas drag.

### 3.1 Vertical Density Profile

In order to find the vertical structure of a non-self-gravitating gaseous disk orbiting a central star, we solve the equation of hydrostatic equilibrium:

$$P=\rho \mathrm{\Phi }$$
(25)
where $`P`$ is the gas pressure of the disk and $`\mathrm{\Phi }`$ is the gravitational potential of the central star. Since we are interested in the vertical structure, we consider only the $`z`$ components and assume that the disk is isothermal in the $`z`$ direction. For a disk without a perturbing protoplanet, $`\mathrm{\Phi }=GM_{}/a`$ where $`M_{}`$ and $`a`$ are the mass and distance of the star, respectively. If $`c_s`$ is the isothermal sound speed, then

$$c_s^2\frac{d\mathrm{ln}\rho }{dz}=\frac{GM_{}}{a^3}z.$$
(26)
If we assume that $`az`$ and define $`hc_sa^{3/2}/(GM_{})^{1/2}=(c_s/v_\varphi )a,`$ then

$$\rho (z)=\rho _0\mathrm{exp}\left(\frac{z^2}{2h^2}\right),$$
(27)
i.e. the density has a Gaussian distribution with a scale height of $`h`$ which is determined by the local temperature.

The insertion of a planet into a passive disk adds an additional term to $`\mathrm{\Phi }`$, representing the potential due to the planet:

$$\mathrm{\Phi }=GM_{}/a+Gm_p/r_p$$

where $`m_p`$ and $`r_p`$ are the mass and distance of the planet. Assuming that the planet is in the midplane ($`z=0`$), Eq. (26) becomes

$$\frac{d\mathrm{ln}\rho }{dz}=\frac{z}{h^2}\left(\frac{M_pa^3}{M_{}h^2}\right)\frac{z}{r_p^3}.$$
(28)
The Hill radius is defined as $`r_\mathrm{H}(M_p/3M_{})^{1/3}a,`$ and this equation integrates to

$$\rho (z)=\rho _1\mathrm{exp}\left[\frac{z^2}{2h^2}+\frac{3r_\mathrm{H}^3}{h^2r_p}\right]$$
(29)
where $`r_p=\sqrt{x^2+y^2+z^2}`$, with the planet as the coordinate origin. So long as $`r_\mathrm{H}h`$ and $`|z|>r_\mathrm{H}`$, the effect of the planet is a small perturbation on the density profile.

We normalize this equation by matching the density at $`z=0`$ to the unperturbed density, thus

$$\rho _p(z)=\rho _0\mathrm{exp}\left[\frac{z^2}{2h^2}+\frac{3r_\mathrm{H}^2}{h^2}\left(\frac{r_\mathrm{H}}{\sqrt{x^2+y^2+z^2}}\frac{r_\mathrm{H}}{\sqrt{x^2+y^2}}\right)\right]$$
(30)
Note that there is a singularity ar $`r_p=0`$. Figure 4 shows the shape of this density profile. In Figure 4a, we hold $`r_\mathrm{H}/h`$ fixed at $`0.5`$ and vary $`\sqrt{x^2+y^2}`$. For a disk with $`h=0.04\mathrm{AU}`$ at $`1\mathrm{AU}`$ around a central star of $`0.5\mathrm{M}_{\mathrm{}}`$, this corresponds to a protoplanetary mass of $`4\mathrm{M}_{}`$. For $`(x^2+y^2)r_\mathrm{H}^2`$, the perturbed density profile does not deviate significantly from the unperturbed density profile. For $`(x^2+y^2)r_\mathrm{H}^2`$, the density profile develops a sharper profile at $`z<h`$. In Figure 4b, we hold $`\sqrt{x^2+y^2}`$ fixed at the Hill radius and vary the mass of the protoplanet, effectively varying $`r_\mathrm{H}`$. Larger protoplanets have the effect of decreasing the overall density profile, whereas position controls the shape of the profile, particularly close to the protoplanet. Since the flow of gas within the Roche lobe cannot be adequately described hydrostatically, we exclude the region $`\sqrt{x^2+y^2}r_\mathrm{H}`$ from our analysis.

The normalization on equation (30) is somewhat arbitrary, but at the disk heights that are of interest to us, the overall normalization has small effect. We demonstrate this as follows: as $`zr_\mathrm{H}`$, equation (30) becomes

$$\rho _p(zr_\mathrm{H})\rho _0\mathrm{exp}\left(\frac{z^2}{2h^2}\frac{3r_\mathrm{H}^3}{h^2\sqrt{x^2+y^2}}\right),$$
(31)
that is, a Gaussian with respect to $`z`$ with some normalization that depends on $`r_\mathrm{H}`$ and $`\sqrt{x^2+y^2}`$. The height of the photosphere, $`H`$, is generally several times the pressure scale height. If we take $`H/h=5`$, then the change in density with vertical distance is

$$\frac{\mathrm{\Delta }\rho }{\rho }=\left|\frac{\mathrm{\Delta }z}{\rho }\frac{d\rho }{dz}(z=H)\right|=\frac{25\mathrm{\Delta }z}{H}.$$
(32)
So, a density change by a factor of $`1/2`$ corresponds to a shift in vertical height of 2%. In other words, the density gradient in the photosphere is so high that changes in the normalization of the density will not significantly affect our results.

### 3.2 Photosphere Height

We calculate the height of the perturbed photosphere from the perturbed density profile as the isodensity contour at the density of the unperturbed photosphere. That is, if $`H_0`$ is the unperturbed photosphere height and $`H=H_0(1ϵ)`$ is the new photosphere height, then we set

$$\rho _d(H_0)=\rho _p(H)$$

or

$$\frac{H_0^2}{2}=\frac{H^2}{2}+3r_\mathrm{H}^3\left(\frac{1}{\sqrt{x^2+y^2+H^2}}\frac{1}{\sqrt{x^2+y^2}}\right).$$
(33)
Setting $`\sqrt{x^2+y^2}=\gamma H_0`$, then to first order in $`ϵ`$

$$HH_0\left\{13\left(\frac{r_\mathrm{H}}{H_0}\right)^3\left[\frac{1}{\gamma }\frac{1}{\sqrt{\gamma ^2+1}}\right]\right\}$$
(34)
so that the depth of the perturbation scales as $`(r_\mathrm{H}/H_0)^3`$. For typical disk models, the ratio of the photosphere height to the thermal scale height, $`H_0/h`$, is between 3 and 5, depending on the distance from the star and is often taken to be a constant (Kenyon & Hartmann, 1987; Chiang & Goldreich, 1997; D’Alessio et al., 1998). Restating the scaling of the depth of the perturbation as $`(r_\mathrm{H}/h)^3`$, we see that this agrees with the estimate made in Sasselov & Lecar (2000).

The shape of this perturbed photosphere is shown in Figure 5 for varying values of $`r_\mathrm{H}`$. The perturbation is small for $`r>r_\mathrm{H}`$, and there is a singularity at $`r=0`$. As the Hill radius increases, the depth of the perturbation also increases. Note that if $`H_0/h=5`$, then $`r_\mathrm{H}/H_0=0.2`$ means that the Hill radius is comparable to the disk pressure scale height and the perturbation is no longer small. Protoplanets this massive are likely to open a gap in the disk rather than simply inducing a perturbation to the photosphere.

## 4 Disk Shear

The differential rotation of the disk causes the disk material near the protoplanet to move through the perturbation at a rate that depends on the velocity with respect to the protoplanet’s orbit. Material moving along a given streamline in the disk may pass through the perturbation too quickly to experience significant heating and cooling, thereby diminishing the effect of the perturbation on the temperature structure of the disk.

The velocity of the disk material with respect to the protoplanet is $`va(\mathrm{\Omega }\mathrm{\Omega }_p)`$ where $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }_p`$ are the orbital angular velocities of the disk material and the protoplanet, respectively, and $`a`$ is the orbital radius of the protoplanet. The orbital angular velocity of a gaseous disk is

$$r^2\mathrm{\Omega }^2=\frac{GM_{}}{r}+\frac{1}{\rho }\frac{dP}{dr}.$$
(35)
In thin disks, such as we are investigating, the pressure gradient term is small so the orbital velocity is nearly equal to the Keplerian velocity.

The heating/cooling rate of the disk material can be expressed as

$$C\frac{T}{t}=F\sigma T^4$$
(36)
where we take $`F=\sigma (T_v^4+T_r^4)`$ and $`C`$ is the specific heat per unit surface area of the disk. We shall adopt a specific heat of $`C=k\mathrm{\Sigma }/m`$ where $`k`$ is the Boltzmann constant, $`\mathrm{\Sigma }`$ is the total surface density, and $`m`$ is the mean molecular mass.

We shall assume that the streamlines of the disk material follow equation (35), and that they are not perturbed by the protoplanet’s graviational potential. This is an adequate approximation for disk material outside the protoplanet’s Hill radius. Along each streamline, we calculate the total radiative flux, $`F`$, at a given position, and using the velocity of the streamline with respect to the protoplanet along with heating/cooling rate, we can calculate the steady state temperature at each position in the disk.

## 5 Model Parameters

Dust properties can be parametrized by $`\alpha `$, $`\kappa _s`$, and $`q`$. Then, the temperature profile in a plane-parallel disk is determined completely by $`F_v`$, $`E_0`$, $`\mu _0`$, and $`\tau _d`$, where $`\tau _d`$ is measured normal to the surface of the photosphere. Stellar properties determine $`F_v`$ and $`E_0`$, and the disk profile determines $`\mu _0`$. If $`H_0r^\xi `$, then the incident angle is determined by

$$\mu _0=\frac{(\xi 1)H_0/a}{[1+(H_0/a)^2][1+\xi ^2(H_0/a)^2]}$$
(37)
A perturbation changes the local angle of incidence and optical depth. In particular, without plane parallel symmetry, there is no single optical depth that parametrizes the distance below the surface. Instead, we sum over different lines of sight, using the optical depth along each line of sight to determine the contribution to the disk flux as in equation (23). The local angle of incidence depends only on the geometry of the perturbed surface, while the optical depth also depends on the density structure along the line of sight. The shape of the perturbed surface depends on $`r_\mathrm{H}`$ and $`H_0`$, and the density structure depends on the disk pressure scale height, $`h`$.

### 5.1 Fiducial Model

Motivated by observations of T Tauri stars, we will assume for our fiducial model that the central star has mass $`M_{}=0.5\mathrm{M}_{\mathrm{}}`$, radius $`R_{}=2R_{}`$, and effective temperature $`T_{}=4000`$ K, and that the protoplanet is at $`a=1\mathrm{AU}`$ from the star. We assume an accretion rate of $`\dot{M}_a=10^8\mathrm{M}_{\mathrm{}}\text{yr}^1`$, so that $`T_v=72`$ K at $`\tau _d=2/3`$ at 1 AU. With these parameters, heating in the photosphere will be dominated by stellar irradiation rather than viscous heating. Near the midplane, however, viscous heating will dominate. Since we only consider the photosphere in our model, the calculation of the midplane temperature is outside the scope of this paper.

We will assume that the fraction of scattered radiation is $`\alpha =0.28`$, with opacity in the optical of $`\kappa _s=400\text{cm}^2\text{g}^1`$, and ratio of opacities $`q=20`$.

We shall assume that $`H_0=0.16\mathrm{AU}(r/1\mathrm{AU})^{9/7}`$, consistent with e.g. D’Alessio et al. (1999a). In calculating the detailed density structure in the photosphere, we shall assume that $`h=3.4\times 10^2\mathrm{AU}`$ so that $`H_0/h=5`$. A protoplanet with mass $`m_p=1\mathrm{M}_{}`$ ($`10\mathrm{M}_{}`$) will have a Hill radius of $`1.26\times 10^2\mathrm{AU}`$ ($`2.72\times 10^2\mathrm{AU}`$). As shown by Bate et al. (2003), planets with masses $`30\mathrm{M}_{}`$ are not massive enough to open a gap, and our perturbative treatment is valid. We have calculated models with $`\mathrm{\hspace{0.25em}1}\mathrm{M}_{}`$, $`10\mathrm{M}_{}`$, and $`20\mathrm{M}_{}`$ planets, and present the results below.

### 5.2 Numerical details

We define our coordinate axes so that $`x`$ is aligned with the radial direction, $`y`$ is aligned with $`\varphi `$ (i.e. the direction of planet’s orbit), and $`z`$ is perpendicular to the orbital plane. The coordinate origin is set at the surface of the disk above the planet position, so that the planet’s coordinates are $`(0,0,H_0)`$.

We calculate the shape of the disk surface in the $`x`$ and $`y`$ directions over a range $`10r_\mathrm{H}x<10r_\mathrm{H}`$ and $`10r_\mathrm{H}y<10r_\mathrm{H}`$, and use this surface to numerically integrate temperatures within the photosphere, to a depth of $`z(\tau _d=2/3)=0.011\mathrm{AU}`$. The temperature is calculated at each point within a $`64\times 64\times 64`$ grid.

When summing over surface elements to find $`B_{\mathrm{tot}}`$ as in equation (23), we can approximate $`\delta \mathrm{\Omega }\nu \delta A/\mathrm{}^2`$ for $`\delta A\mathrm{}^2`$, where $`\mathrm{}`$ is the distance between the point $`P`$ and $`\delta A`$. However, for points close to the surface, this approximation breaks down. In addition, $`\nu `$ may have large excursions over $`\delta A`$. For this reason, we assume that each $`\delta A`$ is approximately planar, so that $`B(\tau _d,\mu )`$ does not vary, but we calculate $`_{\delta \mathrm{\Omega }}\nu 𝑑\mathrm{\Omega }`$ analytically from the limits of $`\delta A`$. In other words, we calculate

$$B_{\mathrm{tot}}=\frac{1}{\pi }B(\tau _d,\mu )\underset{\delta \mathrm{\Omega }}{}\nu 𝑑\mathrm{\Omega }$$
(38)
and solve for the temperature.

## 6 Results

Figures 6, 7, and 8 show the spatial distribution of temperatures in the photosphere of the fiducial model for $`1\mathrm{M}_{}`$ $`10\mathrm{M}_{}`$, and $`20\mathrm{M}_{}`$ planets, respectively. The colors represent temperatures in Kelvin, as indicated by the colorbars. The horizontal axis indicates increasing radial distance from the star so that the photosphere is illuminated from the left. The cross-section is taken at the $`\tau _d=2/3`$ surface. The bottommost row and leftmost column show the unperturbed temperature of $`133.6\mathrm{K}`$ for comparison. This temperature is also indicated in the colorbar by a dotted white line.

The direction of the planet’s orbit is upward on the plots, so the interior disk material moves faster, and the exterior disk material moves slower than the planet. This causes the heated and cooled regions to shear out. The positions of minimum and maximum temperatures are offset from the planet’s position since it takes some time for the heated/cooled material to come back to the equilibrium temperature.

For a $`1\mathrm{M}_{}`$ protoplanet, the shadowing and illumination at the surface results in a slight temperature change, with the area to the right of the planet reaching a maximum temperature of $`135.5\mathrm{K}`$, and the area to the left cooled to a minimun of $`130.7\mathrm{K}`$ – a change of $`2.2\%`$ to $`+1.4\%`$. Thus, we see that an Earth-mass planet is too small significantly affect the interior temperature of the disk.

A larger planet will create a larger perturbation, which will lead to a greater temperature change. Figure 7 and Figure 8 show the temperature distribution in the photosphere for a $`10\mathrm{M}_{}`$ planet and $`20\mathrm{M}_{}`$ planet, respectively. The plots are scaled to the Hill radius, which goes as $`m_p^{1/3}`$, so the physical range in Figures 7 and 8 is actually larger than that in Figure 6. The high temperatures close to the Hill radius are due to the low density in this region, which gives a low optical depth to the illuminated surface. However, since this region is not in hydrostatic equilibrium, we exclude it from our analysis. Outside the Hill radius, the temperatures at $`\tau _d=2/3`$ for a $`10\mathrm{M}_{}`$ protoplanet range from $`123.7\mathrm{K}`$ to $`142.3\mathrm{K}`$, or by $`7.4\%`$ to $`+6.5\%`$. For a $`20\mathrm{M}_{}`$ protoplanet, which is close to but not at the limit where a gap in the disk can form, the temperature range increase slightly, with a minimum of $`122.0\mathrm{K}`$ and maximum of $`146.6\mathrm{K}`$, or by $`9.7\%`$ to $`+8.7\%`$. Thus, a protoplanet with mass 10 $`\mathrm{M}_{}`$ or larger may have a significant effect on the local temperature structure of a protoplanetary disk.

## 7 Discussion

As planet-forming disks are likely to be passive (Calvet et al., 2002), we have explored stellar irradiation effects due to a disk surface perturbation induced by an embedded protoplanet. We find that even a relatively small protoplanet (in the range from 1 to 20 Earth mass) produces enough of a perturbation to cause shadowing. In particular, we find that in a standard passive disk accreting at $`10^8\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$, the temperature variation induced by a $`1\mathrm{M}_{}`$ planet at 1 AU is $`5\mathrm{K}`$, or only about $`4\%`$. However, we find that even though a protoplanet may not be massive enough to open a gap in a disk, it may still be able to significantly change the temperature structure in a the disk material near it. A $`10\mathrm{M}_{}`$ ($`20\mathrm{M}_{}`$) planet can induce a temperature variation of $`19\mathrm{K}`$ ($`25\mathrm{K}`$), or $`14\%`$ ($`18\%`$).

This results in a local temperature perturbation which is large enough to affect the local properties of the disk and the accretion rate, and possibly the migration rate of the protoplanet. Although the perturbation in temperature may be too subtle to directly observe, it can have serious consequences for planet growth and migration.

Ward (1997) has shown that the rate of Type I migration is strongly dependent upon the temperature gradient of the disk. If the temperature decreases with distance from the star, then the total net torques result in inward migration of the secondary. However, if $`k=d\mathrm{ln}T/d\mathrm{ln}r1`$, then the opposite is true and the secondary will migrate outward. The temperature perturbation near a protoplanet locally decreases the value of $`k`$, and may result in halting or even reversing Type I migration. This may resolve the conundrum that Type I migration timescales are typically much less than observed disk lifetimes of $`10^7`$ years (Hartmann et al., 1998). Yet to date, no simulation of planet-disk interactions accounts for local temperature pertubations as described here.

We have not yet considered the behavior of the dust (as a function of size) in the perturbed region. We expect the dust to remain coupled to the gas and compress its scale height accordingly. In fact crude estimates of dust sedimentation in the vicinity of the planet confirm that. However, a detailed estimate of the sedimentation timescale is necessary in order to evaluate the amount of residual dust, e.g. swept by stellar radiation from the inner edge, and its total opacity. These are issues we plan to address in future work.

In a forthcoming paper, we will undertake a parameter study of our model for radiative transfer on a perturbed disk. By varying attributes of the disk, star, and protoplanet, we will investigate the possibility that a relatively small protoplanet which is not large enough to open a gap in a disk might still be able to affect disk structure.

We thank Nuria Calvet for numerous helpful comments in the preparation of this paper.


text:
# Oscillations do not distinguish between massive and tachyonic neutrinos

## Abstract

It is shown that the hypothesis of tachyonic neutrinos leads to the same oscillations effect as if they were usual massive particles. Therefore, the experimental evidence of neutrino oscillations does not distinguish between massive and tachyonic neutrinos.

neutrino oscillations, tachyonic neutrinos

In the last two decades one of the most interesting problems in the particle physics is the issue of neutrino masses Gelmini and Roulet (1995). The flavor oscillations Altmann et al. (2001); Beuthe (2003) are the indirect evidence for non-zero neutrino masses. On the other hand, the experiments devoted to direct measurement of the neutrino mass permanently yield the negative mass squared for the electron neutrino Lobashev et al. (1999); Weinheimer et al. (1999) and the muon neutrino Assamagan et al. (1996) (for a review, see also Hagiwara et al. (2002)). In particular, the most sensitive neutrino mass measurement involving electron antineutrino, is based on fitting the shape of beta spectrum. An as yet not understood event excess near the spectrum endpoint can be explained Ciborowski and Rembieliński (1997, 1999) on the ground of the hypothesis of the tachyonic neutrinos Chodos et al. (1985). The above hypothesis can be also helpful in explanation of some effects in cosmic ray spectrum Ehrlich (1999a, b, 2000). From the theoretical point of view, this hypothesis can be formulated consistently only in the framework of the absolute synchronization scheme Rembieliński (1997); Caban and Rembieliński (1999) which needs the notion of preferred frame.

The aim of this report is to compare the flavor change effect for the massive and tachyonic neutrinos. Calculations will be performed in the preferred frame, in which the metric tensor has the standard Minkowskian form.

Let us consider three flavors of neutrinos. We shall denote the neutrino flavor states by

$$|\nu _e=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)|\stackrel{}{p},|\nu _\mu =\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)|\stackrel{}{p},|\nu _\tau =\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)|\stackrel{}{p}.$$
(1)
The $`|\stackrel{}{p}`$ denotes the vector from the representation space of the Poincaré group, corresponding to the momentum $`\stackrel{}{p}`$.

In the basis (1) the free Hamiltonian describing neutrinos can be written as

$$H=U^{}H_0U,$$
(2)
where $`U`$ is a $`3\times 3`$ unitary mixing matrix, acting in each subspace of the Hilbert space determined by momentum $`\stackrel{}{p}`$, and the Hamiltonian $`H_0`$ in the mass eigenstates basis is of the form

$$H_0=\left(\begin{array}{ccc}\sqrt{\stackrel{}{p}^2\pm m_1^2}& 0& 0\\ 0& \sqrt{\stackrel{}{p}^2\pm m_2^2}& 0\\ 0& 0& \sqrt{\stackrel{}{p}^2\pm m_3^2}\end{array}\right).$$
(3)
Here $`m_j^2`$ ($`j=1,2,3`$) are the absolute values of the squares of neutrino four-momenta and the choice of the sign corresponds to the massive ($`+`$) and tachyonic ($``$) case, moreover, $`|\stackrel{}{p}|>\mathrm{max}_jm_j`$ in the tachyonic case.

Following the standard procedure, we consider the time evolution of the neutrino state. Assume that the initial state $`|\nu _i`$ ($`i=e,\mu ,\tau `$) is an eigenstate of the fractional lepton number operator. Thus, after time $`t`$ the probability that we get the neutrino $`\nu _k`$ is given by

$$P_{\nu _i\nu _k}(t)=\left|\nu _k|U^{}e^{itH_0}U|\nu _i\right|^2.$$
(4)
Now, taking (4) in the limit $`|\stackrel{}{p}|m_j`$, we can expand the matrix elements of the Hamiltonian as follows

$$\sqrt{\stackrel{}{p}^2\pm m_i^2}p\pm \frac{m_i^2}{2p}.$$
(5)
Using (5) we obtain from (4)

$$P_{\nu _i\nu _k}(t)=\left|\underset{j=1}{\overset{3}{}}u_{jk}^{}u_{ji}e^{\frac{it}{2p}m_j^2}\right|^2.$$
(6)
On introducing the polar representation such that

$`u_{2k}^{}u_{2i}u_{3k}u_{3i}^{}`$ $`=r_1(i,k)e^{i\theta _1(i,k)},`$ (7a)
$`u_{3k}^{}u_{3i}u_{1k}u_{1i}^{}`$ $`=r_2(i,k)e^{i\theta _2(i,k)},`$ (7b)
$`u_{1k}^{}u_{1i}u_{2k}u_{2i}^{}`$ $`=r_3(i,k)e^{i\theta _3(i,k)},`$ (7c)
and

$$\omega _1=\frac{m_2^2m_3^2}{2p},\omega _2=\frac{m_3^2m_1^2}{2p},\omega _3=\frac{m_1^2m_2^2}{2p},$$
(8)
we can write (6) in the following form:

$$P_{\nu _i\nu _k}(t)=\underset{j=1}{\overset{3}{}}|u_{jk}|^2|u_{ji}|^2+2\underset{j=1}{\overset{3}{}}r_j(i,k)\mathrm{cos}[\omega _jt\pm \theta _j(i,k)],$$
(9)
where, as before, the signs $`+`$ and $``$ correspond to massive and tachyonic neutrinos, respectively. In particular, if $`k=i`$ then $`\theta _j(i,i)=0`$ ($`j=1,2,3`$) and

$$P_{\nu _i\nu _i}(t)=\underset{j=1}{\overset{3}{}}|u_{ji}|^4+2\underset{j=1}{\overset{3}{}}r_j(i,i)\mathrm{cos}(\omega _jt),$$
(10)
so this probability is the same for massive and tachyonic neutrinos.

We have shown that the only difference between the neutrino oscillations in the massive and tachyonic case lies in the initial phase of oscillations $`\theta _j(i,k)`$. Moreover, the initial phases for tachyonic case are obtained by taking the complex conjugations of the elements of mixing matrix for the massive case (see Eqs. (7)). In the oscillation experiments we cannot decide wheter we should take the mixing matrix $`U`$ or its complex conjugation $`U^{}`$ without taking into account the nature of the neutrinos. More precisely, we can explain the oscillations for massive or tachyonic neutrinos simply by the different choice of the mixing matrix ($`U`$ or $`U^{}`$). Therefore, the experimental evidence of neutrino oscillations does not distinguish between massive and tachyonic neutrinos which, together with other experiments Lobashev et al. (1999); Weinheimer et al. (1999); Assamagan et al. (1996) leave the question of tachyonic character of the neutrinos still open.