text:
# On the Location and Composition of the Dust in the MCG–6-30-15 Warm Absorber

## 1 Introduction

Soft X-ray absorption by photoionized gas was first used by Halpern (1984) and Pan et al. (1990) to explain the unusual spectrum of the quasar MR 2251–178. This “warm absorber” was characterized by ROSAT spectra of Seyfert 1s which showed absorption edges due to O vii and O viii (Nandra & Pounds, 1992; Nandra et al., 1993; Fiore et al., 1993; Turner et al., 1993). Warm absorber studies matured with the launch of the more sensitive ASCA observatory, which allowed more detailed spectral modeling (e.g., Fabian et al., 1994; George et al., 1995; Reynolds & Fabian, 1995). Analysis of a large sample of Seyfert 1s observed by ASCA showed that approximately 50% exhibit absorption features from photoionized gas in their soft X-ray spectra (Reynolds, 1997; George et al., 1998). Currently, the study of the warm absorber is being revolutionized with observations from the dispersion gratings onboard Chandra and XMM-Newton which have been able to resolve individual absorption lines from a myriad of metals and charge states in the warm gas (e.g., Kaastra et al., 2000; Kaspi et al., 2001; Collinge et al., 2001; Lee et al., 2001; Kaspi et al., 2002). Photoionization modeling has then shown that, in many cases, more than one ionization parameter111There are some objects which can still be modeled adequately with only one ionization parameter, e.g., Mrk 509 (Yaqoob et al., 2003). is needed to describe the observed spectrum (e.g, Morales et al., 2000; Kaspi et al., 2001; Kaastra et al., 2002). Also, the positions of the lines point toward outflow velocities on the order of a few hundred to a few thousand km s-1. Thus, this warm absorbing gas seems to be in the form of an outflowing wind (cf., Bottorff et al., 2000; Elvis, 2000). However, there remains a major uncertainty in the location of the gas, with many models considering the broad-line region (e.g., Reynolds & Fabian, 1995; Netzer, 1996; George et al., 1998) or the putative obscuring torus (e.g., Krolik & Kriss, 1995, 2001) as the most likely origin for the warm absorber.

A possible constraint on the location of at least some of this warm gas may be provided if it contains any dust, which sublimates at the radius of the broad-line region (Barvainis, 1987) for a typical active galactic nucleus (AGN). Dusty warm absorbers (hereafter, DWA) were first considered for the quasar IRAS 13349+2438 (Brandt et al., 1996; Siebert, Komossa & Brinkmann, 1999) and the Seyfert 1 galaxy MCG–6-30-15 (Reynolds et al., 1997). In both of these AGN, the column of neutral H inferred from the reddening is significantly larger than that inferred from the neutral absorption in the soft X-ray band, but is of the same order as the column of ionized gas inferred from the warm absorber (Reynolds et al., 1997). This suggests that the dust is similar to Galactic dust and resides within the warm ionized gas, and may significantly affect the observed soft X-ray spectrum (Komossa & Fink (1997) and references therein; Komossa & Bade 1998). Spectroscopic evidence for dust in X-ray warm absorbers has now been found by Lee et al. (2001) in the Chandra gratings observation of MCG–6-30-15 (this seems to have been confirmed by the very recent XMM-Newton data of Turner et al. 2003). This high resolution spectrum exhibited a sharp drop222An alternative explanation for this drop (based on XMM-Newton data) is that it is the blue edge of a relativistically broadened O viii Ly$`\alpha `$ line (Branduardi-Raymont et al., 2001; Sako et al., 2003). However, this interpretation has been challenged from both observational (Lee et al., 2001) and theoretical (Ballantyne et al., 2002) points of view. at $``$0.7 keV consistent with the L3 absorption edge from neutral Fe. A similar feature was also found in a Chandra spectrum of the Galactic X-ray binary Cyg X-1 (Schulz et al., 2002). The column implied by the depth of the Fe edge in MCG–6-30-15 is of the right order to explain the observed reddening ($`A_V\tau _V3E(BV)=1.8`$; Reynolds et al. 1997), assuming Galactic-type dust located within the warm absorber.

Rather than place the DWA near the central engine of the AGN, Kraemer et al. (2000) (see also Crenshaw & Kraemer 2001) argued for the existence of a “lukewarm absorber” outside the narrow-line region. This gas would have sufficient column to explain the observed reddening, and has been ionized to the point where hydrogen is fully stripped, but the metals would only be moderately ionized and would exhibit strong UV absorption lines rather than O vii or O viii edges. Thus, this model requires an inner warm absorber to account for the highly ionized oxygen features (Kraemer et al., 2000). The lukewarm absorber has been shown to be consistent with the X-ray (Kraemer et al., 2000) and UV (Crenshaw et al., 2001) data of NGC 3227, as well as the UV spectrum of Ark 564 (Crenshaw et al., 2002).

In the case of MCG–6-30-15 ($`z=0.008`$, $`L_{210\mathrm{keV}}10^{43}`$ erg s-1; see Table 1 for a summary of the absorbing columns), it was clear from early ASCA variability studies that a multi-zone warm absorber was needed (e.g., Otani et al., 1996; Morales et al., 2000).

In particular, the O viii edge was found to anticorrelate with the source luminosity while the O vii edge seemed to remain constant (Otani et al., 1996; Orr et al., 1997). The Chandra observations found a strong O vii edge ($`\tau _{\mathrm{OVII}}0.7`$) and a series of O vii absorption lines within 200 km s-1 of the systemic velocity of MCG–6-30-15 (Lee et al., 2001). An updated analysis presented by Lee et al. (2002) showed that many of the other low-ionization absorption lines also have very small outflow velocities. These data point to a distribution (in space/velocity) of ionized absorbers, consistent with the earlier ASCA results.

A Hubble Space Telescope (HST) image of MCG–6-30-15 shows a distinct dust lane that cuts across the southern part of the galactic disk (Fig. 1; Malkan et al. 1998; see also Ferruit et al. 2000).

This material is likely intrinsic to the galaxy and not a foreground effect such as Galactic cirrus (e.g., M81; Abolins & Rice 1987). Therefore, in Sect. 2, we follow the ideas of Kraemer et al. (2000), but argue that at least part of the low ionization warm absorber and all of the observed reddening in MCG–6-30-15 originates within this dust lane. It is not our intention to construct a detailed spectral model of the DWA in MCG–6-30-15; rather we use this system as an example of how material in the dusty interstellar medium (ISM) of the host galaxy can have a significant observational impact on the study of absorption in AGN spectra.

The composition of the dust in the MCG–6-30-15 DWA is also a puzzle. For Cyg X-1, where the high signal-to-noise allowed a comparison to laboratory data, Schulz et al. (2002) claimed that pure metallic Fe was the best fit to the edge structure. This is surprising, since the expected microwave emission (at $``$ 90 GHz) from such grains is not observed in the Galaxy (Draine & Lazarian, 1999). Lee et al. (2001) argued that the Cyg X-1 iron edge was similar to the one observed in MCG–6-30-15, and thus pure iron grains may be present in the DWA. In Sect. 3 we discuss the constraints on a pure iron component in astronomical dust, and consider if it is a plausible composition for the dust in the MCG–6-30-15 DWA.

## 2 The Location of the DWA

The dust lane seen in Fig. 1 passes just south of the nucleus, so, judging from the Galactic distribution (Dickey & Lockman, 1990), our line of sight likely passes through a column of a few times 1021 cm-2. This is of the same order as that needed to provide the observed reddening in MCG–6-30-15 assuming Galactic-type grains ($`N_\mathrm{H}4\times 10^{21}`$ cm-2; Reynolds et al. 1997). Further evidence for a low density origin of the DWA is the recombination timescale argument of Otani et al. (1996). They noted that the constancy of the O vii edge in the long 1994 ASCA observation of MCG–6-30-15 implies $`n_\mathrm{H}<2\times 10^5\mathrm{cm}^3`$ in the outer absorber. ISM gas with a fixed hydrogen number density of $`n_\mathrm{H}=10`$ cm-3 will be used as the canonical absorber for the remainder of the paper.

A simple photoionization argument can be made to constrain the distance of the O vii absorber. We define the following ionization parameter

$$U_\mathrm{H}\frac{L_>}{4\pi r^2cn_\mathrm{H}h\nu },$$
(1)
where $`L_>`$ is the AGN luminosity beyond the Lyman edge, $`n_\mathrm{H}`$ is the total hydrogen number density, and $`h\nu `$ is an average ionizing photon energy. Photoionization equilibrium gives $`\xi _\mathrm{H}=n_{\mathrm{H}^0}/n_\mathrm{H}=(\alpha _\mathrm{H}(T)/a_\mathrm{H}c)U_\mathrm{H}^1`$, where $`n_{\mathrm{H}^0}`$ is the number density of neutral H, $`\alpha _\mathrm{H}`$ is the recombination coefficient of H at temperature $`T`$ and $`a_\mathrm{H}`$ is the photoionization cross-section for H. This expression is valid only when $`\xi 1`$. Generalizing for a hydrogenic ion of an element with atomic number $`Z`$ gives

$$\xi _Z=\frac{\alpha _\mathrm{H}(T)}{a_\mathrm{H}c}U_\mathrm{H}^1Z^{4+2\beta }=(1.4\times 10^6)U_\mathrm{H}^1Z^{4+2\beta },$$
(2)
where we have assumed $`L_\nu \nu ^\beta `$ for the ionizing luminosity ($`\beta 1.6`$ for radio-quiet quasars; Telfer et al. 2002), and used $`a_\mathrm{H}=6\times 10^{18}`$ cm2 and $`\alpha _\mathrm{H}`$(20,000 K)=2.5$`\times 10^{13}`$ cm3 s-1 (Osterbrock, 1989). With a typical Seyfert 1 Lyman luminosity of $`L_>=10^{44}`$ erg s-1, and assuming $`h\nu =45`$ eV, the ionization parameter of a $`n_\mathrm{H}=10`$ cm-3 cloud at 1 kpc is 0.04. Equation 2 then gives $`\xi _\mathrm{H}10^5`$ and $`\xi _\mathrm{O}100`$, implying that although hydrogen is still entirely ionized at this distance, oxygen has started to recombine. If we now consider gas only 10 pc from the AGN, the values of $`\xi `$ will drop by 104, and O is predominantly fully-stripped. Therefore, the O vii absorber must lie between 10 and 1000 pc. This calculation assumes photoionization equilibrium and is therefore valid only on timescales longer than the recombination time, $`t_{\mathrm{rec}}1/n_e\alpha _\mathrm{H}(T)10^4`$ years, where $`n_e`$ is the free electron density. The radial range derived above will be reasonable as long as the luminosity of the central engine did not change by factors greater than $`10`$ over $`t_{\mathrm{rec}}`$. To more accurately determine the properties of an ISM warm absorber, numerical photoionization models must be employed.

We used Cloudy 96Beta4 (Ferland, 2002) to predict the ionization structure of a $`n_\mathrm{H}=10`$ cm-3 cloud of ISM gas at various distances from an AGN. The neutral hydrogen column density was fixed at $`N_\mathrm{H}=4\times 10^{21}`$ cm-2, which, for a Galactic gas-to-dust ratio, is the minimum needed to provide the reddening toward MCG–6-30-15 ($`N_\mathrm{H}=4`$–7$`\times 10^{21}`$ cm-2; Reynolds et al. 1997). The dust (silicate plus graphite) and gas-phase metal abundances were fixed at the ISM values described by Ferland (2002). The cloud was illuminated with a “standard” AGN continuum (see pg. 34 in Ferland 2002) with $`\alpha _{\mathrm{ox}}=1.4`$ and a X-ray power-law photon-index $`\mathrm{\Gamma }=2`$. The “big blue bump” was characterized by a temperature of $`1.4\times 10^5`$ K (the maximum Shakura & Sunyaev (1973) accretion disc temperature for a 107 M black hole accreting at 0.1 of its Eddington rate) and a UV slope of $`\alpha (\mathrm{UV})=0.5`$. This spectral energy distribution (SED) is not intended to be a realistic model of the MCG–6-30-15 continuum (which is unknown because of the large reddening), but rather representative of a generic AGN. The normalization of the SED was set by defining the 2–10 keV luminosity to be 1043 erg s-1, typical of many Seyfert 1s (Reynolds, 1997).

