text:
# 1 Introduction

## 1 Introduction

The dynamics of the Universe in the standard cosmological model is completely determined by the present values (1)(1)(1)A subsript $`0`$ will always denote the value taken today by any quantity. of the Hubble parameter, $`H_0`$, and the ratios $`\mathrm{\Omega }_{0i}`$ between the energy density of the different “matter” components and the critical density. The task assigned to cosmologists by Sandage is the precise determination of these numbers. Unfortunately, they are not directly measurable and their determination relies on the measurement of quantities indirectly dependent on them. Several bounds on the $`\mathrm{\Omega }_{0i}`$ follow from the study of the anisotropies of the cosmic microwave background radiation (CMB), the formation of large scale structure, the age of the Universe, etc. The Hubble diagrams (luminosity distances vs. redshifts) of standard candles have traditionally been a useful tool for determining constraints on some of the cosmological parameters. One of the most astonishing physical results of the past few years  is that not only is the Hubble diagram of high redshift supernovae (SNe) compatible with previous independent constraints on the $`\mathrm{\Omega }_{0i}`$, but it is also complementary to them and in principle allows a complete determination ($`H_066\mathrm{km}.\mathrm{s}^1.\mathrm{Mpc}^1`$, $`\mathrm{\Omega }_{0m}.3`$ and $`\mathrm{\Omega }_{0\mathrm{\Lambda }}.7`$). In particular this result reveals the presence of a so-called “dark energy” component accounting for the dimming of high-$`z`$ SNe as a result of a current period of acceleration in the expansion of the Universe. However, both the existence as well as the size of this Dark Energy component jeopardize our current understanding of the fundamental laws of particle physics  and string theory/quantum gravity .

Very recently Csáki, Kaloper and Terning (CKT) have proposed an interesting alternative explanation of combined Type Ia SNe data from the Calán-Tololo  survey, the Supernova Cosmology Project (SCP) , and the High-$`z`$ Supernova Search Team (HZT)  without requiring a period of accelerating expansion. By assuming the existence of a light axion $`\varphi `$ ($`m_\varphi 10^{16}`$ eV$`/c^2`$) and an intergalactic magnetic field of amplitude $`|𝐁|10^{13}`$ Tesla, the authors of  found that the magnitude-redshift curve corresponding to the best fit flat universe with a cosmological constant  ($`\mathrm{\Omega }_{0m}0.28,\mathrm{\Omega }_{0\mathrm{\Lambda }}0.72`$) can be closely reproduced with cosmic strings replacing the cosmological constant. The axion interacts with photons via a term in the effective Lagrangian of the form, $`_{int}=\frac{1}{M_Lc^2}\varphi 𝐄𝐁`$, where $`\varphi `$ is the axion field and $`M_L`$ is the axion coupling mass scale, which CKT fit to $`M_L410^{11}`$ GeV$`/c^2`$. As a result of this coupling, photons and axions will oscillate in a background magnetic field via a Primakoff effect. Assuming a randomly oriented magnetic field with a coherence length of the order of 1 Mpc, almost one third of the highest redshift photons would decay by the time they reach us.

In this letter we study the statistical significance of the CKT results, and discuss possible generalizations of the CKT mechanism. In order to understand the model dependence of generalized photon loss mechanisms we use the SNe data to identify the allowed regions in the parameter space of such models, obtaining fits at least as good as the cosmological constant accelerating universe. The anisotropy of the CMB still requires a Dark Energy component to the energy-momentum tensor. The large scale structure of the Universe then imposes constraints on the equation of state of the invisible Dark Energy. The difference from previous analyses is that the SNe data now allows for a significantly larger range of the equation of state parameter $`w_X`$ (in $`p=w_X\rho `$) for the Dark Energy component. It is interesting that the preferred region for $`w_X`$ remains in the range $`w_X<0`$, in agreement with requirements of large scale structure. The key point is that instead of indicating an accelerating expansion of the Universe, the dimming of SNe at high redshift would reveal the presence of a light axion which mixes with photons in the intergalactic medium.

Because this subject unites several branches of theory and experiment which some readers may not be familiar with, in Section 2 we briefly review the relevant experiments and current bounds on the relevant cosmological and particle physics parameters. In Section 3 we review the theoretical basis of the CKT photon loss mechanism and its generalizations. This mechanism is intrinsically quantum mechanical, as opposed to the essentially classical models of photon absorption by dust . We then describe in Section 4 the constraints placed on the cosmological parameters by the Hubble diagrams of Type Ia SNe and how they are modified in the presence of a photon loss mechanism.

In Section 5 we perform a statistical analysis of the data and our fits, and demonstrate that the CKT model with a best fit photon decay parameter (with $`\mathrm{\Omega }_{0m}=.3,\mathrm{\Omega }_{0X}=.7,w_X=1/3`$) is as good as the $`\mathrm{\Omega }_{0m}=.28,\mathrm{\Omega }_{0\mathrm{\Lambda }}=.72`$ cosmological constant accelerating universe scenario. We present some exclusion regions of the parameter space, but there is a large range of allowed parameters for a generalized CKT mechanism to accommodate the data.

We conclude with proposals for future experiments to better pin down the Sandage numbers and ferret out the nature of the Dark Energy, and perhaps reveal the existence of an axion coupled to the particles of the Standard Model or other new physics beyond the standard model.

## 2 Experiments and bounds

In this section we summarize experimental bounds on and estimates of the relevant cosmological, astrophysical and particle physics parameters.

### 2.1 Cosmological Constraints

Although the energy density ratios $`\mathrm{\Omega }_{0i}`$ are not directly measurable, they affect the time evolution of the scale factor of the Universe which in turn controls at least three types of measurable quantities: (i) the growth of density perturbations responsible for the formation of structure at large scales; (ii) the age of the Universe; (iii) the propagation of light between its emission by some known astrophysical sources and its detection today.

The sequence of structure formation starting from smaller fluctuations and ending with clusters of galaxies is well described by cold dark matter models, and the structure observed today could well have evolved from the density perturbations revealed by the CMB anisotropy at $`z1100`$. The presence of Dark Energy should not inhibit the growth of density perturbations, which implies that the Dark Energy should have started to dominate only recently. This requires an upper bound on $`w_X=p_X/\rho _X`$ at least as strong as $`w_X<w_m=0`$ .

The fluctuations in the CMB can also be used to predict the distribution of matter, baryon and Dark Energy in the Universe. The position of the first peak in the power spectrum as a function of multipole moment depends on the total density ($`\mathrm{\Omega }_{0m}+\mathrm{\Omega }_{0X}`$) and also less sensitively on $`\mathrm{\Omega }_{0m}`$ and $`\mathrm{\Omega }_{0\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{o}\mathrm{n}}`$. The latest determination by the BOOMERANG experiment is about $`\mathrm{\Omega }_0=1.0\pm 0.1`$ and $`\mathrm{\Omega }_{0m}.33\pm .03`$ (and only a small fraction of baryons) . The flat universe is heavily favored by this observation, as any deviation from flatness tends to grow quickly in Big Bang cosmology. This observational result as well as a nearly scale-invariant spectrum of primeval density perturbations is in accordance with the predictions of inflationary models. The estimate of $`\mathrm{\Omega }_{0m}`$ seems also in agreement, at least for its order of magnitude, with the measurements of mass-to-light ratios in clusters. Even if the CMB cannot directly teach us about the Dark Energy (the CMB is an image of the Universe at $`z1000`$ where the Dark Energy was completely washed out), it still points out a missing component of the energy-momentum tensor.

Some other cosmological constraints, for example from the age of the Universe or gravitational lensing, are less predictive due to high uncertainties and limited data.

Finally, we can determine some cosmological constraints from the Hubble diagrams of some “standard candle” sources: if we know the absolute luminosity of the source, by measuring the flux we receive today, we can deduce the distance travelled by the light, and this distance depends on the $`\mathrm{\Omega }_{0i}`$ ratios. At very low redshift ($`z<.1`$), this dependence is universal while at large redshift ($`z5`$) all the components with $`w_i<0`$ have been washed out by the matter energy density. The intermediate interval is thus the best suited to bring some information on the Dark Energy. Among the astrophysical objects observed in this redshift interval are the Type Ia SNe. And by chance, Type Ia SNe have been found to be decent standard candles by means of matching luminosity vs. time curves to the phenomenological Phillips curve . As we will explain in detail later, the information on the cosmological parameters inferred from the Hubble diagrams relies on the nature of the propagation of light. Assuming a free propagation of the photons in the cosmological background, and assuming a flat universe, the best fit of the SNe data gives  $`\mathrm{\Omega }_{0m}=0.28\pm 0.09`$. This is in good agreement with several other determinations, which is perhaps why the relatively poor current statistics in high-$`z`$ SNe data is sometimes overlooked. The aim of this work is to study how this conclusion is modified when the photons are no longer freely propagating but rather interact with an axion in such a way that some photons are lost between their emission by the SNe and their detection today. The following section will review the CKT photon loss mechanism and we will analyze the SNe data taking into account such a photon loss mechanism in Section 5.

### 2.2 Astro-particle Constraints

The axion is a proposed ultralight pseudoscalar field that would couple weakly with ordinary stable matter. The relevant coupling for the present study is that between axions and photons. The coupling is through a term in the Lagrangian of the form, $`\frac{1}{M_Lc^2}\varphi 𝐄𝐁`$, where $`\varphi `$ is the axion and $`M_L`$ is the axion coupling mass scale. Several searches for axions have turned up negative, placing constraints on the axion mass and couplings. In particular, a study of conversion of solar axions to X-rays in a strong magnetic field give $`M_L>1.7\times 10^9`$ GeV$`/c^2`$ . Based on mixing between axions and photons in an external magnetic field, precisely the effect needed for the CKT mechanism, studies of globular-cluster stars provide a rough limit $`M_L>1.7\times 10^{10}`$ GeV$`/c^2`$ . For more complete reviews, see . Bounds on the axion mass are closely tied to the axion coupling, but the absence of intergalactic line emissions and other experiments place a rough bound of $`m_\varphi <10^3`$ eV$`/c^2`$ . The CKT axion of mass $`m_\varphi 10^{16}`$ eV$`/c^2`$ and coupling of $`M_L410^{11}`$ GeV$`/c^2`$ is within current bounds. Some further studies and bounds can be found in .

