text:
# I Introduction

## I Introduction

The successful start of the RHIC spin program marks a new era in spin physics. Very inelastic $`pp`$ collisions at high energies open up unequaled possibilities to thoroughly address many interesting questions, most importantly perhaps those concerning the still unknown role of polarized gluons . QCD analyses of upcoming high-$`p_T`$ prompt photon, jet, heavy flavor, inclusive-hadron, and $`W^\pm `$ boson production data will greatly deepen and enrich the understanding of the nucleon and of QCD that we had gained previously from polarized deeply-inelastic scattering (DIS) of leptons off nucleons.

Investigations of interactions between polarized leptons and nucleons will, however, continue to play a vital role in spin physics. Besides the many observables of DIS, photoproduction processes are particularly interesting. They give access to the parton content of circularly polarized photons, due to the presence of contributions to the cross section where the photon “resolves” into its parton content, for instance, by fluctuating into a vector meson of the same quantum numbers, prior to the hard QCD interaction. Such contributions compete with the “direct” part, for which the photon simply interacts as an elementary pointlike particle.

It is well-known  that at lower (fixed-target) energies, the direct component tends to dominate over the resolved one. Only at collider energies do resolved contributions become important. Extensive studies in the unpolarized case at the HERA and LEP colliders have firmly established the concept of photonic parton densities  and shown that resolved photon processes contribute a sizable fraction of photoproduction events in certain kinematic regions. The fact that $`e^+e^{}`$ and $`ep`$ data can be universally described by the same set of photonic parton densities marks an important test and success of QCD.

It has not yet been possible to observe resolved contributions for polarized photons, due to the fact that only fixed-target lepton-nucleon experiments have been performed so far in the polarized case. In this context, an exciting possibility would be to have a polarized lepton-proton collider such as the planned eRHIC project at BNL , which is currently under discussion.

Motivated by this prospect, we present in this paper a study for photoproduction of single-inclusive high-$`p_T`$ hadrons at a future lepton-proton collider with longitudinally polarized beams, and investigate the sensitivity of this reaction to the contribution from resolved photons. We had made first exploratory studies of this, and related, reactions in . Our present study goes beyond the results of  in that we now consistently include both the direct and the resolved photon processes in the next-to-leading order (NLO) of QCD. This is crucial for reliable quantitative analyses in the future as theoretical uncertainties are expected to be under much better control beyond the lowest order (LO) approximation of QCD. For the resolved contribution, the NLO corrections to polarized QCD hard-scattering have become available only recently in . For the polarized direct component, the corresponding cross sections were already calculated some time ago in . We have rederived the expressions and agree with the published results. In addition to extending our study to NLO, we also tailor our phenomenological results to the eRHIC situation.

In the next section we sketch the technical framework for the description of polarized photoproduction reactions. Section III is devoted to numerical studies for the eRHIC. We conclude in Section IV.

## II Technical Framework

We consider the spin-dependent photoproduction cross section for the reaction $`lpl^{}HX`$, where the hadron $`H`$ is at high transverse momentum $`p_T`$, ensuring a large momentum transfer. We may then write the differential cross section in a factorized form:

$`{\displaystyle \frac{d\mathrm{\Delta }\sigma }{dp_Td\eta }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{d\sigma _{++}}{dp_Td\eta }}{\displaystyle \frac{d\sigma _+}{dp_Td\eta }}\right]`$ (1)
$`=`$ $`{\displaystyle \frac{2p_T}{S}}{\displaystyle \underset{a,b,c}{}}{\displaystyle _{1V+VW}^1}{\displaystyle \frac{dz}{z^2}}{\displaystyle _{VW/z}^{1(1V)/z}}{\displaystyle \frac{dv}{v(1v)}}{\displaystyle _{VW/vz}^1}{\displaystyle \frac{dw}{w}}\mathrm{\Delta }f_a^l(x_l,\mu _f)\mathrm{\Delta }f_b^p(x_p,\mu _f)D_c^H(z,\mu _f^{})`$
$`\times `$ $`\left[{\displaystyle \frac{d\mathrm{\Delta }\widehat{\sigma }_{abcX}^{(0)}(s,v)}{dv}}\delta (1w)+{\displaystyle \frac{\alpha _s(\mu _r)}{\pi }}{\displaystyle \frac{d\mathrm{\Delta }\widehat{\sigma }_{abcX}^{(1)}(s,v,w,\mu _r,\mu _f,\mu _f^{})}{dvdw}}\right].`$
This generic form applies to both the direct and the resolved cases. The subscripts “$`++`$” and “$`+`$” denote the helicities of the colliding lepton and proton, $`\mu _r`$ the renormalization scale, and $`S(P_l+P_p)^2`$ the available c.m.s. energy squared. Neglecting all masses, $`V`$ and $`W`$ in Eq. (1) can be expressed in terms of the c.m.s. pseudorapidity $`\eta `$ and transverse momentum $`p_T`$ of the observed hadron:

$$V=1\frac{p_T}{\sqrt{S}}e^\eta \mathrm{and}W=\frac{p_T^2}{SV(1V)}.$$
(2)
We count positive rapidity in the proton’s forward direction.

The sum in Eq. (1) is over all partonic channels $`a+bc+X`$ contributing to the single-inclusive cross sections, with $`d\mathrm{\Delta }\widehat{\sigma }_{abcX}^{(0)}`$ and $`d\mathrm{\Delta }\widehat{\sigma }_{abcX}^{(1)}`$ the associated polarized LO and NLO partonic cross sections, respectively. For the direct contribution to the cross section, depicted in Fig. 1 (a), we have $`a=\gamma `$, while for the resolved ones \[see Fig. 1 (b)\] $`a`$ stands for the parton emerging from the photon. For convenience, we introduce effective polarized parton densities, $`\mathrm{\Delta }f_a^l(x_l,\mu _f)`$, in the lepton, with $`x_l`$ the fraction of the lepton momentum carried by parton $`a`$ and $`\mu _f`$ the factorization scale. The $`\mathrm{\Delta }f_a^l(y,\mu _f)`$ are defined via the convolution

$$\mathrm{\Delta }f_a^l(x_l,\mu _f)=_{x_l}^1\frac{dy}{y}\mathrm{\Delta }P_{\gamma l}(y)\mathrm{\Delta }f_a^\gamma \left(x_\gamma =\frac{x_l}{y},\mu _f\right),$$
(3)
with

$$\mathrm{\Delta }P_{\gamma l}(y)=\frac{\alpha _{em}}{2\pi }\left\{\left[\frac{1(1y)^2}{y}\right]\mathrm{ln}\frac{Q_{\mathrm{max}}^2(1y)}{m_l^2y^2}+2m_l^2y^2\left(\frac{1}{Q_{\mathrm{max}}^2}\frac{1y}{m_l^2y^2}\right)\right\}$$
(4)
being the spin-dependent Weizsäcker-Williams “equivalent photon” spectrum including its non-logarithmic contributions111The latter turn out to be numerically unimportant in our study, but we keep them for completeness. . Here, $`m_l`$ is the lepton mass and $`Q_{\mathrm{max}}`$ the allowed upper limit on the radiated photon’s virtuality, to be fixed by the experimental conditions. For the direct case, one simply replaces

$$\mathrm{\Delta }f_a^\gamma (x_\gamma ,\mu _f)\delta (1x_\gamma ),$$
(5)
while in the resolved contribution the yet unmeasured parton distributions of circularly polarized photons occur that are the focus of this paper. Likewise, $`\mathrm{\Delta }f_b^p(x_p,\mu _f)`$ denotes the density of parton $`b`$ in the polarized proton, and $`D_c^H(z,\mu _f^{})`$ is the usual unpolarized fragmentation function for parton $`c`$ going into the observed hadron $`H`$, at scale $`\mu _f^{}`$. The relations between parton-level and hadron-level variables are

$$s(p_a+p_b)^2=x_lx_pS,x_l=\frac{VW}{vwz},x_p=\frac{1V}{z(1v)}.$$
(6)
As we mentioned earlier, it is a main feature of our paper that for the first time we include both the direct and the resolved contributions to the polarized cross section at NLO, as indicated by the first-order corrections $`d\mathrm{\Delta }\widehat{\sigma }_{abcX}^{(1)}`$ in Eq. (1). The actual calculation of these corrections is generally a formidable task. Since later on we will mainly consider spin asymmetries, we also need the corresponding spin-averaged partonic cross sections. Let us go through the various ingredients that we need in a little more detail.

The resolved contribution is technically equivalent to a hadroproduction cross section such as $`ppHX`$. The corresponding $`d\mathrm{\Delta }\widehat{\sigma }_{abcX}^{(1)}`$ may therefore be immediately adopted from our recent calculation  for this case. In that reference, we also re-derived the (previously known ) expressions for the unpolarized case. There are in total 16 different partonic channels.

The NLO contributions to the direct-photon component have been calculated in  and  for the unpolarized and polarized cases, respectively. Having the techniques available that we used in , it was a straightforward, though tedious, work to check the old results from the scratch. For the polarized case we fully agree with , at an analytical level. In the unpolarized case, we were able to identify some minor mistakes (beyond those noted in ) in the computer code of ; apart from those we agree with that calculation as well.

The technical details of the NLO calculation for the direct part may be found in . We will not repeat these here, but will only recall an important aspect of photoproduction cross sections beyond the LO: writing the cross section as the sum of the direct and the resolved photon contributions,

$$d\mathrm{\Delta }\sigma =d\mathrm{\Delta }\sigma _{\mathrm{dir}}+d\mathrm{\Delta }\sigma _{\mathrm{res}},$$
(7)
we emphasize that neither $`d\mathrm{\Delta }\sigma _{\mathrm{dir}}`$ nor $`d\mathrm{\Delta }\sigma _{\mathrm{res}}`$ are individually physical, i.e., experimentally measurable. This can be easily seen by noticing that the $`23`$ contributions to the direct cross section have singular configurations arising from a collinear splitting of the incoming photon into a quark-antiquark pair. These contributions are absorbed, according to the factorization theorem, into the so-called “pointlike” part of the NLO photonic parton densities and thus into the resolved component. The subtraction of the singular part is unique only up to finite pieces, and the ensuing ambiguity (i.e., the choice of “factorization scheme”) cancels only when the direct and resolved parts are combined.