Models were calculated with the inner edge of the gas cloud at various distances $`r`$ from the continuum source, and the computed O vii column density was compared with the result from the Chandra observation of MCG–6-30-15, $`N_{\mathrm{OVII}}2.5\times 10^{18}`$ cm-2 (Lee et al., 2001). The maximum O vii column found in the Cloudy runs was $`6.1\times 10^{17}`$ cm-2 for a cloud distance of 150–175 pc. Models with the gas closer in were too ionized, and if the cloud was further out, it was not ionized enough. Of course, lowering the density by a factor $`f`$ would allow a distance $`\sqrt{f}`$ larger. However, if the distance exceeds $``$200–250 pc, then, even with $`n_\mathrm{H}=10\mathrm{cm}^3`$, the physical length of the column exceeds the distance to the AGN. A similar calculation with a power-law SED of energy index $`1.17`$ between 0.013 and 100 keV resulted in $`N_{\mathrm{OVII}}=6.2\times 10^{17}`$ cm-2 at $`r=150`$ pc. While these models cannot account for the entire O vii column in MCG–6-30-15 (it is likely that some fraction of the O vii edge originates within the inner warm absorber), we have shown that this dusty ISM cloud will have an detectable impact on the observed spectrum (see also Komossa & Bade 1998).

To illustrate the extent of this impact, the incident and transmitted continua for the AGN models with $`r=100`$, $`175`$, and $`250`$ pc are shown in Fig. 2.

These runs have $`U_\mathrm{H}=3.7`$, $`1.2`$ and $`0.6`$, respectively. The transmitted spectra show significant warm absorption in the soft X-ray band with strong edges at 0.28, 0.53 and 0.74 keV. These first two are due to C i and O i in the grains (we include both graphite and silicate dust), while the latter is the O vii edge. Of the three models shown here, the depth of the O vii edge is greatest when $`r=175`$ pc, with $`\tau _{\mathrm{OVII}}=0.15`$. The Fe i edge from the dust is noticeable as a small notch just redward of the O vii edge, as was observed by Lee et al. (2001). The majority of the opacity below 0.3 keV is due to absorption by He ii, as was found in the lukewarm absorber of Kraemer et al. (2000). Therefore, such absorption will need to be included in absorption models for sources such as MCG–6-30-15 where distant material is contributing to the warm absorber. Indeed, in their multi-zone DWA fit to the new XMM-Newton data, Turner et al. (2003) require a low-ionization component to account for additional soft X-ray absorption that is inconsistent with cold absorption. It is quite likely that this excess attenuation is caused by He ii absorption.

The fact that there is a significant column of He ii ($`5.9\times 10^{18}`$ cm-2) in the models shows that the gas is very close to recombining to a more neutral configuration. This is in large part due to the presence of dust which, for a Galactic size distribution, dominates the continuum opacity in the extreme UV (EUV) and soft X-ray bands (see Fig. 3). Indeed, our Cloudy models show that neutral oxygen begins to be present in the gas when $`r=175`$ pc. Thus, our results seem to be fine-tuned in the sense that the greatest O vii column is found when the gas is on the verge of recombining. This argues that it may be difficult, in general, for ISM gas to contribute to the warm absorber.

In this section we have argued that ISM material associated with dust lanes in the host galaxy of an AGN is a natural contributor to the observed DWAs. Cloudy models suggest that this gas, when illuminated by a toy AGN SED, can produce significant warm absorption (including a significant O vii edge) when placed $``$150 pc from the AGN. This may be particularly relevant to MCG–6-30-15, where a dust lane does pass just below the nucleus. Using Fig. 1 and the plate scale from Malkan et al. (1998), we can estimate the inner edge of the dust lane from the HST image. Interestingly, this is also of the order of 150–450 pc, assuming $`H_0=50`$ km s-1 Mpc-1.

It may be difficult applying this idea to other sources which harbor DWAs (e.g., IRAS 13349+2438). The dust lane in MCG–6-30-15 lies nearly along our line of sight and may be a relatively unusual alignment. Nevertheless, the strong ionizing power of an AGN will have a great impact on its local ISM. The observational consequences of these effects will depend both on the viewing angle into the source (Crenshaw & Kraemer, 2001) and on the SED of the AGN. Indeed, if our line-of-sight into MCG–6-30-15 were slightly different and passed through the dust lane, rather than just above it, the obscuration could be high enough that the AGN would appear as a Seyfert 2 (see, e.g., Matt, 2000).

## 3 Dust Properties

### 3.1 Composition

As the EUV/soft X-ray opacity in the DWA is dominated by dust, the absorption properties depend on the grain composition. The top panel of Fig. 3 shows the transmitted spectrum for the $`r=150`$ pc model with Galactic silicate and graphite dust abundances (Mathis, Rumpl & Nordsieck, 1977).

In the bottom panel, the spectrum is displayed for the case where the elements in the dust have been returned to the gas. We assumed that the total (gas + dust) abundances are equal to solar (Grevesse & Sauval, 1998) except for O, for which we took O/H=$`4.45\times 10^4`$, a value measured from local young F and G stars (Sofia & Meyer, 2001). This oxygen abundance has the advantage that it does not require depletions greater than what can be accounted for by dust. The DWA model exhibits a C i edge at 0.28 keV of depth $`\tau =0.62`$ from the graphite dust, an O i edge at 0.53 keV of depth $`\tau =0.15`$, and an Fe i edge at 0.72 keV of depth $`\tau =0.09`$, both from the silicate dust. The dust-free model also shows absorption edges, the largest from O vii at 0.74 keV of depth $`\tau =0.26`$.

Unfortunately, the High Energy Transmission Gratings (HETG) Chandra spectra only extend to an energy of 0.45 keV, so we have no X-ray evidence for or against the presence of graphite; however the Low Energy Transmission Gratings (LETG) can reach these low energies although the instrumental C edge will need to be well calibrated to measure the intrinsic edge due to graphite. The Reflection Grating Spectrometer (RGS) onboard XMM-Newton also does not go below 0.3 keV, although the low-resolution European Photon Imaging Cameras (EPIC) are able to do so.

As mentioned in § 1, the strength of the measured Fe i edge from the Chandra observation of MCG–6-30-15 implies $`N_{\mathrm{Fe}\mathrm{I}}4\times 10^{17}\mathrm{cm}^2`$, consistent with the observed reddening for Galactic-type dust. If the Fe really is incorporated into silicate grains, then we expect $`N_{\mathrm{O}\mathrm{I}}>N_{\mathrm{Fe}\mathrm{I}}`$, since a silicate structural unit contains more O atoms than Fe atoms. For example, the Cloudy models discussed above assume a MgFeSiO4 stoichiometry and predict an O i edge that is roughly the same depth as the Fe i edge. Lee et al. (2002) claimed that silicates with an O:Fe number ratio of 2:1 rather than 4:1 were consistent with the Chandra data, but also mentioned that the data were limited by statistics at the position of the O i edge. Therefore, if the silicate dust had an O:Fe ratio of 4:1, it could have been detected by Chandra. However, Sako et al. (2003) quote an upper limit of $`N_{\mathrm{O}\mathrm{I}}=10^{16}`$ cm-2 from their XMM-Newton RGS observation of MCG–6-30-15, which has higher signal-to-noise, but lower resolution, than the Chandra spectrum. These authors also find $`N_{\mathrm{Fe}\mathrm{I}}=7\times 10^{16}`$ cm-2 for the column of neutral iron, which implies a $`V`$-band optical depth of 0.07 (assuming a silicate grain radius of 0.1 $`\mu `$m). Therefore, this column could not account for the observed reddening of MCG–6-30-15, unless it is provided predominantly by iron-poor grains.

Further constraints on the mineralogical structure of the dust can, in principle, be provided by the detailed structure of the Fe i edge. Lee et al. (2001) found that the edge structure in MCG–6-30-15 matches very well that found in Cyg X-1. Schulz et al. (2002) showed that the Cyg X-1 edge could be fairly accurately reproduced if the Fe resides in pure Fe grains (taking experimental data from Kortright & Kim 2000), but not if the Fe is in oxides (taking data from Crocombette et al. 1995). A measurement of the oxygen column in the dust could not be made for Cyg X-1, as the large O i edge ($`\tau >2`$) was dominated by neutral oxygen in the ISM along the line-of-sight. Therefore, Schulz et al. (2002) concluded that pure iron grains were most consistent with the data, and Lee et al. (2001) present this as a possibility for the dust in MCG–6-30-15. However, for Cyg X-1 metallic Fe dust is probably ruled out since it would produce excessive thermal magnetic dipole emission at 90 GHz (Draine & Lazarian, 1999). Draine & Lazarian suggest that at most $``$5% of the Galactic interstellar Fe can be in metallic iron grains. Thus, it seems unlikely that there is a significant column of metallic iron grains toward Cyg X-1. Although we have not performed a detailed spectral fitting, it appears that the observed edge structure in Cyg X-1 could be due to Fe in silicates, based on the spectrum of an olivine sample from Gloter et al. (2000).

The above argument cannot be directly applied to MCG–6-30-15 as there are no microwave observations of this galaxy, so the key to uncovering the dust composition in the DWA is an accurate measurement of other edges such as O i, C i, or Si i at 1.85 keV ($`\tau _{\mathrm{Si}\mathrm{I}}<0.01`$ in the Cloudy models). The Chandra data may imply wüstite (FeO) or hematite ($`\alpha `$Fe2O3) dust, but the XMM-Newton limit by Sako et al. (2003) points to little or no oxygen in the grains.

Perhaps the lack of O i can be explained by selective destruction of O-bearing grains (e.g., silicates). Depletion patterns in the Galaxy suggest that Fe and Si are incorporated into different dust populations, with the Fe-bearing component significantly longer-lived than the Si-bearing component (e.g., Tielens, 1998). If the dust in the DWA in MCG–6-30-15 were subjected to repeated episodes of destruction (in shocks) and re-growth, then the Si-bearing dust component could be largely removed while the Fe-bearing component remains largely intact. Such a process would lead to a reduced $`E(BV)/N_\mathrm{H}`$. However, the metallicity at $`r150`$ pc from the nucleus of MCG–6-30-15 is probably higher than at the location of the Sun in the Galaxy, so that $`E(BV)/N_\mathrm{H}`$ in the DWA in this scenario could be comparable to its value in the local ISM of the Galaxy. One potential difficulty for this scenario (and for silicate dust models in general) is that the destruction timescale inferred by Tielens (1998) for the Si-bearing dust is substantially shorter than predicted for silicate grains (Jones et al., 1996). Counter-intuitively, Jones et al. also find the lifetime of pure Fe dust to be shorter than the lifetimes of graphite or silicate dust, because Fe grains are accelerated to higher speeds in the shock. Weingartner & Draine (1999) suggested the Fe may largely be incorporated into the carbonaceous dust population rather than the silicate component. In this picture, gas-phase Fe atoms rapidly accrete onto polycylic aromatic hydrocarbon molecules, forming organometallic “sandwich” molecules, which may ultimately coagulate to form larger grains. Laboratory Fe L-edge spectra of such compounds would be very useful.

### 3.2 Other implications

The Chandra data also place a constraint on the inner warm absorbing gas which is responsible for the O 8 edge in the ASCA data. From the depth of the high-ionization absorption lines, Lee et al. (2002) found $`N_\mathrm{H}10^{23}\mathrm{cm}^2`$ for this material. If the Fe i edge observation were not available one might argue that the inner absorber contains large grains that do not redden the optical nuclear emission. However, since all of the neutral Fe can be accounted for by the outer absorber, no such Fe-bearing grains exist in the inner absorber, suggesting that the inner absorber is dust-free. This, in turn, suggests than the inner absorber is located within the dust sublimation radius, unless an alternative mechanism can be found for depleting the gas of dust at larger distances.