The mixing between the axion and the photon requires the existence of the magnetic field. Magnetic fields appear to permeate the observable universe over a large range of scales, from relatively strong fields in the interstellar galactic medium, to weaker fields in the outer envelopes of X-ray clusters. The intergalactic magnetic field can be deduced by observations of radiation from relativistic electrons around X-ray clusters. These observations suggest an intergalactic magnetic field of around $`B10^{13}`$ Tesla. (See  for recent summaries.) The magnitude of the intergalactic magnetic field can also be estimated by other means, including studies of radio emission of distant quasars which, assuming a coherence length for the magnetic field of about 1 Mpc, yields roughly the same estimate . The source of these magnetic fields is not well understood, but there have been several suggestions made, from inhomogeneous cosmological lepton number  to expanding quasar outflows . Based on these results, CKT assumed a bulk magnetic field of $`10^{13}`$ Tesla, and assumed a coherence length of 1 Mpc. The distribution and variation of the magnetic fields in the intergalactic medium is still not well understood. For a recent review of the relevant experiments and the significance of magnetic fields in the early universe, and for a more complete list of references, see .

## 3 CKT photon loss mechanism

Following the pioneering works  of Sikivie, the CKT model considers an axion–photon coupling in the inter-galactic magnetic field (2)(2)(2) Throughout the paper, we will use electromagnetic units such that both the electric and magnetic fields have the dimension of squared energy while the axion scalar field has the dimension of energy. This means that the magnetic field in SI units is given by $`B_{SI}=|𝐁|/\sqrt{4\pi ϵ_0\mathrm{}^3c^5}`$. The speed of light will be denoted by $`c`$, Newton’s constant by $`𝒢_N`$, and the vacuum permittivity by $`ϵ_0`$.

$$_{int}=\frac{1}{M_Lc^2}\varphi 𝐄𝐁$$
(3.1)
where the scale $`M_L`$ characterizes the strength of the axion–photon interactions. In a region where the magnetic field is approximately constant, the polarization of the photon whose electric field is parallel to the magnetic field mixes with the axion to form an oscillating system, much like that of the massive neutrino system. The mass eigenstates in the presence of the background magnetic field are not the same as the mass eigenstates in vacuum, and thus the photons oscillate into axions and vice et versa. A quantum mechanical computation gives the probability that a photon has not yet oscillated into an axion over a distance $`l`$ in a uniform magnetic field as ,

$$P_{\gamma \gamma }(l)=1\frac{4\mu ^2ϵ^2}{m_\varphi ^4c^8+4\mu ^2ϵ^2}\mathrm{sin}^2\left((\sqrt{ϵ^2\lambda _+}\sqrt{ϵ^2\lambda _{}})l/(2\mathrm{}c)\right),$$
(3.2)
where $`m_\varphi `$ is the mass of the axion, $`ϵ`$ is the energy of the photon, $`\mu =|𝐁|/(M_Lc^2)`$ and $`\lambda _\pm =(m_\varphi ^2c^4\pm \sqrt{m_\varphi ^4c^8+4\mu ^2ϵ^2})/2`$ are the squared-mass eigenvalues. Note that despite the fact that $`\lambda _{}<0`$, that mode has positive energy when the energy $`ϵ`$ is taken into account.

If $`\mu m_\varphi c^2`$ then there are two distinct energy regimes:

* High energy photons ($`ϵm_\varphi ^2c^4/\mu `$), for which the oscillation is maximal and achromatic:

$$P_{\gamma \gamma }(l)=1\mathrm{sin}^2\left(\mu l/(2\mathrm{}c)\right)$$
(3.3)
* Low energy photons ($`m_\varphi c^2ϵm_\varphi ^2c^4/\mu `$), for which the oscillation is small and energy-dependent:

$$P_{\gamma \gamma }(l)=14\frac{\mu ^2ϵ^2}{m_\varphi ^4c^8}\mathrm{sin}^2\left(m_\varphi ^2c^3l/(4\mathrm{}ϵ)\right)$$
(3.4)

In the Universe, the intergalactic magnetic field is not uniform but is assumed by CKT to vary in domains of size $`L_{\text{dom}}1\text{Mpc}`$, which is much smaller than the distance to high redshift SNe ($`z1`$ is equivalent to distances about $`3\times 10^3`$ Mpc). CKT assume a uniform distribution of randomly oriented magnetic fields of constant magnitude over the 1 Mpc domains(3)(3)(3)As mentioned earlier, these magnetic fields are poorly understood and may be distributed much more irregularly than is assumed here. We thank Hui Li for a discussion of this point.. In that case high energy unpolarized photons (as defined above) approach an equilibrium population evenly divided between the two photon polarizations and the axion at large distances. CKT found that the approach to equilibrium is through an exponential damping of the number of high energy photons :

$$P_{\gamma \gamma }(l)\frac{2}{3}+\frac{1}{3}e^{l/L_{\mathrm{dec}}},$$
(3.5)
where the decay length is approximatively given by:

$$L_{\text{dec}}=\frac{8}{3}\frac{\mathrm{}^2c^6M_L^2}{L_{\mathrm{dom}}|𝐁|^2}.$$
(3.6)
Before studying the implications of this photon loss mechanism on the SNe data, it is worth mentioning how to evade any unwanted effects on the CMB anisotropy. While the photons of the CMB were emitted at energies of the atomic sizes, when the magnetic field appears they have already redshifted down to much lower energies about $`10^4`$ eV. In this low energy regime, the oscillation mixing is energy suppressed and requiring $`1P_{\gamma \gamma }`$ to be less than the CMB anisotropy ($`10^5`$) gives a lower bound on the axion mass around $`m_\varphi >10^{16}`$ eV . (4)(4)(4) It has been noted that the photon loss probability and the above conclusions may be modified in the presence of the intergalactic plasma , depending on the plasma density. It was suggested in  that a small free electron density can in fact weaken the constraints on the axion parameters. In any case, we have not included such modifications in this analysis.

## 4 Cosmography: luminosity distance vs. redshift with and without photon loss

Let us now study the consequences of this photon loss mechanism on the measurement of luminosity distances. A nice review of the cosmological parameters and their determination for non-specialists can be found in .

A homogeneous and isotropic space, as our Universe appears to be at large scales, is described by the Robertson–Walker metric

$$ds^2=c^2dt^2+R_{}^2a^2(t)\left(\frac{dr^2}{1kr^2}+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2\right)$$
(4.7)
where $`a(t)`$ is the dimensionless scale factor (by convention $`a(\text{today})=1`$), $`R_{}`$ characterizes the absolute size of the Universe, while with an appropriate rescaling of the coordinates the curvature parameter $`k`$ takes values $`+1,0,1`$ for closed, flat and open universes.

For energy-momentum sources made of a superposition of perfect fluids specified by their energy densities $`\rho _i`$ and their pressures $`p_i`$, the Einstein equations reduce to the Friedmann equation (a dot stands for a derivative with respect to the time coordinate):

$$H^2=\frac{8\pi 𝒢_N}{3c^2}\underset{i}{}\rho _i\frac{kc^2}{a^2R_{}^2},\text{where }H=\frac{\dot{a}}{a}$$
(4.8)
together with the conservation equation:

$$\dot{\rho }_i+3(\rho _i+p_i)H=0$$
(4.9)
It is common to introduce the ratio of the energy densities to the critical energy density for a flat universe:

$$\mathrm{\Omega }_i=\frac{8\pi 𝒢_N\rho _i}{3c^2H^2}\text{and }\mathrm{\Omega }_{\text{curv.}}=\frac{kc^2}{a^2R_{}^2H^2}$$
(4.10)
Note that $`_i\mathrm{\Omega }_i+\mathrm{\Omega }_{\text{curv.}}=1`$ by (4.8). In terms of the redshift parameter $`z=1/a1`$ (by convention, since $`a_0=1`$, $`z_0=0`$) and assuming a constant equation of state $`p_i=w_i\rho _i`$ for each component of the stress-energy tensor, the conservation equation (4.9) can be integrated with respect to time, giving a relation between the densities $`\rho _i`$, the scale factor $`a`$ and the equation of state $`w_i`$. This gives,

$$\rho _i=\rho _{0i}(1+z)^{3(1+w_i)}.$$
(4.11)
The Friedmann equation (4.8) can then be rewritten (5)(5)(5)From now we will omit the subscript $`0`$ on the energy density ratios $`\mathrm{\Omega }_i`$ evaluating those quantities only at the present time.:

$$H^2(z)=H_0^2(1+z)^2\left(1+\underset{i}{}\mathrm{\Omega }_i((1+z)^{1+3w_i}1)\right)\text{and }1\underset{i}{}\mathrm{\Omega }_i=\frac{kc^2}{R_{}^2H_0^2}$$
(4.12)
A crucial test of cosmological models comes from Hubble diagrams for standard candles by plotting their luminosity distances versus their redshifts. The luminosity distance is defined by considering a source located at the comoving distance $`r_s`$, which has emitted some light at a time $`t_s`$ with an absolute luminosity $``$ (energy per time produced by the source, normalized to that of a standard candle). This light is detected today by a detector located at $`r=0`$ which measures a flux $``$ (energy per time per area) and a redshift $`z_s`$. Taking into account the energy loss due to the redshift of the light and the increase of time intervals between the emission and the reception, the conservation of energy ensures that,

$$=4\pi r_s^2R_{}^2(1+z_s)^2.$$
(4.13)
The luminosity distance between the source and the detector is defined by,

$$=4\pi d_L^2.$$
(4.14)
The quantities characterizing the source are not independent but can be related to each other through the ratios $`\mathrm{\Omega }_i`$ governing the dynamical evolution of the Universe. This is why the independent measurement of two SNe observables can tell us about Sandage’s numbers. The trajectory of light is given by the null geodesic equation, $`ds^2=0`$, which when integrated becomes,

$$_0^{r_s}\frac{dr}{\sqrt{1kr^2}}=_{t_s}^{\mathrm{today}}\frac{cdt}{R_{}a(t)}=_0^{z_s}\frac{cdz}{R_{}H(z)}$$
(4.15)
Using (4.12) this integral equation can be solved to relate the comoving distance $`r_s`$ of the source to its redshift(6)(6)(6)The definition of the function ‘sinn’ depends on the topology of the Universe: $`\text{sinn}x=\mathrm{sinh}x,x,\mathrm{sin}x`$ for an open, flat, closed Universe, respectively. See for instance .:

$$r_s(z_s)=\text{sinn}\left(\frac{c}{R_{}H_0}_0^{z_s}\frac{dz}{(1+z)\sqrt{1+_i\mathrm{\Omega }_i((1+z)^{1+3w_i}1)}}\right)$$
(4.16)
This relation allows one to predict the luminosity distance of the source in terms of its redshift:

$$d_L=(1+z_s)r_s(z_s)R_{}.$$
(4.17)
So far this expression has been derived assuming that once emitted the light propagates freely in all spatial directions. The CKT proposal that photons can disappear into axions in the intergalactic magnetic field changes the relation (4.17) for $`d_L(z_s)`$. Imagine indeed that over the distance $`r_sR_{}`$ a fraction $`(1P_{\gamma \gamma }(r_sR_{}))`$ of the number of photons has decayed into axions, so that only a fraction $`P_{\gamma \gamma }(r_sR_{})`$ of the energy emitted by the source will be detected and the conservation equation (4.13) has to be modified accordingly:

$$=4\pi r_s^2R_{}^2(1+z_s)^2/P_{\gamma \gamma }(r_sR_{})$$
(4.18)
Then the expression for the luminosity distance as a function of the redshift becomes:

$$d_L=\frac{(1+z_s)r_s(z_s)R_{}}{P_{\gamma \gamma }^{1/2}(r_s(z_s)R_{})}$$
(4.19)
where $`r_s(z_s)`$ is still given by (4.16). Usually, the data are not presented in terms of $`d_L(z_s)`$ but rather in terms of the magnitude defined as

$$m=m_0+5\mathrm{log}_{10}(H_0d_L/c)$$
(4.20)
where $`m_0`$ is the magnitude associated to the absolute luminosity of the source (after corrections deduced from the fitting to the Phillips curve, $`m_0`$ is the common quantity to all SNe and it is fit to the data).

We would like to comment briefly on the effect of photon loss on the measured flux. The luminosities of distant SNe fluctuate too much in general for such SNe to make good standard candles. However, it has been noticed that the absolute luminosities of Type Ia SNe can be determined by fitting the time evolution of the measured SNe luminosity to the phenomenological Phillips curve . Roughly speaking, the broader the light curves of a Type Ia SN, the brighter it is intrinsically. This relation is specified phenomenologically by a stretching factor relating the width of the light curve to its intrinsic luminosity. Given the absolute luminosity, together with the measured luminosity, one determines a normalized flux $``$. The effect of photon loss is to reduce the measured flux, but it is important that the Phillips curve not change as a result. The first point is that the Phillips curve measures brightness on a logarithmic scale as a fuction of the age of the SN. Hence, if an energy-independent fraction of photons is lost by a particular source, the light curve will simply be shifted vertically. This is the effect we have been discussing. The only reason for concern is then whether or not the broadness of the observed light curve should change as a result of photon loss, as that is how the absolute luminosity of a Type Ia SN is determined. The luminosity signal is generally high enough above background that we believe this will be a negligeable effect. A more rigorous analysis should be done, however.

The CKT model described above provides a physical example of a photon loss mechanism in which the $`P_{\gamma \gamma }`$ probability is directly related to particle physics parameters that can therefore be constrained by the SNe data.

What the SNe data indicate (Figure 1) is that SNe with $`z1`$ are dimmer than can be explained by a universe with only ordinary matter in it ($`w_M=0`$), assuming no loss of photons. The universe undergoes a period of acceleration if the Dark Energy has an equation of state $`w_X<1/3`$ (the condition for having an accelerating period of acceleration today is $`w_X<1/(33\mathrm{\Omega }_m)`$ and if $`1/(33\mathrm{\Omega }_m)<w_X<1/3`$ the universe will accelerate in the future). Suppose we choose to assume that there is no period of acceleration. Then we can study the general requirements of a photon loss mechanism to explain the SNe data. We think it is worth studying such generalizations because the authors of  make many assumptions regarding observational astrophysics which are not well understood. For example, as mentioned earlier, the structure of the intergalactic magnetic field may be significantly different than the CKT assumption of 1 Mpc domains of constant magnitude. Furthermore, there may be other acceptable photon loss mechanisms, either via new particle physics or via uniform absorption by an unknown intergalactic medium.

For $`z<1`$ we require that the photon loss be so as to force a positive slope for the magnitude-redshift curve. For example, in the absence of a photon loss mechanism, the dotted curve in Figure 1 must be forced upward. The SNe magnitude is proportional to the logarithm of its luminosity, so by (4.18) this magnitude is shifted by the addition of a term proportional to $`\mathrm{log}(P_{\gamma \gamma })`$. If we assume that for low $`z`$, $`P_{\gamma \gamma }1\delta P_{\gamma \gamma }(z)`$ where $`\delta P_{\gamma \gamma }(z)1`$, then we require that $`\delta P_{\gamma \gamma }(z)`$ have a linear term which can control the slope of the magnitude-redshift curve. Otherwise a fit will likely become more difficult. The other requirement is that if we want to match onto the single Type Ia observation at $`z1.7`$ , then the magnitude-redshift curve must be allowed to slope downward by then. The easiest way to accomplish this is to give $`\mathrm{log}(P_{\gamma \gamma }(z))`$ a small negative slope and positive second derivative near $`z1`$. These requirements are satisfied by the CKT mechanism. We have also checked that in a variety of phenomenological photon loss models (without theoretical basis) for which these two requirements are satisfied, it is rather easy to fit the SNe data at least as well as the accelerating universe by a cosmological constant.

## 5 Statistical fits

In this section we describe the statistics behind our $`\chi ^2`$ fits for the parameters in the generalized CKT model. We have included a total of 18 SNe from the Calán-Tololo Supernova Survey  and 36 from the Supernova Cosmology Project . Following  we removed six outlying data points of their 42 total, corresponding to their Fit C. Furthermore, in the fits presented here we have tentatively dropped the $`z1.7`$ point  because of its large error bars compared to the remaining data. We have checked that the results are nearly unchanged when that point is included.