As an example, when using dimensional regularization ($`d=42\epsilon `$), the appropriate “counter cross section” that cancels the $`\gamma q\overline{q}`$ collinear singularity in the subprocess $`\gamma \overline{q}gX`$ is schematically given by

$$\frac{d\mathrm{\Delta }\widehat{\sigma }_{\gamma \overline{q}gX}^{\mathrm{counter}}}{dvdw}\frac{\alpha _{em}}{2\pi }\mathrm{\Delta }H_{q\gamma }d\mathrm{\Delta }\widehat{\sigma }_{q\overline{q}gg}$$
(8)
where $`d\mathrm{\Delta }\widehat{\sigma }_{q\overline{q}gg}`$ is the polarized Born cross section for the reaction $`q\overline{q}gg`$ in $`d`$ dimensions, and ($`e_q`$ denotes the quark charge)

$$\mathrm{\Delta }H_{q\gamma }(x,\mu _f)=\left(\frac{1}{\epsilon }+\gamma _E\mathrm{ln}4\pi \right)\mathrm{\hspace{0.17em}3}e_q^2\left[2x1\right]\left(\frac{s}{\mu _f^2}\right)^\epsilon +\mathrm{\Delta }h_{q\gamma }(x).$$
(9)
The freedom in the choice of the factorization prescription is reflected by the arbitrary finite piece $`\mathrm{\Delta }h_{q\gamma }(x)`$ which may be subtracted alongside the singular $`1/\epsilon `$ contribution. For the $`\overline{\mathrm{MS}}`$ convention, $`\mathrm{\Delta }h_{q\gamma }(x)=0`$. An alternative factorization scheme ($`\mathrm{DIS}_\gamma `$) was proposed in the unpolarized case  to facilitate the analysis of data on the DIS photon structure function $`F_2^\gamma `$. Here the photonic Wilson coefficient $`C_{2,\gamma }`$ in the NLO expression for $`F_2^\gamma `$ is absorbed into the definition of the quark distributions of the photon. A similar scheme can be set up in the polarized case as well  by subtracting the corresponding polarized coefficient function $`\mathrm{\Delta }C_{1,\gamma }`$ in the structure function $`g_1^\gamma `$. The photonic parton densities in the $`\overline{\mathrm{MS}}`$ and $`\mathrm{DIS}_\gamma `$ schemes are then related in the following way :

$$\mathrm{\Delta }f_a^{\gamma ,\overline{\mathrm{MS}}}(x,\mu _f)=\mathrm{\Delta }f_a^{\gamma ,\mathrm{DIS}_\gamma }(x,\mu _f)+\delta \mathrm{\Delta }f_a^\gamma (x),$$
(10)
where

$`\delta \mathrm{\Delta }f_q^\gamma (x)`$ $`=`$ $`\delta \mathrm{\Delta }f_{\overline{q}}^\gamma (x)=2N_Ce_q^2{\displaystyle \frac{\alpha _{em}}{4\pi }}\left[(2x1)\left(\mathrm{ln}{\displaystyle \frac{1x}{x}}1\right)+2(1x)\right],`$
$`\delta \mathrm{\Delta }f_g^\gamma (x)`$ $`=`$ $`0.`$ (11)
The relation between the partonic cross sections for the direct contribution reads accordingly:

$$d\mathrm{\Delta }\widehat{\sigma }_{\gamma bcX}^{\mathrm{DIS}_\gamma }=d\mathrm{\Delta }\widehat{\sigma }_{\gamma bcX}^{\overline{\mathrm{MS}}}+\underset{a}{}\delta \mathrm{\Delta }f_a^\gamma d\mathrm{\Delta }\widehat{\sigma }_{abcX}$$
(12)
with $`\delta \mathrm{\Delta }f_a^\gamma `$ given above and the symbol $``$ denoting a standard convolution. The $`d\mathrm{\Delta }\widehat{\sigma }_{abcX}`$ are the appropriate LO partonic cross sections. For the resolved cross section all partonic cross sections remain unaffected by this transformation. It is straightforward to verify that the scheme transformation indeed leaves the physical cross section in Eq. (7) invariant.

Our explicit analytical expressions for the polarized NLO subprocess cross sections for the direct and resolved cases are too lengthy to be given here, but can be found in our computer code which is available upon request.

## III Numerical Results

Let us now turn to a first numerical application of our analytical results. Instead of presenting a full-fledged phenomenological study of single-inclusive hadron production in polarized $`lp`$ collisions, we only focus on the most interesting questions: the importance of the NLO corrections, the dependence on unphysical scales, and the sensitivity to the yet unmeasured parton densities of circularly polarized photons.

For all our calculations we choose a c.m.s. energy of $`\sqrt{S}=100\mathrm{GeV}`$ as planned for polarized electron-proton collisions at eRHIC. We limit ourselves to the case of neutral pions ($`\pi ^0`$); extension to other hadrons would be straightforward. We will always perform the NLO (LO) calculations using NLO (LO) parton distribution functions, fragmentation functions, and the two-loop (one-loop) expression for $`\alpha _s`$. The value for $`\alpha _s`$ is taken according to our choice of proton parton distributions, which is CTEQ5  for the unpolarized case and GRSV  for the polarized one.

In the equivalent photon spectrum in Eq. (4) we use similar parameters as the H1 and ZEUS experiments at HERA. We employ $`Q_{\mathrm{max}}^2=1\mathrm{GeV}^2`$ and take the photon momentum fraction to be within $`0.2y0.85`$.

The fragmentation functions for neutral pions, $`D_c^\pi `$, have been determined with some accuracy from analyses of data for pion production in $`e^+e^{}`$ collisions . Even though our knowledge of the $`D_c^\pi `$ certainly still needs improvement, we assume for this study that the fragmentation functions will be well known by the time experiments at a future polarized lepton-proton collider are carried out. For our analysis we take the LO and NLO sets of Ref.  which also yield good agreement with recent high-$`p_T`$ data for $`pp\pi ^0X`$ at $`\sqrt{S}=200\mathrm{GeV}`$ by the PHENIX collaboration .

Since the parton distributions $`\mathrm{\Delta }f^\gamma `$ of circularly polarized photons are completely unmeasured so far, we have to invoke some model for them. Two extreme scenarios were considered in , with “maximal” \[$`\mathrm{\Delta }f^\gamma (x,\mu _0)=f^\gamma (x,\mu _0)`$\] and “minimal” \[$`\mathrm{\Delta }f^\gamma (x,\mu _0)=0`$\] saturation of the “positivity constraint” $`|\mathrm{\Delta }f^\gamma (x,\mu _0)|f^\gamma (x,\mu _0)`$. $`\mu _0`$ denotes the scale where the boundary conditions  for the evolution are specified. For the spin-averaged densities $`f^\gamma `$ the phenomenologically successful GRV photon distributions  were used. We employ these two extreme sets for the $`\mathrm{\Delta }f^\gamma `$ in Eq. (3), along with appropriate LO  and NLO  boundary conditions and evolution equations. Our default choice for the numerical studies will be the “maximal” set.

As we already pointed out in previous LO studies , for single-inclusive measurements rapidity differential cross sections are particularly suited for extracting the resolved contributions and hence the polarized photonic parton distributions. This can be easily understood by looking at the momentum fractions of the photon and proton distributions as a function of rapidity. As we count positive rapidity in the proton’s forward direction, large $`x_\gamma 1`$ are probed at large negative rapidities. Here one expects no difference between the two photon scenarios since the cross section is dominated by the direct photon contribution sitting strictly at $`x_\gamma =1`$. In addition, in the region $`x_\gamma 1`$ the $`\mathrm{\Delta }f^\gamma `$ are dominated by the “pointlike” part which only depends on the starting scale $`\mu _0`$ but not on the details of the hadronic input. On the other hand, for large and positive rapidities small $`x_\gamma `$ values are probed, and the resolved contribution is expected to dominate more and more.

Figure 2 shows our results for the unpolarized and polarized photoproduction cross sections for neutral pions at LO and NLO at $`\sqrt{S}=100\mathrm{GeV}`$ as functions of the laboratory-frame pseudorapidity $`\eta _{\mathrm{lab}}`$. For the asymmetric eRHIC kinematics ($`E_p=250\mathrm{GeV}`$, $`E_e=10\mathrm{GeV}`$), $`\eta _{\mathrm{lab}}`$ is related to the c.m.s. frame rapidity via

$$\eta _{\mathrm{lab}}=\eta +\frac{1}{2}\mathrm{ln}\frac{E_p}{E_e}.$$
(13)
We have integrated over the transverse momentum $`p_T4\mathrm{GeV}`$ of the produced pion. All scales have been set equal: $`\mu _r=\mu _f=\mu _f^{}=p_T`$.

The lower part of the figure displays the so-called “$`K`$-factor”

$$K=\frac{d(\mathrm{\Delta })\sigma ^{\mathrm{NLO}}}{d(\mathrm{\Delta })\sigma ^{\mathrm{LO}}}.$$
(14)
One can see that in the unpolarized case the corrections are roughly constant, of size $`60\%80\%`$ over most of the $`\eta _{\mathrm{lab}}`$ region considered. In the polarized case, we find generally much smaller corrections, at least for the above choice of scales. Only for negative $`\eta _{\mathrm{lab}}`$ where the cross section changes sign are the corrections large.