Indeed, if grains are to be found close to the central engine of an AGN, they must be shielded from the outpouring radiation. The radiation pressure opacity of dust is $`10^3`$ times that of electrons (Laor & Draine, 1993), so within a certain radius the dust will see a significantly super-Eddington source. The gravitational and radiation forces will balance when the dust particle is outside a mass $`1000M_{\mathrm{BH}}`$, where $`M_{\mathrm{BH}}`$ is the black hole mass of the AGN. As an example, the inner velocity profiles of S0 galaxies are approximately linear (Seifert & Scorza, 1996; Loyer et al., 1998), implying that $`M(r)r^3`$. Thus, if at $`r=1`$ pc $`M(r)=M_{\mathrm{BH}}`$, then the balancing radius for the dust will be $``$10 pc, independent of the black hole mass. Therefore, if a DWA is at 150 pc it will not feel a significant radiation force, and should not have a large outflow velocity. In the case of MCG–6-30-15, this is consistent with the preliminary results from the Chandra data (Lee et al., 2001, 2002), and the low velocity component inferred by Sako et al. (2003) from the XMM-Newton spectrum.

## 4 Summary

The DWA in MCG–6-30-15 has been the subject of two different X-ray gratings observations, but the interpretation of the spectra has been controversial (Branduardi-Raymont et al., 2001; Lee et al., 2001). However, the optical/UV reddening requires dust along the line-of-sight to the AGN. This dust will imprint observable Fe, C, O and Si edges on an X-ray spectrum, and indeed an Fe i edge was likely detected in the Chandra (Lee et al., 2001) and XMM-Newton (Turner et al., 2003) data. The gas associated with this dust has a column of a few times 1021 cm-2 and may be responsible for much of the low-ionization component of the warm absorber in MCG–6-30-15 (that is, it produces the O vii edge).

Motivated by the the dust lane seen in the HST image of MCG–6-30-15, we have argued in this paper that DWAs may reside in the ISM of the host galaxy. Photoionization modeling shows that a detectable O vii edge can be produced by such material if it is placed $`100`$$`200`$ pc away from the AGN and therefore can contribute to the warm absorption features. Furthermore, the dust contained within this gas is dynamically and thermally stable and can account for any observed reddening. Supporting evidence for this interpretation for MCG–6-30-15 is given by the very small outflow velocity of the low ionization absorption lines, the observed constancy of the O vii edge, and the additional absorption required at low energies.

The composition of the dust in DWAs may be constrained by the relative strength of the Fe i and O i edges, or by the detailed structure around the Fe i edge, as was done in Cyg X-1. The present MCG–6-30-15 data is consistent with silicate dust having a O:Fe ratio of 2:1 or less, which could result from a process of repeated grain destruction and growth in the ISM. The presence of pure iron grains seems unlikely given the limit on their abundance from our own Galaxy. High signal-to-noise data around the O i edge, or a measurement of the C i edge, is needed to fully exploit the observations in determining the dust composition in distant AGN.

###### Acknowledgements.

We thank Gary Ferland for helpful discussions regarding the Cloudy models, and Julia Lee for comments on a draft of the manuscript. This research was supported by the Natural Sciences and Engineering Research Council of Canada and by the Canada Research Chair Program. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.


format: text

text:
# Abrupt field-induced transition triggered by magnetocaloric effect in phase-separated manganites

## I Figure Caption

Fig. 1. Temperature dependence of the magnetization, measured with $`H`$ = 1 T with zero-field-cooling (ZFC), field-cooled-cooling (FCC), and field-cooled-warming (FCW) procedures. The inset shows the time evolution of the normalized magnetization after ZFC to $`T`$ = 8, 40, 60 and 100K.

Fig. 2. (a) Field dependence of the magnetization and (b) specific heat, measured at $`T`$ = 2.5 K. Both measurements show an abrupt change at the same critical field $`H_C`$ $``$ 2.2 T.

Fig. 3. Field dependence of the sample’s temperature showing an abrupt warming from 2.5 to $``$ 30 K at $`H_C`$ $``$ 2.2 T. The inset shows a spontaneous magnetization jump, measured with a fixed magnetic field.

Fig. 4. Time dependence of the magnetization during a magnetic field sweep, for different waiting times between consecutive field increments: $`t_w`$ = 7.5 (squares), 15 (circles), 30 (up triangles) and 60 min. (down triangles), at $`T`$ = 2.5 K. The inset shows an enlarged portion of the region just before the magnetization jump.


text:
# Period doubling in glow discharges: local versus global differential conductivity

## Abstract

Short planar glow discharges coupled to a resistive layer exhibit a wealth of spontaneous spatio-temporal patterns. Several authors have suggested effective reaction-diffusion-models to explore similarities with other pattern forming systems. To test these effective models, we here investigate the temporal oscillations of a glow discharge layer coupled to a linear resistor. We find an unexpected cascade of period doubling events. This shows that the inner structure of the discharge is more complex than can be described by a reaction-diffusion-model with negative differential conductivity.

Glow discharges are part of our daily environment in conventional and energy saving lamps, beamers, flat TV screens, car and street lamps as well as in various industrial applications. While applications typically try to avoid any instabilities, experiments actually exhibit a realm of spontaneous pattern formation, see, e.g.,.

An interesting series of experiments has been performed on short planar dc driven glow discharges with wide lateral extension  where the formed patterns were explored very systematically. The observed structures resemble those observed in Rayleigh-Benard convection in flat cells, in electroconvection in nematic liquid crystals, or in various chemical or biological pattern forming systems. All these systems show the formation of stripes, spots, spirals etc. In comparison to the other systems, the glow discharge system has the advantage of particular convenient experimental handling and time scales . Besides structures known from other physical systems, it continues to exhibit new structures that might be specific for this system . We will focus on the experiment in , where a complete phase diagram of different patterns was identified: homogeneous stationary and homogeneous oscillating modes, patterns with spatial and spatio-temporal structures etc.

The theory for these systems has largely focussed on effective reaction-diffusion models in the two transversal directions, and on the negative differential conductivity of the glow discharge as the driving force of pattern formation. Such models actually have been developed independently by a number of authors . On the other hand, the observation of unconventional patterns like zigzag-destabilized spirals raises doubts whether reaction-diffusion patterns are sufficient to understand the observed patterns.

In the present paper, we examine the concepts of reaction-diffusion models and negative differential conductivity on the particular case of a short DC driven glow discharge in a parameter range that exhibits spontaneous temporal oscillations but no spatial structures transverse to the current . In short, we find $`(i)`$ that a discharge on the transition from Townsend to glow discharge can combine a positive local differential conductivity with a negative global differential conductivity; $`(ii)`$ that a glow discharge in a simple electric circuit shows more complex behavior than can be expected from the proposed reaction-diffusion models  for voltage $`U`$ and current $`J`$ with (global) negative differential conductivity $`dU/dJ<0`$; $`(iii)`$ in particular, that the system can show period doubling bifurcations. Period doubling actually has been observed experimentally in glow discharges, but in more complex geometries and in longer systems . $`(iv)`$ Finally, we derive a new effective dynamical model in terms of a parameter and a function by adiabatic elimination of the electrons. There is no systematic way to reduce this model to a simpler one  with two scalar parameters like voltage $`U`$ and current $`J`$. We draw this conclusion both from direct analysis and from the occurence of period doubling in the numerical solutions.

To be precise, in the experiments of , a planar glow discharge layer with short length in the forward direction and wide lateral dimensions is coupled to a semiconductor layer with low conductivity. The whole structure is sandwiched between two planar electrodes to which a DC voltage $`U_t`$ is applied. Theoretical predictions on how the different spatio-temporal patterns depend on the parameters of the gas discharge, hardly exist. In , an effective reaction-diffusion model in the two dimensions transversal to the current is proposed. Roughly, it consists of two nonlinear partial differential equations for the current $`J`$ and the voltage $`U`$ of the form

$`_tU(x,y,t)=(U,J),_tJ(x,y,t)=𝒢(U,J),`$ (1)
where the nonlinear operators $``$ and $`𝒢`$ contain spatial derivatives $`_x`$, $`_y`$ and possibly also integral kernels. The model is of reactor-inhibitor form as studied extensively in the context of chemical and biological systems in the past decades. If applicable to gas discharges, this identification lays a connection to a realm of analytical and numerical results on reaction-diffusion systems.

To test whether a model like (1) is applicable to the gas discharge system, we will focus on its temporal oscillations that can occur in a spatially completely homogeneous mode ; hence a one-dimensional approximation is appropriate. Similar oscillations have been observed in , and similar effective models for current $`J`$ and voltage $`U`$ of the general form (1) have been proposed in .

Why have different authors come up with the same type of model? The equation for $`U`$ directly results from the simplest form of an external electric circuit: a semiconductor layer of thickness $`d_s`$, linear conductivity $`\sigma _s`$ and dielectricity constant $`ϵ_s`$ will evolve as

$`_tU={\displaystyle \frac{U_tUR_sJ}{T_s}}`$ (2)
where $`U_t`$ is the voltage on the total system, $`J`$ is the total current, and $`U=_0^{d_g}E𝑑z`$ is the voltage over the gas discharge which is the electric field $`E`$ integrated in the $`z`$ direction over the height $`d_g`$ of the discharge. For the experiments in , $`R_s=d_s/\sigma _s`$ is the resistance of the whole semiconductor layer, and $`T_s=ϵ_sϵ_0/\sigma _s=C_sR_s`$ is the Maxwell relaxation time of the semiconductor. In other experimental systems, the quantities $`R_s`$ and $`T_s`$ can have different realizations. Hence the form of the equation for $`U`$ in a reaction-diffusion model (1) is clear.

However, the equation for $`J`$ in a reaction-diffusion model as (1) is based on guesses and plausibility. Different choices have been suggested by different authors, but one thing is clear: to be physically meaningful, the current-voltage characteristics of the glow discharge has to be a stationary solution, so $`𝒢(U,J)=0`$ on the characteristics. Beyond that, there are different suggestions for the functional form of $`𝒢`$ and the intrinsic time scale.

If a model like (1) is applicable to oscillations in glow discharge systems, then the following predictions apply: 1) an oscillation can only occur in a region of negative differential conductivity of the glow discharge characteristics, 2) only a single period can be formed, period doubling is not possible, since this requires at least three independent parameters, 3) in a phase space plot in $`U`$ and $`J`$, the trajectory of an oscillation can intersect the load line $`U=U_tR_sJ`$ only parallel to the $`J`$-axis (since $`_tU=0`$ and $`_tJ0`$), and it can intersect the characteristics of the glow discharge $`U=U(J)`$ only parallel to the $`U`$-axis (since $`_tU0`$ and $`_tJ=0`$).

We now introduce the simplest classical model for a glow discharge , solve it numerically and confront its results with the predictions above.

A discharge between Townsend and glow regime consists of a gas with Ohmic conductivity for the rare charged particles, electrostatic space charge effects and two ionization mechanisms, namely impact ionization by accelerated electrons in the bulk of the discharge (the so-called $`\alpha `$-process) and secondary emission from the cathode (the $`\gamma `$-process). In its simplest form, it can be modelled by continuity equations for electron particle density $`n_e`$ and ion particle density $`n_+`$

$`_tn_e+𝐉_e=𝒮,_tn_++𝐉_+=𝒮,`$ (3)
and the Poisson equation for the electric field E in electrostatic approximation,

$`𝐄={\displaystyle \frac{\mathrm{e}}{\epsilon _0}}(n_+n_e),𝐄=\mathrm{\Phi }.`$ (4)
The particle currents are approximated as purely Ohmic

$`𝐉_e=\mu _en_e𝐄,𝐉_+=\mu _+n_+𝐄.`$ (5)
The source of particles in the continuity equation (3) is written as a sum of generation by impact ionization in Townsend approximation and recombination

