text:
# The response of the Fe K𝛼 line to changes in the X-ray illumination of accretion discs

## 1 Introduction

The discovery of the iron K$`\alpha `$ line and Compton reflection in the X-ray spectra of accreting black holes was an important step in the understanding of the central engine (Mushotzky, Done & Pounds, 1993). These features were predicted to result from the reprocessing of X-rays by an optically thick and cold medium (Lightman & White, 1988; Guilbert & Rees, 1988; George & Fabian, 1991; Matt, Perola & Piro, 1991), quite likely the accretion disc itself. This was spectacularly confirmed by the ASCA observation of a broad Fe K$`\alpha `$ line in the X-ray spectrum of the bright Seyfert 1 galaxy MCG–6-30-15 (Tanaka et al., 1995). The profile of this line was well fit by a model of line emission from relativistically moving material within 10 Schwarszchild radii of a supermassive black hole (Fabian et al., 1989). Alternative explanations for such a broad line suffer from physical inconsistencies and/or fine-tuning of model parameters (Fabian et al., 1995; Reynolds & Wilms, 2000). Suddenly a potentially powerful probe of accretion and black hole physics was observationally accessible (Fabian et al., 2000). Many other active galactic nuclei (AGN) were subsequently observed by ASCA in search of a broad Fe K$`\alpha `$ line, and, although most of the detections were of far less quality than the one of MCG–6-30-15, the mean line profile of a sample of Seyferts seemed to be broadened in a similar manner (Nandra et al., 1997a; Yaqoob et al., 2002).

Since AGN themselves have quite variable X-ray continua, it was anticipated that the Fe K$`\alpha `$ line should also change over an observation. The largest dataset has come from long observations of MCG–6-30-15 where the line is seen to sometimes drastically change energy (Iwasawa et al., 1996, 1999), and its flux varies independently of the continuum (Vaughan & Edelson, 2001). This last property, also seen in other objects (Wang et al., 1999; Nandra et al., 1999; Chiang et al., 2000; Wang, Wang & Zhou, 2001; Weaver, Gelbord & Yaqoob, 2001), is contrary to the predictions of the simple reflection scenario, which predicts the line variations should track the continuum. The fact that this is not observed suggests a more complicated and dynamic pattern of illumination on the disc (Reynolds, 2000) such as in the model of irradiation by magnetic flares in a patchy corona (Galeev, Rosner & Vaiana, 1979; Haardt & Maraschi, 1991; Haardt, Maraschi & Ghisellini, 1994; Merloni & Fabian, 2001). These flares may be rotating (Ruszkowski, 2000), outflowing (Beloborodov, 1999; Malzac, Beloborodov & Poutanen, 2001) and/or temporally connected (Poutanen & Fabian, 1999; Merloni & Fabian, 2001), and so they will clearly have an impact on the observed Fe K$`\alpha `$ line due to the changing pattern of radiation on the surface of the disc.

Many models of X-ray reflection from AGN-like accretion discs have been published (e.g. Ross & Fabian, 1993; Życki et al., 1994; Magdziarz & Zdziarski, 1995; Ross, Fabian & Young, 1999), but they have not progressed to the same dynamic level as the models of accretion disc coronae. The best current models compute the reflection spectrum from a photoionized layer on the surface of an accretion disc in hydrostatic equilibrium (Nayakshin, Kazanas & Kallman, 2000; Ballantyne, Ross & Fabian, 2001; Różańska et al., 2001). If discs are anything like what standard theory predicts, then these calculations should be able to make specific predictions about reflection spectra (e.g. Nayakshin, 2000; Nayakshin & Kallman, 2001).

Here, we compute such spectra from an AGN accretion disc over a wide range of illumination conditions, and examine the behaviour of the energy, equivalent width (EW) and flux of the Fe K$`\alpha `$ line. Our approach differs from that of Nayakshin & Kallman (2001) or Życki & Różańska (2001) by concentrating on the observationally accessible changes to the Fe K$`\alpha `$ line. The parameters are chosen so that conditions appropriate for both flaring and quiescent regions of the disc are explored. Comparison of the model predictions with time-averaged data (to avoid any non-equilibrium effects) may allow constraints to be placed on how discs are irradiated.

The paper is structured as follows. Section 2 describes the model, the assumptions, and the range of parameters that are treated. The results of the calculations and a discussion of the evolution of the Fe K$`\alpha `$ line profile, flux and EW are given in Section 3. A discussion of the results is given in Section 4 before Section 5 summarizes the primary conclusions.

## 2 Computations

We employ the code of Ross & Fabian (1993), recently extended by Ballantyne et al. (2001), to calculate the X-ray reflection spectrum from the ionized layer on the surface of an AGN accretion disc. The details of the calculations are described in the above two papers (see also Ross (1978) and Ross, Weaver & McCray 1978), and so only a brief description is given here.

The model calculates the thermal, ionization and density structure of the top five111When $`F_x/F_{\mathrm{disc}}`$$`=270`$ a total Thomson depth of 10 was used to extend the computational domain beyond the layer that is highly ionized by the illuminating radiation. Thomson depths of an accretion disc that is irradiated by an external source of X-rays from above, and a soft blackbody from the disc below. The calculation is one-dimensional, with the height of the bottom of the atmosphere given by the disc equations of Merloni, Fabian & Ross (2000) in the gas-pressure dominated regime and assuming $`\alpha =0.1`$. The flux of soft radiation entering the layer from below is assumed to result from viscous heating deep in the disc and is calculated using the standard equation from Shakura & Sunyaev (1973). The impinging hard X-rays are in the form of a power-law with photon index $`\mathrm{\Gamma }`$ (so that photon flux $`E^\mathrm{\Gamma }`$), which strikes the surface with net flux $`F_x`$ and at an angle of $`i`$ degrees to the normal. The incident X-rays extend from 1 eV and terminate with a sharp cutoff at 100 keV.

The incident radiation is transferred analytically using a one-stream approximation (e.g., Rybicki & Lightman, 1979). Diffuse radiation within the layer (i.e., incident X-rays that have undergone Compton scattering, other emission from the gas itself, and the soft photons working their way outward from the disc below) is treated via the Fokker-Planck/diffusion approximation (Ross et al., 1978), while the Fe K$`\alpha `$ line is treated in detail via escape and destruction probabilities assuming a Lorentz line profile (see Ross 1978, Ross & Fabian 1993). Hydrogen and helium are assumed to be fully ionized everywhere within the layer, and the following ionization stages of the most abundant metals are treated: C v–vii, O v–ix, Mg ix–xiii, Si xi–xv, and Fe xvi–xxvii. The elemental abundances are given by Morrison & McCammon (1983).

The major differences between this code and the one presented by Nayakshin et al. (2000) are in the method of radiative transfer (we use Fokker-Planck/diffusion, as opposed to a variable Eddington factor scheme) and in the level of detail that atomic physics is treated (Nayakshin et al. (2000) employ xstar for the photoionization calculations and so includes many more transitions than our code). However, comparisons between the two programs yield only minor qualitative differences in the output spectra, particularly around the Fe K$`\alpha `$ line (Péquignot et al., 2002).

To isolate the Fe K$`\alpha `$ variations due to a changing illumination pattern, the accretion disc parameters were set to representative values. A typical black hole mass of 108 M (Gebhardt et al., 2000) was used in the calculations, as was an accretion rate of 0.001 of the Eddington rate, which seems typical of some X-ray bright Seyfert 1 galaxies (Morales & Fabian, 2002). A radius of 7 Schwarszchild radii was selected for where the reflection takes place, but we do not take into account relativistic blurring. With these parameters, the soft flux entering the atmosphere from below is $`F_{\mathrm{disc}}=1.2\times 10^{13}`$ erg cm-2 s-1, and it is assumed that there is no energy dissipation in the corona. Other choices for these system parameters would result in different Fe K$`\alpha `$ EWs and fluxes, but the resulting trends with the illumination parameters should not be affected (Sect. 3.3).

Models were calculated for $`1.5\mathrm{\Gamma }2.1`$, which is the range most often observed for Seyfert 1 galaxies (Nandra et al., 1997a), and with the incident flux covering $`F_x/F_{\mathrm{disc}}`$$``$ 0.5 – 270. To test for geometrical effects (which may be important for magnetic flares where the size of the hard X-ray emitting region is much smaller than the disc), for each value of $`\mathrm{\Gamma }`$ and $`F_x`$ models were calculated for X-ray incidence angles $`20i80`$ degrees in steps of 20 degrees. Finally, we neglect the effects of inclination along the line of sight, but since this does not depend on the disc or illumination model, earlier results (e.g., Matt, Fabian & Ross, 1996) will still apply.

## 3 Results

### 3.1 Fe K$`\alpha `$ line profiles

The properties of the Fe K$`\alpha `$ line depend greatly on the ionization state of the gas (e.g., Matt, Fabian & Ross 1993, 1996; Ross et al. 1999). When the gas is only weakly ionized or neutral (so that Fe i–xvi are the dominant ionization states), a strong line appears at 6.4 keV and the associated absorption edge is close to 7.1 keV. As the gas becomes more ionized, intermediate Fe states such as Fe xvii–xxiii are dominant. However, the K$`\alpha `$ lines from these ions are preferentially suppressed due to Auger destruction (Ross, Fabian & Brandt, 1996). It is not until the gas becomes so ionized that He-like iron (Fe xxv) is the dominant species that the observed K$`\alpha `$ line shifts in energy. This ion has a large fluorescent yield resulting in a strong line and absorption edge at 6.7 keV and $`8.5`$ keV, respectively222Oelgoetz & Pradhan (2001) showed that under certain conditions the dielectronic satellite lines of Fe xxv could move the energy centroid of the line closer to 6.6 keV.. The line from hydrogenic Fe at 6.97 keV is also present in these situations, but is weaker than the helium-like line because of resonant scattering. Finally, it is possible for Fe to be fully stripped and neither a line nor an edge will be imprinted on the reflection spectrum.

In variable density models it is not possible to describe the spectrum by assigning an ionization parameter to it. However, the general evolution of the Fe K$`\alpha `$ line that was described above is seen in these new calculations, but the addition of a realistic density structure can create significant departures from this trend.