An important reason why inclusion of the NLO corrections is generally vital is that they are expected to considerably reduce the dependence of the cross sections on the unphysical factorization and renormalization scales. This improvement in scale dependence when going from LO to NLO is a measure of the impact of the NLO corrections, the residual scale dependence at NLO perhaps providing a rough estimate of the relevance of even higher order QCD corrections. Figure 3 shows the scale dependence of the spin-dependent cross section at LO and NLO. In each case the shaded bands indicate the uncertainties from varying the unphysical scales in the range $`p_T/2\mu _r=\mu _f=\mu _f^{}2p_T`$. The solid lines are for the choice where all scales are set to $`p_T`$. One can see that the scale dependence indeed becomes much smaller at NLO and that the lowest order approximation may capture the main features of the process but will not provide a satisfactory quantitative understanding.

Figure 4 demonstrates, as discussed in Section II, that the direct and resolved contributions individually depend on the choice of the factorization scheme, but that the physical, i.e., experimentally accessible, cross section does not. To make this point, we transform the polarized cross sections from the $`\overline{\mathrm{MS}}`$ to the $`\mathrm{DIS}_\gamma `$ scheme as described in Eqs. (10)-(12).

A first determination of the parton content of circularly polarized photons could be one of the key physics goals of the spin physics program at a future lepton-proton collider. Figure 5 shows the experimentally relevant spin asymmetry

$$A_{LL}^\pi =\frac{d\mathrm{\Delta }\sigma }{d\sigma }=\frac{d\sigma _{++}d\sigma _+}{d\sigma _{++}+d\sigma _+}$$
(15)
for $`\pi ^0`$ photoproduction at eRHIC. As before we have integrated over $`p_T4\mathrm{GeV}`$ and have chosen all scales to be $`p_T`$. At positive rapidities the spin asymmetry shows the expected dependence on the choice of the polarized photon densities. In order to be able to judge whether a future measurement can resolve such a difference we also give in Figure 5 an estimate for the expected statistical errors:

$$\delta A_{LL}^\pi \frac{1}{𝒫_e𝒫_p\sqrt{\sigma _{\mathrm{bin}}}},$$
(16)
where $`𝒫_e`$ and $`𝒫_p`$ are the polarization of the lepton and proton beam, respectively, $``$ denotes the integrated luminosity of the collisions, and $`\sigma _{\mathrm{bin}}`$ is the unpolarized cross section integrated over bins in $`\eta _{\mathrm{lab}}`$. We have used $`𝒫_{e,p}=0.7`$ and $`=1`$/fb, which are targets for eRHIC. It should be noted that an integrated luminosity of $`=1`$/fb is expected to be accumulated after only a few weeks of running so that the statistical accuracy may eventually be much better than the one shown in Fig. 5.

To demonstrate that the sensitivity to $`\mathrm{\Delta }f^\gamma `$ is not masked by the dependence on the spin-dependent proton densities, we show in Fig. 5 also the results obtained for the GRSV “max. $`\mathrm{\Delta }g`$” set , which is characterized by a much larger gluon density $`\mathrm{\Delta }g`$ than the one in the GRSV “standard” set. It is evident that at large positive rapidities the asymmetry is indeed mainly determined by the polarized photon structure. In the region of $`\eta _{\mathrm{lab}}0`$ a measurement of $`A_{LL}^\pi `$ would be also an excellent source of information on the polarized gluon density in the nucleon. We note that, in any case, by the time eRHIC is commissioned we expect the spin structure of the nucleon to be known with much better accuracy than at present.

Figure 6 illustrates how different partonic scatterings $`abcX`$ contribute to the NLO ($`\overline{\mathrm{MS}}`$) polarized resolved cross section, versus rapidity. Recall that parton “$`a`$” emerges from the photon and parton “$`b`$” from the proton. Qualitatively similar results are obtained for the “minimal” photon scenario and hence not shown here.

## IV Conclusions

In this paper we have presented the first complete NLO calculation for the spin-dependent photoproduction of inclusive high-$`p_T`$ hadrons, consistently including the direct and resolved photon contributions at NLO. Our results are relevant for future polarized lepton-proton collider projects like eRHIC, to which we have tailored our numerical results. The spin asymmetry for this process is a promising tool for a first determination of the spin structure of circularly polarized photons. We found that the QCD corrections to the polarized cross section are under good control. The polarized cross section shows a significant reduction of scale dependence when going from LO to NLO.

We finally note that a single process alone, like pion photoproduction, will allow to establish the very existence of a polarized resolved-photon contribution and may perhaps point towards one of the two scenarios discussed above. It will, however, not be sufficient to fully disentangle the various quark flavor distributions and the gluon density of the polarized photon. This will require, as for the nucleon (spin) structure, a multitude of processes, combined in a global QCD analysis. Of course, at eRHIC other photoproduction processes, like jets, di-jets, or heavy flavors, will add further valuable information. Measurements at a future polarized linear $`e^+e^{}`$ collider would be also extremely helpful in this respect.

## Acknowledgments

We are grateful to A. Deshpande for many valuable discussions about the eRHIC project. We thank A. Schäfer for discussions. B.J. and M.S. thank the RIKEN-BNL Research Center for hospitality and support during the final stages of this work, and they acknowledge support from the Program Development Fund PD02-98 of Brookhaven Science Associates. B.J. was supported by the European Commission IHP program under contract number HPRN-CT-2000-00130. W.V. is grateful to RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy (contract number DE-AC02-98CH10886) for providing the facilities essential for the completion of this work. This work was supported in part by the “Bundesministerium für Bildung und Forschung (BMBF)” and the “Deutsche Forschungsgemeinschaft (DFG)”.


format: text

text:
# Extended X-ray emission in the high redshift quasar GB 1508+5714 at 𝑧=4.3

## 1 Introduction

High redshift active galactic nuclei (AGN) are powerful tools with which to study the evolution of super-massive black holes (SMBH) and their host galaxies at an early epoch. One important issue concerns the feedback from SMBH to the host galaxy by means of kinetic energy input carried in energetic particle flows in the form of jets. Relativistic jets are inferred to be present in high-redshift radio-loud AGN, particularly in blazar-like objects (i.e. with a jet oriented close to the line of sight), from their extremely high luminosities ($`10^{47}`$$`\mathrm{erg}\mathrm{s}^1`$) and variability in both the X-ray and radio bands (e.g. Fabian et al. 1997, 1999).

In the X-ray band, thanks to the sub-arcsec spatial resolution of the Chandra X-ray Observatory, direct evidence of resolved X-ray jets and extended emission has been obtained in dozens of radio galaxies and quasars (e.g. Chartas et al. 2000, Schwartz et al. 2000, Worrall et al. 2001, Fabian, Celotti & Johnstone 2003b, Harris et al. 2002, Siemiginowska et al. 2003). Models have been proposed to explain the X-ray emission on such scales involving synchrotron, inverse Compton scattering of photons of the synchrotron, Cosmic Microwave Background (CMB) and radiation from the AGN nucleus (e.g. Celotti et al. 2001, Tavecchio et al. 2000, Harris & Krawczynski 2002). Inverse Compton scattering of the CMB is particularly relevant at high redshifts since, as pointed out by Schwartz (2002a), the X-ray surface brightness of jets is roughly constant with redshift as a consequence of the steep $`(1+z)^4`$ increase in the energy density of the CMB. Another factor to consider at high redshifts (above $`2`$) is that the angular size actually increases with redshift. While there have been efforts to search for X-ray jets at high redshifts e.g. up to $`z6`$ (Schwartz 2002b), the highest redshift at which a convincing detection has been obtained so far is $`z`$=3.8 for the radio galaxy 4C 41.17 as reported by Scharf et al. (2003) recently.

The quasar GB 1508+5714 at $`z`$=4.301 (Hook et al. 1995) is one of a dozen quasars at redshifts above 4 with strong radio emission ($`f_{5\mathrm{G}\mathrm{H}\mathrm{z}}>10\mathrm{m}\mathrm{J}\mathrm{y}`$). Variability in the X-ray and radio bands, as well as a high X-ray luminosity of $`10^{47}`$$`\mathrm{erg}\mathrm{s}^1`$ (Moran & Helfand 1997) suggests relativistically beamed X-ray and radio emission of a blazar type. In several high redshift quasars soft X-ray spectral flattening has been found (e.g. Boller et al. 2000, Yuan et al. 2000, Fabian et al. 2001). However, in our XMM observation of GB 1508+5714 (paper in preparation) no evidence for such an effect was found (see also Moran & Helfand 1997).

We report here on the discovery of extended X-ray emission in GB 1508+5714 from a long Chandra ACIS-S observation (Sect. 2) and discuss its possible origin (Sect. 4). We used $`H_0`$=71 km s-1 Mpc-1, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.73, and $`\mathrm{\Omega }_\mathrm{m}`$=0.27 to calculate luminosities and distances throughout the paper. Errors quoted are at the $`1\sigma `$ level unless stated otherwise.

## 2 Observations and data analysis

The quasar GB 1508+5714 was observed with Chandra with its back-illuminated CCD ACIS-S3 in the faint mode for a duration of 91.8 ksec starting from 2001 June 10. The data were taken from the Chandra archive (PI: Paerels, Paerels et al. 2002). The target was placed close to the nominal aim-point with a small offset of 20 arcsec. The data analysis was performed using the standard CIAO tools (version 2.3, caldb version 2.21). The good exposure time is 88.97 ksec after filtering the data. A level 2 events file was created from level 1 data following the standard procedure. Events with grade of 0,2,3,4,6 were selected. Only events with energy in the 0.3–8 keV range were included.

### 2.1 Extended X-ray emission feature

The X-ray image of GB 1508+5714 in the 0.3–2 keV band is shown in Fig. 1, with the bin size as 1 original pixel of the CCD, i.e. 0.5 arcsec. We used the 0.3–2 keV band for spatial analysis as the Chandra point spread function (PSF) is best understood within this energy range. The core of the X-ray source is unresolved (the source profile being consistent with that of the PSF) and coincident with the optical/radio position of the quasar (15h10m02.9s +57d02m43s J2000). In addition, to the south-west of the core the X-ray emission is apparently resolved and extended. This extended part of the emission is peaked at 15h10m02.6s +57d02m42s (indicated by a smaller contour of 22 cts pixel-1 level in Fig. 1), which is $`3`$ arcsec away from the position of the quasar core and has a position angle of $`250^{}`$ relative to the core. At a redshift of 4.3, this angular distance corresponds to a projected distance of 21 kpc. The feature extends up to $`4`$ arcsec (28 kpc) from the core along the direction defined by this position angle, and spans $`2`$ arcsec (14 kpc) in the transverse direction.