$`𝒮=|n_e\mu _eE|\alpha _0\text{e}^{E_0/\left|E\right|}\beta n_en_+.`$ (6)
Finally, the secondary emission from the cathode enters as a boundary condition at the position $`d_g`$ of the cathode

$`\mu _en_e(d_g,t)=\gamma \mu _+n_+(d_g,t).`$ (7)
This is the classical glow discharge model .

We reduce the problem to one spatial dimension $`z`$ transverse to the layers which is an excellent approximation for the experiment . Furthermore, we introduce dimensionless quantities as in  by rescaling all parameters and fields as $`z=r_z/X_0`$, $`\tau =t/t_0`$, $`L=d_g/X_0`$, $`\sigma (z,\tau )=n_e(r_z,t)/n_0`$, $`\rho =n_+/n_0`$, $`=E_z/E_0`$ with the scales $`X_0=\alpha _0^1`$, $`t_0=(\alpha _0\mu _eE_0)^1`$ and $`n_0=ϵ_0\alpha _0E_0/\mathrm{e}`$. A key role is played by the small parameter $`\mu =\mu _+/\mu _e`$, which is the mobility ratio of ions and electrons.

The gas discharge layer is now modelled by

$`_\tau \sigma `$ $`=`$ $`_z(\sigma )+\sigma \alpha (),\alpha ()=e^{1/||},`$ (8)
$`_\tau \rho `$ $`=`$ $`\mu _z(\rho )+\sigma \alpha (),`$ (10)
$`\sigma (L,\tau )=\gamma \mu \rho (L,\tau ),`$
$`\rho \sigma `$ $`=`$ $`_z,`$ (11)
where recombination was neglected \[$`\beta =0`$ in (6), a discussion of this approximation follows below\], while the external circuit is described by

$`_\tau 𝒰={\displaystyle \frac{𝒰_t𝒰_sj}{\tau _s}},𝒰(\tau )={\displaystyle _0^L}(z,\tau )𝑑z`$ (12)
with the dimensionless voltage $`𝒰=U/(E_0X_0)`$, time scale $`\tau _s=T_s/t_0`$ and resistance $`_s=R_s/R_0`$, $`R_0=X_0/(\mathrm{e}\mu _en_0)`$ and with a spatially conserved total current

$`j(\tau )=_\tau +\mu \rho +\sigma ,_zj(\tau )=0,`$ (13)
where $`_zj=0`$ follows from (8), (10) and (11) as usual.

As a result, the gas discharge is parametrized by the three dimensionless parameters of system length over ionization length $`L`$, secondary emission coefficient $`\gamma `$ and mobility ratio $`\mu `$ (as discussed in ), and the external circuit is parametrized by relative resistance $`_s`$, ratio of time scales $`\tau _s`$ and dimensionless applied voltage $`𝒰_t`$.

For calculational purposes, the ion density $`\rho `$ can be completely eliminated from the one-dimensional gas discharge equations (8)–(11) with the help of the Poisson equation (11) and the total current $`j`$, see . The result are two equations of motion for $`_\tau \sigma `$ and $`_\tau `$. In our numerical calculations, the system was implemented in this form. Our choice of parameters was guided by the experiments in : we chose the secondary emission coefficient $`\gamma =0.08`$, the mobility ratio $`\mu =0.0035`$ for nitrogen and the dimensionless system size $`L=50`$ which amounts to 1.4 mm at a pressure of 40 mbar. The external circuit has $`_s=30597`$, $`\tau _s=7435`$ and a dimensionless total voltage $`𝒰_t`$ in the range between 18 and 20. This corresponds to a GaAs layer with $`ϵ_s=13.1`$, conductivity $`\sigma _s=(2.610^5\mathrm{\Omega }\mathrm{cm})^1`$ and thickness $`d_s=1.5\mathrm{mm}`$, and a voltage range between 513 and 570 V.

Fig. 1 shows electric current $`j`$ and voltage on the gas discharge $`𝒰`$ as a function of time for a total stationary voltage $`𝒰_t=19`$ applied to the complete system of gas discharge and semiconductor layer. The system exhibits spontaneous oscillations with sharp current peaks: when the voltage $`𝒰`$ on the gas layer becomes high enough, the discharge ignites. The conductivity of the gas increases rapidly and produces a current pulse that deposits a surface charge on the gas-semiconductor interface. Therefore the voltage $`𝒰`$ over the gas layer breaks down. Due to the low conductivity of the semiconductor, the voltage $`𝒰`$ recovers only slowly. Eventually the gas dicharge ignites again, and the cycle is repeated.

Note that the oscillations in Fig. 1 are not quite periodic. This is not due initial transients since the system is observed after the long relaxation time $`\tau =4.74510^6`$. The nature of this temporal structure becomes clear when the trajectory is plotted in the plane spanned by current $`j`$ and voltage $`𝒰`$ in Fig. 2(b). The figure contains the data of the time span from $`\tau =310^6`$ to $`610^6`$ which amounts to approximately 90 current pulses. The phase space plot shows that the system is actually periodic, with a period of 8 current pulses. Fig. 1 shows precisely one period.

This discovery raises the question whether our system actually follows the well-known scenario of period doubling. Indeed, it does. Fig. 2(a) for $`𝒰_t=18`$ shows an oscillation where one current pulse is repeated periodically as observed experimentally in . For $`𝒰_t=18.5`$, a period consists of two current pulses (not shown). For $`𝒰_t=19`$, the period is 8 pulses as in Fig. 1 and Fig. 2(b). For $`𝒰_t=20`$, the systems seems to have reached the chaotic state as can be seen in Fig. 2(c).

A detailed comparison of the experiments in  with simple oscillations as in Fig. 2(a) will be given elsewhere, and we only state here that there is semi-quantitative agreement of several features. Here we emphasize that period doubling events in glow discharges have been observed experimentally in other systems . However, this was always in systems with more complicated geometries like long narrow tubes, and the authors allude to general knowledge on nonlinear dynamics rather than to solutions of explicit models. We state that period doubling can be a generic feature of a simple, strictly one-dimensional glow discharge when coupled to the simple circuit (2). We propose to search experimentally for a period doubling route to chaos in such simple systems which would then allow quantitative comparison with theory.

Let us return to the initial question: is a 2-component reaction diffusion model like (1) with negative differential conductivity appropriate for the present system? Above Eq. (3), we gave a list of predictions for the reaction diffusion model (1) to be applicable. Prediction 2 is falsified by the observation of period doubling. Prediction 3 is also falsified by a simple check of either of the three figures in Fig. 2: the trajectories definitely do not intersect with the characteristics or the load line with the angle prescribed by (1), in particular not in the upper part of the figures that represent the rapid current pulses.

There rests prediction 1: is negative differential conductivity required for the oscillations to occur? We have not found a numerical counterexample where oscillations would occur while the current voltage characteristics of the gas discharge shows a positive differential conductivity, but we have found no reason to exclude its existence. Furthermore we note that the characteristics is a global property of the whole discharge layer with its boundary conditions  while the local differential conductivity in our model is always positive: the field dependent stationary ionization is $`n_+=|\mu _eE|\alpha _0e^{E_0/|E|}/\beta `$ according to (6); hence the local conductivity increases monotonically with the applied field $`|E|`$. The global negative differential conductivity is due to electrode effects being much stronger than bulk recombination $`\beta `$.

Last but not least, we have derived an analytical approximation of the model (8)–(13) that can be confronted with the suggested form (1). Electron and ion current in the gas are of the same order of magnitude. Since the electrons are much more mobile, their density is appropriately lower. Rescaling this density like $`s=\sigma /\mu `$ and time like $`\overline{\tau }=\mu \tau `$, the electrons can be eliminated adiabatically in the limit of $`\mu 0`$. Space charges in the gas discharge are then due to the ions only $`\rho =_z`$, and $`\rho `$ can be expressed by $``$. Splitting the field $`(z,t)=_L(t)+ϵ(z,t)`$ into the field on the cathode $`_L`$ and a correction $`ϵ`$ with $`ϵ(L,t)=0`$, the complete system for $`\mu 0`$ can be expressed by two dynamic equations

$`_{\overline{\tau }}_L(t)=F(_L,ϵ),_{\overline{\tau }}ϵ(z,t)=G(_L,ϵ),`$ (14)
details will be given elsewhere. While the equation for the time dependent parameter $`_L`$ corresponds to the equation for $`U`$ in (1), the space dependent field correction $`ϵ`$ within the gas layer cannot be reduced to a single component like the current $`J`$ in (1). E.g., for the ion density on the cathode $`\rho _L=_zϵ|_L`$, we can derive the equation of motion

$`_{\overline{\tau }}\rho _L=\rho _L^2_L(_z\rho )|_L+\gamma \rho _L_L\alpha (_L)`$ (15)
A two component reaction-diffusion equation for $`\rho _L`$ and $`_L`$ could result from the completely unsystematic approximation $`(_z\rho )|_L=0`$. Rather the transport of ions $`\rho `$ from the bulk of the gas towards the cathode is a central feature of the system. The field $`ϵ(z,t)`$ in (14) indeed accounts for the ion distribution within the gas gap with its intricate dynamics.

We acknowledge support of D.S. by the Dutch physics funding agency FOM and of I.R. by ERCIM.


text:
# Zeros of the Partition Function and Pseudospinodals in Long-Range Ising Models

## I Introduction

Mean-field treatments of fluids and Ising models yield metastable and unstable regions, separated by a well-defined line known as the spinodal kac . As the spinodal is approached, the system shows phenomena similar to that at mean-field critical points. In particular, the isothermal susceptibility $`\chi `$ diverges as a power law as the spinodal value of the magnetic field $`h_s`$ is approached from the metastable state uk ,

$$\chi (h_sh)^{1/2},$$
(1)
where $`h`$ is the (dimensionless) magnetic field. A mean-field Ising system can be realized in a well defined way by assuming an infinite-range interaction between the spins kac .

If the interaction range is long but finite, the system is no longer mean-field, but can be described as near-mean-field, and the spinodal singularity is replaced by a pseudospinodal, which still has physical effects. As with apparent critical points in finite size systems, the susceptibility for a finite interaction range can be fit to a power law over a limited range of the scaling field $`(h_sh)`$. For example, the susceptibility in the metastable state of a long but finite-range interaction Ising model appears to diverge over several decades in $`(h_sh)`$ as the pseudospinodal is approached stauffer . However, the divergence is suppressed if the pseudospinodal is approached too closely, indicating that the spinodal singularity has been smeared out stauffer . It also is found that the properties of the pseudospinodal converge rapidly with increasing interaction range to those predicted for the spinodal in mean-field theory novotny .

Because the spinodal is a line of critical points, we expect that the Ising spinodal has properties similar to those of the Ising critical point. In particular, we expect that the spinodal is related to the zeros of the partition function. Lee and Yang leeyang  showed that the singularity in the free energy for $`T`$ less than $`T_c`$, the critical temperature, arises from the presence of a real positive zero of the partition function in the thermodynamic limit. In finite systems there is no real positive zero for $`T<T_c`$.

The zeros of the partition function for real temperature $`T`$ lie on the imaginary axis of the complex magnetic field plane latticegas . For finite systems there is a gap in the distribution of the zeros around $`h=0`$. This gap shrinks to zero for $`T<T_c`$ as $`N\mathrm{}`$. The Lee-Yang relation between the zeros of the partition function and critical points is valid for all interaction ranges at the Ising critical point. This ideas were extended to complex temperature by Fisher fisher .

In order to generalize these ideas to spinodals in Ising models, we consider an Ising model in a magnetic field $`h`$. We will consider the “infinite-range” Ising model in which each spin interacts with all other spins milotti , and the Domb-Dalton domb  version of the Ising model in which each spin interacts with its neighbors within a given interaction range $`R`$ with a constant interaction. These models can be described by mean-field theory in the limit $`N\mathrm{}`$ and $`R\mathrm{}`$ (for $`N`$ infinite), respectively. Our main result is that the pseudospinodal is related to the zeros of the partition function in four-dimensional complex temperature/magnetic field space. In addition, the zeros approach the real temperature/magnetic field plane as the system becomes more mean-field.