Figure 1 shows how the Fe K$`\alpha `$ region of the reflection spectrum can vary for two different values of $`F_x/F_{\mathrm{disc}}`$ when the photon index of the illuminating spectrum changes. At the lower value of $`F_x/F_{\mathrm{disc}}`$, the Fe K$`\alpha `$ emission lines become more ionized (i.e., has a significant He-like component) with larger values of $`\mathrm{\Gamma }`$ than with smaller values. This is possibly counter-intuitive because the harder irradiating spectra have more Fe-ionizing power than the softer ones. However, when $`F_x/F_{\mathrm{disc}}`$ is increased (the bottom panel of Fig. 1) the situation is reversed, and the Fe K$`\alpha `$ line becomes highly ionized and weak at the lowest values of $`\mathrm{\Gamma }`$.

The explanation for this behaviour is illustrated

in Figure 2. The plots show that when $`\mathrm{\Gamma }=1.5`$ the temperature of the photoionized gas falls abruptly from a very hot ($``$ Compton temperature) fully ionized phase to a cold, recombined phase. Nayakshin et al. (2000) showed that such behaviour is due to the well-known thermal ionization instability of gases in pressure balance (e.g., Krolik, McKee & Tarter 1981). There is a negligible fraction of gas with intermediate ionization states, so, when $`F_x/F_{\mathrm{disc}}`$ is relatively small, only a neutral Fe K$`\alpha `$ line is emitted with any real strength333Thermal conduction between the two temperature zones may be important under some conditions, allowing enough gas at intermediate temperatures to imprint spectral features on the reflection spectrum (Różańska, 1999; Li, Gu & Kahn, 2001).. At higher values of $`F_x/F_{\mathrm{disc}}`$ the incident radiation ionizes further into the atmosphere and the transition occurs at a non-negligible Thomson depth. Few K-shell ionizing photons can traverse this layer to the neutral material below without being scattered outwards or absorbed by a higher ionized species of iron. Therefore, at higher values of $`F_x/F_{\mathrm{disc}}`$, the line weakens and disappears at the lowest values of $`\mathrm{\Gamma }`$. As the photon-index is increased, the maximum temperature of the gas decreases (Guilbert, Fabian & McCray 1983 and Sect. 4.1) and the instability cannot operate. Other ionized species of Fe can exist near the surface and a He-like line appears in the spectrum. From Fig. 1 we see that one effect of this instability is to suppress ionized Fe K$`\alpha `$ lines when $`\mathrm{\Gamma }1.6`$ (but see Sect. 4.1).

A typical example of how the Fe K$`\alpha `$ line evolves with increasing $`F_x/F_{\mathrm{disc}}`$ is shown

in Figure 3. When $`F_x/F_{\mathrm{disc}}`$=0.5 or 1.0 then the surface of the accretion disc is not highly ionized and a neutral line at 6.4 keV dominates the line profile. As $`F_x/F_{\mathrm{disc}}`$ is increased the He-like line strengthens and then weakens as the atmosphere is further ionized. Also, the line profile gradually broadens due to Compton scattering. This trend is similar to the predictions of constant density reflection models (Ross et al., 1999).

Finally, we will briefly discuss the impact of the incidence angle of the incoming radiation.

When the angle of the beam is increased, the radiation encounters any specific line-of-sight Thomson depth at a smaller physical depth into the atmosphere. Therefore, the energy is dumped into a smaller physical region, which increases the temperature of the gas (see also Nayakshin et al., 2000). When $`F_x/F_{\mathrm{disc}}`$ is relatively small this can enhance the thermal instability, as is shown in the top panel of Figure 4. In that case, the He-like component in the line profile disappears when $`i=80`$ degrees. However, when $`F_x/F_{\mathrm{disc}}`$ is greater (the bottom panel of Fig. 4) the line photons are emitted at a significant Thomson depth, and the hot outer layers result in a weaker overall line profile.

### 3.2 Equivalent width and line flux

The two most important observational diagnostics of the strength of the Fe K$`\alpha `$ line are its equivalent width (EW) and the flux in the line. Therefore, it is important to investigate how they vary as the illumination of the disc is changed. To facilitate comparison with observations, the illuminating power-law was added to the reflection spectrum before the measurements were made. This can be thought of as a reflection fraction of unity, where the reflection fraction measures the strength of the reflected component in an observed spectrum. In the disc-corona model of an AGN, reflection fractions close to one are expected.

The EWs of the Fe K$`\alpha `$ lines were calculated by performing the following integral on the total reflected+incident spectrum:

$$\mathrm{EW}=_{5.50\mathrm{keV}}^{7.11\mathrm{keV}}\frac{EF_s(E)EF_c(E)}{EF_c(E)}𝑑E,$$
(1)
where $`EF_s(E)`$ is the spectral flux of the sum, and $`EF_c(E)`$ is the estimated spectral flux of the continuum at energy $`E`$. The continuum was estimated by fitting a straight line to the reflection spectrum between 5.50 and 7.11 keV (see Figure 5).

The line flux was computed by an analogous integral over the same limits.

The lower limit of the integration was chosen to take into account any Compton broadening of the line profile, and the upper limit was set to the energy of the photoelectric absorption edge of neutral Fe. Using the models where $`F_x/F_{\mathrm{disc}}`$$`=33`$, $`i=20`$ degrees and $`\mathrm{\Gamma }=`$1.5–2.1 a comparison was made between the EWs calculated by this method and the ones calculated by xspec using a simple power-law plus Gaussian model. An ASCA SIS-0 response matrix was used to fake a 500 ks observation for each model spectrum, and the fits were performed on the simulated 3–10 keV data with the width of the Gaussian ($`\sigma `$) fixed at 0.25 keV. The EWs calculated by xspec (43–264 eV) were only $`50`$% larger than those obtained from the above method (26–156 eV). Larger EWs would result if the reflection fraction was greater than unity.

Figure 6 shows how the flux in the Fe K$`\alpha `$ line depends on the parameters of the illuminating radiation. Each panel plots the line flux (normalized to $`F_{\mathrm{disc}}`$) versus $`F_x/F_{\mathrm{disc}}`$ for the four different incident angles at a specific photon index.

We find that the line flux is directly proportional to the flux incident on the disc since the amount of energy emitted by the atmosphere must follow the amount that it absorbs. However, the correlation flattens as the line becomes highly ionized and very weak. This is noticeable for $`\mathrm{\Gamma }1.7`$. An interesting feature about this plot is that it illustrates roughly at which $`F_x/F_{\mathrm{disc}}`$ optical depth effects become important. At low $`F_x/F_{\mathrm{disc}}`$ the line flux increases with increasing incidence angle because line photons are produced closer to the surface and can easily escape (George & Fabian, 1991). But this trend reverses at higher values of $`F_x/F_{\mathrm{disc}}`$ where there is a thicker ionized skin and ionizing photons are easily scattered out of the layer before reaching any neutral material. The amount of photons “lost” to scattering increases with $`i`$ and so the line flux drops as the incidence angle increases. The value of $`F_x/F_{\mathrm{disc}}`$ at which the crossover occurs increases with $`\mathrm{\Gamma }`$ because greater fluxes are needed to produce a deep enough ionized layer for the softer incident spectra.

Turning to the measurements of the equivalent width, Figure 7 plots the Fe K$`\alpha `$ EW versus $`F_x/F_{\mathrm{disc}}`$.

At low $`F_x/F_{\mathrm{disc}}`$ a strong neutral Fe K$`\alpha `$ line is present in all the reflection spectra, so the EWs are high. As $`F_x/F_{\mathrm{disc}}`$ is increased the atmosphere becomes ionized, and where the ionization instability operates ($`\mathrm{\Gamma }=1.5`$ & 1.6) no other line replaces the neutral line so the EW rapidly becomes very small. For softer illuminating spectra, a He-like line will replace the neutral line, and the EW can remain high, but, as $`F_x/F_{\mathrm{disc}}`$ grows, this line will ultimately weaken and the EW will decrease. As seen in Fig. 6 the incidence angle crossover marks the onset of an optically thick ionized layer at the surface of the atmosphere. This point moves to higher $`F_x/F_{\mathrm{disc}}`$ with larger $`\mathrm{\Gamma }`$. Finally, it is worth emphasizing that for $`\mathrm{\Gamma }1.8`$ the Fe K$`\alpha `$ EW can stay roughly constant over a factor of about 10 in $`F_x/F_{\mathrm{disc}}`$.

One other interesting correlation that can be shown with this dataset is illustrated in Figure 8.

In this plot the models with $`F_x/F_{\mathrm{disc}}`$$`2.08`$ are shown with thin solid lines and crosses. In these cases the Fe K$`\alpha `$ line complex is dominated by the neutral line for most values of $`\mathrm{\Gamma }`$, and there is only a weak dependence of the EW on the photon index. However, when there is a strong ionized component to the line (the curves with solid points), there is a significant positive correlation between the Fe K$`\alpha `$ EW and $`\mathrm{\Gamma }`$. This is a possible explanation for the tentative EW-$`\mathrm{\Gamma }`$ correlation in the sample of Lubiński & Zdziarski (2001). Furthermore, it may be possible to use these variations as an independent observational diagnostic to determine the ionization state of the accretion disc when other effects (e.g., relativistic broadening) makes it difficult to estimate the energy centroid of the line.

### 3.3 On the self-consistency of the results

These computations make the zeroth-order assumption that the changes in the irradiating flux, $`F_x`$, will not affect the accretion disc structure. In reality, changes in $`F_x`$ are a result of the changes in disc structure, which will also likely affect $`F_{\mathrm{disc}}`$. Here, we consider the impact of our assumption on the evolution of the Fe K$`\alpha `$ line.