The source counts of the extended emission were extracted from a circle with a radius of 2.5 pixel (1.25 arcsec), as shown in Fig. 2. It covers the extended emission region and lies outside the circle of 95 per cent encircled energy of the central point source with a radius of 1.4 arcsec (the circle around the core in Fig. 2). To extract background counts, we have to take into account the contamination of photons from the nearby core. We define three circles as background regions which have the same sizes and are located at the same radial distances from the quasar core as the source region (see Fig. 2). Since the PSF of the point source is nearly symmetric (the object is almost on–axis), the contamination from the quasar core is expected to be eliminated to a large extent with this background subtraction. In the 0.3–8 keV band, the total number of counts in the extended source region is 158. The averaged number of counts in each of the three background regions is 22.3, consistent with what is expected from the contribution of the core for the Chandra PSF. The net source counts of the extended emission (with possible contamination from the quasar core excluded) is estimated to be 136$`\pm `$13, yielding a highly significant detection (11 $`\sigma `$). Its count rate is 1.50$`\pm `$0.14$`\times 10^3`$$`\mathrm{cts}\mathrm{s}^1`$ (0.3–8 keV), which is $``$2.5 per cent of that of the quasar (0.06 $`\mathrm{cts}\mathrm{s}^1`$).

We do not consider this extended emission feature as an artifact. In the data validation and verification report the observation and data processing were noted as having a good aspect solution and no identified problem. Furthermore, the count rates of the extended emission show no sign of time variability on any timescale during the observation.

We attempted to distinguish whether the extended emission is point-like or diffuse by modelling the 2-D quasar core emission using the Chandra PSF and then subtracting the model from the image data. However, the result is inconclusive, given the small number of net source counts of the extended emission in contrast to the nearby bright quasar core, and the uncertainty in modelling the PSF. Future deep X-ray observations with higher spatial resolution are needed to resolve this emission.

### 2.2 X-ray spectrum of the extended emission

The spectra of the extended emission and the background were extracted from the regions described above. The RMF and ARF files were created at the source position. We performed spectral fitting using XSPEC (v.11.2). The effect of degradation in the ACIS quantum efficiency was taken into account by applying the apply\_acisabs provided in CIAO (Chartas & Getman 2003). Given the small number of source counts, the spectrum was binned to have at least 15 counts per energy bin, and the C-statistic is adopted in the spectral fitting. Absorption with Galactic column density (1.6$`\times 10^{20}`$$`\mathrm{cm}^2`$) was included and the value was fixed in the fit. The background-subtracted spectrum is best fitted with a power law model, (C-stat=4.2 for 9 bins), yielding a photon index of $`1.92_{0.33}^{+0.38}`$ (90 per cent confidence level). The spectrum and fitted residuals are shown in Fig.3. A redshifted thermal plasma model (Raymond-Smith) gave a slightly worse but acceptable fit (C-stat=5.2 for 9 bins) with a rest frame temperature $`T=17.2_{7.4}^{+21.1}`$ keV.