The structure of the paper is as follows. In Sec. II we consider the infinite-range Ising model and show both analytically and numerically that the zeros of the partition function approach the real $`\beta `$ and $`h`$ plane as $`N`$ increases. We find a similar result in Sec. III by estimating the partition function in the metastable state by using the Metropolis algorithm and the single histogram method. In Sec. IV we consider the Domb-Dalton version of the Ising model and estimate the partition function in the metastable state in the same way.

## II The infinite-range Ising Model: Analytical Approach

We first consider an Ising model in which every spin interacts with every other spin. We will refer to this model as the infinite-range Ising model milotti , although the interaction range becomes infinite only in the limit $`N\mathrm{}`$. The Hamiltonian is

$$H=J_N\underset{ij=1}{\overset{N}{}}\sigma _i\sigma _jh\underset{i}{}\sigma _i.$$
(2)
We need to rescale the interaction so that the total interaction energy seen by a given spin remains the same as $`N`$ is increased kac . We will take

$$J_N=\frac{4J}{N1}.$$
(3)
This choice of $`J_N`$ yields the Ising mean-field critical temperature $`T_c=4`$ when $`N\mathrm{}`$ domb , where we have chosen units such that $`J/k=1`$.

The exact density of states is easily calculated for this model and is given by

$$g(M)=\frac{N!}{n!(Nn)!},$$
(4)
where $`n`$ is the number of up spins. We have $`M=2nN`$ and the total energy $`E=J_N(NM^2)/2`$, where $`M`$ is the magnetization. (In general, the density of states depends on both $`E`$ and $`M`$, but because $`E`$ is a unique function of $`M`$ in the infinite-range Ising model, we need only write $`g(M)`$.) Hence the partition function $`Z`$ can be expressed analytically for arbitrary complex inverse temperatures $`\beta `$ and magnetic fields $`h`$:

$$Z(\beta ,\beta h)=\underset{\mathrm{all}\mathrm{M}}{}g(M)e^{\beta E}e^{\beta hM}.$$
(5)
To understand the nature of $`Z`$ in the metastable state, imagine a simulation of an Ising model in equilibrium with a heat bath at inverse temperature $`\beta `$ in the magnetic field $`h=h_0>0`$. Because $`h_0>0`$, the magnetization values will be positive. Then we let $`hh_0`$. If $`h_0`$ is not too large, the system will remain in the metastable state for a reasonable amount of time and sample positive values of $`M`$. Hence, to determine the zeros of $`Z`$ associated with the pseudospinodal, we need to restrict the sum in Eq. (5) to magnetization values, $`M`$, that are representative of the metastable state. The following examples will illustrate the need for this restriction and the procedure for determining the zeros of $`Z`$.

The notion of using a restricted partition function sum to describe the metastable state has a long history. Penrose and Lebowitz penleb ; lebowitz  review such restricted partition functions and their properties. A more physical approach can be found in the discussion of the non-interacting droplet model by Langer lang . In this model fluctuations are restricted to non-interacting compact droplets of the stable phase occurring in the metastable phase. The partition function sum is restricted to droplets less than the critical size. This approximation is reasonable for low temperatures close to the coexistence curve. Langer showed that this restriction gives the same metastable state free energy as the analytic continuation of the stable state free energy in the same model. However, there are additional properties of the analytic continuation that do not appear in the restricted sum which are related to the decay of the metastable state rather than the description of the metastable state itself.

We first consider $`N=4`$ and retain only the terms in the partition function sum that correspond to the two positive values of $`M`$. From Eqs. (4) and (5), we have

$$Z_r(\beta ,\beta h)=\underset{M=2,4}{}g(M)e^{\beta E}e^{\beta hM}=e^{8\beta }e^{4\beta h}+4e^{2\beta h}.$$
(6)
The subscript $`r`$ denotes that the sum over $`M`$ is restricted. If we let $`x=e^{2\beta h}`$, the equation $`Z_r=0`$ is equivalent to,

$$e^{8\beta }x+4=0,$$
(7)
and has the solution

$$\beta h=\mathrm{ln}2i\frac{\pi }{2}+4\beta .$$
(8)
In general, we have four unknowns (the real and imaginary parts of $`\beta `$ and $`h`$); Eq. (8) yields two conditions. In the following we will fix

$$\mathrm{Re}\beta =9/16,$$
(9)
which is equivalent to a temperature of $`T=\frac{4}{9}T_c`$. For this value of $`T`$, the value of the spinodal magnetic field is known to be $`h_s1.2704`$ monette . Equation (7) then gives a line of zeros in complex $`(\beta ,\beta h)`$ space. However, if we are interested only in the zero closest to the real $`\beta `$, $`\beta h`$ plane, we need a fourth condition. This condition is found by requiring that the quantity,

$$D^2=(\mathrm{Im}\beta )^2+(\mathrm{Im}\beta h)^2,$$
(10)
be a minimum, which is equivalent to requiring that the leading zero of $`Z_r`$, the zero closest to the real $`\beta `$ and $`\beta h`$ plane, be as close to this plane as possible. If we let $`y=\mathrm{Im}\beta `$ and use Eq. (8), we can rewrite $`D^2`$ as

$$D^2=y^2+(\frac{\pi }{2}+4y)^2.$$
(11)
Because we want $`D^2`$ to be a minimum, we require

$$\frac{dD^2}{dy}=2y+2(\frac{\pi }{2}+4y)(4)=0.$$
(12)
The solution is

$$y=\mathrm{Im}\beta =\frac{2\pi }{17}0.3696,$$
(13)
and $`D=0.38097`$. Note that $`\mathrm{Re}T=\mathrm{Re}\beta /(\mathrm{Re}\beta ^2+\mathrm{Im}\beta ^2)=1.2417`$. We finally use Eq. (8) to obtain the value of complex $`h`$. The result is summarized in the first row of Table 1.

The solutions for $`N=9,16,100`$, and 1000, keeping all the positive $`M`$ contributions to the partition function, also are shown in Table 1. We see that although $`D`$ becomes smaller as $`N`$ is increased, $`|\mathrm{Re}h|`$ overshoots the mean-field value of $`h_s1.27`$ (for $`\beta =9/16`$). Hence, retaining all the positive $`M`$ terms in the partition function allows the system to explore more than the metastable state, and we need to further restrict the sum over values of $`M`$. Physically we want to exclude values of $`M`$ that would drive the system to the stable phase.

What are the appropriate values of the magnetization that will keep the system in the metastable state? One way to determine these values is to look at $`P(M)`$, the probability that the system has magnetization $`M`$ for a particular value of $`\beta `$ and $`h`$:

$$P(M)=g(M)e^{\beta E}e^{\beta hM}.$$
(14)
Figure 1 shows $`P(M)`$ for a system of $`N=400`$ spins for $`h=1.0`$ and $`\beta =9/16`$. The negative magnetization values have a relatively high probability (because $`h<0`$) and correspond to the stable phase. The positive values of magnetization have a much lower probability, and the peak at $`M=360`$ corresponds to the most probable value of $`M`$ in the metastable state. In between the peak at $`M=400`$ and the peak at $`M=360`$, $`P(M)`$ has a minimum at $`M_{\mathrm{min}}=192`$ for this value of $`h`$ and an inflection point at $`M_I=298`$. We will only include values of the magnetization in the partition function sum that are greater than $`M_I`$.

The reason for this choice of the cutoff has to do with the nature of metastability. We expect that $`P(M)\mathrm{exp}(\beta F(M))`$, where $`F(M)`$ is the metastable state free energy. We want to exclude from the partition function values of $`M`$ that correspond to states which are not characteristic of equilibrium. In the infinite range model the configurations with $`M_{\mathrm{min}}<M<M_I`$ are unstable in that the initial evolution of a fluctuation does not monotonically decay to the metastable well. This behavior translates into an initial growth of fluctuations rather than the monotonic decay expected in equilibrium.

This behavior is a consequence of the fact that the free energy, which is proportional to $`\mathrm{log}(P(M))`$, is not convex for these values of $`M`$. Obviously, this behavior holds for a range of values of $`M<M_{\mathrm{min}}`$ as well. However, we can exclude all configurations with $`M<M_{\mathrm{min}}`$ because they are in the stable free energy well and do not occur in the metastable state. Our particular choice of the cutoff is well defined, but is arbitrary to some extent as long as $`M`$ is greater than $`M_I`$.

To determine $`M_I`$, we calculate the second derivative of $`P(M)`$ as given in Eq. (14). We find the value of $`M`$ that satisfies

$$\frac{^2P(M)}{M^2}=\frac{N}{(NM)(N+M)}+\beta J=0.$$
(15)
Clearly there will be two inflection points (see Fig. 1). We choose the one closest to the metastable state maximum of $`P(M)`$. We find that the value of $`M_I`$ is independent of $`h`$, which is consistent with the idea that the free energy for this system can be written in the Landau-Ginzburg form where the magnetic field appears only in a term linear in $`M`$.

We now write the restricted partition function $`Z_r(\beta ,\beta h)`$ as

$$Z_r(\beta ,\beta h)=\underset{M=M_I}{\overset{N}{}}C_Mx^{M/2},$$
(16)
The coefficients, $`C_M`$, extend over a wide range of values and are as large as $`10^{200}`$ for the values of $`N`$ that we considered. For this reason we computed $`C_M`$ to arbitrary precision so as not to lose accuracy. The zeros of $`Z_r`$, which is a polynomial in $`x`$, were found using MPSolve bini ; mpsolve .

For a given value of $`\mathrm{Im}\beta `$, we solve for the zeros of $`Z_r`$ in Eq. (16) and find the value of $`x`$ that corresponds to the leading zero, the zero that minimizes $`D`$ in Eq. (10). We repeat this step for a range of values of $`\mathrm{Im}\beta `$ and determine numerically the value of $`\mathrm{Im}\beta `$ that yields the minimum value of $`D`$. The typical dependence of $`D`$ on $`\mathrm{Im}\beta `$ is shown in Fig. 2. From Fig. 2 we see that for $`N=400`$, $`D`$ is a minimum for $`\mathrm{Im}\beta 0.035`$. Once we know this value of $`\mathrm{Im}\beta `$, we solve for $`h`$ from the relation $`h=\mathrm{log}x/(2\beta )`$. (The value of $`x`$ was determined from the solution of $`Z_r=0`$.)

We repeat the above steps for a range of values of $`N`$ and obtain $`D`$, $`\mathrm{Im}\beta `$, $`\mathrm{Re}h`$, and $`\mathrm{Im}h`$. Our results are summarized in Table 2. Note that $`\mathrm{Im}h`$, $`\mathrm{Im}\beta `$, and $`D`$ decrease as $`N`$ increases and $`|\mathrm{Re}h|`$ approaches $`h_s=1.27`$. A plot of the zeros of $`Z`$ for the infinite-range Ising model in the $`\mathrm{Im}x`$, $`\mathrm{Re}x`$ plane is shown in Fig. 3. The values of $`D`$ listed in Table 2 are plotted as a function of $`N`$ in Fig. 4. Because this log-log plot indicates a power-law dependence, we write

$$DN^a.$$
(17)
A least squares fit gives $`a=0.659\pm 0.003`$. The estimate of the error is only statistical.

Our numerical result for the exponent $`a`$ can be understood by a simple scaling argument. In order for a mean-field approach, including the idea of a spinodal, to be a reasonable approximation, the system must satisfy the Ginzburg criterion, that is, the Ginzburg parameter $`G`$ must be much greater than unity. For the infinite-range Ising model, the Ginzburg criterion can be written as landau :

$$\frac{\xi ^d\chi }{\xi ^{2d}\varphi ^2}<<1,$$
(18)
where $`\xi `$ is the correlation length, $`\varphi `$ is the order parameter, and $`d`$ is the spatial dimension. We have $`\chi \mathrm{\Delta }h^{1/2}`$ and $`\varphi \mathrm{\Delta }h^{1/2}`$ uk , where $`\mathrm{\Delta }h=h_sh`$, to obtain