An increase in $`F_x`$ would most probably result from an increase in $`f`$, the fraction of the viscous energy dissipated in the corona. However, $`f`$ is very likely to be a function of the accretion rate $`\dot{m}`$. The strength of the dependence of these three parameters on each other is subject to the unknown details of how energy is dissipated and transported in accretion discs, but the dependence of the disc structure on $`\dot{m}`$ and $`f`$ is well known (e.g., Svensson & Zdziarski, 1994; Merloni et al., 2000). For the low accretion rate assumed for these calculations, the disc is gas-pressure dominated and the scale height of the disc $`H\dot{m}^{1/5}(1f)^{1/10}`$. A drop in $`H`$ will increase the density of the disc, making it more difficult to produce ionized lines. However, this will only occur if $`F_x`$ is strongly dependent on $`\dot{m}`$ and/or $`f`$. It is likely that as $`F_x`$ is made larger by an increase in $`\dot{m}`$ the disc moves from a gas-pressure dominated phase to a radiation pressure dominated one. In this case $`H\dot{m}(1f)`$, but is much larger in amplitude (by a factor of $`10^3`$) than in the gas-pressure dominated regime. Therefore, the density is much lower and ionized Fe K$`\alpha `$ lines are easily produced over a wide range of $`F_x/F_{\mathrm{disc}}`$ (e.g., $`F_x/F_{\mathrm{disc}}`$=0.48–2.88, Ballantyne, Ross & Fabian, 2002).

The disc flux $`F_{\mathrm{disc}}\dot{m}(1f)`$, and we have seen that ionized lines are common when $`F_x/F_{\mathrm{disc}}`$$`2`$. Therefore, if $`F_x`$ has an equal or greater dependence on $`\dot{m}(1f)`$ then ionized lines should remain common. However, as commented above, there are very few constraints on the relations between $`F_x`$, $`\dot{m}`$, and $`f`$. Detailed models of energy transport and dissipation in a disc-corona system are required to fully test the validity of these calculations.

## 4 Discussion

### 4.1 The role of the Compton temperature

The results presented in the previous section have shown that the thermal ionization instability can influence the evolution of the Fe K$`\alpha `$ line in certain circumstances. This instability has been investigated previously in the context of broad-line region clouds (Krolik et al., 1981) and irradiated accretion discs (Różańska, 1999; Nayakshin et al., 2000). Here, we combine the previous work with the insight gained by our large dataset to determine, with the minimum of assumptions, the relative importance of the instability on the Fe K$`\alpha `$ line.

As discussed by Krolik et al. (1981) and Guilbert et al. (1983), the ionization instability is triggered only when the gas temperature is very high (as seen in Figs. 12). In a photoionized gas in thermal equilibrium the maximum temperature obtainable444In reality, the maximum temperature is less than this limit due to cooling by bremsstrahlung radiation. is the Compton temperature $`T_\mathrm{C}`$:

$$4kT_\mathrm{C}_{E_l}^{E_h}u_E𝑑E=_{E_l}^{E_h}u_EE𝑑E,$$
(2)
where $`u_E`$ is the spectral energy density of the incident radiation field which extends from $`E_l`$ to $`E_h`$. However, this equation overestimates $`T_\mathrm{C}`$ when $`E_h30`$ keV because it does not include the Klein-Nishina reduction in Compton heating at higher energies. For radiation in the form of a single power-law $`u_EE^{1\mathrm{\Gamma }}`$, making the first order Klein-Nishina correction (e.g., Ross et al., 1999) results in

$$kT_\mathrm{C}=\frac{(2\mathrm{\Gamma })E_l}{4(\epsilon ^{2\mathrm{\Gamma }}1)}\left[\frac{\epsilon ^{3\mathrm{\Gamma }}}{(3\mathrm{\Gamma })}\frac{21E_l\epsilon ^{4\mathrm{\Gamma }}}{5m_ec^2(4\mathrm{\Gamma })}\right],$$
(3)
where $`\epsilon E_h/E_l`$ and $`\mathrm{\Gamma }2.5`$. This equation can safely determine $`T_\mathrm{C}`$ for $`E_h<100`$ keV.

The dominant factors in determining $`T_\mathrm{C}`$ are $`\epsilon ^{3\mathrm{\Gamma }}`$ and $`\epsilon ^{4\mathrm{\Gamma }}`$, thus a larger value of $`\mathrm{\Gamma }`$ results in a lower Compton temperature. Alternatively, $`T_\mathrm{C}`$ rises if the irradiating spectrum is cut off at larger values of $`E_l`$ due to the reduction in inverse Compton cooling. For example, with $`E_h=100`$ keV and $`\mathrm{\Gamma }=2`$, $`kT_\mathrm{C}=1.3`$, 1.6 & 2.1 keV at $`E_l=1`$, 10 & 100 eV. This result and Fig. 2 implies that the lower-limit to the X-ray power-law would have to cut off in the hard X-rays, contrary to observations, in order to increase the temperature to high enough values to initiate the instability when $`\mathrm{\Gamma }=2`$.

Other effects not described by Eq. 3 can also change $`T_\mathrm{C}`$. As already noted in Sect. 3.1, $`T_\mathrm{C}`$ rises with an increase in the incidence angle (Nayakshin et al., 2000), although this can only suppress ionized Fe lines when $`\mathrm{\Gamma }2`$ if $`F_x/F_{\mathrm{disc}}`$$`1`$ where they are unlikely to occur for most conditions. Also, $`T_\mathrm{C}`$ can decrease if the disc flux increases due to, say, an increase in accretion rate (Nayakshin & Kallman, 2001). The rise in disc flux increases the amount of Compton cooling and is, in a sense, equivalent to steepening the incident spectrum.

Observations of a small sample of Seyfert 1 galaxies by BeppoSAX have detected high-energy cutoffs in the power-law continuum between 70 and 300 keV (Matt, 2000). Increasing the cutoff energy of the incident power-law will increase the Compton temperature at any value of $`\mathrm{\Gamma }`$ (possibly to the point of initiating the ionization instability; Nayakshin & Kazanas 2002), but it is difficult to estimate the value of $`T_\mathrm{C}`$ for $`E_h>100`$ keV as higher-order Klein-Nishina corrections need to be included. We are therefore placed in the interesting position where the strength of an object at 200 keV can affect its spectral features at 6–7 keV. Models of thermal Comptonization show that the cutoff energy ($`2kT_e`$ of the corona) decreases with larger $`\mathrm{\Gamma }`$, but the exact form of the relationship depends on the geometry (e.g., Svensson, 1996). If AGN with $`\mathrm{\Gamma }`$=1.8–2.1 are found to have spectral cutoffs $`<200`$ keV then it is likely that these objects will exhibit ionized Fe K$`\alpha `$ lines. However, it is clear that the calculations need to be extended to higher energies to make more general predictions, and we intend to do this in future work. It is also necessary to obtain more measurements of the overall spectral shape of Seyfert 1 galaxies.

### 4.2 Comparison with observations

In principle, it is possible to compare the observed long time-scale variations of the broad Fe K$`\alpha `$ line in AGN with the model predictions. There have been many recent studies on Fe K$`\alpha `$ variability based on ASCA or RXTE data (see references listed in Sect. 1 or, for a recent review, Reynolds, 2001). However, recent XMM-Newton and Chandra results have shown that many AGN have narrow Fe K$`\alpha `$ components in their spectra (Kaspi et al., 2001; Yaqoob et al., 2001; Reeves et al., 2001; Pounds et al., 2001). These lines are expected in AGN unification scenarios where reflection of X-rays may occur from the inner edge of a molecular torus (e.g., Krolik, Madau & Życki, 1994). It thus seems likely that narrow Fe K$`\alpha `$ lines may be common in Seyfert galaxies. These components would have been unresolvable in ASCA, and so the published variations in any broad line would be contaminated by this feature.

Nevertheless, it may be possible to compare our predictions to the observations of MCG–6-30-15. Iwasawa et al. (1999) reported that the broad Fe K$`\alpha `$ line profile of this Seyfert 1 galaxy shifted dramatically to the red during a flare event. There was no noticeable narrow component ($`EW<60`$ eV) remaining at 6.4 keV, suggesting that MCG–6-30-15 may harbour only a very weak narrow component. A subsequent long ($`400`$ ks) observation of MCG–6-30-15 by RXTE reported by Lee et al. (2000) has resulted in a comprehensive study of the Fe K$`\alpha `$ variability. Lee et al. (2000) split their data into 4 different continuum flux levels and found that the EW of the Fe K$`\alpha `$ line (as fitted by a Gaussian model) dropped as the flux increased, but that the line flux remained the same. However, if more complex continuum models that included reflection were used, the line flux decreases with increasing continuum flux. A re-analysis of this RXTE data by Vaughan & Edelson (2001) found that the line flux varied with time, but it was not correlated with the continuum. As with Lee et al. (2000), Vaughan & Edelson (2001) found that the photon index was correlated with the 2–10 keV flux. Finally, Reynolds (2000) used sophisticated numerical techniques on the RXTE data to search for a correlation between the Fe K$`\alpha `$ line and the continuum over a wider range of timescales than covered by Lee et al. (2000). The search was unsuccessful which suggests that the line flux is uncorrelated with the continuum down to timescales of 0.5 ks.

Figure 6 predicts that the Fe K$`\alpha `$ line flux should correlate with the flux striking the disc. Clearly, this does not agree with the observations. The EW was observed to drop with the continuum flux, which does agree with the predictions of Fig. 7. However, the continuum flux in MCG–6-30-15 is observed to correlate with the photon index, which has the opposite effect on the EW (Fig. 8). To give an idea of which effect might dominate the evolution of the Fe K$`\alpha `$ EW, we plot in Figure 9 how the EW varies with $`\mathrm{\Gamma }`$, but tying $`F_x/F_{\mathrm{disc}}`$ to the photon-index.

The solid points shows the EW variations when $`F_x/F_{\mathrm{disc}}`$$`=0.52`$ at $`\mathrm{\Gamma }=1.5`$, but with $`F_x/F_{\mathrm{disc}}`$ doubling for an increase in $`\mathrm{\Gamma }`$ of 0.1. In this case the line is mainly neutral and the EW initially drops due to the increase in $`\mathrm{\Gamma }`$ (recall that increasing $`F_x/F_{\mathrm{disc}}`$ does not change the EW by much when the line is dominated by a 6.4 keV line), before slowly increasing as the ionized component comes in. The decrease continues the longest when $`i=80`$ degrees. The open stars begin at $`F_x/F_{\mathrm{disc}}`$$`=4.16`$ when $`\mathrm{\Gamma }=1.5`$ and trace the ionized evolution. All four panels show the same decrease-increase behaviour with the absolute value decreasing with incidence angle. From these curves and the photon-index range of MCG–6-30-15 (1.8–2.1) we can conclude that the EW evolution of MCG–6-30-15 is consistent with that of a single purely neutral line, or of a hybrid line that is excited by radiation at a large incidence angle. However, this cannot explain the strange variation in line flux which may be indicative of more than one emission component. Indeed, both types of models can fit the data (Ballantyne & Fabian, 2001).