The unabsorbed flux density in the 0.5–8keV band is $`9.2\times 10^{15}`$$`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, using the best-fit power law model, which corresponds to a 2.7–42.4 keV luminosity of $`1.6\times 10^{45}`$$`\mathrm{erg}\mathrm{s}^1`$ in the quasar rest frame.

## 3 Constraints from other wavebands

In the radio band, the quasar was detected as a radio source at low and high frequencies, e.g. 222 mJy at 365 MHz (Douglas et al. 1996), 279 mJy at 4.85 GHz (Becker et al. 1991). The radio emission GB 1508+5714 was unresolved with VLA A-configuration at 1.4 (234 mJy) and 8.4 GHz (152 mJy) with 5-min observations (Moran & Helfand 1997). On mas scales the quasar remains unresolved with VLBI observations (the EVN, Frey et al. 1997). The compactness of the radio emission suggests that the postulated radio jet is closely aligned to the line of sight.

The optical imaging data for GB 1508+5714 taken with the Wide-Field-Camera at the Isaac Newton Telescope show no emission at the position of the extended X-ray emission at $`5\sigma `$ limiting magnitudes of 23.20 and 23.57 in the $`i^{}`$ and $`g^{}`$-band, respectively. These data also suggest that it is unlikely that the extended X-ray emission is coming from a foreground or background AGN. This is because at this X-ray flux level the majority of AGNs have I-band magnitudes brighter than $`I`$23.3 (corresponds to $`f_\mathrm{x}`$/$`f_{\mathrm{opt}}`$10), as found in the Chandra deep surveys (e.g. Alexander et al. 2001). In the following section of this paper we consider only the most likely scenario that the extended X-ray emission is associated with the quasar.

## 4 The emission mechanism and nature of the extended emission

We have found a blob of extended X-ray emission about 3 arcsec away from the powerful blazar GB 1508+5714. This corresponds to a projected distance of about 21 kpc, but could be 100 kpc or more if it is the outer part of the blazar jet. Interpretation of this X-ray emission is not straightforward given the paucity of information from other wavebands.

The most obvious emission mechanism is inverse Compton scattering of Cosmic Microwave Background photons. At the quasar redshift of 4.3 the energy density in the CMB is 790 times that now. The relativistic electron population scattering the CMB photons is plausibly due to the blazar jet. It can either be a quasi-static diffuse cloud fed by the jet, or an outer extension of the jet with a bulk Lorentz factor $`\mathrm{\Gamma }10`$. This last possibility has been proposed for the X-ray jets seen around several lower redshift quasars (e.g. Celotti et al. 2001, Tavecchio et al. 2000, Siemiginowska et al. 2003).

Infrared emission from the blazar itself may contribute to the seed photons for inverse Compton scattering (e.g. Brunetti et al. 1997) but the emission would have to be very luminous ($`10^{48}`$$`\mathrm{erg}\mathrm{s}^1`$) to compete with the CMB at that redshift. We have found no published far infrared measurements of this object with which to make an estimate of any such contribution.

An important issue is to explain why the X-ray blob is not seen at radio wavelengths. High resolution radio imaging by the EVN (Frey et al. 1997) shows that the source is unresolved on the 10 milliarcsec scale at 5 GHz. The flux in that component is within 2 mJy of the total flux seen with the VLA. Although the source may vary between the observations, we shall assume that the radio emission from this component is less than a few mJy.

The lack of a detection of the extended component could be a consequence of Compton losses on the electron population, or of a low magnetic field. If the electrons have no relativistic bulk motion, then the X-ray emission requires electrons with Lorentz factor $`\gamma 10^3\gamma _3`$ with $`\gamma _31`$. The magnetic field in the same region required to produce 5 GHz radio emission in our rest frame is $`B10^2\gamma _3^2`$ G. The energy density of this magnetic field is about $`10^4`$ times that in the CMB, which would be reflected in the radio luminosity being $`10^4`$ times that in the X-ray, rather than the observed limit of $`0.1`$. Consequently the magnetic field must be much lower and $`\gamma `$ much higher in order to produce 5 GHz radio emission.

The lifetime of the electrons provides a strong constraint, $`t_\mathrm{c}3.2\times 10^6\gamma _3^1`$ yr. Radio-emitting electrons with $`\gamma 10^5`$, and a plausible galactic field of 1$`\mu `$ G, thus have a lifetime of only $`3\times 10^4`$ yr. If these particles diffuse out (without further acceleration) at a velocity $`v0.1c`$ then they will have gone only 1 kpc while emitting in the radio. X-ray emitting electrons with $`\gamma 10^3`$ can have travelled 100 times further. So the lack of any corresponding strong radio emission from the X-ray blob could just be due to Compton losses on the highest energy electrons required to make that radio emission.

We note that the optical limits may require that the electron population has a low energy turnoff at $`\gamma <100`$.

If the X-ray emission is produced instead within a relativistic jet with $`\mathrm{\Gamma }10`$, the required $`\gamma `$ from the X-rays drops by 10 and the lifetime of the electrons $`t_\mathrm{c}`$ by a similar factor. The magnetic field can now be higher without any conflict with the radio limit. Deep radio and optical imaging will provide valuable information with which to discriminate between the possibilities.

It is thus likely that the luminous X-ray blob is a by-product of the powerful inner blazar jet which must have been in operation for at least $`10^5`$ yr (the light crossing time to the blob). It only represents about one per cent of the power of the blazar but is emitted at considerable distance from the nucleus. There should be energetic protons and nuclei associated with the electrons which mean that significant energy is being deposited into the halo of the host galaxy. The fact that it is best seen in the X-ray band is a consequence of the steep $`(1+z)^4`$ increase in the energy density of the CMB target photons for inverse Compton scattering. The predominance of extended inverse Compton X-ray emission in distant radio sources is becoming apparent from the results on 3C 294 at $`z`$=1.786 (Fabian et al. 2003a), 3C9 (Fabian, Celotti & Johnstone 2003b) and 4C 41.17 (Scharf et al. 2003). Further high resolution X-ray observations of distant radio sources should give more insight into this energy channel, which could influence the formation and growth of galaxies, groups and clusters.

## Acknowledgements

ACF thanks the Royal Society for their support. AC acknowledges the MIUR (COFIN) and ASI for financial support. The Chandra data presented in this paper were obtained from the CXC Chandra archive. We thank all the members of the Chandra team for making the observation and data analysis available. This work has made use of the INT-WFC optical imaging data obtained by the XID imaging programme of the XMM Survey Science Centre (XMMSSC). This research has made use of the NASA/IPAC Extragalactic Database (NED).


text:
# Environmental Effects on the Kinematics of Virgo Cluster Galaxies

## 1. Introduction

Environmental effects play an important role on the evolution of galaxies in clusters. Different mechanisms are invoked to explain the morphology segregation observed in clusters like tidal interactions between galaxies and ram pressure stripping exerted by the intra-cluster medium on disc galaxies (see Combes 2003, van Gorkom 2003 and references therein).

We are completing a survey of Virgo cluster spiral galaxies by FP observations (Chemin 2003) in order to determine the role of tidal interactions and ram-pressure stripping (ICM-ISM interactions) on their ionized gas kinematics. Observed data-cubes will be compared with the results of 3D simulations of ICM-ISM and/or tidal interactions (eg Vollmer 2003).

We will also investigate the influence of the environmental effects on the dark matter distribution by comparing mass models with the results of GHASP, a FP survey of local isolated galaxies (Garrido et al. 2002).

## 2. Observations

We observe galaxies for which exist deep H$`\alpha `$ images (Koopmann et al. 2001, Gavazzi et al. 2002) and which present perturbations in their rotation curves as derived from long-slit measurements (Rubin et al. 1999). The survey spans all morphological types from S0a to Irregulars. We have already observed 25 objects which are distributed over the whole cluster, with a concentration around M87 in the cluster core where the environmental effects are expected to be stronger. The observations take place at the 1.6m Observatoire du mont Mégantic (Canada), 1.93m Observatoire de Haute-Provence (France), 3.6m Canada-France-Hawaii and 3.6m ESO telescopes, equipped with a scanning FP interferometer coupled with a photon-counting camera of very high sensitivity (FaNTOmM and GHASP instruments, Gach et al. 2002).

The observation of NGC 4438 (Fig. 1) has been done on April 2002 at the 3.6m ESO telescope. The field of view is 3.6 and the spatial sampling is 0.42<sup>′′</sup>. The spectral step and the total exposure time were 16 km/s and 9360s respectively.

## 3. Results

Figure 1 (left) presents the H$`\alpha `$ emission map of NGC 4438, a highly inclined disc located near the cluster core. It shows bright HII regions which are associated to the galactic disc and to the nuclear outflow recently evidenced by Kenney & Yale (2002). Diffuse filaments along which are observed a few brighter clumps are seen on the western side of the disc. The right panel of the figure presents the complex velocity field of NGC 4438 (Chemin et al., in prep.). Both fields indicate strong signs of perturbations, probably induced by an ICM-ISM interaction and/or a tidal interaction with a companion.

## References

Chemin, L. 2003, PhD Thesis, Université Paris VI

Combes, F. 2003, astro-ph/0308293

Gach, J.-L., et al. 2002, PASP, 114, 1043

Garrido, O., et al. 2002, A&A, 387, 821

Gavazzi, G., et al. 2002, A&A, 386, 114

Kenney, J. D. P., & Yale E. E. 2002, ApJ, 567, 865

Koopmann, R. A., Kenney, J. D. P., & Young, J. 2001, ApJS, 135, 125

Rubin, V. C., Waterman, A. H., & Kenney, J. D. P. 1999, AJ, 118, 236

van Gorkom, J. H. 2003, astro-ph/0308209

Vollmer, B. 2003, A&A, 398, 525


text:
# Crystal nucleation for a model of globular proteins

## I Introduction

In recent years there has been an enormous increase in the number of proteins that can be isolated, due to the rapid advances in biotechnology. However, the determination of the function of these proteins has been slowed by the difficulty of determining their crystal structure by standard X-ray crystallography. A major problem is that it is difficult to grow good quality protein crystals. Experiments have clearly shown that this crystallization depends sensitively on the physical factors of the initial solution of proteins. An important observation was made by George and Wilson , who showed that x-ray quality globular protein crystals only result when the second virial coefficient, $`B_2`$, of the osmotic pressure of the protein in solution lies within a narrow range. This corresponds to a rather narrow temperature window. For large positive $`B_2`$, crystallization does not occur on observable time scales, whereas for large negative $`B_2`$, amorphous precipitation occurs. Rosenbaum, Zamora and Zukoski then showed  the crystallization of globular proteins could be explained as arising from attractive interactions whose range is small compared with the molecule’s diameter (corresponding to the narrow window of $`B_2`$. In this case the gas-fluid coexistence curve is in a metastable region below the liquidus-solidus coexistence lines, terminating in a metastable critical point. In a seminal study, ten Wolde and Frenkel then showed  in a simulation of a modified Lennard-Jones model with short range attractions that crystal nucleation is significantly enhanced in the vicinity of the metastable critical point, due to a reduction in the nucleation free energy barrier. Thus the conditions for optimal crystallization correspond to being near the metastable critical point. The reduction in the nucleation barrier results from the large density fluctuations that exist near the metastable critical point, which also affect the structure of the critical nucleus. Namely, the initial step toward formation of the critical droplet is the creation of a dense, liquid-like droplet which then forms a crystal nucleus in its interior. This process of first increasing the density and then forming a crystalline order is the opposite of what normally occurs in the liquid-solid phase transition at high temperatures and close to the liquidus line. The authors also showed that classical nucleation theory is invalid in the vicinity of the critical point.
The idea that protein crystal nucleation could be enhanced by critical density fluctuations was subsequently reinforced in a density functional study by Talanquer and Oxtoby. They studied a phase field model with two order parameters, the local density $`\rho (r)`$ and a structural order parameter $`m(r)`$. The phenomenological free energy functional was assumed to have a van der Waals fluid branch and a van der Waals-like solid branch. Their numerical solution of the Euler-Lagrange equations for the nucleating critical droplet and free energy barrier confirmed that the nature of crystal nucleation changes qualitatively in the vicinity of the metastable critical point. They also found that the nucleation barrier is significantly lower in this region than elsewhere in the phase diagram and that classical theory is invalid. Subsequently Sear studied a continuum model in which the structural order parameter was omitted . He carried out a mean field theory near the metastable critical point and showed that the liquid ”tail” of the density profile of the critical droplet is described by an Ornstein-Zernike like exponential behavior, whose range is given by the fluid correlation length. Since this length diverges as one approaches the critical point, so does the critical size of the droplet and the excess number of particles within it. In a subsequent paper Sear  extended his arguments via a scaling theory using the critical exponents for the Ising universality class, so as to be valid in the vicinity of the critical point. In particular, the singular part of the free energy barrier associated with the tail was shown to satisfy the same scaling as the order parameter associated with the metastable critical point (the density difference) and vanishes as one approaches $`T_c`$. Sear did not calculate the shape of the droplet in the core and interface region and hence was unable to calculate the value of the barrier to nucleation.

In this paper we extend the work of Talanquer and Oxtoby, to obtain a better understanding of homogeneous crystal nucleation. We present additional numerical results for the model as well as an approximate analytical description of the critical droplet, including its crystal core and tail (the latter is based on Sear’s work ). From this we compute the free-energy barrier to nucleation in the region between the fluid-solid coexistence curve and the metastable fluid-fluid coexistence curve. As expected, this barrier is smallest in the vicinity of the metastable critical point. We show that classical nucleation theory is incorrect in the region between the liquidus and solid coexistence lines, except presumably very close to the liquidus line.

The outline of the paper is the following. In section 2 we provide a brief review of the model. In section 3 we extend the numerical analysis of  for the nucleation free energy barrier along various thermodynamic paths. We also calculate the profile of the critical droplet as a function of the density and structural order parameter. This yields the excess number of molecules, N , and the number of molecules in the crystalline phase, $`N_c`$ ,of the critical droplet. In the next section we present a brief discussion of the range of validity of classical nucleation theory and show that it is incorrect in the entire metastable region of the phase diagram that we explore. We also show results for the surface tension along different paths in the phase diagram. In section 5 we present an approximate theory for the shape of the critical droplet and from this obtain N and $`N_c`$. The theory accurately describes the critical droplet except in the interface region, and is in qualitative agreement with the numerical results for $`N/N_c`$. Section 6 contains a brief conclusion. Some details of the numerical analysis are given in an Appendix.

## II  Model

We use a model due to Talanquer and Oxtoby  to describe globular protein crystallization. Their phase field model is based on the following grand canonical free-energy functional:

$$\mathrm{\Omega }[\rho ,m]=𝑑\stackrel{}{r}\left[f(\rho ,m)\mu \rho +\frac{1}{2}K_\rho (\rho )^2+\frac{1}{2}K_m\rho _s^2(m)^2\right]$$
(1)
The free energy depends on two order parameters: the (conserved) local density $`\rho (\text{r},t)`$ and a (non-conserved) local structural order parameter that shows whether the system is in a solid or fluid phase $`m(\text{r},t)`$. Here $`f(\rho ,m)`$ is the Helmholtz free-energy density and $`\mu `$ is the chemical potential. Talanquer and Oxtoby  use the van der Waals free energy density for the fluid branch:

$$f_f(\rho ,m)=k_BT\rho \left[\mathrm{ln}\rho 1\mathrm{ln}(1\rho b)\right]a\rho ^2+k_bT\alpha _1m^2$$
(2)
and a corresponding van der Waals free energy for the solid phase:

$$f_s(\rho ,m)=k_BT\rho \left[\mathrm{ln}\rho 1\mathrm{ln}(1\overline{\rho }b)\right](a+a_mm_s^2)\rho ^2+k_bT(\alpha _1m^2+\alpha _2)$$
(3)
Here $`\overline{\rho }`$ plays a role of a weighted coarse grained density $`\overline{\rho }=\rho (1\alpha _3m(2m))`$ where the term in square brackets implements a difference between the fluid and solid close-packing limits. For the fluid close-packing limit, $`\rho b=1`$, while for the solid close-packing limit $`\overline{\rho }b=1`$. The quantity $`m_s(\rho )`$ is the equilibrium value of the order parameter in the solid phase and hence is a solution of the equation $`(f_s/m)_{m_s(\rho )}=0`$. Thus in the solid close-packing limit ($`b\overline{\rho }=1`$) $`m_s=1`$. The parameter $`a_m`$ in equation (3) has been introduced in order to change the range of the attractive interactions between molecules in the solid phase from that in the liquid phase (as described by the parameter a in equation (2)).

The chemical potential in the solid and fluid phases is the first derivative of the free energy with respect to density: $`\mu =(f/\rho )_T`$. In order to get coexisting densities $`\rho _\alpha `$ and $`\rho _\beta `$ we have to solve the equations $`\mu (\rho _\alpha )=\mu (\rho _\beta )`$ and $`\omega (\rho _\alpha )=\omega (\rho _\beta )`$, where $`\omega =f\mu \rho `$. Graphically this gives the well-known ”common tangent” rule for coexistence. By repeating this for different temperatures we can obtain the entire phase diagram.

In order to obtain a critical droplet profile we must solve the Euler-Lagrange equations with appropriate boundary conditions:

$$\frac{\delta \mathrm{\Omega }}{\delta \rho }=0\text{ }\text{ }\text{ and }\text{ }\text{ }\frac{\delta \mathrm{\Omega }}{\delta m}=0$$

i.e.,

$`K_\rho ^2\rho +{\displaystyle \frac{f}{\rho }}\mu =0`$ (4)
$`K_m^2m+{\displaystyle \frac{f}{m}}=0`$
Using the solutions of (4) for $`\overline{\rho }(\text{r})`$ and $`\overline{m}(\text{r})`$ we can obtain such properties of the inhomogeneous system as the free energy barrier $`\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Omega }\{\overline{\rho }(\text{r}),\overline{m}(\text{r})\}\mathrm{\Omega }\{\rho _0,0\}`$ and the surface tension $`\sigma =1/2(_0^{\mathrm{}}K_\rho \left|\rho \right|^2𝑑r+_0^{\mathrm{}}K_m\left|m\right|^2𝑑r)`$.

## III Numerical Results

In Appendix A we describe the numerical methods that were used to solving equations (4). Here we present some results, extending the work of , for the following choice of parameters: $`a=1`$, $`b=1`$, $`\alpha _1=0.25`$, $`\alpha _2=2`$, $`\alpha _3=0.3`$, $`a_m=1`$. The phase diagram for these values is shown in Figure 1, which shows in particular a metastable fluid-fluid coexistence curve. The existence of the metastable critical point affects the nucleation and growth processes in the vicinity of this point.

The main quantity characterizing the nucleation process is a nucleation rate $`I`$, given by

$$I=I_0e^{\mathrm{\Delta }\mathrm{\Omega }/kT}$$
(5)
where $`I_0`$ is a prefactor given by the product of dynamical and statistical parts . In order to understand how this metastable critical point affects nucleation, we calculate the free energy barrier for different thermodynamic paths. (Some results for this barrier were presented in .) The barrier dependence on temperature and density in the metastable region between the liquidus (solubility curve) and solidus lines is shown in Figures 2 and 3. Figure 2 shows that as the temperature decreases, the barrier also decreases, while Figure 3 shows that as the density increases, the barrier decreases. This is the normal behavior of a free energy barrier, since at any point on the solubility curve the free energy is infinite, decreases as one moves away from it and vanishes at the solid coexistence curve.

However, the existence of the metastable fluid-fluid coexistence curve changes the pathways of the constant free energy barrier and constant supersaturation lines (Fig.4) as compared with the case in which the coexistence curve is not metastable. We can also see that along the lines with constant supersaturation the free energy barrier decreases as one approaches the critical point (Fig. 5), which is consistent with . Figure 5 shows that the free energy barrier has a minimum near, but not at, the critical point. The location of this minimum changes from above the critical point for low supersaturation to below the critical point for high supersaturation. It also can be seen that the increase of the free energy barrier below the critical point becomes very sharp for the cases where the constant supersaturation lines intersect the fluid-fluid coexistence curve . In order to understand the behavior of the free energy barrier along the constant supersaturation lines, we first consider their shape (Fig. 4). At high temperatures these lines are more vertical (constant density lines) and the free energy barrier decreases with temperature, as shown in Fig. 2. Near the critical point the curves become almost horizontal (isothermal lines) and the barrier starts to increase (Fig. 3). Thus the shape of the constant supersaturation curve in some sense determines the behavior of the free energy barrier. To understand this in more detail, we note that the derivative of temperature with respect to density at constant supersaturation is

$$\left(\frac{T}{\rho _0}\right)_{\mathrm{\Delta }\mu }=\frac{\left(\frac{\mathrm{\Delta }\mu }{\rho _0}\right)}{\left(\frac{\mathrm{\Delta }\mu }{T}\right)}=\left[\frac{f_{\rho \rho }^{(f)}}{\left(\frac{\mathrm{\Delta }\mu }{T}\right)}\right]_{\rho =\rho _0}$$
(6)
where $`f_{\rho \rho }`$ is evaluated at the background fluid density. This derivative vanishes at the critical point, so that the constant supersaturation lines become horizontal in the vicinity of the critical point. On the other hand, at large temperatures the denominator vanishes, so that the constant supersaturation lines become vertical. Thus we see that the presence of the critical point changes the behavior of the free energy barrier from decreasing as we lower the temperature far from the critical point (the vertical part of the $`\mathrm{\Delta }\mu =const`$ lines) to increasing as we lower the temperature to its critical value (the horizontal part of the $`\mathrm{\Delta }\mu =const`$ lines). Therefore somewhere in between there is a minimum of the free energy barrier. Thus we can conclude that the free energy barrier is relatively unaffected by the existence of the critical point along paths of constant temperature or constant density , whereas it plays a crucial role along paths of constant supersaturation

One quantity of interest is the excess number of molecules, defined as the number of molecules in the presence of the droplet relative to the number of molecules in the spatially homogeneous metastable state:

$$N=_V(\rho (r)\rho _0)𝑑\stackrel{}{r}$$
(7)
We can calculate the number of molecules in the crystalline state using $`m(r)`$, since this phase field shows whether a point is in the solid or fluid state. Thus:

$$N_c=_Vm(r)(\rho (r)\rho _0)𝑑\stackrel{}{r}$$
(8)
Mean field theory yields a divergence in the excess number of particles at the metastable critical point, whereas the number of particles in the crystalline state remains finite. Figure 6 shows $`N_c`$ versus $`N`$ for different supersaturations. We can check our numerical results with the nucleation theorem :

$$\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Delta }\mu }=N$$
(9)
Thus we can take a derivative of the free energy barrier with respect to supersaturation and compare this with our results for the dependence of the excess number of particles on the supersaturation (Fig.7). We see that the free energy barrier decreases rapidly near the critical point as a function of the supersaturation. This happens because the background fluid density as a function of the supersaturation becomes flat in the vicinity of the critical point.

## IV Breakdown of classical nucleation theory

Classical nucleation theory (CNT) assumes that the critical droplet is large and that its interface is very sharp. It predicts that the free energy barrier is

$$\mathrm{\Delta }\mathrm{\Omega }_{cl}=\frac{16\pi \sigma ^3}{3(\rho _s\rho _0)^2\mathrm{\Delta }\mu ^2}$$
(10)
where $`\mathrm{\Delta }\mu `$ is $`\mu \mu _{coex}`$. However, classical nucleation theory becomes increasingly inaccurate as one approaches the metastable critical point . This is because the critical droplet has a diffuse interface in this region, in contradiction with the assumption of CNT. Our goal in this section is to compare CNT with our numerical results for the free energy barrier and to show that CNT has at best a very limited region of validity. It is useful to divide the contributions to the free energy barrier into two parts, i.e. the bulk and surface contributions, and compare our numerical results with the predictions of CNT for each of these separately. First, define

$$\mathrm{\Delta }\mathrm{\Omega }_{bulk}=4\pi \mathrm{\Delta }\omega (r)r^2𝑑r,$$
(11)
$$\mathrm{\Delta }\mathrm{\Omega }_\sigma =4\pi \frac{1}{2}\left[K_\rho |\rho |^2+K_m|m|^2\right]r^2𝑑r$$
(12)
so the free energy barrier can be written as

$$\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Delta }\mathrm{\Omega }_{bulk}+\mathrm{\Delta }\mathrm{\Omega }_\sigma .$$
(13)
The classical approximation for this barrier is

$$\mathrm{\Delta }\mathrm{\Omega }_{cl}=\frac{4\pi R^3}{3}\mathrm{\Delta }\omega _{cl}+4\pi R^2\sigma $$
(14)
In the above $`\mathrm{\Delta }\omega (r)=[f(\rho (r),m(r))\mu \rho (r)][f(\rho _0,0)\mu \rho _0]`$ is the difference between the value of the grand canonical free energy at the points $`r`$ and infinity, $`\mathrm{\Delta }\omega _{cl}=[f(\rho _s,m_s)\mu \rho _s][f(\rho _0,0)\mu \rho _0]`$ and $`\sigma `$ is the surface tension. This division also allows us to determine which of these contributions to the free energy barrier is more important in the breakdown of CNT. We first note that if CNT is correct at some point ($`\rho `$, $`T`$) of the phase diagram, then the first terms of (13) and (14) should be equal; likewise, the second terms also should be equal. Therefore the values of the radii that result from comparison of the first terms, i.e. $`R_1=(3\mathrm{\Delta }\mathrm{\Omega }_{bulk}/|\mathrm{\Delta }\omega _{cl}|)^{1/3}`$, as well as the second terms, $`R_2=(\mathrm{\Delta }\mathrm{\Omega }_\sigma /\sigma )^{1/2}`$ should be equal to each other and to the classical radius of the droplet, $`R_{cl}=2\sigma /|\mathrm{\Delta }\omega _{cl}|`$.

Using this observation, we have checked the accuracy of CNT close to the liquidus line at low temperatures, as well as close to the critical point. The point close to the critical point is $`T/T_c=1.01`$ and $`\rho b=0.33`$, and those close to the liquidus line are $`T/T_c=1.5,\rho b=0.33`$ and $`T/T_c=1.4,\rho =0.2`$. For all these points CNT is incorrect. (For reasons described at the end of the Appendix we have been unable to explore the region very close to the liquidus line, where CNT presumably becomes correct.) Figure 8 shows the density profiles for these cases, while Table $`1`$ shows the values of the corresponding quantities, $`\mathrm{\Delta }\mathrm{\Omega }`$, $`\mathrm{\Delta }\mathrm{\Omega }_{cl}`$, $`R_1`$, $`R_2`$, $`R_{cl}`$ and the actual radius of the droplet $`R_a`$. We define the latter as the radius of a spherical droplet which has an infinitely sharp interface and the same solid density and number of particles as our droplet. As can be seen from Table $`1`$, CNT underestimates the nucleation barrier by a factor of ten as compared with the actual value for the model. It also underestimates the critical droplet size by a factor of two. As one would expect, the ratio of the actual and classical nucleation barriers decreases as we approach the liquidus line. However, it is important to note that even in the region which is ’close’ to the liquidus line the difference between the actual and classical values of the barrier is still large. It seems, therefore, that CNT is accurate only very close to the liquidus line.

Next we try to determine how CNT breaks down. From Table 1 we see that the radius $`R_2`$ extracted from the surface contribution to the free energy barrier is closer to the actual value of the droplet radius than the radius $`R_1`$ extracted from the bulk part. In order to understand why $`R_1`$ is significantly smaller than the actual radius, we examine the dependence of $`\mathrm{\Delta }\omega `$ on $`r`$ (Fig. 9). At the solid core $`\mathrm{\Delta }\omega `$ is negative and equal to the classical $`\mathrm{\Delta }\omega _{cl}=\omega _{solid}\omega _{fluid}`$ (dashed line in Fig. 9). As we reach the interface region of the droplet, $`\mathrm{\Delta }\omega `$ starts to increase, becomes positive and reaches its maximum at the inflection point of the density profile . As $`r`$ increases beyond the interface region, $`\mathrm{\Delta }\omega `$ decays and approaches zero as we reach the background value of the fluid density. In order to obtain the bulk contribution to the barrier one must multiply $`\mathrm{\Delta }\omega `$ by $`r^2`$ and then integrate over $`r`$. The dependencies of $`\mathrm{\Delta }\omega (r)r^2`$ and the ’incomplete’ barrier $`\mathrm{\Delta }\mathrm{\Omega }(r)=_0^r\mathrm{\Delta }\omega (r)r^2𝑑r`$ on $`r`$ are shown in Figures 10 and 11, respectively. The dashed lines in these figures are $`\mathrm{\Delta }\omega _{cl}r^2`$ and $`4\pi /3r^3\mathrm{\Delta }\omega _{cl}`$. The dotted line in Figure 11 is the value of the bulk contribution to the free energy barrier (which is equal to the value of the incomplete barrier as $`r\mathrm{}`$). The value of $`R_1`$ corresponds to the radius at which the bulk contribution to the classical free energy is equal to the bulk contribution of the actual free energy barrier (given by the intersection of the dotted and dashed lines in Fig. 11).

Therefore, one can see that since the interface of the droplet is not sufficiently sharp, the radius $`R_1`$, and therefore the classical radius $`R_{cl}`$ , is significantly different from the actual radius, $`R_a`$, of the droplet. This is responsible for the large difference between the classical and nonclassical values of the free energy barrier. The main contribution to the interface is its ”tail”, which is given by the correlation length and hence becomes very large close to the critical point. At temperatures close to the critical temperature the system should be in the very narrow region near the fluid coexistence curve in order for the radius of the droplet to be much larger than the interface thickness, as given by the correlation length.

## V Theoretical analysis

### V.1 Approximations for the critical droplet.

Equations (4) with appropriate boundary conditions define the saddle point of the free energy landscape in functional space. The boundary conditions are a) that the density and structural order parameter at infinity are given by the metastable values and b) their spatial derivatives are zero at the origin, i.e. $`\rho (\mathrm{})=\rho _0`$, $`m(\mathrm{})=0`$, $`d\rho /dr(0)=0`$ and $`dm/dr(0)=0`$. There is a trivial, spatially uniform solution, satisfying these boundary conditions, given by $`\rho (r)=\rho _0`$, $`m(r)=0`$. This uniform solution corresponds to the supercooled liquid, which is a metastable ”point” in functional space. The critical droplet is a spatially inhomogeneous, unstable solution satisfying the same boundary conditions, which we find by first solving the Euler-Lagrange equations around the droplet core ($`r=0`$). We assume that at the origin the density has a value close to $`\rho _s`$ and the structural order parameter has a value close to $`m_s`$.  (In general the values of the density and structural order parameter at the droplet center are not equal to their equilibrium values, in contrast to the assumption of classical nucleation theory.) First, we represent the solution in the following form:

$`\rho (r)=\rho _s+{\displaystyle \frac{\chi _\rho (r)}{r}}`$ (15)
$`m(r)=m_s+{\displaystyle \frac{\chi _m(r)}{r}}`$ (16)
We substitute (15) and (16) into (4), rewrite these equations in terms of the $`\chi `$’s and then linearize them to obtain

$`K_\rho {\displaystyle \frac{d^2\chi _\rho }{dr^2}}+f_{\rho \rho }\chi _\rho +f_{\rho m}\chi _m=0`$ (17)
$`K_m{\displaystyle \frac{d^2\chi _m}{dr^2}}+f_{m\rho }\chi _\rho +f_{mm}\chi _m=0`$ (18)
Here $`f_{\rho \rho }`$, $`f_{\rho m}`$ and $`f_m`$ are the second derivatives of f at $`\rho _s`$ and $`m_s`$. Because this is an approximation near the droplet core, we obviously cannot use the boundary conditions at infinity. So, in addition to the zero derivatives at $`r=0`$, we add two arbitrary boundary conditions at the origin. The boundary conditions for the $`\chi `$’s therefore become the following: $`\chi _\rho (0)=0`$, $`\chi _m(0)=0`$, $`\chi _\rho ^{}(0)=\mathrm{\Delta }\rho `$ and $`\chi _m^{}(0)=\mathrm{\Delta }m`$, where $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$ are the (unknown) shifts in the values of density and order parameter from $`\rho _s`$ and $`m_s`$. The values of $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$ cannot be equal to zero, because in this case the solution is spatially uniform (see Appendix A). Solving equations (17) and (18), we obtain

$$\chi _\rho =A\left(\frac{sinh(q_2r)}{q_2}\frac{sinh(q_1r)}{q_1}\right)+\frac{\mathrm{\Delta }\rho }{q_2^2q_1^2}\left(\frac{q_2^2}{q_1}sinh(q_1r)\frac{q_1^2}{q_2}sinh(q_2r)\right)$$
(19)
and

$$\chi _m=B\left(\frac{sinh(q_2r)}{q_2}\frac{sinh(q_1r)}{q_1}\right)+\frac{\mathrm{\Delta }m}{q_2^2q_1^2}\left(\frac{q_2^2}{q_1}sinh(q_1r)\frac{q_1^2}{q_2}sinh(q_2r)\right)$$
(20)
where

$`A={\displaystyle \frac{f_{\rho \rho }\mathrm{\Delta }\rho +f_{\rho m}\mathrm{\Delta }m}{K_\rho (q_2^2q_1^2)}}\text{}B={\displaystyle \frac{f_{m\rho }\mathrm{\Delta }\rho +f_{mm}\mathrm{\Delta }m}{K_m(q_2^2q_1^2)}}`$
$`q^2={\displaystyle \frac{1}{2}}\left[\left({\displaystyle \frac{f_{\rho \rho }}{K_\rho }}+{\displaystyle \frac{f_{mm}}{K_m}}\right)\pm \sqrt{\left({\displaystyle \frac{f_{\rho \rho }}{K_\rho }}{\displaystyle \frac{f_{mm}}{K_m}}\right)^2+{\displaystyle \frac{4f_{\rho m}^2}{K_\rho K_m}}}\right]`$ (21)
Thus we obtain the solution at the center of the droplet in terms of two arbitrary parameters, $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$. In principle these parameters could be obtained from the exact solution by matching this form to the boundary conditions at infinity.

Next we find the solution for the tail of the droplet  . As $`r\mathrm{}`$ we have metastable boundary conditions, so we choose for the solution small deviations from the metastable values:

$`\rho (r)=\rho _0+{\displaystyle \frac{\delta _\rho }{r}}e^{q_\rho r}`$
$`m(r)={\displaystyle \frac{\delta _m}{r}}e^{q_mr}`$ (22)
where $`\delta _\rho `$ and $`\delta _m`$ are arbitrary parameters and $`q_\rho =1/\xi _\rho ^{(f)}=(K_\rho /f_{\rho \rho }^{(f)})^{1/2}`$ and $`q_m=1/\xi _m^{(f)}=(K_m/f_{mm}^{(f)})^{1/2}`$. The tail solution for $`m`$ is in fact exact. We now have solutions in both the core and tail regions, but not in the interface region. Figure 12 shows the comparison of the numerically obtained density profile at $`T/T_c=1.2`$ and $`\rho _0=0.33`$ with our approximations for the tail and the core. The tail approximation is rather good, but the core approximation becomes incorrect as one approaches the interface region. For comparison purposes, we used the values of $`\mathrm{\Delta }\rho `$, $`\mathrm{\Delta }m`$, $`\delta \rho `$ and $`\delta m`$ obtained numerically.

### V.2 Approximation for $`N_c/N`$.

$`N`$ and $`N_c`$ can each be considered as the sum of three terms: the number in the core, the tail and the interface. We can calculate the contributions to each of these from the core and tail, using the solutions for these two regions from the above section. Hence we can write

$$N^{(core)}=_0^{R_{core}}(\rho _s+\frac{\chi _\rho (r)}{r}\rho _0)r^2𝑑r$$
(23)
$$N^{(tail)}=_{R_{tail}}^{\mathrm{}}\frac{\delta _\rho e^{q_\rho r}}{r}r^2𝑑r$$
(24)
where $`R_{core}`$ is the maximum value of $`r`$ for which we can still use the approximation near the center of the droplet, and $`R_{tail}`$ is the minimum value of $`r`$ for which the tail approximation can be considered as correct. The interface region is $`R_{core}<r<R_{tail}`$. In the core region, $`\frac{\chi (r)}{r}<<\rho _s`$ so we can drop the second term in (23). Integrating (23) and (24) we find

$`N^{(core)}=(\rho _s\rho _0)V^{(core)}`$ (25)
$`N^{(tail)}=\delta _\rho \xi _\rho ^{(f)}(\xi _\rho ^{(f)}+R_{tail})e^{R_{tail}/\xi _\rho ^{(f)}}`$ (26)
and similarly

$$N_c^{(core)}=_0^{R_{core}}(m_s+\frac{\chi _m(r)}{r})(\rho _s+\frac{\chi _\rho (r)}{r}\rho _0)r^2𝑑r$$
(27)
$$N_c^{(tail)}=_{R_{tail}}^{\mathrm{}}\frac{\delta _\rho e^{q_\rho r}}{r}\frac{\delta _me^{q_mr}}{r}r^2𝑑r$$
(28)
which gives us

$`N_c^{(core)}=m_s(\rho _s\rho _0)V^{(core)}`$ (29)
$`N_c^{(tail)}=\delta _\rho \delta _m{\displaystyle \frac{\xi _\rho ^{(f)}\xi _m^{(f)}}{\xi _\rho ^{(f)}+\xi _m^{(f)}}}e^{R_{tail}(1/\xi _\rho ^{(f)}+1/\xi _m^{(f)})}`$ (30)
Next, we assume we can neglect contributions from the interface region, i.e.

$$\frac{N_c^{(core)}+N_c^{(tail)}}{N^{(core)}+N^{(tail)}}\frac{N_c}{N}$$
(31)
In order to check this approximation, we compare its results with the numerical results. However, to do this we must discuss how we choose to determine the droplet profile, in order to estimate the ratio of $`N_c/N`$.

We first note that the LHS of (31) depends on two parameters, $`R_{core}`$ and $`R_{tail}`$ (for simplicity we didn’t introduce different radii for $`\rho `$ and $`m`$). $`R_{core}`$ depends on the amplitudes $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$, while $`R_{tail}`$ depends on the amplitudes $`\delta _\rho `$ and $`\delta _m`$. Therefore the LHS of (31) depends on four parameters. If we knew the interface solution analytically, we would be able to determine these four parameters by matching these solutions. As the approximation we could match the core and the tail solutions (so $`R_{core}=R_{tail}`$). This is the easiest to do but the least accurate. By setting $`R_{core}=R_{tail}=R_{match}`$ we must take into account that the core solution is obtained for the solid branch and the tail solution is found for the fluid branch of the free energy. Therefore $`R_{match}`$ is a point where the solid branch of the free energy intersects the fluid branch.

Figure 13 shows the dependence of the ratio of the excess number of particles and number of crystalline particles on the background densities for the case in which we match the core solution directly with the tail solution. As can be seen, this method gives the best estimate for this ratio near the critical point, since the tail plays a more important role in that region and the approximation for the tail of $`\rho `$ is more accurate there. However, the actual values of $`N`$ and $`N_c`$ are each underestimated by a factor of ten .

As an alternative method we use the numerically determined $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$ to determine the core solution. Then we match the tail solution with this core solution at the point of intersection of the fluid and solid free energies. In this method the derivatives of $`\rho `$ and $`m`$ with respect to $`r`$ are not equal at the fitting point (Fig. 14). As can be seen from Figures 15 and 16 this approximation is qualitatively good close to the critical point, but becomes worse as $`T/T_c`$ becomes larger. We still underestimate $`N`$ and $`N_c`$, but not by as large a factor as in the above case.

However, close to the metastable critical point $`N^{(tail)}`$ (and therefore $`N`$) diverges and becomes the leading term in the ratio (31), while $`N_c`$ stays finite. Thus we can obtain the asymptotic behavior of the ratio (31) close to the metastable critical point. In this region we can rewrite (31) as

$$\frac{N_c}{N}\frac{N_c}{N^{(tail)}}\frac{1}{\left(\xi _\rho ^{(f)}\right)^2}f_{\rho \rho }^{(f)}$$
(32)
Therefore close to the metastable critical point this ratio vanishes linearly with $`\tau `$ and quadratically with $`\mathrm{\Delta }\rho _c`$:

$$\left(\frac{N_c}{N}\right)_{metastable}\left(\tau +\frac{27}{8}(\mathrm{\Delta }\rho _c)^2\right)$$
(33)
where $`\tau =(TT_c)/T_c`$ and $`\mathrm{\Delta }\rho _c=\rho _0\rho _c`$. Another asymptotic limit is the approach to the liquidus line. In this case the radius of the solid core diverges and therefore $`N^{(core)}`$ and $`N_c^{(core)}`$ become leading terms in (27) so at the liquidus line the ratio becomes

$$\left(\frac{N_c}{N}\right)_{solubility}=m_s$$
(34)
We can see that there are two cases where the number of molecules in the critical droplet diverges, namely close to the critical point and close to the liquidus line, respectively. However, the number of molecules in the solid core of the droplet diverges only in the latter case. Thus close to the critical point most of the molecules are in the tail of the droplet, due to the divergence of the fluid correlation length, as first pointed out by Sear . At the liquidus line most of the molecules are in the solid core, whose radius diverges. Thus instead of one length scale (the radius of the droplet) we have two: the radius of the solid core and the fluid correlation length. The solid core of the globular protein crystal droplet can be observed experimentally , , where it has been shown that close to the metastable critical point the core consists of only a few molecules( about 4-10) . However, the existence of the long tail in the vicinity of the metastable critical point has not been directly observed.

## VI Conclusion

In this paper we have extended the numerical results  to obtain a better understanding of nucleation for their model. In particular, we have calculated the density and structure order parameter profiles of the critical nucleating droplet for different temperatures. This solution of the saddle point equations then yields the free energy barrier to nucleation, which we have shown for a variety of different paths in the phase diagram. We also show in detail how classical nucleation theory is invalid everywhere in the region we studied, except presumably very close to the liquidus line. An approximate theory for the shape and properties of the critical droplet is developed that is in reasonable agreement with the numerical results, although there is clearly a need for a better description of the interface region that separates the core and tail of the droplet. Finally, we note that experiments show the existence of a gel state for globular proteins, in part of the metastable region of the phase diagram . As a consequence, our discussion of the free energy barrier to nucleation for globular proteins would only apply to the region of the phase diagram in which there is no gelation, which is (roughly speaking) for densities $`\rho \rho _c`$. The model studied here does not exhibit gelation.

## VII Acknowledgments:

We wish to thank D. Oxtoby, V. Talanquer and P. Vekilov for helpful discussions. This work was supported by a grant from the National Science Foundation, DMR-0302598.

## VIII Appendix A. Numerical solution of saddle point equations.

Solving the equations (4) is a nonlinear two-boundary problem. One of the standard methods of solving this type of equation is the shooting method. Because is it very difficult to use the conditions at both boundaries simultaneously, we implement the boundary conditions on one side, and choose the number of free parameters equal to the number of boundary conditions on the other side. Then we can solve the resulting initial value problem. Depending on the final values (on the second boundary), we can decide how to change these free parameters. In our case we have two boundary conditions for $`r=0`$, which we can use directly, and two boundary conditions at $`r=\mathrm{}`$, which we cannot use directly. However, we can use $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$ as free parameters.

We still have several difficulties. First of all, the second boundary is at infinity, so we can’t explicitly get the values on this boundary. We could use some large finite value instead of infinity and avoid this problem, as long as the solution doesn’t diverge at infinity for an incorrect choice of the parameters $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$. Unfortunately it does diverge, so the linear shooting method doesn’t work in this simple way. The reasons for the divergence are the error involved with the numerical method and the absence of solutions for most choices of $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$. The error involved with the numerical method increases as we approach a cusp in free energy landscape; this is one of the reasons for the divergence. By varying the parameters $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$ we can push the beginning of this divergence to larger distances, so we can obtain part of the solution (including part of the tail) before the divergence occurs. However, this gives large errors in estimating such quantities as the free energy barrier, surface tension, etc. We could decrease the numerical error of the shooting method in the interface region by using Runge-Kutta methods with adaptive step size, but if we desire reasonably small values of the errors, the step size becomes too small for simulations .

To avoid these difficulties we use a more advanced shooting technique - shooting with a fitting point. The fitting point was chosen at $`r=R_{core}`$ where the solid and fluid branches of the free energy intersect. The value of $`R_{core}`$ depends on the parameters $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$. Next, we know the exact solution for the fluid part of $`m(r)`$ (22), where $`\delta _\rho `$ is a parameter. By matching the fluid and solid parts of $`m(r)`$ and their derivatives at $`r=R_{core}`$ we determine $`\mathrm{\Delta }m`$ for a fixed values of $`\mathrm{\Delta }\rho `$ and $`\delta \rho `$. The next step is to obtain $`\mathrm{\Delta }\rho `$ and $`\delta \rho `$ by matching the values and derivatives of the fluid and solid parts of $`\rho (r)`$. Note that the value of $`\mathrm{\Delta }m`$ depends on the value of $`\mathrm{\Delta }\rho `$, so when we shoot with a particular value of $`\mathrm{\Delta }\rho `$ we obtain a value for $`\mathrm{\Delta }m`$ each time. This shooting with a fitting point minimizes the numerical error as compared with standard shotting method.

There still is one remaining problem, even if we use shooting with a fitting point. Namely, as we approach the fluid coexistence curve, the values of $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }m`$ become very small, so for shooting and fitting we need to vary them very slightly. Thus it becomes beyond the numerical limit of precision of the calculations to obtain a matching at the fitting point as we change these parameters. This so far prevented us from exploring the free energy barrier as one approaches the coexistence curve and hence we have not been able to determine where the classical theory becomes valid.