$$G=\xi ^d\mathrm{\Delta }h^{3/2}>>1.$$
(19)
For the infinite-range model, $`N\xi ^d`$. Hence,

$$G=N\mathrm{\Delta }h^{3/2},$$
(20)
up to a constant of order unity. The Ginzburg parameter $`G`$ is a measure of how mean-field the system is for finite $`N`$; the larger the value of $`G`$, the more mean-field the system is. If we keep $`G`$ constant as we approach the spinodal, we see that

$$\mathrm{\Delta }hN^{2/3}.$$
(21)
We present an argument for why $`G`$ should be held constant in Sec. V.

From Eq. (21) we see that $`\mathrm{\Delta }h`$ and the distance $`D`$ approach zero with an exponent $`a=2/3`$, in good agreement with our numerical result. Because $`\mathrm{\Delta }h`$ in Eq. (21) could be associated with $`\mathrm{Im}h`$, $`(h_s|\mathrm{Re}h|)`$, or $`D`$, we note that $`\mathrm{Im}h`$ in Table 2 goes to zero with the same exponent as in Eq. (17). We also find that $`(h_s|\mathrm{Re}h|)N^{0.61}`$.

## III The infinite-range Ising Model: Monte Carlo Approach

In general, the partition function is not known analytically. However, we can use a Monte Carlo (MC) method to determine the density of states from which we can determine an estimate of the partition function. Such an approach has been used to find the density of states for the nearest-neighbor Ising model bhanot ; alves . In this way the leading partition function zeros at the critical point have been computed and the critical exponent $`\nu `$ and corrections to scaling have been found with high precision alves . Our goal is not to obtain precise estimates of the critical exponents near the pseudospinodal, but to show that the same simulations that show an apparent divergence in the susceptibility also yield an estimate for the leading zero of the partition function which behaves as expected as the system becomes more mean-field.

To this end we use the Metropolis algorithm to equilibrate the system at temperature $`T=16/9`$ and applied magnetic field $`h=h_0`$ for about 100 MC steps per spin. (The system equilibriates as quickly as 10 MC steps per spin depending on the strength of the field.) Then we flip the magnetic field and compute the histogram $`H(E,M)`$ from which we determine the densify of states $`g(E,M)`$ and the partition function for complex $`\beta `$ and $`h`$. We save the values of $`M`$ after $`hh_0`$ and run until $`M`$ changes sign or for 5000 MC steps per spin, whichever comes first. We then throw away the first 20% of the data to ensure that the system is in metastable equilibrium and the last 20% of the data to ensure that we do not retain values of $`M`$ that are too close to the stable state. The remaining 60% of the run is used to obtain $`H(E,M)`$. We also omit any run whose lifetime in the metastable state is less than 100 MC steps per spin. Our results for $`H(E,M)`$ are not sensitive to the choice of the minimum lifetime nor the percentage of each run that we use to estimate $`H(E,M)`$. (Because $`E`$ is a function of $`M`$ for the infinite-range Ising model, we need only to compute $`H(M)`$. However, we need to compute $`H(E,M)`$ in Sec. IV.) We averaged $`H(E,M)`$ over approximately 5000 runs for each value of $`h_0`$ for a total of approximately $`1.5\times 10^7`$ MC steps per spin for a given value of $`h`$ and $`N`$. Our results for the susceptibility $`\chi `$ are given in Fig. 5. As mentioned in Sec. I, $`\chi `$ shows an apparent divergence with an mean-field exponent of 1/2 until the pseudospinodal is approached too closely.

Given the histogram $`H(E,M)`$ at $`\beta _0=9/16`$ and $`h=h_0`$, we use the usual single histogram method ferrenberg  and express the partition function for arbitrary (complex) $`\beta `$ and $`h`$ as:

$$Z_m(\beta ,\beta h)=\underset{E,M}{}H(E,M)e^{(\beta _0\beta )E}e^{(\beta _0h_0\beta h)M}.$$
(22)
Note that we do not have to determine the lower cutoff for $`M`$ because the Monte Carlo simulation only samples values of $`M`$ while the system is in a metastable state (noted by the subscript $`m`$ in Eq. (22)). As $`h_0`$ is increased for fixed $`N`$, the distance $`D`$ initially decreases, but then begins to increase as $`h_s`$ is approached too closely. That is, $`D`$ shows a minimum as a function of $`h_0`$. For each of value of $`N`$ we choose the value of $`h_0`$ for which $`D`$ is a minimum. A comparison of the histogram that was determined analytically in Sec. II and estimated by the Metropolis Monte Carlo algorithm shows similar qualitative behavior, except that the latter is approximately a Gaussian and extends to lower values of $`M`$ (see Table 3), but with a smaller amplitude. \[xx note change xx\]

Table 3 shows our results for $`h`$, $`\mathrm{Im}\beta `$, and $`D`$ for a range of values of $`N`$ at the values of $`h_0`$ that minimize the distance $`D`$. A log-log plot of $`D`$ versus $`N`$ for the data in Table 3 is shown in Fig. 6; a least squares fit gives $`a=0.60\pm 0.03`$, which is consistent with the result obtained using the exact density of states (with a cutoff). Note that $`\tau `$, the average lifetime of the metastable state, for each value of $`N`$ is approximately a constant at the value of $`h_0`$ that was chosen to minimize $`D`$ (see Fig. 7). We will use this fact to choose the value of $`h_0`$ for the long-range Ising model in Sec. IV.

## IV Long-range two-dimensional Ising model

As discussed in Sec. I, the susceptibility of long-range Ising models in the metastable state shows an apparent divergence as the applied magnetic field is increased. In the following we show that this effect of a pseudospinodal in long-range Ising models is reflected in the behavior of the zeros of the partition function as a function of the complex temperature and magnetic field. We will show that as the interaction range $`R`$ increases, the leading zero moves closer to the real plane.

Following Refs. kac,  and domb, , we consider an Ising model such that each spin interacts with its neighbors within a given interaction range $`R`$ with a constant interaction $`J=4/q`$, where $`q`$ is the number of interaction neighbors. (The factor 4 is included so that $`J=1`$ for the usual Ising model on the square lattice.) If the thermodynamic limit is taken first kac , the system is mean-field in the $`R\mathrm{}`$ limit, and the system is described by Curie-Weiss theory lebowitz . In this limit the metastable state ends at a spinodal point. The spinodal is a critical point and the susceptibility $`\chi `$ diverges as in Eq. (1).

We consider the Ising system on a square lattice with linear dimension $`L=240`$ and $`N=57600`$. The interaction range $`R`$ is defined such that a given spin interacts with any spin that is within a circle of radius $`R`$. The number of neighbors of a given spin is shown in the second column of Table 4. This system is large enough ($`L`$ is 10 times larger than the maximum value of $`R`$) for finite size effects to be minimal, but the finite size of the system implies that the zeros of the partition function must be complex for any range $`R`$.

As in Sec. III, we equilibrate the system at inverse temperature $`\beta =9/16`$ and applied magnetic field $`h_0`$ for 100 MC steps per spin. Then we flip the field and compute the histogram $`H(E,M)`$ from which we determine the density of states $`g(E,M)`$ and the partition function for complex $`\beta `$ and $`h`$. The most time consuming part of our procedure is determining the appropriate value of $`h_0`$. Although we could determine $`h_0`$ by minimizing $`D`$ as in Sec. III, we instead determined the mean lifetime $`\tau `$ of the metastable state as a function of $`h_0`$. For small $`h_0`$ (away from the spinodal), $`\tau `$ is greater than the duration of our runs which are $`2\times 10^4`$ Monte Carlo steps per spin. As $`h_0`$ approaches $`h_s`$, the lifetime begins to decrease. As mentioned in Sec. III, we found that the values of $`h_0`$ that minimize $`D`$ are near $`h_s`$ and also yield approximately constant lifetimes for different values of $`N`$. (This behavior was found for both the infinite- and long-range Ising models.) We choose the largest value of $`h_0`$ for which the mean lifetime of the metastable state just begins to decrease below $`2\times 10^4`$ Monte Carlo steps per spin. This criterion for choosing $`h_0`$ is not as sensitive as minimizing $`D`$, but is much quicker although some additional error is introduced when $`h_0`$ is chosen this way. For a given $`h_0`$, $`10^3`$ runs of $`2\times 10^4`$ MC steps per spin were done to determine the histogram $`H(E,M)`$.

The resulting values of the distance $`D`$ as a function of the interaction range $`R`$ are listed in Table 4 and are plotted in Fig. 8. We find

$$DR^{1.8\pm 0.3}.$$
(23)
The largest contribution to the above error estimate arises from the uncertainty in the values of the applied fields $`h_0`$. Similarly, we find that the difference $`(h_sh_0)`$ varies with $`R`$ as (see Fig. 9)

$$h_sh_0R^{1.97\pm 0.06}.$$
(24)
The scaling behavior of $`(h_s|\mathrm{Re}h|)`$ is similar. The systematic error due to the uncertainty in $`h_0`$ is the largest contribution to the error estimates.

The scaling behavior of $`D`$ and $`h_sh_0`$ can be understood by a scaling argument similar to that given in Sec. II. For a finite range system, the Ginzburg parameter is monette

$$G=R^d\mathrm{\Delta }h^{3/2d/4}.$$
(25)
For two dimensions, Eq. (25) becomes

$$G=R^2\mathrm{\Delta }h.$$
(26)
Hence, if we keep the Ginzburg parameter fixed, we conclude that

$$\mathrm{\Delta }hR^2,$$
(27)
which is consistent with Eqs. (23) and (24).

## V Discussion

We have shown that the pseudospinodal in Ising models has a well defined thermodynamic interpretation and can be associated with the leading zero of the partition function in complex temperature/magnetic field space in analogy with the behavior of the Ising critical point in finite systems. Our results for the nature of the approach of the leading zero of the partition function to the real temperature and magnetic field plane are consistent with simple scaling arguments.

An essential ingredient in the scaling arguments was the condition that the Ginzburg parameter $`G`$ was held constant as the system approached the pseudospinodal. As was seen in Sec. III, choosing the value of $`h_0`$ that minimizes the distance of the leading zero to the real temperature/magnetic field plane also leads to a metastable state lifetime that is found numerically to be constant. From nucleation theory near the spinodal, we know that the lifetime of the metastable state, $`\tau `$, is given by uk ,

$$\tau \frac{\mathrm{\Delta }h^{1/2}e^G}{R^d\mathrm{\Delta }h^{d/4}},$$
(28)
where $`G`$ is defined in Eq. (25). For large $`G`$ and $`R`$ as well as small $`\mathrm{\Delta }h`$ (close to the pseudospinodal), it is easily shown that constant $`\tau `$ implies constant $`G`$. To see this relation we simply replace $`R`$ by $`R+\delta R`$ and $`\mathrm{\Delta }h`$ by $`\mathrm{\Delta }h+\delta h`$ in Eq. (28) and demand that $`\tau `$ remain constant. For $`\delta R`$ and $`\delta h`$ small, constant $`\tau `$ implies constant $`G`$. Because $`G`$ is constant, the scaling arguments of Secs. III and IV follow.

The relation between the zeroes of the partition function and the spinodal provides a mathematical foundation for the notion of a pseudospinodal and clarifies the extent to which spinodals act like critical points. It also provides a possible route by which pseudospinodals in supercooled liquids can be characterized klein .

###### Acknowledgements.

We would like to thank Greg Johnson and Frank Alexander for very useful discussions. This work was carried out in part at Los Alamos National Laboratory (LA-UR 03-5959) under the auspices of the Department of Energy and supported by LDRD-2001501 and was supported in part by the National Science Foundation under grants PHY99-07949 and PHY98-01878.

## Tables

## Figures


text:
# Evolution of rotation in binaries: physical processes

## 1. Introduction