In our fits we have assumed a two component universe with a fractional density $`\mathrm{\Omega }_m`$ of ordinary matter and $`\mathrm{\Omega }_X`$ of Dark Energy with equation of state $`p_X=w_X\rho _X`$. We assume the CKT model for photon loss parametrized by the photon decay length $`L_{\mathrm{dec}}`$. Our theoretical models thus involve five parameters: $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_X`$, $`w_X`$, $`m_0`$ and $`L_{\mathrm{dec}}`$. So far all the statistical analyses of the SNe data have arbitrarly reduced this parameter space by choosing $`L_{\mathrm{dec}}=\mathrm{}`$ a priori. As we will see now, the introduction of this extra parameter to fit the data changes qualitatively and quantitatively the confidence levels in the determination of the cosmological parameters.

Because the absolute magnitude, $`m_0`$, is related to the calibration of the measurements, we have chosen to systematically fit this variable to the SNe data (both low and high redshift) with a uniform prior. Hence the number of parameters in the fits presented here is always one greater than the effective theoretical parameter space. This procedure is slightly different than that which is sometimes performed, in which the additional annoyance parameters are integrated over. Because the distribution in $`m_0`$ is quite highly peaked, the two procedures are nearly equivalent. There will in addition be a small correction to these fits from the error in the stretching factor of the Phillips curve (we used the effective magnitudes as presented by the SCP collaboration, which already involved fitting the stretching factor). In the various fits presented in this letter, some further constraints have been placed on the remaining parameters of the model in order to make use of previous experimental results and to limit the dimensionality of our confidence plots.

The likelihood of a model to explain experimental data is measured by a $`\chi ^2`$ together with the number of data points and model parameters, which are then translated into a confidence level about the best fit:

$$\chi ^2=\underset{i,j}{}(m_i^{\mathrm{exp}}m_i^{\mathrm{th}})\sigma _{ij}^2(m_j^{\mathrm{exp}}m_j^{\mathrm{th}}).$$
(5.21)
The correlated error matrix $`\sigma _{ij}`$ can be found at http://www-supernova.lbl.gov for the SCP data .

The data are displayed together in Figure 1, together with magnitude-redshift curves of various models including the CKT model. The data are displayed in terms of residual magnitudes with respect to the SCP flat accelerated universe. The CKT model differs only slightly from the SCP model for small $`z`$. Note that for $`1<z<5`$ the CKT magnitude dips below the SCP accelerated universe. This fact may allow the models to be distinguished if better high-$`z`$ data is obtained: a precision of .05 in magnitude at $`z1`$ would allow this distinction between the two scenarios to be made. The deviation at $`z>100`$ should not be taken seriously, as the universal magnetic field was likely insignificant before large scale structure formation (but see for example  for an alternative point of view).

Figure 2 demonstrates the maximum likelihood regions for a flat universe with the CMB preferred value $`\mathrm{\Omega }_m=0.3`$. In this first statistical analysis, the effective theoretical parameter space is two dimensional only, being spanned by the photon decay length, $`L_{\mathrm{dec}}`$, and the equation of state of the Dark Energy component, $`w_X`$. Neglecting a photon-axion coupling would have restricted ourselves to the line $`L_{\mathrm{dec}}=\mathrm{}`$ in the plot, with the conclusion that at $`99\%`$ confidence level the Universe is accelerating. Clearly a non-vanishing photon-axion coupling opens up a valley that extends beyond $`w_X<1/3`$ and thus allows the SNe data to accomodate a Universe which does not accelerate, either today or in the future. Also note that even with the enlarged parameter space due to the photon loss mechanism, negative values for $`w_X`$ are strongly favored, which is an independent result consistent with the constraints coming from the CMB anisotropies and structure formation.

It is worth studying what can be learned from current SNe data while relaxing some constraints coming from the other cosmological experiments. Of course we expect a larger acceptable region of parameter space with the additional photon loss parameter, but in particular it is interesting to determine the most likely region in the enlarged parameter space. First, let us vary the matter density ratio, $`\mathrm{\Omega }_m`$, keeping a flat Universe ($`\mathrm{\Omega }_X=1\mathrm{\Omega }_m`$). The parameter space is now three dimensional. To arrive at the 2D confidence levels plotted in Figure 3, for each point $`(w_X,\mathrm{\Omega }_m)`$, we have minimized $`\chi ^2`$ with respect to $`L_{\mathrm{dec}}`$ (a minimization over the range of values for $`L_{\mathrm{dec}}`$ allowed by the different constraints on the photon-axion coupling and the intergalactic magnetic field). The two figures illustrate the bias introduced in the statistical analysis when removing the decay length as a free parameter. Notice again that while the usual interpretation of the SNe data (in the absence of photon loss) prefers a current period of acceleration at the $`90\%`$ level, the plot on the right shows that the CKT mechanism can easily accommodate a cosmological evolution without acceleration.

In the last fit reproduced in Figure 4, we have assumed that the Dark Energy is formed of a vacuum energy ($`w_X=1`$) and we have examined the bias in the determination of the remaining cosmological parameters. Again, for each point $`(\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_m)`$ in the plot on the right of Figure 4, $`\chi ^2`$ is first minimized with respect to $`L_{\mathrm{dec}}`$. Of course, since the dimming of SNe is now a combined effect of acceleration and the CKT photon loss mechanism, the confidence level regions are enlarged compared to the usual SNe data analysis.

## 6 Conclusions

We have presented an analysis of the SNe data taking into account a photon loss mechanism as recently proposed by CKT. This mechanism offers a new parameter, namely the photon decay length, to accomodate the data. The main conclusion of our analysis, as already announced by CKT, is that the dimming of SNe at large redshift no longer requires an accelerating expansion of the Universe, and the photon-axion coupling opens up the window for the equation of state of the Dark Energy towards less negative values. We argue that the accelerating universe and the CKT models may be distinguishable if better high-$`z`$ data becomes available. In this way future experiments like SNAP may reveal new physics beyond the standard models of both particle physics and cosmology. In the new allowed region of the parameter space, the SNe data could not only confirm the existence of a Dark Energy component but may also reveal the existence of a light pseudoscalar that has so far evaded any direct particle physics detection.

## Acknowledgments

We are especially grateful to Csaba Csáki, Nemanja Kaloper and John Terning for sharing their work prior to publication, and for many interesting discussions. We are also happy to thank Tanmoy Bhattacharya, Guillaume Blanc, Salman Habib and Hui Li for useful discussions. C.G. thanks the hospitality of the T-8 group at Los Alamos, where this work was performed. J.E. is supported in part by the US Department of Energy under contract W-7405-ENG-36.


format: text

text:
# The Thermal Sunyaev Zel’dovich effect : a powerful probe for missing baryons

## 1. Introduction

The universe is turning out to be quite complicated. The favorite cosmology postulates that a mysterious negative pressure field or a cosmological constant accounts for about $`2/3`$ of the total energy budget. The remaining $`2/3`$ is believed to be dominantly cold dark matter. The ordinary baryonic matter only accounts for about $`5\%`$ of the total matter in the universe. But the directly observed components such as stars, interstellar medium, intracluster gas, etc., only account for about $`10\%`$ of the total baryon budget (Persic & Salucci 1992; Fukugita, Hogan & Peebles 1998). This is the so called missing baryon problem. These missing baryons are believed to be in the form of the intergalactic medium (IGM) and are elusive to detect. To understand their state, such as density, temperature, peculiar velocity, metalicity, etc., stands as a major challenge to both observation and theory, and is crucial to understanding the thermal history of the universe, galaxy formation, etc.

The IGM has various observable tracers. (1) Neutral hydrogen absorbs the background light of Quasars and produces the Lyman-$`\alpha `$ forest. (2) The ionized electrons have thermal and peculiar motions and are capable of scattering CMB photons and generate secondary CMB temperature fluctuation, which are known as the thermal and kinetic Sunyaev Zel’dovich effects (TSZ and KSZ effect), respectively. Their precision measurements are becoming routinely available with devoted CMB experiments such as AMIBA. (3) The ionized electrons and protons interact with each other and emit X-rays through thermal bremsstrahlung and contribute to the soft X-ray ($`0.52`$ keV) background (XRB).

These effects have different dependences on the IGM state and probe different phases. The Lyman-$`\alpha `$ forest probes the low density neutral IGM phase where the overdensity $`\delta \delta \rho /\rho 11`$. Since our universe is highly ionized, the Lyman-$`\alpha `$ forest only directly probes a tiny fraction of the missing baryons and radiative modelling is required to extrapolate to the state of the bulk of the gas. The X-ray emission is proportional to $`\rho ^2`$ and primarily probes the densest IGM regions. X-ray flux is diluted by distance and further by the expansion of the universe. So the observed IGM XRB is mostly from the nearby universe ($`z0.5`$) and it is difficult to look very deep. X-ray absorption is another possibility (Perna & Loeb 1998), but is strongly temperature and metalicity dependent. In contrast to these two effects, the SZ effect (both TSZ and KSZ) is a universal and powerful probe of baryons. Firstly, it provides an unbiased sampling of the missing baryons. All free electrons participate in Thomson scattering and contribute to the SZ effect. The Thomson optical depth from the epoch of reionization $`z10`$ to the present is $`\tau 0.1`$, which means that a few percent of CMB photons have been scattered by electrons. Secondly, the Compton scattering does not depend on redshift and is not affected by distance or the expansion of the universe. So, the SZ effect can probe the more distant universe. Thirdly, the SZ effect has a weak dependence on the gas density in contrast to the X-ray emission and thus probes a large range of baryon. Fourthly, from a theoretical viewpoint, the SZ effect is easier to model and simulate than the other two effects. In observations, the next generation CMB experiments such as AMIBA are able to measure the SZ effect with high accuracy and cover a large proportion of the sky. All these advantages suggest that the SZ effect is a powerful probe for the missing baryons. Since the thermal SZ effect is about an order of magnitude larger than the kinetic SZ effect, hereafter we focus on the thermal SZ effect.

The biggest shortcoming of the SZ effect is its lack of redshift information. Since scattering is independent of redshift, the SZ effect on one hand allows a direct probe of the IGM to high redshift, on the other hand makes it challenging to disentangle the contributions arising from different redshifts. We have proposed a method to recover the space distribution and time evolution of the IGM (Zhang & Pen 2001) to take full advantage of the SZ effect.

In this paper, we first describe the SZ effect, its observational feasibility and our variational method to extract the 3-D information of the IGM pressure power spectrum in section section 2. We outline our theoretical model of the SZ effect in section section 3. We conclude in section 4.

## 2. Observation

The Cosmic Microwave Background (CMB) scatters off all free electrons through inverse Compton scattering, causing a change in the CMB temperature and allows us to “see” the IGM through scattering, in analogy to absorption. This effect is known as the Sunyaev Zel’dovich effect. The temperature distortion caused by the TSZ effect (Zeldovich & Sunyaev 1969) is:

$$\mathrm{\Theta }(\widehat{n})\frac{\mathrm{\Delta }T_{CMB}(\widehat{n})}{T_{CMB}}=y(\widehat{n})\frac{xe^x}{e^x1}\left[4x/\mathrm{tanh}(x/2)\right]2yS(x),$$
(1)
where $`\widehat{n}`$ is the direction on the sky and $`xh\nu /(kT_{CMB})`$. Since the typical temperature of the IGM $`10^7`$ K is much higher than $`T_{CMB}=2.73`$K, the photons always gain energy and are kicked up to higher frequencies. Thus, the number density of photons on the (less energetic) left side of the CMB peak decreases and results in an apparent cooling of the CMB temperature in the Rayleigh-Jeans regime. From photon number conservation, the number density of photons on the right side increases and results in a heating. With a multiple frequency SZ survey, this frequency dependence enables us to distinguish the SZ effect from the primary and secondary CMB anisotropies (Cooray, Hu & Tegmark 2000).

For a fixed frequency, what matters is the Compton $`y`$ parameter, which is defined as

$$y(\widehat{n})=\frac{\sigma _T}{m_ec^2}_0^{l(z_{cmb})}n_ekT_g𝑑l=\frac{\sigma _T}{m_ec^2}P_e(\widehat{n})𝑑l.$$
(2)
Here, $`T_g`$ and $`n_e`$ are the temperature and number density of free electrons, respectively. $`P_e`$ is the gas pressure. $`dl`$ is the proper distance along the path of CMB protons. The typical value of the temperature distortion is a few $`\mu `$K. The SZ effect has been detected for targeted clusters where $`\mathrm{\Delta }Tm`$K. The next step is to carry out a random sky SZ survey to measure an unbiased statistical sample. The next generation CMB experiments such as AMIBA can realize this goal. In such observations, we observe the sky map of the CMB temperature fluctuation and thus the angular power spectrum, which is defined by $`C_l_{m=l}^la_{lm}a_{lm}^{}/(2l+1)`$ (Here, $`\mathrm{\Theta }(\widehat{q})\frac{\delta T}{T}(\widehat{q})=^{l,m}a_{lm}Y_{lm}(\widehat{q})`$). In the small angle approximation,

$$C_l=4(\frac{\sigma _T}{m_ec^2})^2_0^{x(z_{cmb})}P_{SZ}(k,z)|_{k=l/x}a^2x(z)\mathrm{d}x(z),$$
(3)
where $`x(z)`$ is the comoving distance and $`P_{SZ}(k,z)`$ is the power spectrum of the gas pressure. $`P_{SZ}`$ contains detailed information about the IGM, but is smeared by redshifts. Because the SZ effect is mostly contributed by nonlinear structures at $`z0.52`$ (Zhang & Pen 2001), it should have a strong cross correlation with galaxies in that redshift range. Taking advantage of the cross correlations with a photometric redshift survey such as Sloan Digital Sky Survey, we are able to statistically extract IGM 3-D correlations and evolution. The basic idea is sketched in fig. 1.

We show that the galaxy selection function to maximize the SZ-galaxy cross correlation coefficient is actually the ratio of the 3-D pressure-galaxy cross correlation power spectrum and the galaxy auto correlation power spectrum. With further investigation, we can extract the pressure power spectrum (Zhang & Pen 2001).

## 3. Theoretical estimation

With all these promising upcoming observational capabilities, we need to understand the SZ effect theoretically. One would like to predict the value of the $`y`$ parameter, the shape of the SZ power spectrum, the angular scale where it becomes dominant over the primary CMB, contributions to $`y`$ and the power spectrum from different redshifts, the dependence on the cosmological parameters and the gas parameters and any other information about the intervening IGM that we can infer.

Two analytical models have been considered to answer these questions. The first is the halo model method (see Komatsu & Kitayama 1999 for a review). The contribution to the SZ effect by a single halo is described by the gas distribution in the halo, for example, an isothermal distribution. The collective effect is integrated by the Press-Schechter formalism.

The second approach is to directly deal with the underlying density field as we proposed (Zhang & Pen 2001). The gas temperature depends on the gravitational potential ($`T\varphi `$) and is thus determined by the underlying density field. So roughly $`P_en_e\varphi `$. Then, we can estimate the evolution of the mean gas pressure, which is shown in fig. 2. Furthermore, the pressure correlation function $`P_e(x)P_e(x+r)n_e(x)\varphi (x)n_e(x+r)\varphi (x+r)`$. With the hierarchical model of the density field (Fry 1984, Scoccimarro & Frieman 1999), we can decompose this high order density correlation into the products of the 2-point density correlations. The result is shown in fig. 3. Our calculation shows that, for a typical $`\mathrm{\Lambda }`$CDM model, (1) the mean temperature distortion $`\mathrm{\Delta }T6\times 10^6`$ K. (2) The SZ power spectrum exceeds the primary CMB anisotropy at $`l2000`$ and peaks around $`l3000`$ ($`\theta 10^{}`$) as shown in fig. 4. The position of the peak is an unambiguous signature of the the underlying gas-dark matter correlation. (3) $`C_l`$ has a strong dependence on $`\sigma _8`$ ($`C_l\sigma _8^{69}`$). (4) The dominant contribution to $`C_l`$ depends on the angular scale, varying from $`z0.4`$ at $`l1000`$ to $`z1.5`$ at $`l10000`$. Smaller angular scales probe the more distant universe (fig. 5).

Non-gravitational heating is a key ingredient needed to understand the IGM (Pen 1999). We have estimated its effect and found that it will leave an unambiguous signature in the 3-D power spectrum (Zhang & Pen 2001). If non-gravitational heating is proportional to the gravitational heating, the only signature is an enhancement of the amplitude of the power spectrum. A spatially homogeneous energy injection will produce an uprising tail in the power spectrum towards smaller scales (fig. 6). The non-gravitational heating parameters including epoch and energy can be extracted from a precise measurement of the SZ map in conjunction with a galaxy photometric redshift survey.

## 4. Conclusion

We summarize the benefits of using the SZ effect as the probe for the missing baryons:

* The SZ effect is based on the Thomson scattering, which is sensitive to the majority of baryons: the ionized IGM. The interaction between photons and electrons (about $`10\%`$ CMB photons are scattered) guarantees the SZ effect to be an unbiased sample of the missing baryons.
* Compton scattering is completely understood physics, so the SZ effect can be cleanly computed from first principles.
* Compton scattering does not depend on redshift and thus the SZ effect can readily probe the IGM at high redshift.
* The upcoming experiments such as AMIBA will provide a wealth of information on the SZ effect. The dependence of the SZ effect on frequency allows the separation from primary and other secondary CMB anisotropies.
* With the aid of galaxy photometric redshift surveys and our variational method, we can obtain the full time and space resolved IGM state.
* We can reasonably estimate the SZ effect in theory both analytically and in simulations from first principles.
* The SZ effect can help us resolve some key questions on IGM such as the role of non-gravitational heating, which is of key importance to galaxy formation.

The SZ effect encodes key cosmological information, is observationally feasible and theoretically tractable. It will allow us to uncover the state of the missing baryons.

## References

Persic, M. & Salucci, P., 1992, MNRAS, 258, 14P

Fukugita, M., Hogan, C.J. & Peebles, P.J.E., 1998, ApJ, 503, 518

Array for Microwave Background Anisotropy 2003, http://www.asia a.sinica.edu.tw/amiba/

Perna, R. & Loeb, A., 1998, ApJ, 503, L137

Zhang,P.J. & Pen, U.L., 2001, ApJ, 549, 18

Zel’dovich, Y.B. and Sunyaev, R., 1969, Ap&SS, 4, 301

Cooray, A., Hu, W. & Tegmark, M., 2000, ApJ, 540, 1

Sloan Digital Sky Survey, http://www.sdss.org

Komatsu, E. & Kitayama, T., 1999, ApJ, 526, L1

Fry, J.N., 1984, ApJ, 279, 499

Scoccimarro, R. and Frieman, J., 1999, ApJ, 520, 35

Pen, U.L., 1999, ApJ, 510, L1


text:
# 1 Introduction

## 1 Introduction

While a satisfactory understanding of nucleation and growth in homoepitaxial systems is beginning to emerge , the dramatic effects of tiny traces of adsorbates continue to surprise researchers. Among a wealth of recent examples, we may mention the influence of CO on island shapes and interlayer transport on Pt(111) , the modification of attachment kinetics by H preadsorbed on Si(001)  and the H-induced enhancement of Pt self-diffusion on Pt(110) . In this situation it seems useful to explore different generic scenarios for adsorbate effects on submonolayer growth within a reasonably simple, yet flexible, theoretical model. The present paper reports on an ongoing effort aimed in this direction.

In Section 2 we briefly review our earlier findings and present new results concerning the dependence of the island density and morphology on the temperature and the characteristics of the adsorbate (mobility and the strength of interaction with adatoms). In the Sec. 3 we sketch a rate equation theory which provides a unified description of two-dimensional nucleation in the presence of predeposited mobile and immobile adsorbates, and derive analytic solutions for some special cases.

## 2 Kinetic Monte Carlo Simulations

### 2.1 Model

We employ a recently introduced solid-on-solid growth model with two surface species $`A`$ and $`B`$, representing the adatoms of the growing material and the adsorbates, respectively . The simulation starts on a flat substrate composed only of $`A`$-atoms. The basic microscopic processes are deposition and migration; desorption is not allowed. Two deposition modes can be studied: (i) simultaneous deposition (codeposition) of both species, with generally different fluxes, $`F_A`$ and $`F_B`$, and (ii) predeposition of a certain adsorbate coverage $`\theta _B`$ prior to growth.

Adatoms and adsorbates migrate with a nearest-neighbor hopping rate $`R_D=k_0e^{E_D/k_BT}`$, where $`k_0=10^{13}`$ Hz and the energy barrier for a particle of type $`X=`$ A or $`B`$ is given by

$$E_D^X=\underset{Y=A,B}{}\left(n_0^YE_{\mathrm{sub}}^{XY}+n_1^{XY}E_\mathrm{n}^{XY}\right).$$
(1)
Here $`E_{\mathrm{sub}}^{XY}`$ is the hopping barrier for a free $`X`$ atom on a substrate atom $`Y`$, $`n_0^Y`$ is equal to one if the substrate atom is of type $`Y`$ and zero otherwise, $`n_1^{XY}`$ is the number of nearest-neighbor $`X`$-$`Y`$ pairs, and $`E_\mathrm{n}^{XY}`$ is the corresponding contribution to the barrier (symmetric in $`X`$ and $`Y`$). There are no lateral interactions between adsorbate atoms ($`E_\mathrm{n}^{BB}=0`$). The diffusion rate of a free adatom is $`D=k_0e^{E_{\mathrm{sub}}^{AA}/k_BT}`$.

In previous work , we have shown that in addition a process of place exchange between an adsorbate and an adatom with at most one lateral bond is necessary to achieve decorated island edges. This occurs at a rate $`k_0e^{E_{\mathrm{ex}}/k_BT}`$. When the barrier $`E_{\mathrm{ex}}`$ is sufficiently low, then adsorbates are floating on the edges of growing islands.

Adsorbates can be mobile or immobile and strongly or weakly interacting with the adatoms. Their mobility is controlled by the strength of interaction of adsorbates with the substrate which is described by the energy barrier $`E_{\mathrm{sub}}^{BA}`$. The interaction with adatoms is determined by the energy barrier $`E_\mathrm{n}^{AB}`$.

### 2.2 Summary of previous results

We have investigated previously the case of floating mobile adsorbates in the situation when the concentration of adsorbates is comparable with that of the growing material . We studied the dependence of the island density $`N`$ on the flux $`F_A`$ and the coverage $`\theta _A`$ for both codeposition and predeposition. We found that the adsorbates strongly increase the island density without appreciably changing its power-law dependence on flux, the increase being stronger for larger $`E_\mathrm{n}^{AB}`$. The increase was only slightly higher for predeposition than for codeposition.

We also observed a stronger coverage dependence of the island density in comparison with homoepitaxy. This was interpreted as a delay of the saturation regime, where island density becomes independent of coverage. The coverage dependence is again more pronounced for strongly interacting adsorbates, but it is weaker in the case of predeposition than for codeposition. A further noteworthy feature is the much higher density of free adatoms than in homoepitaxy, which shows an intriguing, oscillatory coverage dependence in the case of codeposition .

Many of these simulation results could be qualitatively explained in terms of a simple rate equation theory, which rests on the assumption that the adsorbates affect the growth process only by slowing down the diffusion of adatoms. If, in addition, the adsorbates are treated as immobile traps (which may be justified at least for sufficiently high adsorbate coverage $`\theta _B`$, see Section 3.4), the adatom migration is described by an effective diffusion coefficient

$$\overline{D}(\theta _B)=\frac{D}{1\theta _B+\theta _Be^{E_\mathrm{n}^{AB}/k_BT}}.$$
(2)
Recently, we have considered a situation when only a low concentration of predeposited adsorbates is present on the surface . We investigated the dependence of island density on flux and on the concentration of adsorbates $`\theta _B`$. Even concentrations as low as $`\theta _B=0.002`$ ML can lead to a severalfold increase of the island density. However, the increase significantly depends on the mobility of the adsorbates. In the case of essentially immobile adsorbates ($`E_{\mathrm{sub}}^{BA}=5`$ eV), a new feature in the flux dependence of the island density appears. Instead of a single power law relationship, there is a plateau where $`N\theta _B`$, reflecting the dominance of heterogeneous nucleation in a flux interval $`F_1>F>F_2`$. Scaling arguments yield the estimates $`F_1D\theta _B^2`$, $`F_2\theta _BDe^{E_\mathrm{n}^{AB}/k_BT}`$ for the plateau boundaries in terms of impurity coverage and strength. For sufficiently small $`\theta _B`$ these were confirmed by scaling plots of $`N/\theta _B`$ vs. $`F/F_1`$ and $`F/F_2`$ .

In this work we further extend the above results. In the following subsection, we discuss effects of adsorbates on the island morphology, and in the subsequent subsection we present results on the temperature dependence of the island density. The energy barriers $`E_{\mathrm{sub}}^{AA}=0.8`$ eV, $`E_{\mathrm{sub}}^{AB}=0.1`$ eV, $`E_{\mathrm{sub}}^{BB}=0.1`$ eV, and $`E_\mathrm{n}^{AA}=0.3`$ eV were fixed in the following. Other energy barriers as well as the ratio $`D/F`$ were varied.

### 2.3 Island morphology

The island shape in homoepitaxy depends on the rates of the various kinetic processes involved. The size of islands varies with the ratio $`D/F`$. To be specific, we set in this subsection the temperature to $`T=500`$ K. Then, the island shapes in homoepitaxy with the parameters mentioned above in Section 2.2 are regular (approximately rectangular, close to square shape, nondendritic) for fluxes in the interval $`F=0.000250.64`$ ML/s. The shapes become irregular for higher fluxes. In order to evaluate the effect of adsorbates, we performed several sets of simulation with fixed flux $`F=0.001`$ ML/s and varying adsorbate properties.

Let us consider first the case of floating adsorbates. It is useful to distinguish two situations: (i) well decorated island and (ii) poorly decorated islands. In the former case, the majority of edge sites is decorated by an adsorbate. We found that the shape is very regular, only the islands are smaller; in fact it can be shown that the decoration enhances edge diffusion and thus makes the island edges smoother . The shape is not very sensitive to the mobility of adsorbates. It only has to be sufficiently high so that adsorbates can reach the island edges. The exchange of adatoms approaching the step edge and the adsorbates decorating the step edge is the dominant mechanism.

If the islands are poorly decorated (the number of adsorbates attached to the island edges is considerably less than the number of perimeter sites), then, the shape changes to irregular and regularity is restored again only for very small $`\theta _B`$. Rather than promoting edge diffusion, in this regime the adsorbates bound on the island edges act as sinks for adatoms provided that the adatom-adsorbate interaction (given by $`E_\mathrm{n}^{AB}`$) is sufficiently strong. Another factor affecting the island shape in this situation is the change of the mobility of adsorbates (given by $`E_{\mathrm{sub}}^{BA}`$). Examples of the different island mophologies for various values of parameters $`E_\mathrm{n}^{AB}`$ and $`E_{\mathrm{sub}}^{BA}`$ are shown in Fig. 1. We can see that the increase of both energy barriers leads to more irregular island shapes and also to an increase of island density.

The irregularity is not so apparent for weak adsorbate-adatom interaction (left column) or high mobility of adsorbates (top row). The immobile strongly interacting adsorbate catches more diffusing adatoms and at the same time it is floating on the growing island edge. Hence, the island grows faster in the neighborhood of such adsorbates. This kind of adsorbates also block the motion of adatoms along the edge and thus suppresses the smoothening of the island edges.

Let us discuss briefly the island morphology for non-floating adsorbates. The mechanism of floating becomes less relevant for adsorbates with low mobility and does not act at all in the limit of immobile adsorbates. Non-floating adsorbates are incorporated inside a growing island and beyond a certain island size have no effect on growth. In this context, we want to note that surface defects can be treated in our approach as immobile adsorbates.

### 2.4 Temperature dependent island density

While in our simulations we have so far focused on the flux dependence of the island density \[5-7\] in experiments it is most often measured as a function of temperature at fixed flux. Therefore, we complement our previous results by the calculation of the temperature dependence. We restrict ourselves here to the case of predeposition of immobile (non-floating) adsorbates. In Fig.2 we show data for several adsorbate concentrations. We find, similar to the flux dependence , a plateau interval $`T_1<T<T_2`$, where $`N`$ is almost independent of $`T`$. For low temperatures, $`T<T_1`$, homogeneous nucleation is more important than heterogeneous nucleation by adsorbates, because the adatoms are not sufficiently mobile to reach the adsorbates. In the plateau region nucleation is predominantly heterogeneous. For high temperatures, $`T>T_2`$, adatoms start to detach from adsorbates and homogeneous nucleation is again important. This kind of behavior is observed experimentally in defect nucleation on oxide and halide surfaces  and in growth on surfaces with strain-relief patterns, which act as an ordered array of defects .

For a quantitative estimate of the transition temperatures $`T_1`$ and $`T_2`$, we recall the three relevant time scales in the problem : (i) The diffusion time $`\tau _D1/D\theta _B`$, which is the time required for an atom to explore the “capture zone” of area $`1/\theta _B`$ associated with a single impurity, and hence to get trapped at the impurity. (ii) The deposition time $`\tau _F\theta _B/F`$, which is the time between subsequent arrivals of adatoms within the capture zone. (iii) The trapping time $`\tau (1/D)e^{E_\mathrm{n}^{AB}/k_BT}`$ at the impurity. The plateau regime is characterized by $`\tau >\tau _F>\tau _D`$, which yields the expressions

$$T_1=\frac{E_{\mathrm{sub}}^{AA}}{k_B\mathrm{ln}(k_0\theta _B^2/F)},T_2=\frac{E_{\mathrm{sub}}^{AA}+E_\mathrm{n}^{AB}}{k_B\mathrm{ln}(k_0\theta _B/F)}.$$
(3)
Using $`E_{\mathrm{sub}}^{AA}=0.8`$ eV, $`F=0.001`$ ML/s and $`\theta _B=0.001`$ ML, we get $`T_1=`$ 403 K, while $`T_2=`$ 465 K for $`E_\mathrm{n}^{AB}=0.4`$ eV and $`T_2=`$ 542 K for $`E_\mathrm{n}^{AB}=0.6`$ eV, in good agreement with the simulations. From the expression for $`T_2`$, we directly see how the width of the plateau grows with increasing $`E_\mathrm{n}^{AB}`$.

## 3 Rate Equations

In this section we develop a set of rate equations for the description of two-dimensional nucleation in the presence of mobile and immobile adsorbates. Throughout, the impurities will be assumed to have been predeposited at a coverage $`\theta _B`$, though codeposition is easily incorporated along the lines of our earlier work . To keep matters transparent, we also restrict ourselves to the case of irreversible aggregation (critical nucleus size $`i^{}=1`$). Detachment of impurities from adatom islands will, however, be included, since this process plays a crucial role at low deposition fluxes.

In addition to the adatom density $`n`$ and the island density $`N`$, a minimal description of nucleation in the presence of adsorbates turns out to require the monitoring of two additional quantities : The density of free impurities $`n_0`$, and the density $`n_c`$ of critical clusters bound at impurities; in the present case $`i^{}=1`$ the latter is simply the density of bound adatom-impurity pairs. Taking into account all relevant processes, one arrives at the following set of evolution equations:

$$dn/dt=FDn(2n+N)Dn(n_0+n_c)+n_c/\tau ,$$
(4)
$$dN/dt=Dn(n+n_c),$$
(5)
$$dn_0/dt=Dnn_0\stackrel{~}{D}n_0N+n_c/\tau +N/\tau ^{},$$
(6)
$$dn_c/dt=Dn(n_0n_c)n_c/\tau .$$
(7)
Here $`1/\tau k_0e^{(E_{\mathrm{sub}}^{AA}+E_\mathrm{n}^{AB})/k_BT}`$ denotes the rate of detachment of an adatom from an impurity, $`\stackrel{~}{D}=k_0e^{E_{\mathrm{sub}}^{BA}/k_BT}`$ is the impurity diffusion constant, and $`1/\tau ^{}k_0e^{(E_{\mathrm{sub}}^{BA}+E_\mathrm{n}^{AB})/k_BT}`$ is the detachment rate of impurities from islands. In the expression for $`\tau `$, we have assumed that the impurities are (much) less mobile than the adatoms, so that the relative mobility is dominated by the adatom mobility.

Most of the impurity terms in Eqs.(4-7) should be self-explanatory. We merely point out the term $`n_c/\tau `$ in (4,6,7) which describes the dissociation of adatom-impurity pairs, and the term $`N/\tau ^{}`$ in (7) which describes detachment of impurities from the islands. Strictly speaking, all reaction terms should be adorned with dimensionless capture coefficients. However, since our objective here is a qualitative, rather than a quantitative description, all these coefficients have been set to unity.

### 3.1 Dimensionless formulation

Following Tang  we now rescale the adatom and island densities by $`\sqrt{F/D}`$ and time by $`1/\sqrt{DF}`$. In addition, the impurity densities $`n_0`$ and $`n_c`$ will be rescaled by the initial impurity coverage $`\theta _B`$. Marking dimensionless quantities with a hat, we obtain the rescaled evolution equations

$$d\widehat{n}/d\widehat{t}=1\widehat{n}(2\widehat{n}+\widehat{N})a\widehat{n}(\widehat{n}_0+\widehat{n}_c)+b\widehat{n}_c,$$
(8)
$$d\widehat{N}/d\widehat{t}=\widehat{n}(\widehat{n}+a\widehat{n}_c),$$
(9)
$$d\widehat{n}_0/d\widehat{t}=\widehat{n}\widehat{n}_0c\widehat{n}_0\widehat{N}+(b/a)\widehat{n}_c+(bc/a^2)\widehat{N},$$
(10)
$$d\widehat{n}_c/d\widehat{t}=\widehat{n}(\widehat{n}_0\widehat{n}_c)(b/a)\widehat{n}_c,$$
(11)
where $`a=\theta _B\sqrt{D/F}`$, $`b=\theta _B/F\tau `$ and $`c=\stackrel{~}{D}/D`$. The initial conditions are $`\widehat{n}(0)=\widehat{N}(0)=\widehat{n}_c(0)=0`$, $`\widehat{n}_0(0)=1`$.

The parameters $`a`$ and $`b`$ may be written in terms of the characteristic fluxes, $`F_1`$ and $`F_2`$, which limit the plateau regime in the case of immobile impurities, as $`a=\sqrt{F_1/F}`$ and $`b=F_2/F`$. For impurities to be at all relevant, it is necessary that $`F_1/F_2=D\theta _B\tau =a^2/b1`$. By retaining only the most important terms, we can analytically extract the behavior predicted by Eqs.(8-11) in simple cases. So far this has been achieved only for immobile impurities ($`c=0`$).

### 3.2 Immobile impurities without detachment

Here we consider the plateau regime, where $`a1`$ and $`b1`$. Numerical integration of (8-11) with $`b=c=0`$ and $`a1`$ shows that the adatom density is time independent, taking the value $`n=F/D\theta _B`$ to high accuracy. For early times, this reflects the balance between deposition and capture at impurities, while at late times the capture at islands dominates; since $`N=\theta _B`$ at late times, the resulting adatom density is the same in both regimes. Setting $`nF/D\theta _B`$, the remaining three rate equations become linear and are readily integrated, with the result

$$N=\theta _B(1(1+\theta /\theta _B)e^{\theta /\theta _B}),$$
(12)
$$n_0=\theta _Be^{\theta /\theta _B},$$
(13)
$$n_c=\theta e^{\theta /\theta _B}.$$
(14)
Both impurity species decay exponentially at late times.

### 3.3 Immobile impurities with detachment

The regime $`FF_2`$ is detachment-dominated, in the sense that an adatom typically visits many impurities before being incorporated in an island. The primary effect of the impurities is then to slow down the diffusion by temporarily trapping the adatoms . To see how this emerges from the rate equations, we infer from numerical integration that, in the relevant late time regime, an equilibrium is established between attachment and detachment of adatoms at free impurities. This implies that $`Dnn_0n_c/\tau `$. Equation (6) then yields $`dn_0/dt0`$, so that the impurity concentration essentially retains its initial value $`n_0=\theta _B`$, and $`n_c(D\theta _B\tau )n`$. In this way, the impurity rate equations (6,7) have effectively been eliminated, and we are left with a pair of modified equations for $`n`$ and $`N`$. Taking into account that $`D\theta _B\tau =F_1/F_21`$, and therefore $`n_cn`$, they read

$$dn/dt=FDn[(D\theta _B\tau )n+N],$$
(15)
$$dN/dt=(D\theta _B\tau )Dn^2.$$
(16)
Standard analysis  shows that the asymptotic island density is of the order $`N(F\tau \theta _B)^{1/3}`$, which can be brought into the familiar form $`N(F/\overline{D})^{1/3}`$ by identifying the effective diffusion constant as $`\overline{D}=1/(\tau \theta _B)`$. This is precisely the strong impurity limit of (2).

### 3.4 The effect of impurity mobility

The second term on the right-hand side of Equation (6) describes the depletion of free impurities due to capture at stable islands. This is partly compensated by the detachment of impurities from islands described by the last term, but nevertheless, the net effect is expected to be a decrease of the island density, because the depletion reduces the ability of mobile impurities to act as nucleation centers.

In the absence of an analytic study of the full set of rate equations, in Fig. 3 we provide some sample results obtained by numerical integration. The adatom diffusion constant $`D`$ and the detachment time $`\tau `$ were chosen in accordance with the parameters employed in the KMC simulations. Data for three different choices of the impurity diffusion barrier $`E_{\mathrm{sub}}^{BA}`$ are shown, and the behavior in the absence of impurities is included for reference. The case of low impurity mobility ($`E_{\mathrm{sub}}^{BA}=1.2`$ eV, corresponding to $`\stackrel{~}{D}/D10^4`$) is particularly noteworthy, as it indicates an intermediate scaling regime where the island density scaling exponent, defined by $`NF^\chi `$, is larger than the homoepitaxial value $`\chi =1/3`$ (a fit yields $`\chi 0.42`$). Further exploration of this phenomenon will be presented elsewhere.

#### Acknowledgments.

J.K. would like to thank CAMP and the Department of Physics at DTU for their kind hospitality while this paper was prepared. Support by the COST project P3.130 and by Volkswagenstiftung are gratefully acknowledged.


text:
# Derivatives of the fixed-node energy

## I Introduction

Quantum Monte Carlo simulations are emerging as a promising alternative route towards high-accuracy predictions of molecular and materials properties. The fixed-node approximation fnc  is widely used in quantum Monte Carlo simulations mitas  of many-fermion systems to avoid the sign problem kalos . It consists of solving for the lowest-energy eigenstate, $`\mathrm{\Phi }`$, of the Hamiltonian, $`H`$, with the boundary conditions that the wave function vanishes on the nodal surface of an anti-symmetric trial function, $`\mathrm{\Psi }`$, i.e. at the points $`R`$ in configuration space where $`\mathrm{\Psi }(R)=0`$. The fixed-node energy (FNE)—i.e. the eigenvalue $`E`$ of the Hamiltonian corresponding to $`\mathrm{\Phi }`$—is an upper bound for the ground state energy, and it becomes exact if $`\mathrm{\Psi }`$ has the same nodes as the true ground state wave function, $`\mathrm{\Phi }_{ex}`$. Therefore, for a given trial function, the FNE is more accurate than the variational energy (the expectation value of $`H`$ on $`\mathrm{\Psi }`$), because exactness of the latter requires the stronger condition that $`\mathrm{\Psi }`$ coincide with $`\mathrm{\Phi }_{ex}`$ for all $`R`$. Indeed, the fixed-node approximation offers the best viable numerical solution for several many-fermion problems mitas .

The ability to calculate low-order derivatives of the Born-Oppenheimer, energy with respect to an external parameter, $`\lambda `$, is an important ingredient of the calculation of materials properties. Forces acting on individual atoms or the stress acting within an extended system are examples of first-order derivatives whose vanishing determines the mechanical equilibrium of a system. Similarly, vibrational frequencies—as well as any kind of static generalized susceptibility—are determined by second-order derivatives. The calculation of derivatives of the FNE, $`E/\lambda `$, has the potential of largely extending the scope of applications of quantum Monte Carlo simulations and spawning further methodological developments (besides the examples given above, $`\lambda `$ could also be the strength of any external field, or a variational parameter which modifies the nodal surface of $`\mathrm{\Psi }`$). However technical difficulties hammond  have long hindered routine calculation of FNE derivatives, and despite recent progress zong ; filippi ; assaraf  open questions still remain. In particular, the nodal displacement upon variation of $`\lambda `$ has always been neglected in toto or in part commento , thus adding a further bias on top of the fixed-node approximation: referring to recent literature, $`\lambda `$-independent nodes have been chosen in the path integral Monte Carlo calculation of electronic forces zong ; they have to be chosen so in applications of the Hellmann-Feynman theorem assaraf ; flad ; finally, the drift-diffusion factor of the propagator, together with its nodal contributions, has been dropped from the weights in correlated-sampling quantum Monte Carlo simulations filippi .

In this paper we address the calculation of the derivatives of the FNE without any further approximations. To this end, we examine in detail a simple but significant model that, we believe, is representative of the difficulties which are met in the calculation of these derivatives in more realistic situations. The focus being on the role of nodal displacements, we consider a system where such a role is crucial and a check against exact results is possible, namely a spinless Fermi gas in two dimensions in a weak external field.

## II Methodological aspects and the model

The Hamiltonian of our model system is:

$$H=\underset{i=1}{\overset{N}{}}\frac{1}{2}_i^2+\lambda \underset{i=1}{\overset{N}{}}\mathrm{cos}(𝐪𝐫_i),$$
(1)
where the strength of the external field, $`\lambda `$, is the parameter upon which the FNE is considered to depend, and where periodic boundary conditions have been assumed. The ground state is a Slater determinant, $`D`$, of $`\lambda `$-dependent Mathieu functions gradstein , whereas our trial function is $`\mathrm{\Psi }=J\times D`$ with a symmetric Jastrow factor $`J(R)=\mathrm{exp}(_{i,j}u(r_{ij}))`$ which introduces spurious pair correlations. The fixed-node solution, exact in this case, is recovered as an absolutely non-trivial result of the simulation, since the trial function, wrong on purpose, can be made poor at will by tuning the pseudo-potential, $`u`$.

We use the reptation quantum Monte Carlo (RQMC) method rqmc , in which the calculation of derivatives is simpler to implement than with branching algorithms—just as simple as in path integral Monte Carlo zong . The FNE can be computed from the mixed estimate of the Hamiltonian,

$$E=\frac{\mathrm{\Psi }|H|\mathrm{\Phi }}{\mathrm{\Psi }|\mathrm{\Phi }}=\frac{E_L(R)f(R)𝑑R}{f(R)𝑑R},$$
(2)
by averaging the local energy $`E_L(R)=H\mathrm{\Psi }(R)/\mathrm{\Psi }(R)`$ over the mixed distribution $`f(R)=\mathrm{\Phi }(R)\mathrm{\Psi }(R)`$. The unknown state $`\mathrm{\Phi }`$ is expressed in terms of $`\mathrm{\Psi }`$ through the imaginary-time propagator, $`\mathrm{\Phi }=lim_\tau \mathrm{}e^{H\tau }\mathrm{\Psi }`$. Upon discretization of the imaginary time evolution and introduction of importance sampling, the mixed distribution can be cast into the form:

$$f(R_M)=\underset{M\mathrm{}}{lim}\mathrm{\Pi }_{i=0}^{M1}G(R_i,R_{i+1};ϵ)\mathrm{\Psi }(R_0)^2,$$
(3)
where $`ϵ=\tau /M`$ is the time step. The importance-sampled Green’s function, $`G(R,R^{};ϵ)=\mathrm{\Psi }(R)R|e^{ϵH}|R^{}/\mathrm{\Psi }(R^{})`$, is replaced by a short-time approximation, a common choice being mitas

$$G_F(R,R^{};ϵ)e^{(R^{}RϵF(R)/2)^2/2ϵ}e^{\frac{ϵ}{2}(E_L(R)+E_L(R^{}))},$$
(4)
with $`F=\mathrm{ln}\mathrm{\Psi }^2`$.

In the RQMC method a generalized Metropolis algorithm is used to sample the probability distribution for a path of finite length in imaginary time (a reptile, $`X=\{R_0,\mathrm{},R_M\}`$,

$$P(X)=\mathrm{\Pi }_{i=0}^{M1}G(R_i,R_{i+1};ϵ)\mathrm{\Psi }(R_0)^2,$$
(5)
the limit in Eq. (3) being reached for large enough $`M`$. Replacement of Eq. (3) in (2) shows that the FNE can be evaluated by averaging the local energy on the distribution $`P`$,

$$E=E_L(R_M)\frac{𝑑XE_L(R_M)P(X)}{𝑑XP(X)}.$$
(6)
Since $`P(X)`$ is explicitly known, we can differentiate to get

$`_\lambda E`$ $`=`$ $`A(X),`$ (7)
$`A(X)`$ $`=`$ $`\left[_\lambda E_L(R_M)+(E_L(R_M)E)_\lambda \mathrm{ln}P(X)\right]`$
Finite increments of amplitude $`\mathrm{\Delta }`$ are used for the evaluation of $`_\lambda E_L(R)`$ and $`_\lambda \mathrm{ln}P(X)`$ on the path configurations sampled from $`P(X)`$. Equivalently, finite increments can be used for $`E`$ via the correlated-sampling method:

$`E^{}`$ $`=`$ $`{\displaystyle \frac{𝑑XP^{}(X)E_L^{}(R_M)}{𝑑XP^{}(X)}}`$
$`=`$ $`{\displaystyle \frac{𝑑XP(X)W(X)E_L^{}(R_M)}{𝑑XP(X)W(X)}},`$
with $`W(X)=P^{}(X)/P(X)`$, primed quantities being evaluated at $`\lambda +\mathrm{\Delta }`$.

In the variational Quantum Monte Carlo (VMC) approximation, the probability distribution $`P(X)`$ in Eq. (II) or (II) is replaced by $`\mathrm{\Psi }^2(R_M)`$. Small energy differences are routinely calculated by VMC hammond . The more accurate FNA requires instead the full probability distribution of Eq (5), which entails nontrivial difficulties through the term $`_\lambda \mathrm{ln}P(X)`$ of Eq. (II). Therefore we now examine in some detail the behavior of the probability distribution $`P`$ when a configuration along the path, say $`R_i`$, approaches a node of $`\mathrm{\Psi }`$.

When the distance $`d`$ from the nodal surface is small enough, the dynamics of the system is dominated by the hard wall potential due to the fixed-node boundary conditions nota . The importance-sampled Green’s function should behave like that of a one-dimensional motion close to a hard wall $`[G_0(d,d^{};ϵ)\frac{d^{}}{d}e^{(dd^{})^2/2ϵ}(1e^{2dd^{}/ϵ})]`$, as shown in the left panel of Fig. (1). The short time approximation $`G_F(R,R^{},ϵ)`$, instead, is a very poor approximation near a node, because it assumes $`F=\mathrm{ln}\mathrm{\Psi }^2`$ to be constant during a time step, whereas $`F`$ diverges as $`d^1`$. The approximation (4) for a one-dimensional motion close to a hard wall, shown in the middle panel of Fig. (1), also bears little resemblance to $`G_0`$. First, $`G_F`$ does not vanish at $`d=0`$ where $`G_0`$ does. Sampling from $`G_F`$ would thus imply an unphysical finite density at the nodes; the correct density, which is proportional to $`d^2`$ can be recovered including in the algorithm a Metropolis test based on the square of the trial function fnc ; umrigar . Second, $`G_F`$ vanishes at $`d^{}=0`$ where $`G_0`$ does not. This would imply a divergence proportional to $`d^2`$ of $`\mathrm{ln}P`$ at the nodes, thus making the evaluation of the derivative, Eq. (II), impossible. Such a divergence can be eliminated, leaving the zero-time-step limit unchanged, by a suitable cutoff on the quantum force, $`\mathrm{ln}\mathrm{\Psi }`$. We thus replace $`F`$ in Eq. 4 with $`\overline{F}=4F(1+\sqrt{1+F^2ϵ/2})/(F^2ϵ)`$ umrigar . The resulting Green’s function, $`G_{\overline{F}}`$, is shown in the right panel of Fig. (1). With this replacement the variance of the estimator (II) of the derivative is finite. However we empirically find that the extrapolation of the result at $`ϵ0`$ can be problematic (see also Ref. )

In the simulation of electrons in atoms and molecules, it proves convenient to use an interpolation between Eq. (4) and a locally better Green’s function near the nuclei umrigar . In the same spirit, we use an interpolation $`\stackrel{~}{G}`$ between $`G_{\overline{F}}`$ and a locally better Green’s function $`G_N`$ near the nodes:

$$\stackrel{~}{G}(R,R^{};ϵ)=\alpha (d)G_{\overline{F}}(R,R^{};ϵ)+(1\alpha (d))G_N(R,R^{};ϵ)$$
(9)
where the function $`\alpha (d)`$ smoothly varies from 0 to 1 as $`d`$ goes from 0 to $`ϵ^{1/2}`$. In the new Green’s function

$`G_N(R,R^{};ϵ)`$ $``$ $`{\displaystyle \frac{\mathrm{\Psi }(R^{})}{\mathrm{\Psi }(R)}}e^{(R^{}R)^2/4ϵ}(1e^{d^{}d/ϵ})`$
$`\times `$ $`e^{ϵ(V(R^{})+V(R))/2}`$
the potential $`V`$ is treated in the primitive approximation, whereas the nodal dependence has the correct behavior of $`G_0`$. Although the variance due to $`_\lambda \mathrm{ln}P(X)`$ is finite using $`\stackrel{~}{G}`$, a time-step dependent cutoff on $`d`$ is legitimate and can be useful.

## III Results

We turn now to the discussion of our results for the Fermi gas. We first choose a very poor trial function $`\mathrm{\Psi }_I`$ with a pseudo-potential, $`u=a\mathrm{exp}(br^2)`$ where $`a=0.9`$ and $`b=0.7`$, yielding a variational energy about 20% higher than the exact ground state energy (whereas the FNE and the exact result coincide in this case, since the Slater determinant is the true ground state wave function). The derivative $`_\lambda E|_{x_0}`$, with $`\lambda _0=0.02`$ and $`q=1.585`$, is calculated using $`\stackrel{~}{G}`$. The result, shown in Tab. 1, does indeed recover the correct value.

For comparison, we also perform calculations using the approximation of Ref. filippi , which drops the drift-diffusion part of the propagator in Eq. (II) and introduces an extra factor $`\mathrm{\Psi }^2(R_M)`$ in order to recover the exact result in the limit of perfect importance sampling. As shown in Tab. 1, the approximate algorithm is more efficient (it gives a statistical error 3 times smaller for the same number of path configurations). However its result for the derivative is clearly biased. Using a better trial function $`\mathrm{\Psi }_{II}`$ with $`a=0.45`$ the systematic bias of the approximate algorithm becomes smaller than statistical error. It is interesting to note that the ratio of the statistical error between the two methods does not seem to depend strongly on the quality of $`\mathrm{\Psi }`$. We also have evidence that the efficiency depends on the number of particles only through the extra cost of sampling configurations for a larger system (simulations have been performed for $`5,9`$ and $`21`$ particles).

In Fig. 2 we show the time-step dependence of various derivatives calculated using different Green’s functions. Data obtained with $`G_{\overline{F}}`$ have smaller statistical errors than with $`\stackrel{~}{G}`$, but their extrapolation to zero time-step is harder to guess in advance. The upper panel shows derivatives of the FNE, $`_\lambda E|_{\lambda =0.002}`$. The calculation with $`\stackrel{~}{G}`$ extrapolates to the exact result assuming a quadratic time-step dependence, whereas an $`ϵ^{1/2}`$ term seems necessary in order to account for the $`G_{\overline{F}}`$ data.

The problems with the $`ϵ0`$ limit using $`G_{\overline{F}}`$ are more evident for the derivative of the density fluctuation, $`_\lambda \rho _𝐪|_{\lambda =0}`$, where $`\rho _𝐪(R)=_i\mathrm{exp}(i𝐪𝐫_i)`$: as we can see from the lower panel of Fig. (2), in this case an even smaller power of $`ϵ`$ would be required to extrapolate to the exact result at $`ϵ=0`$.

The physical interest in $`_\lambda \rho _𝐪|_{\lambda =0}`$ stems from its relation to the static dielectric response. The derivative of $`\rho _𝐪`$ is evaluated as $`B(X)`$, with $`B(X)=\rho _𝐪(R_{M/2})+\rho _𝐪(R_{M/2})_\lambda \mathrm{ln}P(X)`$ (cfr. Eq. (II) notaV ). The same physical information can be formally obtained from $`_\lambda ^2E|_{\lambda =0}`$. Second derivatives are in principle straightforward to evaluate, but differentiation of Eq. (II) leads to terms containing $`_\lambda ^2\mathrm{ln}P`$ and $`(_\lambda \mathrm{ln}P)^2`$ which have large fluctuations.

Despite its low efficiency, the calculation of second derivatives of the energy offers an interesting conceptual advantage over the first derivative of the density fluctuation. Suppose we don’t know the exact nodal dependence of $`\mathrm{\Psi }`$ on the strength $`\lambda `$ of the external field. Starting from plane-wave nodes at $`\lambda =0`$, we can explore a family of wave functions, $`\mathrm{\Psi }_\alpha `$, in which the one-particle orbitals are Mathieu functions calculated for an external potential of strength $`\alpha \lambda `$ different from the physical value $`\lambda `$. The nodal dependence of $`\mathrm{\Psi }_\alpha `$ on $`\lambda `$ varies with $`\alpha `$: for $`\alpha =0`$, no nodal displacements are allowed; for $`\alpha =1`$, the exact nodal dependence is recovered. We can compute within the same run (hence in a correlated way) the values of $`_\lambda ^2E|_{\lambda =0}^{(\alpha )}`$ pertaining to several choices of $`\alpha `$. The result is shown in Fig. (3). Because of the variational character of the fixed node approximation, the optimal nodal dependence (within the family $`\mathrm{\Psi }_\alpha `$, and within error bars) is given by the value $`\alpha `$ corresponding to the minimum of $`_\lambda ^2E|_{\lambda =0}^{(\alpha )}`$. Inspection of Fig. (3) leads us to conclude that: (i) suppression of the nodal dependence ($`\alpha =0`$) leads to a small statistical error, but produces a large bias; (ii) the minimum of $`_\lambda ^2E|_{\lambda =0}^{(\alpha )}`$ is consistent with the exact nodal dependence ($`\alpha =1`$); (iii) the dependence of $`_\lambda ^2E|_{\lambda =0}^{(\alpha )}`$ on $`\alpha `$ is quadratic, so that calculations for just 3 values of $`\alpha `$ would suffice (this could be useful for a higher-dimensional optimization).

Unfortunately, the calculation of second derivatives has still very large statistical errors (if nodal displacement is properly included). We hope that a successful implementation of noise reduction techniques caffarel  could eventually turn it into a convenient computational tool. Further perspectives in nodal optimization are offered by the calculation of first derivatives of the FNE, for instance in conjunction with the stochastic gradient approximation method sga .


text:
# Classical nonlocal models for states of a modified quantized Klein-Gordon field

## Abstract

This paper is withdrawn.