### 4.3 Neutral and ionized Fe K$`\alpha `$ lines

The calculations presented here predict that ionized Fe K$`\alpha `$ lines at $``$ 6.7 keV should be common in almost all cases if $`F_x/F_{\mathrm{disc}}`$$`1`$. Evidence for ionized lines have been seen in the X-ray spectra of some Narrow Line Seyfert 1 Galaxies (Comastri et al., 1998, 2001; Turner, George & Nandra, 1998; Vaughan et al., 1999; Ballantyne, Iwasawa & Fabian, 2001; Turner et al., 2001a, b) and high-luminosity Seyferts/low-luminosity quasars (Reeves et al., 2000; Pounds et al., 2001; Orr et al., 2001), both of which are thought to be high accretion rate objects and so might not be directly applicable to the current computations. Most Seyfert galaxies seem to exhibit lines of neutral Fe at 6.4 keV (e.g. Nandra et al., 1997a), although this may be a result of contamination from a narrow component. Taking these results at face value suggests that discs are either weakly illuminated (i.e., $`F_x/F_{\mathrm{disc}}`$$`<1`$) or much denser than previously considered.

If discs are weakly illuminated then it is difficult to reconcile this with the otherwise very successful theory of coronal heating by magnetic flares which necessarily predicts a highly irradiated disc. It is possible that measurements of the centroid energy of the Fe K$`\alpha `$ line can be confused by poor signal-to-noise data, calibration problems, or gravitational redshifting effects. Searching for EW variations with the photon index through relations shown in Fig. 8 may provide an alternative test for determining the ionization state of the Fe K$`\alpha `$ line.

Extensions to the basic magnetic flare model can also reconcile the models with the data. For example, Nayakshin et al. (2000), Nayakshin (2000) and Done & Nayakshin (2001) discuss the idea that a local wind will be driven by a magnetic flare event. These authors argue that as a result of this wind, surface material on the disc will be blown away revealing denser gas deeper into the disc. It is this gas that produces the observed reflection spectrum complete with a neutral Fe K$`\alpha `$ line. However, the effects (or even the existence) of such an X-ray driven wind has not been properly modelled, so, in light of the results of this paper, the following model is proposed as an alternative.

Consider the situation when separate magnetic flares are occurring over the accretion disc (see Figure 10 and Merloni & Fabian, 2001).

Directly underneath each illuminating event, the surface of the disc is highly ionized and therefore produces no Fe K$`\alpha `$ emission. As the angle of incidence of the radiation increases, the flux incident on the gas drops resulting in ionized Fe K$`\alpha `$ lines. As the incidence angle further increases, the local $`F_x/F_{\mathrm{disc}}`$ decreases until eventually neutral lines dominate. Such a gradient in ionization on the disc has also been considered by Reynolds (2000) for AGN and by Done & Życki (1999), Young et al. (2001) and Done & Nayakshin (2001) to explain the low reflection fractions observed in low-state Galactic Black Hole Candidates. Since there is much more surface area on the disc that is weakly illuminated and/or illuminated at grazing incidence, the majority of the line emission detected by a telescope will be at neutral energies. This effect should work for any value of $`\mathrm{\Gamma }`$, as the combination of low $`F_x/F_{\mathrm{disc}}`$ and large $`i`$ will result in a neutral line (Fig. 4). Similarly, the complete ionization below each flare can occur at any photon-index, but its assurance will depend on the physical and geometrical details of the event itself.

Spectral evolution of the illumination events will also cause changes in the ionization state of the Fe K$`\alpha `$ line. The recent RXTE observation of the black hole candidate GX 339-4 reported by Feng et al. (2001) showed that as the X-ray flux decreased from its low state to an “off state”, its Fe K$`\alpha `$ line evolved from neutral at 6.4 keV to ionized at 6.7 keV. More importantly, the photon index of the continuum evolved from about 1.5 to 2.2 as the flux decreased. As discussed earlier (Sect 4.1), when $`\mathrm{\Gamma }=1.5`$ the thermal ionization instability can operate and a neutral line is expected, but it is highly unlikely to work when $`\mathrm{\Gamma }=2.2`$, and a 6.7 keV line would be observed. Perhaps this is an observation of a decaying and expanding flare.

#### 4.3.1 X-ray Baldwin effect

Observations of many Seyfert galaxies and quasars by ASCA have shown that as the X-ray luminosity of the source increases, the EW of the Fe K$`\alpha `$ line decreases until it becomes undetectable at a 2–10 keV luminosity of 1046 erg s-1 (Iwasawa & Taniguchi, 1993; Nandra et al., 1997b; Reeves & Turner, 2000). This trend has been named the X-ray Baldwin effect in analogy to the decrease in EW of the C iv $`\lambda `$1550 line in the UV continuum of QSOs (Baldwin, 1977). In terms of the constant density models, this effect could easily be explained as the ionization parameter of the gas increasing to the point where the outer layers of the disc are fully ionized and therefore would not imprint any line emission on the reflection spectrum. The ionization parameter would increase if the higher luminosity in quasars was due to a higher accretion rate (Matt et al., 1993). Newer hydrostatic disc models also can explain the X-ray Baldwin effect by ionizing a thick skin on the surface of the accretion disk (Nayakshin, 2000; Życki & Różańska, 2001).

As discussed earlier, ionization gradients along the accretion disc by different flare locations may be able to account for the observed energies of the Fe K$`\alpha `$ line. If the power and number of flares above the disk increases as one moves from Seyferts to quasars then this would increase the X-ray luminosity, and also increase the amount of surface area that is subject to a highly ionizing flux (see Fig. 10). In such a situation only a very weak, highly ionized Fe K$`\alpha `$ line would be emitted from the surface. Indeed, recent XMM-Newton observations show that ionized lines are observed in higher-luminosity objects (Reeves et al., 2001; Pounds et al., 2001).

Such a scenario is not too different from that envisaged by the other models, but suggests a more complex geometry. There are a number of details to work out in order for such a model to placed on firm theoretical ground. Clearly, in a situation such as the one presented in the lower part of Fig. 10, feedback effects between the reflection spectrum and the coronal activity would be vitally important and may even assist in erasing spectral features through Comptonization (Petrucci et al., 2001).

## 5 Conclusions

The broad Fe K$`\alpha `$ line which (most likely) results from reflection from an accretion disc is the most important spectral signature in the X-ray study of AGN. The model calculations presented in this paper have shown that the behaviour and properties of the Fe line from an irradiated disc in equilibrium depend on how it is illuminated. Specifically, when the incident flux is larger than the disc flux, we predict that a He-like Fe K$`\alpha `$ line at 6.7 keV will be a common feature in the reflection spectrum. The thermal ionization instability is able to change this and produce a neutral line at 6.4 keV, but only when $`\mathrm{\Gamma }1.6`$ (for a cutoff energy of 100 keV) and the incident flux is not too high.

The Fe K$`\alpha `$ EW will remain roughly constant as $`F_x/F_{\mathrm{disc}}`$ is increased until the atmosphere becomes highly ionized and then it will decrease to small values. The rate of decline is lower for higher values of $`\mathrm{\Gamma }`$. When the line complex is dominated by He-like Fe, the EW will strongly increase with $`\mathrm{\Gamma }`$, but will slightly decline with the photon-index if it is neutral. This effect could be used as an observational diagnostic to determine the ionization state of the line, independent of the line centroid.

Long observations of Seyferts & quasars by XMM-Newton would be required to test the predictions of the models. However, any narrow component to the line profile must be taken into account before a comparison could be made. The results of such a study would be important, as it is clear that the ionization state of the disc depends greatly on how it is irradiated which would greatly constrain coronal models.

The prediction that ionized lines are common can be reconciled with their observational rarity by employing multiple ionization gradients on the accretion disc, which would be a natural result of a magnetically active and patchy corona. The challenge is now extended to modellers to compute observationally testable predictions for such complex scenarios.

## Acknowledgments

The authors thank A. Fabian, S. Vaughan, K. Iwasawa & A. Merloni for valuable discussions, and C. Berger for assistance with the models. DRB acknowledges financial support from the Commonwealth Scholarship and Fellowship Plan and the Natural Sciences and Engineering Research Council of Canada. RRR acknowledges support from the College of the Holy Cross.


format: text

text:
# Simulations of pulsar wind formation

## 1. Introduction

Most of the pulsar wind models assume a relativistic MHD-type flow with various degrees of symmetry (e.g., Kennel and Coroniti 1984, Begelman and Li 1992), and are reasonably successful in explaining the overall morphology and energetics of the surrounding plerion. However, they run into various difficulties under attempts to extrapolate the solution back towards the pulsar. The most famous of these is the problem of the magnetization of the wind ($`\sigma `$-paradox), which is the contradiction between the low ratio of magnetic to kinetic energy ($`\sigma 10^3`$) at the wind-nebula interface inferred from observations, expected strong magnetization near the pulsar ($`\sigma 1`$), and the conservation of $`\sigma `$ in an ideal MHD flow. Recent observations of the Crab and Vela with Chandra and HST have underscored the importance of understanding the origins of the pulsar wind. Observations of pulsar jets and orthogonal equatorial outflows suggest that rotation of the neutron star is imprinted in the wind structure, possibly at a very early stage. How and where acceleration and collimation of the wind occurs can be reliably studied only by constructing a wind model that begins from the pulsar itself. The main difficulty in the way of theoretical progress on this front is the lack of intuition about the phenomena occurring near a pulsar. For a general case of an oblique rotator, and, as we show in this paper, even for an aligned rotator, the problem is intrinsically three dimensional, nonaxisymmetric, and time-dependent, or in other words “too complicated”. With observations unlikely to resolve light cylinder scales ($`10^8`$ cm) any time soon, our only resort is numerical simulations. In this paper we describe our effort to simulate the formation of pulsar winds from the first principles. In particular, we study the behavior of plasma that is either emitted from or injected in the vicinity of a strongly magnetized rotating conducting sphere with arbitrary inclination angles between the magnetic and rotational axes. Despite its implied simplicity, this problem has very interesting solutions that may at first appear strange and unusual. Hence, in describing these preliminary results our emphasis is more on developing intuition, rather than on direct relevance to pulsars and observational implications. Such issues as well as more realistic models are currently being investigated. The paper is organized as follows: in § 2 we describe the particle-in-cell method used in our simulations and the setup of the problem, in § 3 and § 4 we describe numerical experiments with aligned and oblique rotators, and in § 5 we conclude and discuss further work.