The evolution of rotation is usually associated with an evolution of angular momentum; changing the angular momentum of any body requires torques and stars do not escape from this law of physics. In binary stars there is a permanent source of torques: the tides. Hence, understanding the evolution of rotation of stars in a binary system demands the understanding of the action of tides. This compulsory exercise was started more than thirty years ago by Jean-Paul Zahn during his “Thèse d’état”, Les marées dans une étoiles double serrée (Zahn 1966). All the concepts needed to understand tides and their actions in the dynamical evolution of binaries are presented in this work.

### 1.1. Why should we consider the rotation in binary stars ?

Surely, as in isolated stars, rotation is an important ingredient of evolution through the induced mixing processes: turbulence in stably stratified radiative zones, circulations… All these processes will influence the abundances of elements in the atmospheres or the internal profile of angular velocity, for instance.

However, in binary stars new phenomena appear: tides. They make the orbit evolving, force some mixing processes (through eigenmode resonances for instance) or may even generate instabilities leading, to some turbulence (see below).

These new phenomena need also to be understood if one wishes to decipher the observations of binary stars. In this respect binary stars offer more observables than single stars like the parametres of the orbit, masses of the stars, their radii, etc. If the system has not exchanged mass during evolution and if other parameters like luminosity, surface gravity, abundances can also be determined unambiguously, binary stars offer new constrains on the stars which may be useful for our understanding of stellar evolution. Also, a statistical view of orbital parameters may constrain the dissipative processes at work in these stars (Mathieu et al. 1992).

### 1.2. Synchronisation and circularization : how it works.

Let us consider an isolated system made of two stars of mass M1, M2, of moment of inertia I1, I2 and of spin angular velocity $`\stackrel{}{\mathrm{\Omega }}_1`$, $`\stackrel{}{\mathrm{\Omega }}_2`$. The semi-major axis of the orbit is $`a`$ and the eccentricity $`e`$. For simplicity we shall assume that the angular momentum vectors are all aligned. Hence, the total (projected) angular momentum of the system, which is conserved during evolution, reads:

$$L=\frac{M_1M_2}{M_1+M_2}a^2\mathrm{\Omega }_{\mathrm{orb}}\sqrt{1e^2}+I_1\mathrm{\Omega }_1+I_2\mathrm{\Omega }_2$$
(1)
On the other hand, the total energy of the system, namely,

$$E=\frac{GM_1M_2}{2a}+\frac{1}{2}I_1\mathrm{\Omega }_1^2+\frac{1}{2}I_2\mathrm{\Omega }_2^2$$

decreases because of dissipative processes.

To appreciate the natural evolution of such a system, let us consider the even more simplified system where the angular momentum and the energy of the spin of the stars are negligible compared to their orbital equivalent. Using Kepler third law to eliminate the mean angular velocity of the orbital motion $`\mathrm{\Omega }_{\mathrm{orb}}`$, the previous equations lead to

$$L_{\mathrm{system}}a^{1/2}\sqrt{1e^2}\mathrm{and}E_{\mathrm{system}}\frac{GM_1M_2}{2a}$$

During evolution the system loses energy through dissipative mechanisms, thus $`a`$ decreases which implies that $`e`$ also decreases to insure a constant angular momentum. Thus, with time the orbit slowly circularizes.

Once the orbit is circular or nearly circular, the system may continue to evolve if the orbital angular velocity and spin angular velocity are not identical: this is the synchronization process after which the system has reached its minimum mechanical energy state: all the stars rotate at the same rate, i.e. $`\mathrm{\Omega }_{\mathrm{spin}}=\mathrm{\Omega }_{\mathrm{orb}}`$ like the moon on its terrestrial orbit.

### 1.3. Tides

In the foregoing section we presented a global view of the evolution of the system, however the way the orbit or the spin change is controlled by the torques applied to the stars. As we said at the beginning, a permanent source of torques is given by the tides which therefore need to be studied. But what is a tide?

The tide is the fluid flow generated by the tidal potential, i.e. the azimuth dependent part of the gravitational potential inside a star. In short, if you sit on a star, it is the forced flow generated by the celestial bodies orbiting around you. If you sit on Earth you feel the tides of the moon and the sun essentially.

Now let us assume that the tidal forcing is milde enough so that the fluid flow obeys linear equations; formally, we may write the system like

$$\rho \frac{\stackrel{}{v}}{t}+(\stackrel{}{v})=\stackrel{}{F}_T(\stackrel{}{r})f(t)$$

where we assume that the tidal force $`\stackrel{}{F}_T(\stackrel{}{r})f(t)`$ can be separated into its spatial and temporal dependence. Written in this way we immediately see that if the inertia of the fluid can be neglected (i.e. the term $`\rho \frac{\stackrel{}{v}}{t}`$), then the velocity field can be computed with the same temporal dependence as the exciting force. The response is instantaneous. Moreover, if Coriolis acceleration and viscosity are negligible, the only response of the fluid is through a pressure perturbation, i.e. it is purely hydrostatic. This extremely simple, but not unrealistic, case is called the equilibrium tide.

On Earth, this tide is enough to understand the basic properties of terrestrial tides: i.e. that there are two tides a day, that their amplitude is $`50`$ cm or that they are stronger at full or new moon; the hydrostatic perturbation describes the famous tidal bulge which is often representing tides in elementary courses. Such a description is appropriate if you sit on the mediterranean beaches or on those of the Pacific ocean; however, if you sit on the Atlantic shore, like here in Cancun, you will easily notice that the tide is much larger than the expected 50 cm. In the Mont Saint-Michel bay it easily reaches 10 meters !

The difference comes from what we neglected: the inertia of the fluid and the ensuing resonance phenomenon. For the tidal wave, whose wavelength is a few thousand kilometers, the Atlantic ocean is a mere puddle five kilometers deep. Surface gravity waves may thus be studied using the shallow water approximation and their speed is given by

$$V_{\mathrm{wave}}=\sqrt{gh}220\mathrm{m}/\mathrm{s}$$

where $`g`$ is the gravity and $`h`$ the depth. With a mean width of 5000 km, the Atlantic is crossed twice in 12.6 hours; but the tidal forcing is back after 12.4 hours. Obviously, we are close to a resonance and in this case the equilibrium tide is insufficient to describe the tidal response of the fluid. In this case the tide will be qualified as dynamical.

Quite clearly, the equilibrium tide is much easier to handle than the dynamical one; this is why it was first studied (Zahn 1966); the dynamical tide received first serious considerations by Zahn (1970) and a proper treatment by Zahn (1975).

### 1.4. Torques and tidal evolution

In order to compute the dynamical evolution of the system, namely the parameters of the orbit and the spin of the stars it is necessary to evaluate both the torques and the dissipation inside the stars.

The torque suffered by a star results from an unsymmetrical distribution of mass with respect to the tidal potential; mathematically,

$$𝒯_z=_{(V)}\left(\stackrel{}{r}\times \rho \stackrel{}{}\mathrm{\Phi }_T\right)_z𝑑V=_{(V)}\frac{\mathrm{\Phi }_T}{\phi }\rho ^{}𝑑V$$

where $`\rho ^{}`$ is the density perturbation generated by the tidal potential $`\mathrm{\Phi }_T`$; $`\phi `$ is the angular azimuthal variable. We see from this expression that torques can only exist if the excitation $`\mathrm{\Phi }_T`$ and the response $`\rho ^{}`$ are out of phase (or antiphase). The phase lag between these two quantities comes from the dissipative mechanisms at work in the stars and is usually a small but important quantity.

From the expressions of the energy and angular momentum of the orbit, we can derive:

$$\frac{d\mathrm{ln}a}{dt}=\frac{\dot{E}_{\mathrm{orb}}}{|E_{\mathrm{orb}}|}$$

$$\frac{de^2}{dt}=|E_{\mathrm{orb}}|^1\left((1e^2)\dot{E}_{\mathrm{orb}}\mathrm{\Omega }\sqrt{1e^2}\dot{L}_{\mathrm{orb}}\right)$$

where $`\dot{E}_{\mathrm{orb}}`$ and $`\dot{L}_{\mathrm{orb}}`$ are respectively the power dissipated in the stars and the torques exerted on the orbit (the parallelism of angular momentum vectors is still assumed).

The foregoing discussion therefore shows that, provided the stars do not change, all the evolution of the system is controlled by the dissipation of energy by the tidal flow.

In the case stellar evolution is important, for instance before or after the main sequence, changes in the inertia momenta will change the spin of the stars and the applied torques.

## 2. Dissipation mechanisms

### 2.1. Viscosity

The first and most obvious physical mechanism to dissipate mechanical energy is viscosity. However, the viscosity of stellar plasma is far too weak to be efficient and only turbulent viscosity of convection zones can significantly affect the evolution of the system.

But the effectiveness of turbulent viscosity is hampered by the fact that tidal flows are periodic in time: in any region of a convective zone where the lifetime of eddies is longer than the period of the tidal flows, the turbulent eddy viscosity is reduced. The question is, of course, how much it is reduced compared to its usual approximation

$$\nu _T^0=\frac{1}{3}V_{\mathrm{turb}}\mathrm{}_{\mathrm{turb}}$$

where $`V_{\mathrm{turb}}`$ and $`\mathrm{}_{\mathrm{turb}}`$ are respectively the velocity of turbulent eddies and their mean free path.

Presently, two prescriptions coexist: the first by Zahn (1966) says that

$$\nu _T=\nu _T^0\mathrm{min}(1,\frac{V_{\mathrm{turb}}P_{\mathrm{tide}}}{2\mathrm{}_{\mathrm{turb}}})$$

which means that when the period of tides is shorter than the turnover time of the eddies, the mean-free path of the eddies should be reduced to the distance covered by an eddy during half a period.

The second prescription is originally due to Goldreich and Keeley (1977) but adapted by Golman and Mazeh (1991); it says that

$$\nu _T=\nu _T^0\mathrm{min}(1,\left(\frac{V_{\mathrm{turb}}P_{\mathrm{tide}}}{2\pi \mathrm{}_{\mathrm{turb}}}\right)^2)$$

This prescription means that the turbulent viscosity should be such that

$$\nu _T\left(\frac{V_{\mathrm{tide}}}{\mathrm{}_{\mathrm{tide}}}\right)^2=<ϵ>$$

namely that the dissipation by tidal currents fits the dissipation by the turbulent cascade $`<\epsilon >`$ (recall that in the Kolmogorov cascade $`<ϵ>V_{\mathrm{}}^3/\mathrm{}`$).

These two prescriptions yield different exponent in the dependence of the evolution time scale with the period of the system; in principle, they can be discriminated by observations. However, this exercise turns out to be difficult; moreover, using binaries in stellar clusters, Mathieu et al. (1992) found that no reduction of viscosity was necessary to explain observations ! Obviously, more work is needed to clarify this question.

Another viscous damping process was also put forward by Tassoul (1987), namely the dissipation in Ekman boundary layers. However, we have shown that because of the stress-free surfaces of the stars, such layers are rather regions with lower dissipation than the rest of the star (see Rieutord 1992, Rieutord and Zahn 1997).

### 2.2. Radiative damping

Another way to dissipate energy of tidal flows is through radiative damping. This mechanism will affect essentially radiative zones; indeed, from the tidal excitation, one needs to generate temperature fluctuations which are dissipated by radiative diffusion.

The most natural way to generate these temperature fluctuations comes from the excitation of gravity modes since, mechanically, they are associated with the buoyancy force. It is therefore quite clear that the dissipation and the ensuing torques will be most important when the forcing frequency is near that of an eigenmode of the star. From this remark, it is also clear that only the dynamical tide will be relevant in this process.

On general grounds we may observe that the tidal forcing is low-frequency and that low-frequency gravity modes are high order modes which can then be described by an asymptotic theory. This was the way chosen by Zahn (1975) and later by Goldreich and Nicholson (1989).

This mechanism is essentially relevant for early-type stars which own an outer radiative zone. As shown by the previous authors the tidal excitation is most intense near the core-envelope boundary but since gravity waves are only partially reflected at the star surface, they deposit their angular momentum there and these layers are synchronized first. This argument was developed by Golreich and Nicholson (1989) to explain the higher synchonization rate observed in early-type stars compared to the theoretical predictions of Zahn (1977).

These studies are based on an asymptotic approach and ignore the role of the Coriolis force. This force is, however, unavoidable since the stars are always rotating and in most cases the tidal forcing is also in the band of the so-called ’inertial modes’ whose restoring force is precisely the Coriolis force. In fact inertial modes always combine with gravity modes as they are also low-frequency modes; both together, they occupy the band $`[0,\mathrm{max}(2\mathrm{\Omega }_s,N_{\mathrm{max}})]`$, where $`N`$ is the Brunt-Väisälä frequency.

Recently, much progress was achieved in the direction of including the whole spectrum of low-frequency modes in the tidal response of an early-type star. Indeed, Witte and Savonije (1999a,b,2001) computed numerically the response of a 10 M main sequence star and the ensuing evolution of the system with various inital eccentricities. Their results show many interesting features:

* The efficiency of resonance crossing at decreasing the eccentricity when it is initially small (a few percent).
* A new phenomenon, which they called “resonance locking” and by which two resonant modes, one tending to spin-up the star and the other trying to spin it down, yield equilibrating torques which maintain the two stellar modes close to resonance and therefore force a strong evolution of eccentricity (see figure 3 in Witte and Savonije 1999b).

The calculations of Witte and Savonije show that radiative damping is efficient at reducing the eccentricity of the orbit in a fraction of the lifetime of the star: in their examples, the 10 M star has a lifetime of 20 Myrs; an initial eccentricity of 0.02 is erased in 1 Myr, but an initial eccentricity of 0.25 or 0.7 reduces to 0.1 in 3 Myrs; a further evolution apparently operates on a much longer time scale.

Obviously, if the model of Witte and Savonije could be used as a true system, it would mean that orbits with $`e=0.1`$ may be as evolved as $`e=0`$ ones and differ only by initial conditions.

### 2.3. A new possibility: the elliptic instability

Likely, all the mechanisms by which energy can be dissipated have not been examined. We shall now present a new one caused by the elliptic instability and which is potentially a rather strong source of dissipation.

The elliptic instability has been studied mainly in simple geometries and even in these cases it is a difficult problem. For a recent review we refer the reader to the work of Kerswell (2002). Presently, the work closest to astrophysical applications is that of Seyed-Mahmoud et al. (2000) who investigated this instability in an ellipsoidal configuration for the core of the Earth.

The basic result we need from fluid mechanics is the growth rate $`\gamma `$ of this instability, namely

$$\gamma \epsilon \mathrm{\Omega }$$

where $`\epsilon `$ is the ellipticity of streamlines and $`\mathrm{\Omega }`$ the rotation rate of the vortex.

Schematically, this instability may be seen as a parametric instability: the solid body rotating fluid feels a perturbation (the ellipticity of the stream lines) with a frequency $`2\mathrm{\Omega }`$. Such a periodic forcing can destabilize modes at half its frequency namely $`\mathrm{\Omega }`$; this is precisely the frequency of the so-called spin-over modes (or Poincaré modes). Such modes are solid-body rotation around an axis in the equatorial plane. When the instability develops and the amplitude gets sufficiently large, the vortex start to precess and is usually completely destroyed as observed in experiments (Malkus 1989). Hence, this instability may be very important as it is able to dissipate the total kinetic energy of the vortex. But let us consider the astrophysical case.

For the sake of simplicity we consider a main sequence star in a binary system with a circular orbit. The tidal perturbation is provided by a point-mass object. The system is not synchronized and therefore $`\mathrm{\Omega }_{\mathrm{spin}}\mathrm{\Omega }_{\mathrm{orb}}`$. In a reference frame rotating with the tidal potential, the main sequence star is like a strained vortex rotating with the angular velocity $`\mathrm{\Omega }_{\mathrm{spin}}\mathrm{\Omega }_{\mathrm{orb}}`$ and enduring an elongation $`\epsilon `$ of the tidal potential.

In such a case the instability develop on a time scale $`\left(\epsilon (\mathrm{\Omega }_{\mathrm{spin}}\mathrm{\Omega }_{\mathrm{orb}})\right)^1`$. The energy available for dissipation is

$$E\frac{1}{2}I_{EIZ}(\mathrm{\Omega }_{\mathrm{spin}}\mathrm{\Omega }_{\mathrm{orb}})^2$$

where $`I_{EIZ}`$ is the inertia momentum of the zone developing the elliptic instability. In this first attempt, we shall consider a late type stars and restrict, conservatively, the action of this instability to the convection zone so that it does not interfere with the stratification. We note that the spin-over modes which are destabilized are rigid rotations and are thus not affected by the large turbulent viscosity of the fluid.

Hence, the power dissipated is approximately $`\frac{\epsilon }{2}I_{\mathrm{CZ}}(\mathrm{\Omega }_{\mathrm{spin}}\mathrm{\Omega }_{\mathrm{orb}})^3`$

With these two quantities, one can estimate the time scale over which synchronization occurs: this is typically the growth of the elliptical instability. To have orders of magnitude, let us take two solar type stars orbiting their center of mass in 10 days and let their spin be twice faster (i.e. $`\mathrm{\Omega }_{\mathrm{spin}}=2\mathrm{\Omega }_{\mathrm{orb}}`$). The time scale is then 64 years, thus very short.

As far as circularization is concerned, we need to restrict ourselves to weakly eccentric orbits ($`e<0.1`$ say). Using the stars to dissipate energy while the angular momentum of the orbit is assumed constant, the reduction of the eccentricity also occurs on a short time scale, namely a thousand of years.

The numbers thus derived are rather robust as they depend on geometric quantities of the system. In view of the observations, which show that they are systems with a 10 days period which are not circularized nor synchronous, we may wonder why the elliptic instability is not so efficient in binary stars. The answer is likely in the much more complex setup yielded by stars compared to their equivalent in the laboratory. In stars, there are stratification, magnetic fields, time variable ellipticity (on eccentric orbits), rotation etc… All these effects have been explored but most of the time using systems which are not quite similar to a star. Nevertheless, we shall list them and speculate about their effects in a stellar situation.

* Rotation of the frame associated with the orbital motion influences the elliptic instability through the Coriolis force. Craik (1989) has shown, but using an unbounded strained flow, that rotation is either stabilizing or destabilizing. It is stabilizing if, for instance, in an inertial frame, the star does not rotate.
* Stratification: Similarly to rotation, stratification does not act in a unique sense; using the context of an elliptical cylinder Kerswell (1993) has shown that stratification is stabilizing for polar regions of a star but destabilizing for equatorial regions.
* Magnetic fields were found invariably stabilizing (see Kerswell 1994).
* Time-periodic ellipticity is likely destabilizing through parametric instabilities. However, the only studied cases are those of an unbounded strained fluid (Kerswell 2002).
* Nonlinear effects are the most difficult to appreciate. Experiments have shown that the nonlinear development is violent (Aldridge et al. 1997) and rapidly leads to a turbulent state (Malkus 1989). In fact, it seems that the saturated state exists only in a very narrow range of parameters and beyond this range secondary instabilities due to triadic interactions of inertial modes lead to small scale motions and a turbulent state (Mason and Kerswell 1999).

Globally, the way the instability saturates may be thought as a change of the spin axis of the fluid so as to reduce the ellipticity of the streamlines. In stars such an effect is possible if the spin axis of stars is inclined in the plane normal to the orbital plane and passing through the centres of mass. This shows that the apparently inocuous hypothesis of angular momentum vectors alignment may be crucial to the instability. It may also be an observational signature. Indeed, we noticed that the turbulent viscosity of convection zone could hardly inhibit the resonance of spin-over modes; however, it can easily suppress secondaries instabilities which, in laboratory experiments lead to a turbulent state. It may well be that the saturated state is difficult to obtain in the laboratory but much more easily in stars.

The foregoing discussion shows that many points of the elliptic instability in the stellar context remain in the shadows; but because of its great potential dissipative power, this instability deserves more study.

## 3. Concluding remarks

Born together, stars of a binary system share the same age and the same initial metallicity. In most favourable cases, their mass, radii, rotation rates can be determined. These are of course very interesting pieces for the puzzle of stellar evolution but there is a price to pay for these additional informations: the two stars interact during their whole lifetime. In its simplest form this interaction is of gravitational origin and generates tides. As we have shown, tides generate various fluid flows giving rise to transport processes which may be observationnally constrained if a comparison with analogous isolated stars is possible.

The foregoing presentation which sketched out all the known mechanisms by which energy is dissipated also shows that the situation is not simple. Various processes can dissipate energy and we suggest that among them the elliptic instability plays a non negligible part. Quite surprisingly, this instability has been overlooked until now. Despite its strength, shown by laboratory experiments, observations do not show, yet, an evidence of this instability. In the list of the mechanisms which may inhibit this instability, we noticed the misalignement of rotation axis. As such a misalignement may also result from a saturated state of this instability and can be an observational signature, it deserves some study.

To conclude this contribution let us point out some directions of research:

* On the theoretical side more work is obviously needed on the effects of the elliptic instability.
* On models, the integration of stellar evolution combined with dynamical (tidal) evolution following up the works of Claret and Cunha (1997) and Witte and Savonije (1999b) will be useful to understand, for instance, the statistical properties of eccentricities as a function of periods and ages, or the relative importance of pre-main sequence, main sequence, post-main sequence phases.
* On the observational side, much data are needed. First, the elements describing the dynamics are very much desired: these are the elements of the orbit $`(a,e,i,\mathrm{\Omega },\omega )`$, the masses M1, M2, the spins $`\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2`$ and their variations with time. For instance, the motion of the apsidal line $`\dot{\omega }`$ is a quantity constraining the mass distribution of the stars and therefore their internal rotation. It is clear that such data will require a lot of efforts but they will help us much in our understanding of stellar evolution. A pulsar like J0045-7319, which travels around an early-type star on a highly eccentric orbit, offers a good step in this direction (Lai et al. 1995).

## References

Aldridge K., Seyed-Mahmoud B., Henderson G. and van Wijngaarden W., 1997, Phys. Earth Plan. Int., 103, 365

Claret A. and Cunha N., 1997, A&A, 318, 187

Craik A., 1989, J. Fluid Mech., 198, 275

Goldreich, P. and Keeley, D., 1977, ApJ, 212, 243

Kerswell, R., 1993, Geophys. Astrophys. Fluid Dyn., 71, 105

Kerswell, R., 1994, J. Fluid Mech., 274, 219

Kerswell, R., 2002, Ann. Rev. Fluid Mech., 34, 83

Lai D., Bildsten L. and Kaspi V., 1995, ApJ, 452, 819

Malkus W., 1989, Geophys. Astrophys. Fluid Dyn., 48, 123

Mason D. and Kerswell R., 1999, J. Fluid Mech., 396, 73

Mathieu R. D., Duquennoy A., Latham D. W., Mayor M., Mermilliod T., Mazeh, J. C., 1992, Binaries as Tracers of Stellar Formation. Proceedings of a Workshop held in Bettmeralp, Switzerland, Sept. 1991, Edts A. Duquennoy, M. Mayor, p. 278

Rieutord, M., 1992, A&A, 259, 581

Rieutord, M. and Zahn, J.-P., 1997, ApJ, 474, 760

Seyed-Mahmoud B., Henderson G. and Aldridge K.D., 2000, Phys. Earth Plan. Int., 117, 51

Tassoul, J.-L., 1987, ApJ, 322, 856

Witte, M. and Savonije, G., 1999a, A&A, 341, 842

Witte, M. and Savonije, G., 1999b, A&A, 350, 129

Witte, M. and Savonije, G., 2001, A&A, 366, 840

Zahn, J.-P., 1966, Ann. Astrophys., 29, 313, 489, 565.

Zahn, J.-P., 1970, A&A, 4, 452

Zahn, J.-P., 1975, A&A, 41, 329