## 2. Numerical method

We are interested in studying the outer magnetosphere of the pulsar – the region around the light cylinder ($`R_{LC}=c/\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ is the angular frequency of the star). The possible physical conditions in this region put stringent demands on a successful simulation method. Previous theoretical investigations suggest that the method should be able to handle relativistic flows, vacuum gaps, charge-separated configurations, space-charge limited flows and countersteaming flows, as well as MHD flow when it is applicable. For oblique rotators the light cylinder is a region of emission of large-amplitude electromagnetic waves, so strong wave-plasma interaction should be also addressable by the method. The requirement of simulating counterstreaming, i.e. flows with multivalued velocity, dispenses with the option of using fluid grid-based schemes. As our method we therefore chose a particle-based “Particle-In-Cell” (PIC) algorithm. It is a fully relativistic multidimensional method for self-consistent solution of Maxwell’s equations in the presence of plasma. The underlying philosophy of PIC is the following: a plasma is represented as a collection of macroparticles which carry charge and mass; the macroparticles are moved by leap-frog integration of the relativistic equations of motion with Lorentz force; the currents and charges associated with the macroparticles are deposited on an Eulerian grid, on which Maxwell’s equations are also discretized; the fields are updated (with particle currents and charges as sources) and extrapolated to the position of the macroparticles to contribute to the Lorentz force for the next time step.

The main advantage of using the PIC method for the pulsar wind problem is that one does not have to impose restrictions on plasma behavior in advance. The method can easily handle charge-separated flows and counterstreaming. It is intrinsically a kinetic model, so particle acceleration can be studied. The ideal MHD condition $`\stackrel{}{E}+(\stackrel{}{v}/c)\times \stackrel{}{B}=0`$ can be achieved in the appropriate regime, and has been observed in our tests. The disadvantage of using such a fundamental plasma model is that multiple spatial and time scales need to be resolved, which introduces constraints on the time step and spatial resolution (Birdsall and Langdon 1991).

We have modified the publically available 3D relativistic PIC code TRISTAN (Buneman 1993) to simulate plasma dynamics in the magnetosphere of a magnetized conducting rotating sphere. As is well known, electric fields are induced on the surface of such a body, and should serve as the boundary conditions for the field solve. In order to avoid the asymmetries introduced from the discretization of a rotating sphere on a Cartesian grid, we used the linearity of Maxwell’s equations to represent the field $`\stackrel{}{E}_{total}`$, which is used to advance the plasma particles, as a superposition: $`\stackrel{}{E}_{total}=\stackrel{}{E}_{vacuum}+\stackrel{}{E}_{plasma}`$, where $`\stackrel{}{E}_{vacuum}`$ is the field of the “vacuum rotator”, given by an analytic formula (Deutsch, 1955), and $`\stackrel{}{E}_{plasma}`$ is the “plasma” field computed from the currents deposited by macroparticles. The same decomposition is done for the magnetic field. We note that this is not a small amplitude expansion, this decomposition is exact, and two components can be of equal magnitude. A spherical surface at the center of the domain acts a source and sink of macroparticles. This surface can represent the neutron star or a rigidly corotating inner magnetosphere extending to a large fraction of the light cylinder (presently the radius of the sphere $`a>0.1R_{LC}`$). The outer walls of the computational domain in the simulation have radiation boundary condition for the fields and particles.

## 3. Aligned rotator

We begin by considering the case of magnetic axis being parallel to the rotation axis. There are no time-varying magnetic fields in this case, and the plasma in the magnetosphere will know that the star is rotating only through the electric fields of the unipolar inductor. Therefore, it is interesting to see whether the magnetosphere can fill and reach corotation given a source of plasma. It is natural to expect the plasma to be provided from the surface of the star, as the induced electric field is likely to overwhelm any conceivable work function. At the start of the simulation the rotating star has an induced quadrupolar surface charge and a central monopole giving rise to the following electric field components in spherical coordinates (e.g., Michel and Li 1999): $`E_r=\frac{2}{3}\varphi _0\frac{a}{r^2}+\varphi _0\frac{a^3}{r^4}(13\mathrm{cos}^2\theta )`$, $`E_\theta =\varphi _0\frac{a^3}{r^4}\mathrm{sin}2\theta `$, where $`\varphi _0=\mathrm{\Omega }B_0a^2/c`$, with $`B_0`$ magnetic field at the pole. Plasma is introduced as electron-positron pairs at the stellar surface at each time step. The surface field extracts one sign of charge while driving the other sign into the star where it is arrested. This is equivalent to the surface charge being torn off the star, eventually achieving a configuration where the surface charge due to the vacuum fields is equal and opposite to the surface charge induced by extracting the plasma particles, and the $`𝐄_{total}𝐁=0`$ is satisfied on the surface (strictly speaking, however, the charge left behind after extracting the plasma is not distributed over the star; in this sense the star in the simulation is a dielectric, rather than a conductor; see discussion in §5). The vacuum fields are such that electrons are extracted over the poles of the star and positrons in the equatorial region. This creates a charge-separated plasma with a vacuum region in between. The quadrupolar field of this plasma eventually counteracts the induced vacuum field at the surface, emission of new charges is slowed down and the magnetosphere reaches a quasi-stable disk-dome configuration as shown in Fig. 1. So far our results are in agreement with the previos axisymmetric simulations using a different method (Krause-Polstorff and Michel 1984 (KPM1) and 1985 (KPM2); Smith et al 2001 (SMT)). However, the 3D nature of our simulation allows us to address new physics.

The disk-dome configuration of Fig. 1 is not in rigid corotation with the star. The equatorial disk consists of charges of only one sign and has both velocity shear due to $`E\times B`$ drift in the total vacuum+plasma field and a density gradient. These conditions are favorable to the growth of the non-neutral plasma analog of the Kelvin-Helmholtz instability – the diocotron instability. This instability proceeds by amplification of the azimuthal perturbations in the charge density, and is illustrated in Fig. 1. The instability feeds on the energy of shear in the flow, and in the process can move the plasma across the field lines (hence the expanding radius of the disk in fig. 1). Indeed, an azimuthal perturbation in charge density has a nonzero $`E_\theta `$ component of the field, which results in a radial component of the $`E\times B`$ drift.

Such radial drifts represent the solution to the paradox of Holloway (1973), who noted that there are certain regions in the Goldreich-Julian (1969, GJ) solution that, if evacuated, would not be directly replenished from the star with the correct sign of charge. What probably happens in reality is that regions that are in direct contact with the star will lose corotation and transfer charge across the field lines due to the growing azimuthal perturbations in the charge density. They would then be replenished directly from the surface. Our simulations show similar behavior. The bottom panel of Fig. 1 displays the evolution of charge density in the equatorial plane as a function of radius. As the instability develops, the density tends to the Goldreich-Julian value. The agreement is not perfect, however: the GJ solution implies an infinite extent of the charge. In our case, the magnetosphere filled only up to the light cylinder, with occasional filaments of charge being torn off after moving past $`R_{LC}`$. In this sense, we did not get a self-consistent wind with the aligned model. Part of the reason is the absence of the global axial current that would significantly alter the magnetic topology. We believe this can be modified with the more realistic boundary conditions (see §5). However, even now it is clear that the GJ solution, corotation and the filling of the magnetosphere with plasma are likely an end result of a robust dynamical process, which is intrinsically three-dimensional in nature.

## 4. Oblique rotator

With the three dimensional problem set up as described above, studying oblique rotator amounts to just dialing the inclination angle in the formula for the Deutsch field. For the test cases considered here we also modified the plasma injection condition. At every step we injected cold uniform pair plasma within a sphere of radius $`R_{LC}`$, in a deliberately crude imitation of the pair formation process. If the injected plasma density is larger than $`\rho _{GJ}`$, the plasma has a significant dynamical effect on the field, and after an initial transient the system achieves a steady state (number of particles in the domain is constant). The flow consists of two main regions: the closed, quasi-corotating magnetosphere and the outflow (Fig. 2a). The first is the region of the size of $`R_{LC}`$ encompassing the closed magnetic field lines. Here the initially neutral plasma is separated with excess positive charge around the equator and excess negative charge around the poles, similar to the GJ solution. The near-corotation is accomplished due to this charge separation and due to the induced electric field that arises from time variation of the magnetic field. Near the light cylinder the plasma is accelerated due to the radiation pressure of the large amplitude electromagnetic wave and forms an outflow. An interesting feature of this wind is that it is formed at the edge of the closed zone, and is mostly in the sector around the rotational equator, where it has the characteristic spiral wave pattern (Fig. 2b). Other regions of outflow are also possible. In Fig. 2a we introduced an additional plasma source in the polar flux tube for a $`20^{}`$ rotator. The plasma acts to cancel the parallel electric fields and moves along the magnetic field tracing a cone around the rotation axis.

## 5. Conclusion

Our findings can be summarized as follows: pulsar magnetospheres can be filled by surface emission with the help of a 3D diocotron instability; simple models of oblique rotators can form relativistic winds; such winds have variation in properties with latitude, with equatorial flow being formed at the edge of the closed magnetosphere, and polar cap flow probably contributing in the direction of rotational axis. Much work remains to be done to quantify these qualitative results and to understand the dependence of the wind geometry, acceleration and collimation on the parameter regime. While tempting, it would be too early to suggest that the equatorial and polar outflows observed in these simulations are directly responsible for features in the Crab and Vela. The models described in this paper are self-consistent in the sense that the field and particle equations are solved simultaneously and correctly. However, the models are not currently self-consistent with respect to the boundary conditions: the induced quadrupole electric fields are those from a rotating magnetized conductor, but as far as the plasma is concerned the central body is a dielectric with $`ϵ=1`$. Both KPM1, KPM2 and SMT simulations suffer from the same problem as well. This restriction does not allow currents to close inside the star, and instead leads to charging of the central body at the locations where the current is extracted and returned. Consequently, any current flow in the system is temporary. We are currently working on overcoming this deficiency. The more realistic models with true central conductors will be described in detail in an upcoming publication.

## References

Begelman, M. C. and Li, Z.-Y. 1992, ApJ, 397, 187

Birdsall, C. K. and Langdon, A. B. 1991, “Plasma Physics via Computer Simulation” (Adam Hilger: Bristol)

Buneman, O. 1993, in Computer Space Plasma Physics, ed. H. Matsumoto & Y. Okura (Terra Scientific), 67

Goldreich, P. and Julian, W. H. 1969, ApJ, 160, 971

Deutsch, A. 1955, Annales d’Astrophysique, 18, 1

Holloway, N. J., 1973, Nature Phys. Sci., 246, 6

Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 694

Krause-Polstorff, J. and Michel, F. C. 1984, MNRAS, 213, 43 (KPM1)

Krause-Polstorff, J. and Michel, F. C. 1985, A&A, 144, 72 (KPM2)

Michel, F. C. and Li, H. 1999, Phys. Reports, 318, 227

Smith, I. A., Michel, F. C., Thacker, P. D. 2001, MNRAS, 322, 209 (SMT)


text:
# Heterogeneous nucleation near a metastable vapour-liquid transition: the effect of wetting transitions

## 1 Introduction

When water is cooled below $`0^{}`$C at atmospheric pressure it freezes, it turns into ice. This conversion of one phase, water, into another, ice, starts with the nucleation of a microscopic nucleus of ice. This nucleus consists of only of order 10 molecules, its formation costs free energy and occurs not in the bulk of the water but at a surface in contact with the water. The free energy cost provides a barrier to the nucleation of ice. If the free energy cost or barrier is large it will limit the rate of crystallisation. When the barrier is very large the phase which is not the equilibrium one, for example water below $`0^{}`$C, will persist for very long times. The fluid phase is then called metastable. So to determine whether a phase which is not an equilibrium phase is metastable or whether the equilibrium phase nucleates rapidly we need to calculate the free energy barrier to nucleation. We do this here near to and at another phase transition, a phase transition between two phases neither of which is the true equilibrium phase. The process of nucleation at a surface is called heterogeneous nucleation to distinguish it from homogeneous nucleation which occurs in the bulk. See Ref.  for an introduction to nucleation.

So, here we study the rate of heterogeneous nucleation of one phase transition, a fluid-crystal phase transition, near a second phase transition — a phase transition between phases which are both meta- or unstable with respect to crystallisation. This second phase transition is a vapour-liquid transition. At phase transitions the thermodynamic functions, including interfacial tensions, exhibit singular behaviour which is universal in the sense that many different systems show behaviour which is identical up to a few scale factors. Here we show that at the vapour-liquid transition the free energy barrier to heterogeneous nucleation of a crystal shows behaviour which although not truly universal is the same near all vapour-liquid transitions, up to a few scale factors, assuming that the longest-range interactions are dispersion forces. This is due to wetting: the formation of a thick layer of the liquid phase on a surface in contact with the vapour phase . The layer forms as coexistence is approached and causes a drop in the nucleation barrier to the nucleation of a dense phase such as a crystalline phase. The dependence on the nature of longest-range forces makes our findings not-quite universal (unlike homogeneous nucleation near a bulk critical point  which is universal).

Motivation for our study is provided by the fact that some globular proteins have the correct set of phase transitions to observe heterogeneous nucleation of a crystal near a vapour-liquid transition. See Refs.  for the phase diagrams of a number of globular proteins. The crystallisation of globular proteins is of great interest because protein crystals are required in order to determine the all-important three-dimensional structure of a protein . Also, although we will always refer to the nucleus as being crystalline and the other transition as being a vapour-liquid transition, our findings are much more general. They refer to the nucleation of any noncritical phase near another, Ising-type, phase transition. We simply describe the phases as crystal, vapour and liquid for simplicity and because having definite phases in mind is useful for pedagogical purposes.

Earlier work has found universal behaviour of the nucleation barrier for homogeneous nucleation near a bulk critical point, see Refs. . This earlier work, in particular that of ten Wolde and Frenkel , inspired this study of heterogeneous nucleation and the results of Refs.  are in a sense the homogeneous nucleation analogues of the results we will obtain for heterogeneous nucleation. Also, Talanquer and Oxtoby have studied heterogeneous nucleation of a liquid from a vapour phase . They studied the nucleation of the liquid phase at a surface when the liquid phase itself is close to wetting this surface. So, although they studied heterogeneous nucleation in the vicinity of a wetting transition, as we do below, they studied the nucleation of the liquid phase, the phase which is doing the wetting, whereas here we study the heterogeneous nucleation of another phase, the crystalline phase.

In the next section, we introduce both our model of the process of heterogeneous nucleation and the phase behaviour of the systems we are interested in. In section 3 we derive the variation in the rate of heterogeneous nucleation near wetting transitions. Then in section 4 we compare with homogeneous nucleation of a crystal near a bulk critical point. The last section is a conclusion.

## 2 Heterogeneous nucleation

Heterogeneous nucleation is an activated process  and as such occurs at a rate which decreases exponentially with the height of the barrier $`\mathrm{\Delta }F`$, which must be overcome. If $`N_n`$ is the number of nuclei per unit area crossing the barrier per unit time then $`N_n`$ is given by an expression of the form

$$N_n=\sigma \tau ^1\mathrm{exp}(\mathrm{\Delta }F/kT),$$
(1)
where $`\sigma `$ is a surface density, i.e., it has dimensions of inverse area, and $`\tau `$ is a characteristic time. We will refer to $`N_n`$ as the nucleation rate or heterogeneous nucleation rate. The surface is smooth, perfectly planar and chemically homogeneous. Equation (1) comes from the fact that the nucleus is a large, i.e., improbable, fluctuation. As a fluctuation its probability of occurring in unit area is $`\sigma \mathrm{exp}(\mathrm{\Delta }F/kT)`$. The rate at which these fluctuations cross the barrier is then estimated as the number of fluctuations divided by $`\tau `$ which is an estimate of how long it takes the nucleus to acquire one or a few extra molecules, enough for the nucleus to be big enough to grow irreversibly into a crystallite. Equation (1) is far from rigorous but has been found for homogeneous nucleation to be a reasonable estimate. See the book of Debenedetti  for a discussion and  for a detailed comparison of an expression of the form of Eq. (1) with the results of computer simulation (for homogeneous nucleation). For the remainder of this work we will assume that $`\tau `$ and $`\sigma `$ vary weakly with temperature and chemical potential and so the variation of the rate of heterogeneous nucleation is dominated by the variation in the free energy barrier $`\mathrm{\Delta }F`$. Work on homogeneous nucleation has shown this assumption to be most often justified, except near a glass transition where $`\tau `$ increases sharply. However, verifying it requires a detailed calculation for a specific model system, which we do not do here.

Consider the phase diagram in the chemical-potential–temperature plane shown in Fig. 1. It is the phase diagram of a simple model of a globular protein, calculated using an approximate theory. See Ref.  for the precise definition of the model; the model parameters have the same values as they do for Fig. 4 of that reference. Also, see this reference for the same phase diagram in the density-temperature plane which can be seen to be qualitatively like that of a number of globular proteins . At true equilibrium there is only one phase transition, from a fluid phase, the equilibrium phase below the solid curve, to a crystal, the equilibrium phase above the curve. However, if the barrier to formation of the crystal phase is high then the chemical potential can be increased at constant temperature along a path such as that indicated by the arrow in Fig. 1, until the metastable fluid undergoes another phase transition: a transition from a vapour phase to a liquid phase. The vapour is below the dashed curve, the liquid above. We have used the phase diagram of a model protein because it has the correct form: it has a vapour-liquid transition near where we expect the nucleation of the crystalline phase to occur. However, the behaviour we find will apply whenever heterogeneous nucleation occurs at a surface which passes through a wetting transition.

Figure 2 is a schematic of a crystalline nucleus in contact with a surface and immersed in a liquid layer of thickness $`l`$, with a vapour phase on top. This is the situation of interest, there is a bulk vapour phase against a smooth surface which attracts the molecules causing a layer of liquid to form near and at vapour-liquid coexistence. At coexistence and if we are above the wetting temperature, see Refs.  and section 3, there is a very thick layer of liquid covering the surface and separating the vapour phase from the surface. The thickness $`l`$ is then limited only by gravity. This wetting layer will reduce the nucleation barrier greatly if the interfacial tension between the liquid and the nucleus, $`\gamma _{xl}`$ is lower than that between the vapour and the nucleus, $`\gamma _{xv}`$. If the surface area of the nucleus not in contact with the surface is $`S`$, then the free energy reduction will be $`S(\gamma _{xl}\gamma _{xv})`$ when the wetting layer forms. For $`S`$ of order 10 times the area per molecule of the surface of a crystal, and the difference $`\gamma _{xl}\gamma _{xv}`$ of order $`kT`$ divided by the area of one molecule, the reduction in the surface contribution to the free energy barrier is of order $`10kT`$. This reduction will occur on moving in the vapour phase from conditions of chemical potential and temperature far from vapour-liquid coexistence, where there is no wetting layer, to at or very near coexistence. A large reduction which will lead, Eq. (1), to a very large increase in the nucleation rate $`N_N`$, which should be easily large enough to observe in an experiment.

## 3 Variation in the rate of heterogeneous nucleation as a wetting transition is approached

In this section we start with the assumption that the variation of the rate is dominated by that in the free-energy barrier $`\mathrm{\Delta }F`$ and then calculate how $`\mathrm{\Delta }F`$ varies near wetting transitions of different types. Near a wetting transition the qualitative behaviour, in particular the form of the singularities, can be determined without knowing any specific details of the interactions or of the phase which is nucleating, we only require that the longest range interactions be dispersion forces . This is true as the singularities come from long lengthscale phenomena for which the small lengthscale chemical details are irrelevant.

There are a number of different wetting phase transitions, see for example the excellent review of Schick . We will deal with the three most common, starting with complete wetting. In each case we will work very close to the wetting transition, temperatures or chemical potentials very close to their values at the transition. We will determine the leading order singular terms in the temperature or chemical potential variation of the nucleation rate near the transition.

### 3.1 Complete wetting

As heterogeneous nucleation occurs at a surface in contact with the fluid it is rather obviously controlled by what happens at the surface. Now, if the surface attracts the molecules of the fluid then not-too-far from the critical point , we will have wetting. Wetting is where at vapour-liquid coexistence a thick layer of the liquid phase interposes itself between the surface and the vapour phase, replacing the surface-vapour interface by a surface-liquid interface plus a liquid-vapour interface. It is also called complete wetting . At coexistence the thickness of the layer is generally limited only by gravity. Thus, our first result is the, rather obvious, fact that if the surface is wet by a thick layer of liquid the rate of heterogeneous nucleation will be the same at surfaces in contact with the two coexisting phases, the vapour and the liquid. This result is universal, it is true whenever the surface is wet. The free energy cost of forming a nucleus at the surface in contact with the bulk vapour is the same as the cost of forming a nucleus at the surface in contact with the bulk liquid, as both surfaces are covered by liquid.

This is at coexistence, the chemical potential $`\mu =\mu _c`$, where $`\mu _c`$ is the chemical potential at vapour-liquid coexistence at this temperature. In the vapour phase off coexistence, $`\mu <\mu _c`$, there is still a layer of liquid at the surface provided $`h=\mu \mu _c`$ is not too negative. But this layer thins as $`h`$ becomes more negative as we move farther from coexistence. How it thins depends on the nature of the longest-range forces present in the system . Here we assume that the longest range forces are dispersion forces, which is most often the case. Then the free energy per unit area, $`f`$, of a film of liquid between the surface and the vapour of thickness $`l`$ is

$$f=\frac{a}{l^2}h\delta \rho l,$$
(2)
where we have included only those parts which depend on the thickness of the layer. The first term is the contribution of the dispersion forces, $`a`$ is a positive constant and the $`l^2`$ dependence comes from starting with the usual $`l^6`$ dependence, integrating over a volume, for the solid substrate, and then over the thickness of the layer of liquid. This fourfold integration changes the $`l^6`$ dependence to $`l^2`$ . The second term is the increase in free energy due to the formation of a layer of liquid of thickness $`l`$ when the chemical potential of the liquid, $`\mu _c`$, is higher than the chemical potential. In this term $`\delta \rho `$ is the difference in number density between the liquid and vapour phases. Minimising the free energy Eq. (2) we obtain the thickness

$$l=\left(2a/(h\delta \rho )\right)^{1/3},$$
(3)
which diverges at coexistence: this divergence is in practice cutoff at some large thickness by gravity. The divergence at coexistence is a genuine phase transition: $`l^1`$ is analogous to an order parameter and the exponent $`1/3`$ in Eq. (3) is a critical exponent .

So, what is the free energy barrier to heterogeneous nucleation off coexistence: $`h<0`$ but small? We know that at coexistence the liquid layer is very thick and so the free energy of a nucleus in the vapour phase is the same as that in the liquid, call it $`\mathrm{\Delta }F_L(\mu ,T)`$. As we move off coexistence the liquid layer thins and the nucleus will notice this because it will interact with the vapour phase once $`l`$ is not too large. See Fig. 2 for a schematic of a nucleus at a substrate in a liquid layer. For $`l`$ larger than the radius of the nucleus, the nucleus will interact with the vapour as a point object. The interaction of a small object with the vapour across the liquid layer varies as $`l^3`$ and is proportional to the volume of the nucleus, $`v_n`$ . So the nucleus-vapour interaction increases the nucleation barrier to $`\mathrm{\Delta }F`$, where $`\mathrm{\Delta }F`$ is given by

$$\mathrm{\Delta }F(\mu ,T)=\mathrm{\Delta }F_L(\mu ,T)+\frac{Av_n}{l^3}l\text{large},$$
(4)
where $`A`$ is a coefficient for the interaction of the nucleus with the vapour across a slab of liquid. A dense nucleus will generally repel a dilute vapour and so $`A`$ will then be positive . Here, we focus on nucleation of a crystal phase and crystalline phases are rather dense. Also, note that $`\mathrm{\Delta }F_L(\mu ,T)`$ is a function of $`\mu `$ and we will use it for $`\mu <\mu _c`$ where the vapour phase is more stable than the liquid. We assume that we can continue $`\mathrm{\Delta }F_L(\mu ,T)`$ into the region where the liquid is metastable with respect to the vapour, $`h`$ negative but small, and that $`\mathrm{\Delta }F_L(\mu ,T)`$ is analytic at $`h=0`$, at coexistence. Using Eq. (3) for the thickness we obtain

$$\mathrm{\Delta }F(\mu ,T)=\mathrm{\Delta }F_L(\mu ,T)h\frac{A\delta \rho v_n}{2a}h<0,|h|\text{small},$$
(5)
which implies that $`\mathrm{\Delta }F\mathrm{\Delta }F_L`$ as $`h`$ is negative in the vapour. For small $`h`$ the variation in $`\mathrm{\Delta }F`$ is linear — the exponent is 1. This result holds for any system at fluid-fluid coexistence with dispersion forces near a complete wetting transition.

Assuming Eq. (1) holds for the rate and using Eq. (5) for $`\mathrm{\Delta }F`$ we see that if $`\mathrm{\Delta }F_L(\mu ,T)`$ and $`\sigma \tau ^1`$ vary smoothly with through $`\mu `$ around $`\mu _c`$ that near coexistence $`\mathrm{\Delta }F`$ has the following form

$$N_n=\{\begin{array}{cc}N_{nL}+\left[N_{nL}^{}+N_{nL}\frac{A\delta \rho v_n}{2a}\right](\mu \mu _c)\hfill & \mu <\mu _c\hfill \\ N_{nL}+N_{nL}^{}(\mu \mu _c)\hfill & \mu >\mu _c\hfill \end{array},$$
(6)
where $`N_{nL}`$ is the nucleation rate in the liquid at coexistence, $`\mu =\mu _c`$, and $`N_{nL}^{}`$ is the derivative of the nucleation rate, with respect to the chemical potential in the liquid, at $`\mu =\mu _c`$. Equation (6) states that the first derivative of the nucleation rate is discontinuous at the vapour-liquid transition because this derivative contains a contribution from the thickening wetting layer on the vapour side of the coexistence curve but not on the liquid side. Figure 3 is a schematic of the variation of the rate of heterogeneous nucleation near the vapour-liquid transition (the solid curve).

### 3.2 Critical wetting

For complete wetting, the thickness $`l`$ of the layer of liquid diverges as coexistence is approached, as in Eq. (3). Now, on moving away from the critical point of the vapour-liquid transition, if the attraction of the surface for the molecules is not too strong, then the surface-vapour interface may cease to be wet by the liquid phase. Then as coexistence is approached the thickness of the layer of liquid between the surface and the vapour does not diverge, it remains finite. This is called partial wetting . The transition from complete wetting to partial wetting is a phase transition. It can be continuous, called a critical wetting transition, or it can be first order. We will deal with each in turn.

First critical wetting. This occurs when the coefficient of the $`l^2`$ term in the free energy per unit area, $`f`$, changes sign. To deal with this we need the next order term in an expansion in $`l^1`$. This is the $`l^3`$ term and adding such a term to the free energy of Eq. (2) we have

$$f=\frac{t}{l^2}+\frac{b}{l^3}hl,$$
(7)
where $`t`$ is a measure of the distance from the critical wetting transition which occurs at $`t=0`$, and $`b`$ is a positive constant. Above the transition, where there is complete wetting, as both terms are of the same sign we can neglect the $`l^3`$ as being small for thick layers, then with $`t`$ fixed we recover Eq. (2). The transition occurs at coexistence so $`h=0`$. Then minimising $`f`$ we have

$$l=\{\begin{array}{cc}\mathrm{}\hfill & t0\hfill \\ 3b/(2t)\hfill & t<0,\hfill \end{array}$$
(8)
$`l`$ varies as $`(t)^1`$ near and below the transition. The leading order interaction of the nucleus with the vapour phase is still given by the second term in Eq. (4) when the layer is thick, i.e., near the transition. Thus, using Eq. (8) for $`l`$ we have

$$\mathrm{\Delta }F(T)=\{\begin{array}{cc}\mathrm{\Delta }F_L(T)\hfill & t0\hfill \\ \mathrm{\Delta }F_L(T)\left[8Av_n/(27b^3)\right]t^3\hfill & t<0,\hfill \end{array}$$
(9)
the difference between the nucleation barrier in the vapour phase and that in the liquid varies as $`t^3`$ ($`t`$ small) below the transition, they are the same above it. Again this holds for any system at fluid-fluid coexistence with dispersion forces near a critical wetting transition. The exponent of three is rather large, it means that $`\mathrm{\Delta }F`$ and its first, second and third derivatives are all continuous at the transition. The singularity in $`\mathrm{\Delta }F`$ is very weak and so detecting its effect on the variation in the nucleation rate $`N_n`$ near the wetting transition in an experiment may be very difficult.

Below the wetting transition, $`t<0`$, the surface is partially wet: either covered by a thin film of molecules or with only a few molecules on the surface. Then as the coexistence curve is crossed at constant $`t<0`$ we go from a nucleus on a surface which is essentially in direct contact with the vapour (with at most a thin film between them) to a nucleus on a surface in contact with the liquid phase. The rate of heterogeneous nucleation then has a discontinuity at coexistence, $`\mu =\mu _c`$. It will jump upwards, as shown schematically in Fig. 3 (the dashed curve).

### 3.3 First-order wetting

Now for a first-order partial-wetting-to-complete-wetting transition. This occurs when the $`l^3`$ term is negative, $`b<0`$. Then for stability we require a $`l^4`$ term. This is just as in a Landau expansion for the free energy near a phase transition. Adding a $`l^4`$ term with a positive coefficient $`c`$ to the $`f`$ of Eq. (7) we have

$$f=\left(\frac{b^2}{4c}+t\right)\frac{1}{l^2}+\frac{b}{l^3}+\frac{c}{l^4}hl,$$
(10)
where in order to keep the transition at $`t=0`$ we have added a constant, $`b^2/(4c)`$ to the coefficient of $`l^2`$. At coexistence $`h=0`$ and we minimise to obtain $`l`$

$$l=\{\begin{array}{cc}\mathrm{}\hfill & t0\hfill \\ 8c/\left[\left(b^232ct\right)^{1/2}3b\right]\hfill & t0,\hfill \end{array}$$
(11)
giving a jump from $`l=2c/b`$ to $`\mathrm{}`$ at the transition. Putting this jump in Eq. (4) for the free energy of the nucleus we have a jump in the free energy barrier of $`Av_n/(2c/b)^3`$. Near the transition the barrier varies as

$$\mathrm{\Delta }F(T)=\{\begin{array}{cc}\mathrm{\Delta }F_L(T)\hfill & t0\hfill \\ \mathrm{\Delta }F_L(T)+\left[Av_n/(2c/b)^3\right]\left[112ct/b^2+O(t^2)\right]\hfill & t<0,\hfill \end{array}$$
(12)
Above the transition, $`t>0`$, the free energy barriers are the same in the vapour and liquid phases while just below the transition the difference between the two is $`Av_n/(2c/b)^3`$. The jump in the free-energy barrier to heterogeneous nucleation will cause a jump in the nucleation rate $`N_n`$, from Eq. (1). As coexisting vapour and liquid phases are cooled the rate of heterogeneous nucleation at surfaces in contact with the vapour phase will jump downwards when the first-order wetting transition is crossed. Assuming $`\sigma `$ and $`\tau `$ vary smoothly through the transition the ratio of the nucleation rate $`N_n`$ just above the wetting transition to that just below it is $`\mathrm{exp}[Av_n/(2c/b)^3]`$.

In general we expect that the deeper we are into the region where the crystal is the equilibrium phase, the larger is the nucleation rate. Thus we expect that if we cool coexisting vapour and liquid phases, that the rate of heterogeneous nucleation will increase in both, see Fig. 1. However, we have just shown that if there is a first-order wetting transition that the rate of heterogeneous nucleation will jump downwards as we cross this transition. Potentially at least, the nucleation rate may not be a monotonic function of temperature: in the vapour phase it may increase and then jump downwards as the transition is crossed. A non monotonic variation in $`\mathrm{\Delta }F`$ is very rare, we are aware of only one example . The non-monotonic variation in $`\mathrm{\Delta }F`$ found by Auer and Frenkel  is (presumably) not due the presence of another phase transition.

When the wetting transition at coexistence is first order, a prewetting transition branches off from the coexistence curve into the vapour off coexistence . This prewetting transition does not go far into the vapour phase, it ends at a critical point which is at a value of $`h`$ which depends on $`b`$ and $`c`$ but is always small. At the prewetting transition there is a jump in the value of $`l`$, this jump decreases as $`h`$ decreases until the jump reaches zero at the prewetting critical point. At the prewetting transition as $`l`$ jumps then so does the barrier to nucleation, from Eq. (4). This transition including the critical point may be calculated from the free energy Eq. (10). As this free energy is analytic it yields mean-field exponents for the critical point . Thus at the temperature of the prewetting critical point, $`t_{cp}`$, and near the critical point the thickness difference $`ll_{cp}\text{sgn}(hh_{cp})|hh_{cp}|^{1/3}`$, where $`l_{cp}`$ and $`h_{cp}`$ are the liquid layer thickness and value of $`h`$ at the critical point. This corresponds to the critical exponent $`\delta =3`$ – its mean-field value. Putting this variation of $`l`$ into our expression for the interaction of the nucleus with the vapour phase across the liquid layer, we obtain

$$\mathrm{\Delta }F\mathrm{\Delta }F_{cp}\text{sgn}(hh_{cp})|hh_{cp}|^{1/3}\text{mean-field},$$
(13)
the nucleation barrier varies with the chemical potential minus that at the critical point to the one third power. $`\mathrm{\Delta }F_{cp}`$ is the free energy of the nucleus at the prewetting critical point, and Eq. (13) holds for $`t=t_{cp}`$ and $`|hh_{cp}|`$ small.

Finally, we note that unlike complete and critical wetting, a first-order wetting transition proceeds via nucleation and growth, see Refs. . So for example on cooling below the transition, a metastable thick wetting layer may persist, where by metastable we mean the that the thickness of layer is not the thickness which occurs at the absolute minimum of the free energy Eq. (10). Note that this layer is then doubly metastable, its free energy is higher than that of a thinner film and the system with either of these two layer thicknesses has of course a higher free energy than at true equilibrium where there is dilute-fluid–crystal coexistence. Above we have, for simplicity, neglected the time taken to reach equilibrium thickness and assumed that the thickness of the layer is always that at the minimum in the free energy Eq. (10). In reality the transition from a thick to thin layer will follow some set of dynamics which will complicate the analysis. See Refs.  and references therein for the dynamics at first-order wetting transitions. See Ref.  for an explicit study of the analogous problem of homogeneous nucleation near another transition where nucleation of the equilibrium and a metastable phase compete.

## 4 Comparison with homogeneous nucleation near a bulk critical point

Essentially by definition, phase transitions are where the thermodynamic functions of an equilibrium system have singularities. In earlier work  we showed that the rate of homogeneous nucleation has a singularity at a Ising-type phase transition in the bulk and here we have shown that the rate of heterogeneous nucleation has a singularity at surface phase transitions associated with an Ising-type phase transition in the bulk. The presence of a singularity at a phase transition is common to both homogeneous and heterogeneous nucleation.

Our findings here for heterogeneous nucleation near a prewetting critical point and our earlier findings for homogeneous nucleation are particularly closely related. In both cases we have a nucleus which is a small (point-like) perturbation which couples to the order parameter of the transition . The order parameter is the density for the bulk transition and the thickness in the prewetting surface transition. In earlier work  we used scaling arguments to obtain the correct exponents for homogeneous nucleation near the critical point of an Ising-type transition. In heterogeneous nucleation the nucleus will couple to the order parameter of the prewetting transition, which is an Ising-type transition in two dimensions. In Ref.  we showed that near an Ising-type critical point the behaviour is fixed and universal providing only that the nucleus couples to the order parameter. Thus heterogeneous nucleation near a prewetting critical point is completely analogous to homogeneous nucleation in a two-dimensional system near a bulk critical point. We can apply the scaling approach of Ref.  to heterogeneous nucleation near a prewetting critical point. We then obtain the correct and universal exponents for the variation of the free-energy barrier near the critical point. For example along the prewetting critical isotherm the free-energy barrier scales with distance to the critical point as

$$\mathrm{\Delta }F\mathrm{\Delta }F_{cp}\text{sgn}(hh_{cp})|hh_{cp}|^{1/15}\text{universal},$$
(14)
where the exponent, which is $`1/\delta `$, is obtained from the exact value $`\delta =15`$ for the Ising model in two dimensions. In two dimensions the mean-field predictions for the critical exponents like $`\delta `$ are very poor: the mean-field prediction, Eq. (13), has an exponent which is five times too large . See Ref.  for definitions of the critical exponents. Equation (14) is just Eq. (14) of Ref.  in two dimensions. See that reference for a derivation.

## 5 Conclusion

Almost invariably, a first-order phase transformation starts with heterogeneous nucleation. The nucleus of the new phase forms at a surface, see Fig. 2. Thus, the free energy barrier to the formation of the nucleus and therefore the rate of nucleation in a phase depend sensitively on anything which happens at the interface between the surface and the phase. If the phase in contact with the surface is a vapour phase close to a second, vapour-liquid, phase transition, then if the surface attracts the molecules a wetting layer may form at the surface. This is a layer of liquid at the surface, separating the surface from the vapour. The wetting layer will reduce the nucleation barrier greatly if the interfacial tension between the liquid and the nucleus, $`\gamma _{xl}`$ is lower than that between the vapour and the nucleus, $`\gamma _{xv}`$. We estimated in section 2 that the reduction in the surface contribution to the free energy barrier is of order $`10kT`$. A large reduction which should be easily large enough to observe in an experiment. The formation of a wetting layer, either as coexistence is approached (complete wetting), or along the coexistence curve (critical wetting or a first-order wetting transition) is a phase transition. We found that at complete wetting the derivative of the barrier as a function of chemical potential was discontinuous while the barrier itself has a discontinuity as the temperature is varied through a first-order wetting transition. Thus the rate of heterogeneous nucleation has a discontinuity in its slope as the coexistence curve is crossed at constant temperature, when there is complete wetting. The rate of change of the nucleation rate in the vapour phase just below coexistence is not the same as its rate of change in the liquid phase just above coexistence. Above and below mean at values of the chemical potential above and below that at coexistence. The rate of heterogeneous nucleation has a discontinuity as a first-order wetting transition is crossed.

Our model system is highly idealised, the surface is assumed perfectly smooth and homogeneous, and the dynamics of the formation of wetting layers have been neglected: the thickness was always taken to be at equilibrium. Future work should address how the dynamics of formation of wetting layers can effect nucleation; near a first-order wetting transition this will presumably be analogous to homogeneous nucleation near a metastable first-order bulk transition . Also, an understanding of the effects of chemical heterogeneity and of curvature would be useful as in practice surfaces will not be perfectly homogeneous or smooth and heterogeneous nucleation can occur on particles whose surfaces are inherently curved. But the most urgent requirement is for experiments on heterogeneous nucleation on simple, well characterised surfaces.

It is a pleasure to acknowledge helpful discussions with A. Parry. Work supported by EPSRC (GR/N36981).