text:
# Quantum Extension of the Jarzynski Relation; Analogy with Stochastic Dephasing

## I Introduction

Jarzynski 1 ; 3  had established a remarkably simple and general relationship between the distribution of work performed on a classical system by an external force and the free energy difference between the initial and final states 2 ; 4 ; 6 . This interesting prediction, connecting equilibrium quantities (free energies) to nonequilibrium trajectories corresponding to measurements performed on an arbitrary timescale, has recently been verified in single molecule measurements 5 .

In this letter we utilize the adiabatic representation 9  to prove that this relation holds equally for quantum systems. For an isolated driven quantum system (i.e., not coupled to a bath) the work does not depend on the path, and the proof follows directly from a sum rule of nonadiabatic couplings. We further develop a path integral representation of systems coupled to a bath described by a master equation in the adiabatic basis using “trajectories” over adiabatic state populations, and generalize Jarzynski’s relation to obtain the joint distribution of work and energy in a nonequillibrium measurement. A close resemblance is established between the physical picture underlying the distribution of work and the distribution of phase acquired by a quantum coherence due to coupling to a bath.

Consider a quantum system coupled to an external force which changes its Hamiltonian $`H(t)`$ by switching a parameter $`\lambda (t)`$ from $`\lambda =0`$ at $`t=0`$ to $`\lambda =1`$ at time $`t`$. Initially the system is in thermal equilibrium and the instantaneous energy levels (eigenvalues of $`H(t)`$) are denoted $`ϵ_n(t)`$ with the corresponding eigenstates $`|\phi _n(t)`$. The canonical partition function for a given value of $`\lambda `$ is

$$Z_\lambda =\underset{n}{}exp[\beta ϵ_n(t)]exp(\beta F_\lambda )$$
(1)
where $`F_\lambda `$ is the Helmholtz free energy.

Let us consider an externally driven, but otherwise isolated, system. In the adiabatic representation 9  we expand its density matrix as

$$\rho (t)=\underset{nm}{}\rho _{nm}(t)|\phi _n(t)><\phi _m(t)|,$$
(2)
where the coefficients $`\rho _{nm}(t)`$ satisfy the Liouville equation

$$\frac{d}{dt}\rho _{kl}(t)=i\omega _{kl}(t)\rho _{kl}(t)\underset{mn}{}S_{kl,mn}(t)\rho _{mn}(t)$$
(3)
with $`\omega _{kl}(t)ϵ_k(t)ϵ_l(t)`$ and

$$S_{kl,mn}(t)=\phi _k(t)|\dot{\phi }_m(t)\delta _{ln}+\dot{\phi }_n(t)|\phi _l(t)\delta _{km}$$
(4)
The nonadiabatic coupling $`S`$ is of purely quantum origin and has no classical analogue.

For a given realization of $`\lambda (t)`$, the solution of Eq. (3) allows us to compute the conditional probability $`K_{mn}(t)`$ of the system to be in the state $`|\phi _m(t)`$ at time $`t`$ given that it started at state $`|\phi _n(0)`$, so that $`\rho _{mm}(t)=K_{mn}(t)\rho _{nn}(0)`$. This probability is normalized by summing over *final* states $`_mK_{mn}(t)=1`$. A second, less obvious, sum rule is obtained by summation over *initial* states $`_nK_{mn}(t)=1`$. When the external force varies slowly, it follows from the adiabatic theorem 9  that no nonadiabatic transitions take place, $`K_{mn}(t)=\delta _{mn}`$, and this sum rule trivially holds. For finite switching timescales there are nonadiabatic transitions and the instantaneous eigenstates $`|\phi _n(t)`$ are no longer the solutions of the time dependent Schrödinger equation. Nevertheless, the second sum rule follows from the relation

$$\underset{n}{}S_{kl,nn}(t)=\frac{d}{dt}\phi _k(t)|\phi _l(t)=0.$$
(5)
A direct consequence of Eq. (5) is that the following uniform distribution in adiabatic population space, obtained by setting $`\rho _{nm}(t)=\delta _{nm}`$ at all times in Eq. (2)

$$\rho (t)=\underset{n}{}|\phi _n(t)\phi _n(t)|,$$
(6)
satisfies the Liouville equation. This interesting identity gives rise to the second sum rule, which holds despite the nonadiabatic transitions and is the key for proving the quantum Jarzynski identity. Alternatively, Eq. (6) can be viewed as a unitary transformation of the unit operator between adiabatic basis sets at different times. This equation then simply states that the unit operator is invariant to a unitary transformation.

We shall now compute the ensemble average of $`exp(\beta W)`$ where $`W`$ is the cumulative work performed up to time $`t`$. If the system starts in state $`|\phi _n(0)`$ and ends up at time $`t`$ in the state $`|\phi _m(t)`$, then the work made by the external force is $`W=ϵ_m(t)ϵ_n(0)`$, and we have

$$exp(\beta W)=\frac{1}{Z_0}\underset{mn}{}exp[\beta ϵ_n(0)]K_{mn}(t)exp[\beta (ϵ_m(t)ϵ_n(0))].$$
(7)
This gives

$$exp(\beta W)=\frac{1}{Z_0}\underset{mn}{}K_{mn}(t)exp[\beta ϵ_m(t)].$$
(8)
Using the second sum rule we can carry out the $`n`$ summation and obtain Jarzynski relation,

$$exp(\beta W)=Z_1/Z_0$$
(9)
where $`Z_1`$ is the final state partition function for $`\lambda =1`$. Eq. (8) implies that it is generally possible to compute the free energy change $`F_1F_0\beta ^1\mathrm{log}(Z_1/Z_0)`$ using the distribution of work.

We next turn to a system coupled to a bath and described by a master equation in the adiabatic basis 3 ; 6 . We denote the transition rate from state $`n^{}`$ to $`n`$ by $`R_{nn^{}}(t)`$. $`R`$ depends on $`t`$ since it is recast in the time dependent, adiabatic, basis. It further satisfies the detailed balance condition 13

$$R_{nn^{}}(t)/R_{n^{}n}(t)=exp(\beta \omega _{nn^{}}(t))$$
(10)
For an isolated system the work equals the energy change and is therefore a state function. When the system is not isolated, the work becomes path dependent and may not be computed from the initial and final states alone, as was done in Eq. (7). For our model, the change in energy along the path stems from the variation of the eigenvalues $`ϵ_n(t)`$ which is induced by the driving force, as well as the changes in their occupations, which are induced by the coupling to the bath (Fig. 1). The former is the work whereas the latter is the heat. Both work and heat are path dependent whereas the overall energy change $`ϵ_m(t)ϵ_n(0)`$ is a state function.

We shall compute the joint distribution of work performed $`(W)`$ and the change in system’s energy $`(E)`$ in the process by defining the joint energy-work generating function

$$S(\gamma ,\delta ;t)=exp[\gamma W(t)\delta E(t)],$$
(11)
where $`\mathrm{}.`$ denotes an ensemble average over trajectories for a fixed realization of $`\lambda (t)`$. The moments of $`W`$ and $`E`$ can then be computed as derivatives with respect to the parameters $`\gamma `$ and $`\delta `$

$$W(t)^pE^q(t)=(1)^{p+q}\frac{^p}{\gamma ^p}\frac{^q}{\delta ^q}S(\gamma ,\delta ;t)|_{\gamma =\delta =0}$$
(12)
By expanding the solution of the master equation perturbatively in the off diagonal elements of $`R_{nn^{}}(nn^{})`$ it can be represented as a path integral in the adiabatic population space (see Fig. 1): Consider a family of trajectories which start at state $`|\phi _n(0)`$ at $`t=0`$ and end in state $`|\phi _m(t)`$ at time $`t`$. In the course of a given trajectory the system had jumped $`j`$ times between different adiabatic states with $`j=0,1,2\mathrm{}.`$. The conditional probability of the system to start at state $`|\phi _n(0)`$ and end in the state $`|\phi _m(t)`$ is then given by a sum over all possible values of $`j`$. The work can be computed by integrating $`ϵ_n/t`$ along the continuous segments of each trajectory. The contribution of this family to the generating function can be calculated by solving the equation of motion for the following Green function

$$\frac{dG_{mn}(t)}{dt}=\underset{n^{}}{}R_{mn^{}}(t)G_{nn}(t)\gamma \frac{ϵ_m}{t}G_{mn}(t)$$
(13)
with the initial condition $`G_{mn}(0)=\delta _{mn}`$. It immediately follows from the path integral representation of the solution of the master equation that the generating function is given by averaging the Green function over the equilibrium distribution of initial states and a weighted summation over final states

$$S(\gamma ,\delta ;t)=\underset{nm}{}exp[\delta [ϵ_m(t)ϵ_n(0)]]G_{mn}(t;\gamma )\frac{exp(\beta ϵ_n(0))}{Z_0}.$$
(14)
Eq. (14) provides the following physical picture for the process: We consider an ensemble of trajectories over adiabatic state populations. In the classical case 1  the ensemble is generated by the distribution of initial states, whereas here the ensemble originates from the random nature of the stochastic jumps among adiabatic states. The Green function represents an open system where population can flow to (or from) state $`m`$ with a rate $`\gamma ϵ_m/t`$ which depends on the variation of the adiabatic energies with time and on the parameter $`\gamma `$. As a result, the total population $`_mG_{mn}(t)`$ is not conserved and its value determines the generating function for the work and energy. The adiabatic representation thus provides a “classical trajectory” picture for the joint distribution of work and energy. To recover the Jarzynski relation we set $`\gamma =\beta `$ and $`\delta =0`$ in Eq. (14) which then assumes the form

$$S(\beta ;0,t)=\underset{m}{}\text{G}_m(\beta ,t)$$
(15)
where

$$\text{G}_m(\beta ,t)=\underset{n}{}G_{mn}(\beta ,t)\frac{exp(\beta ϵ_n(0))}{Z_0}$$
(16)
Let us consider the (unnormalized) distribution of states $`m`$

$$\text{G}_m(\beta ,t)=\frac{exp(\beta ϵ_m(t))}{Z_0}$$
(17)
It can be easily verified that Eq. (17) is a solution to Eq. (13); this follows from the detailed balance of $`R`$ that guarantees that the first term in the r.h.s. of Eq. (13) vanishes when acting on Eq. (17). Substituting Eq. (17) in Eq. (15) immediately gives the Jarzynski relation (Eq. (9)). It will be interesting to extend these results to systems described by more elaborate kinetic schemes such as fractional kinetics 10  rather than the ordinary master equation.

The validity of the Jarzynski relation for a system coupled to a bath at constant temperature has been established by Jarzynski 3  and Crooks 6 . The present derivation and the definition of work closely resemble Jarzynski’s master equation approach. However, the path integral picture of continuous segments with discrete jumps provides a new insight for quantum systems and helps establish an interesting analogy with the theory of pure dephasing of quantum coherence, as will be shown next.

Eq. (13) closely resembles the stochastic Liouville equation of Kubo 7 ; 8 ; 12 , which has been widely used in the theory of NMR and optical lineshapes as well as in single molecule spectroscopy 14 . Computing the work generating function is then formally equivalent to computing the dephasing of a two level system ($`g`$ and $`e`$) coupled to a bath. In that problem we consider equilibrium fluctuations of a quantum coherence $`\sigma `$ between the two levels which are coupled to a classical bath represented by a phase space $`\mathrm{\Gamma }`$. $`\sigma `$ satisfies the stochastic Liouville equation 7

$$\frac{d\sigma (\mathrm{\Gamma },t)}{dt}=\underset{\mathrm{\Gamma }^{}}{}W(\mathrm{\Gamma },\mathrm{\Gamma }^{})\sigma (\mathrm{\Gamma }^{},t)iU(\mathrm{\Gamma })\sigma (\mathrm{\Gamma },t).$$
(18)
Here $`W`$ is a Markovian rate matrix and $`U(\mathrm{\Gamma })`$ is the fluctuating energy gap between the two levels. $`\sigma `$ serves as a generating function for the absorption lineshape which is given by the Fourier transform of its zero’th moment $`J(t)=\sigma (\mathrm{\Gamma },t)𝑑\mathrm{\Gamma }`$. The normalization of $`\sigma `$ is thus the physically interesting quantity whose Fourier transform gives the spectral lineshape.

Eq. (13) and Eq. (18) describe very different physical phenomena: The former corresponds to an externally driven nonequilibrium system whereas the latter represents a non driven system at equilibrium. However, the two equations have a close formal connection. To better see the analogy we recast the solution of Eq. (18) as

$$\sigma (\mathrm{\Gamma },t)=exp[i_0^t𝑑\tau U(\tau )]$$
(19)
where the averaging is over a subensemble of trajectories which assume the same value $`\mathrm{\Gamma }`$ at time $`t`$. For comparison we write the solution of Eq. (13) as

$$G_m(\gamma ,t)=exp\left[\gamma _0^t𝑑\tau ϵ/\tau \right].$$
(20)
The average here is over a subensemble of trajectories (Fig. 1) which start at thermal equilibrium at $`t=0`$ and end up in state $`m`$ at time $`t`$. We reiterate that $`ϵ/\tau `$ only corresponds to the continuous part of the trajectory; the jumps should be excluded.

In both equations we compute a phase space distribution whose time dependent normalization gives the generating function for the quantity of interest (either work or absorption spectrum). The first term in the r.h.s. of Eq. (13)) or (18) represents a regular dynamics which conserves the normalization, whereas the second term is responsible for varying the normalization, which is the physically interesting property. In Eq. (13) or Eq. (20) this comes from the flow of adiabatic state populations in and out of the system, originating from the time dependence of adiabatic energies. In the dephasing problem (Eq. (18) or (19)), the coherence $`\sigma `$ acquires a random phase $`(U)`$ which is different for different trajectories and it thus constitutes path function. This translates to loss of magnitude (dephasing and decoherence) upon ensemble averaging. The work in the nonadiabatic dynamics is related to an integral over the continuous part of the energy trajectory in the same way that the phase in the lineshape problem is obtained by an integral over the fluctuating two level frequency. This analogy suggests that it should be possible to carry out nonequilibrium lineshape measurements by e.g. switching on an electric field for a solute in a polar solvent and observing the spectroscopic analogues of the Jarzynski relations by looking, for example, at the time dependent Stokes Shift.

Acknowledgement:

The support of the National Science Foundation grant no. CHE-0132571 is gratefully acknowledged. I wish to thank the referee for most useful comments that helped clarify the generality of these results.


format: text

text:
# 1 Introduction

## 1 Introduction

Calculations of ultra-high energy (UHE) neutrino-nucleon cross sections $`\sigma ^{\nu N}`$  have attracted attention for many years since, on the one hand, UHE neutrino fluxes are predicted by cosmological models but have not yet been measured experimentally and, on the other hand, theoretical models exist for the evaluation of these cross sections.

Neutrino telescopes are the most promising tools for probing the distant stars and galaxies. Due to their high stability and neutrality, neutrinos arrive at a detector on a direct line from their sources, undeflected by intervening magnetic fields, with expected energy up to $`10^{20}`$ eV. Whereas high-energy photons are completely absorbed by a few hundred grams/cm2 of material, the interaction length of a 1 TeV neutrino is about 250 kilotonnes/cm2 .

An important scientific goal of large-scale neutrino telescopes is the detection of UHE cosmic neutrinos produced outside the atmosphere, i.e. neutrinos produced by galactic cosmic rays interacting with interstellar gas, and extra-galactic neutrinos. High neutrino energy brings a number of advantages . Firstly, the charged-current cross section increases as $`\sigma E_\nu `$ for $`E_\nu <10^{12}`$ eV and is believed to rise as $`\sigma E_\nu ^{0.4}`$ for higher $`E_\nu `$. (In this paper we question the last point.) Secondly, the background of atmospheric neutrinos decreases as compared to the flux from extra-galactic sources, approximately as $`E_\nu ^{1.6}.`$

Another important and promising aspect of UHE neutrino physics is the possibility to study the production of exotic objects, such as supersymmetric particles  or black holes .

The phenomenology of UHE neutrinos and their detection depend on the neutrino-nucleon total cross section, which has been calculated in the standard model (see e.g. ). A striking feature of these calculations is the continued power-law growth of the cross sections with $`E_\nu `$ at highest energies. The rise with $`E_\nu `$, on the other hand, is directly related to the very low-$`x`$ behavior of the nucleon parton distributions derived from e.g. the HERA data. Problems related to the unitarity bound have been discussed in .

At ultra-high neutrino energies the total cross section is completely dominated by its deep inelastic component. The energy dependence of neutrino-nucleon cross section, as suggested from the small-$`x`$ nucleon structure functions (SF) measured at HERA, should show a rapid increase with energy, like $`s^\lambda `$ , where $`\lambda `$ scales as $`\mathrm{ln}(Q^2/\mathrm{\Lambda }^2)`$ and ranges approximately in the interval $`0.1÷0.4`$ for $`1`$ GeV$`{}_{}{}^{2}<Q^2<100`$ GeV2 (see Ref. ).

The highest available $`ep`$ energy at HERA, $`\sqrt{s_{ep}}=314`$ GeV, is still too low as compared to the expected UHE neutrino-nucleon collision energies. For the neutrino energy $`E_\nu 10^{12}`$ GeV, the relevant $`x`$-region is $`x10^8`$, while the HERA measurements extend now to $`x10^6`$; however such low values of $`x`$ are measured at $`Q^2<0.1`$ GeV2. For the interesting region of $`Q^2M_W^2`$, SF are measured only up to $`x10^3`$ in the D0 experiment of inclusive jets . Since studies of the UHE $`\nu N`$ cross sections began nearly a decade ago  our knowledge about quark distributions and nucleon SF improved significantly .

In the present paper we study the UHE behavior of the neutrino-nucleon charged-current cross section starting from models where nucleon SF have a milder increase. Our analysis relies strongly on recent progress and our previous experience in the studies of the nucleon SF , based on phenomenological fits, extended by QCD evolution.

The paper is organized as follows: in the next Section we briefly review the steps involved in the determination of the neutrino-nucleon cross section and illustrate our approach; in Section 3 we study the asymptotic behavior of the total neutrino-nucleon cross section; in Section 4 we discuss the implications of our results for the observation of UHE neutrinos in future experiments and finally in Section 5 we draw some conclusions.

## 2 Nucleon structure functions and neutrino-nucleon cross sections

In order to compute the cross section for neutrino-nucleon interactions at high energies, one needs a detailed description of the quark structure of the nucleon. The special challenge of the UHE neutrino-nucleon cross section is that the $`W`$-boson propagator emphasizes smaller and smaller values of $`x`$ as the neutrino energy $`E_\nu `$ increases. In the UHE domain, the most important contribution to the $`\nu N`$ cross section comes from $`xM_W^2/(2ME_\nu )`$ and hence the greatest uncertainty in the predictions comes from the small-$`x`$ extrapolation.

Deep inelastic scattering of UHE neutrinos is unique in the sense that it explores extreme regions of the $`(Q^2,x)`$ phase space where no accelerators data exist. One thus must rely on predictions or assumptions regarding the nature of these distributions, in particular in the region of ultra-low $`x`$ and low $`Q^2`$, as well as in the region of low $`x`$ and low $`Q^2.`$

Let us consider the double differential cross section for the deep inelastic scattering of neutrinos on an isoscalar target for charged-current processes. We will limit ourselves to the leading order corrections to the simple parton model and hence all parton model formulas remain unchanged except that the parton distributions now depend on $`x`$ and $`Q^2`$ and not only on $`x`$. In particular, the Callan-Gross relation, $`F_2=2xF_1`$ or $`F_L=0`$, holds in leading order . With the usual notation we can then write for charged-current neutrino interactions

$`\left({\displaystyle \frac{d\sigma }{dxdy}}\right)^\nu `$ $`=`$ $`{\displaystyle \frac{G_F^2ME_\nu }{\pi }}\left({\displaystyle \frac{M_W^2}{Q^2+M_W^2}}\right)^2`$ (1)
$`\times `$ $`\left\{{\displaystyle \frac{1+(1y)^2}{2}}F_2^\nu (x,Q^2)+\left(1{\displaystyle \frac{y}{2}}\right)yxF_3^\nu (x,Q^2)\right\}.`$
For anti-neutrino charged-current processes one must change the sign in front of $`F_3^\nu (x,Q^2)`$. (The changes necessary in order to describe neutral-current neutrino interactions are given in many textbooks and review articles, for example in ). In Eq. (1), $`G_F`$ is the Fermi constant, $`M`$ is the nucleon mass and the variable $`y`$ is related to the Bjorken $`x`$ through the relation

$$y=\frac{Q^2}{x(sM^2)}\frac{Q^2}{xs},$$
(2)
where $`s`$ is the square of the c.m. energy for the neutrino-nucleon scattering. The laboratory neutrino energy $`E_\nu =(sM^2)/(2M)`$ is approximately $`s/(2M)`$ in the energy region of interest.

The main purpose of this paper is to evaluate the asymptotic behavior of the total cross section for charged-current neutrino-nucleon process

$$\sigma ^{\nu N(CC)}=_0^1𝑑x_0^1𝑑y\left(\frac{d\sigma }{dxdy}\right)^\nu $$
(3)
in a model for the structure function $`F_2`$ based on a soft non-perturbative input complemented by an approximate QCD evolution. The contribution of $`xF_3`$ can in fact be neglected at ultra-high energies. Unfortunately, experimental data for the neutrino SF do not extend to the “deep sea” region, i.e. to $`x`$ smaller than 0.001, and to large $`Q^2`$. Both these limits, for the electromagnetic structure functions $`F^{ep}`$, are subject of debate regarding the presence of possible saturation effects at small $`x`$ and discrepancies in the $`Q^2`$ evolution for large $`Q^2`$. $`F_2^\nu `$ differs from $`F_2^{ep}`$ at small $`Q^2`$ as can be easily seen by formally integrating the double differential cross section (1) over $`x`$ at fixed value of the product $`(E_\nu y)`$

$$\underset{y0}{lim}\frac{1}{E_\nu }\frac{d\sigma }{dy}=\frac{G_F^2M}{\pi }_0^1𝑑xF_2^\nu (x,Q^2=0),$$
(4)
formula that has been used for the normalization of the neutrino flux . The integral in Eq. (4), when evaluated for $`F_2^{ep}`$, would vanish by gauge invariance. At increasing $`Q^2`$, this difference will disappear but, at the moment, no smooth transition has been found in all the range of the variables appearing in Eq. (3).

### 2.1 The non-perturbative input

In a number of papers  a diffraction model has been developed, describing equally well hadronic reactions at high energies and DIS at small $`x`$. The main idea behind the model is that SF at small and moderate values of $`Q^2`$ are Regge-behaved. The $`Q^2`$-dependence was introduced in a phenomenological way.

It was shown  that data on deep inelastic $`ep`$ scattering from HERA and elsewhere can be fitted perfectly well by a “soft” Pomeron alone. This Pomeron could be a dipole with a $`\mathrm{ln}(1/x)`$ dependence or a squared logarithm, saturating asymptotically the Froissart bound. For the electromagnetic SF, parametrizations exist  that allow for a simplified evolution until very large $`Q^2`$, while maintaining a soft dependence on the Bjorken $`x`$. As an example we will consider a structure function $`F_2^\nu (x,Q^2)`$ of the form

$$F_2^\nu (x,Q^2)=\left[a(Q^2)+b(Q^2)\mathrm{ln}\left(\frac{1}{x}\right)\right](1x)^{\nu (Q^2)}$$
(5)
and derive for it an explicit approximate evolution starting from $`Q^2=Q_0^2`$. The main motivation for our choice of the model is that at small $`Q^2`$ there is similarity between the small-$`x`$ behavior of the structure function and the $`s`$-dependence, $`xQ^2/s`$, of the hadronic total cross section at high energies. The alternative approach, with a power behaved dependence, has been discussed in the literature in several publications. The differences between the two approaches can be significant and will be discussed in Section 3.

We notice that the coefficients $`a(0)`$, $`b(0)`$ and $`\nu (0)`$ are constrained by Eq. (4) whose left hand side can be determined from the experiment, Ref. . The constraint is

$$\underset{y0}{lim}\frac{1}{E_\nu }\frac{d\sigma }{dy}=\frac{1}{\nu (0)+1}\left\{a(0)+b(0)\left(\gamma +\psi [\nu (0)+2]\right)\right\},$$

where $`\gamma `$ denotes the Euler’s constant and $`\psi [z]`$ is the logarithmic derivative of the $`\mathrm{\Gamma }`$-function. When low-$`Q^2`$ data  will extend to lower $`x`$ values, it will be possible to check the choice (5).

By imposing a lower cut in $`Q^2`$, $`Q^2=Q_0^2`$, and rewriting Eq. (3) as

$$\sigma ^{\nu N(CC)}\frac{1}{2ME_\nu }_{Q_0^2}^s𝑑Q^2_{Q^2/s}^1\frac{dx}{x}\left(\frac{d\sigma }{dxdy}\right)^\nu ,$$

where in the differential cross section $`y=Q^2/(xs)`$, the $`F_2`$ contribution to the total cross section can be written in the form

$$\overline{\sigma }^{\nu N}\frac{\sigma ^{\nu N}+\sigma ^{\overline{\nu }N}}{2}$$

$$=\frac{G_F^2}{2\pi }_{Q_0^2}^s𝑑Q^2\left(\frac{M_W^2}{Q^2+M_W^2}\right)^2_{Q^2/s}^1\frac{dx}{x}\frac{1+(1Q^2/(xs))^2}{2}F_2^\nu (x,Q^2).$$
(6)
### 2.2 $`Q^2`$ evolution

Let us consider a simplified approach to the DGLAP  evolution, where the parton distribution functions satisfy the relation $`q(x)=\overline{q}(x)=u(x)=d(x)=s(x)=2c(x)=2b(x)`$ and $`F_2^\nu `$, for five flavors, is given by $`F_2^\nu =8xq(x)`$. This assumption gives $`F_2^{ep}=17xq(x)/9`$ and, as shown in , with this simple form for the SF a high quality fit to the HERA data can be obtained.

For each parton distribution we assume the same function (5) chosen for $`F_2^\nu (x,Q^2)`$

$$xf_i(x,Q^2)=\left[a_i(Q^2)+b_i(Q^2)\mathrm{ln}\left(\frac{1}{x}\right)\right](1x)^{\nu (Q^2)},(i=q,g)$$
(7)
with the same functional form for quark and gluon but with different coefficients. Experimental data for quarks and global analysis, for example the parton analysis of MRST , for the gluon could determine the coefficients in the expressions for $`xf_i(x,Q^2)`$, where $`i=q,g`$ stands for quark or gluon, respectively.

In the leading order of perturbation theory, an approximate evolution is feasible by using the results of Ref. . The only condition is that $`a_ib_i`$ $`(i=q,g)`$, where $`a_ia_i(Q_0^2)`$ and $`b_ib_i(Q_0^2)`$, so that no interference appears in the $`Q^2`$ evolution of the coefficients multiplying the different powers of the logarithm111The requirement $`a_ib_i`$ is needed for all $`Q^2`$ values. It drastically simplifies our formulas since it avoids a mixing between $`a_i`$ and $`b_i`$ parameters during the $`Q^2`$ evolution.. We define $`\beta _0=112f/3`$, where $`f=5`$ is the number of flavors in our approach,

$$t=\mathrm{ln}\left[\frac{\alpha (Q_0^2)}{\alpha (Q^2)}\right]=\mathrm{ln}\left[\frac{\mathrm{ln}(Q^2/\mathrm{\Lambda }^2)}{\mathrm{ln}(Q_0^2/\mathrm{\Lambda }^2)}\right]$$

and recall the form of the scaling variables,

$$\sigma =2\sqrt{\widehat{d}_{gg}t\mathrm{ln}(1/x)},\rho =\sqrt{\frac{\widehat{d}_{gg}t}{\mathrm{ln}(1/x)}}=\frac{\sigma }{2\mathrm{ln}(1/x)},$$
(8)
where $`\widehat{d}_{gg}=12/\beta _0`$ 222We follow here the notation of : $`d`$ is the ratio between the anomalous dimension and twice the QCD $`\beta `$-function in leading order, while $`\widehat{d}`$ and $`\overline{d}`$ are, respectively, the singular and regular part of $`d`$ when $`n1`$. . Then, from , we get for small $`x`$ the evolution of the structure function $`F_2(x,Q^2)=8xf_q(x,Q^2)=8x[f_q^{}(x,Q^2)+f_q^+(x,Q^2)]`$ as

$`xf_q^{}(x,Q^2)`$ $`=`$ $`\left(a_q+b_q\left[\mathrm{ln}\left({\displaystyle \frac{1}{x}}\right)C_1(\nu )\right]\right)e^{d_{}(1)t}+O(x)`$ (9)
$`xf_q^+(x,Q^2)`$ $`=`$ $`\left(a_q^+\rho I_1(\sigma )+b_q^+I_0(\sigma )\right)e^{\overline{d}_+(1)t}+O(\rho ),`$ (10)
where $`d_{}(1)=16f/(27\beta _0),\overline{d}_+(1)=1+20f/(27\beta _0)`$, $`a_q^+=f(a_g+4a_q/9)/9`$ and an analogous expression for $`b_q^+`$. $`C_1(\nu )`$ is related to the logarithmic derivative of the $`\mathrm{\Gamma }`$-function: $`C_1(\nu )=\psi [1+\nu ]\psi [1]`$. $`I_0(z)`$ and $`I_1(z)`$ are modified Bessel functions of the first kind.

We notice that the presence of higher powers of the logarithm of $`1/x`$ can be easily handled in this formalism provided that $`a_ib_i\mathrm{}`$. While the evolution given by Eqs. (9)-(10) is correct in the limit of very large $`\sigma `$ and $`Q^2Q_0^2`$, the requirement that the initial input must be reproduced for $`Q^2=Q_0^2`$ introduces a corrective term in $`f_q^+`$ that does not affect the evolution for large $`Q^2`$. Eq. (10) should then be replaced by

$$xf_q^+(x,Q^2)=\left(a_q^+\rho I_1(\sigma )+b_q^+\left[I_0(\sigma )+\mathrm{\Delta }\right]\right)e^{\overline{d}_+(1)t}+O(\rho ),$$
(11)
where $`\mathrm{\Delta }=b_qC_1(\nu )/b_q^+1`$.

## 3 Asymptotic behavior of $`\sigma ^{\nu N}`$

We rewrite Eq. (6) in the form

$$\overline{\sigma }^{\nu N}=\frac{2G_F^2}{\pi }_{Q_0^2}^s𝑑Q^2\left(\frac{M_W^2}{Q^2+M_W^2}\right)^2$$

$$\times _{Q^2/s}^1\frac{dx}{x}[1+(1\frac{Q^2}{xs})^2]\{xf_q^{}(x,Q^2)+xf_q^+(x,Q^2)\},$$
(12)
where $`xf_q^{}(x,Q^2)`$ are given in Eqs. (9)-(10) respectively. The integral over $`x`$ can be done exactly for $`f_q^{}(x,Q^2)`$. Let us denote this integral with the symbol $`J_{}(s,Q^2)`$. The result is presented in Eq. (A.1) of the Appendix, and the leading behavior, for large $`s`$, of the $`x`$-integral is given by $`\mathrm{ln}^2(s/Q^2)`$.

The presence of many scales $`(\mathrm{\Lambda }^2,Q_0^2,M_W^2\text{and}s)`$ makes difficult any precise estimate of the asymptotic behavior of the $`xf_q^{}(x,Q^2)`$ contribution to $`\overline{\sigma }^{\nu N}`$. From the Ostrowski’s inequality  for integrals it follows, however, that this contribution has an upper bound for $`s\mathrm{}`$:

$$\left|_{Q_0^2}^s𝑑Q^2\left(\frac{M_W^2}{Q^2+M_W^2}\right)^2J_{}(s,Q^2)\right|K\mathrm{ln}^2\left(\frac{s}{Q_0^2}\right),$$
(13)
where $`K`$ is a suitable constant.

Only the small-$`x`$ evolution of $`f_q^+(x,Q^2)`$ is known from (10) and the condition $`\rho 1`$ requires a cutoff at the upper limit of the $`x`$-integral. Let us call it $`x_M`$ with $`x_M<1`$, while the $`Q^2`$-integral has now $`x_Ms`$ as the upper limit of integration. The first integral we meet, when considering the term with $`f_q^+(x,Q^2)`$, is

$$_{Q^2/s}^{x_M}\frac{dx}{x}xf_q^+(x,Q^2)=e^{\overline{d}_+(1)t}\left(a_q^+_1+b_q^+_2\right)$$
(14)
and, with the change of variable $`x=\mathrm{exp}[z^2/(4\widehat{d}_{gg}t)]`$, the solution can be easily found. Setting $`\sigma _u=2\sqrt{\widehat{d}_{gg}t\mathrm{ln}(s/Q^2)}`$ and $`\sigma _{\mathrm{}}=2\sqrt{\widehat{d}_{gg}t\mathrm{ln}(1/x_M)}`$, we obtain:

$$_1(s,Q^2)=_\sigma _{\mathrm{}}^{\sigma _u}𝑑zI_1(z)=I_0(\sigma _u)I_0(\sigma _{\mathrm{}})$$
(15)
and

$$_2(s,Q^2)=\frac{1}{2\widehat{d}_{gg}t}_\sigma _{\mathrm{}}^{\sigma _u}𝑑zzI_0(z)=\frac{1}{2\widehat{d}_{gg}t}\left(\sigma _uI_1(\sigma _u)\sigma _{\mathrm{}}I_1(\sigma _{\mathrm{}})\right).$$
(16)
Other integrals, coming from powers of $`Q^2/(xs)`$ in Eq. (12), are $`O(\rho )`$ with respect to those in Eqs. (15) and (16) and it is consistent with our approximation in Eq. (10) to neglect them. The proof can be found in the Appendix.

We define now the integrals

$$𝒥_k(s)=_{Q_0^2}^{x_Ms}𝑑Q^2\left(\frac{M_W^2}{Q^2+M_W^2}\right)^2e^{\overline{d}_+(1)t}_k(s,Q^2),(k=1,2)$$
(17)
and notice that their large-$`s`$ behavior determines completely the $`f_q^+`$ contribution to the mean total cross section in this limit. Of the two integrals $`𝒥_1(s),𝒥_2(s)`$ only the first one needs to be evaluated explicitly since

$$\frac{d𝒥_2(s)}{d\mathrm{ln}s}=𝒥_1(s)+O\left(\frac{1}{\mathrm{ln}s}\right).$$

We consider then the integral $`𝒥_1(s)`$. The function under the integral sign in $`𝒥_1(s)`$ vanishes at both integration limits and has a bell-shaped form with a well defined maximum. The maximum moves at larger values of $`Q^2`$, when $`s`$ increases, and its position is defined, in the limit $`s\mathrm{}`$, by the solution of the transcendental equation

$$\beta _0\alpha _s(Q_0^2)\mathrm{ln}(s/Q^2)\sqrt{\frac{t\mathrm{ln}(s/Q^2)}{\widehat{d}_{gg}}}$$

$$\times (\frac{8\pi Q^2}{M_W^2+Q^2}e^t+\beta _0\alpha _s(Q_0^2)\overline{d}_+(1))4\pi e^t=0,$$
(18)
where the variable $`t`$ has been rewritten as

$$t=\mathrm{ln}\left[1+\frac{\beta _0\alpha _s(Q_0^2)}{4\pi }\mathrm{ln}(Q^2/Q_0^2)\right].$$

Eq. (18) can be solved for $`\mathrm{ln}(s/Q^2)`$, obtaining an explicit expression for $`s`$ as a complicated function of the value of $`Q^2`$ at the maximum $`Q_{\text{max}}^2`$. If we rewrite Eq. (18) as $`f(s,Q^2)=0`$, we get the slope

$$\frac{dQ_{\text{max}}^2}{ds}=\frac{\frac{df(s,Q^2)}{ds}}{\frac{df(s,Q^2)}{dQ^2}}|_{Q^2=Q_{\text{max}}^2}$$

that, as can be easily shown, vanishes when $`s\mathrm{}`$.

In order to estimate the position of the maximum as a function of $`s`$, in the range $`10^7`$ GeV $`s10^{15}`$ GeV, we have chosen the parameters in intervals that include their most probable value and found that the position of the maximum can be roughly reproduced by a term linear in $`(\mathrm{ln}s)^{1.5}`$, while the value of the integrand, at the maximum, is near the value of the Bessel function there. The width at half height of the peak goes like $`\mathrm{ln}s`$. By multiplying the value, at the maximum, of the integrand in $`𝒥_1(s)`$ by the width of the peak, we get the right order of magnitude for the corresponding integral in Eq. (17), as can be verified by a numerical computation. In this limit the mass of the $`W`$-boson plays little role since $`Q_{\text{max}}^2M_W^2`$ until $`s10^{15}`$ GeV2. It is then possible to conclude that, in this approach, the total cross section for neutrino-nucleon scattering increases with $`s`$ faster than any power of $`\mathrm{ln}s`$ but slower than any power of $`s`$.

Most of the existing models of $`F_2`$ assume a fast, power-like increase in $`1/x`$ that transforms in a like increase of $`\sigma ^{\nu N}(E)`$ . In an alternative approach  the HERA data are shown to fit a model of the SF, at low and intermediate $`Q^2`$ with a moderate, logarithmic, increase in $`1/x`$, the observed “HERA effect” being attributed to the decrease of the non-leading contributions, relevant at larger $`x`$. DGLAP evolution transforms the logarithmic behavior in the way we have shown in Eqs. (9)-(10). It is important, however, to realize that this asymptotic behavior is quite different from the one obtained starting with a power, $`x^\mu `$. According to , if the starting distributions have a power-law form in $`x`$, the leading term keeps the same $`x`$-dependence under DGLAP evolution also at the next-to-leading order. The component $`xf_q^+(x,Q^2)`$ in Eq. (10) can be represented, in the leading order of perturbation theory and in the small-$`x`$ limit, by the approximated form

$$xf_q^+(x,Q^2)\frac{f}{9}\mu xf_g(x,Q_0^2)e^{d^+t},$$
(19)
if $`xf_i(x,Q_0^2)=h_i(Q_0^2)x^\mu ,(i=q,g)`$. In Eq. (19), $`d^+=\gamma ^+/(2\beta _0)`$ is the largest eigenvalue of the anomalous dimension matrix evaluated at $`(1+\mu )`$. 333Note that the power-like behavior $`x^{\widehat{d}_{gg}t}`$ of parton distributions is in agreement also with the DGLAP dynamics in the $`Q^2`$ range near the starting point $`Q_0^2`$ (see Refs. ). However, the power-like behavior cannot be combined together with (19) to some unified Regge-like $`Q^2`$-evolution of parton distributions. For $`\mu 0.4`$, as in Ref. , $`d^+`$ is negative and $`d^+3.4`$. In this approach it is easy to show that $`\overline{\sigma }^{\nu N}(s/M_W^2)^\mu +O(s^0)`$.

It is interesting to compare our result based on Eqs. (15) and (16) with a model where the cross section increases like a power of $`s`$ . We notice that Eq. (12) can be written as

$$\overline{\sigma }^{\nu N}\frac{4G_F^2}{\pi }\left(a_q^+𝒥_1+b_q^+𝒥_2\right)$$
(20)
and that all the integrals can be evaluated numerically, choosing for example $`Q_0^2=10`$ GeV2, $`\alpha _s(Q_0^2)=0.23`$ and the values $`x_M=0.001,0.01`$. The choice of $`Q_0^2=10`$ GeV2 is reasonable because in the region of $`Q^2`$ near this value the rise in $`x`$ of SF is well described by a $`\mathrm{ln}(1/x)`$ and/or $`\mathrm{ln}^2(1/x)`$ behavior 444For larger $`Q^2`$ values, the rise of $`F_2`$ is stronger and is better described by a power law $`x^\lambda `$, with $`\lambda 0.3`$ at $`Q^2100`$ GeV2. Evidently, our initial conditions (7) cannot be applied here. However, after $`Q^2`$ evolution in our approach, we get a behavior with a $`Q^2`$-dependent effective slope $`\lambda _{\text{eff}}(Q^2)`$ in agreement with experimental data ..

By evaluating the integrals numerically, a comparison can be done with the values of $`\overline{\sigma }^{\nu N}`$ obtained from Table I of Ref. . We plot in Fig. 1 the values of $`P_i(4G_F^2/\pi )𝒥_i`$ and the result of Ref. . As can be seen from Fig. 1, the power-behaved cross section increases with $`s`$ faster than the components of $`\overline{\sigma }^{\nu N}`$ in the relation (20). The explicit values of the parameters $`a_q^+`$ and $`b_q^+`$ do not enter this comparison, but it turns out that their order of magnitude can be made consistent with the parametrizations of the HERA data studied in Ref. .

The Bessel-like asymptotic behavior of the total cross section in Eq. (20) violates the Froissart bound since it is steeper than any power of $`\mathrm{ln}s`$, though being flatter than any power of $`s`$. The comparison with a power-behaved cross section, strongly violating the Froissart bound, makes us confident on the fact that additional contributions needed to restore unitarity should be less relevant in our approach.

## 4 Observation of ultra-high energy neutrinos

The evolution in leading order considered in the present paper overestimates the cross sections. Next-to-leading calculations show (see, for example, Ref. ) that the ratio $`\sigma _{\text{NLO}}/\sigma _{\text{LO}}`$ decreases when the neutrino energy increases and it is reasonable to assume that our estimates constitute an upper bound for the true predictions of the model. Hence the NLO evolution of the soft non-perturbative input could lead to a more pronounced difference in the predictions with respect to a power-like input. Another possible cause for the decrease of the true cross sections, that has not been taken into account here, resides in the gluon recombination at small $`x`$. 555The results of numerous analyses of the higher-twist contributions in DIS  and of recent papers discussing screening effects in neutrino-nucleon cross-section  are still not conclusive. This recombination can be mimicked by considering a large, “effective” gluon anomalous dimension, with a consequent increase of the value of $`\overline{d}_+(1)`$ that appears in the exponential $`\mathrm{exp}(\overline{d}_+(1)t)`$ in Eq. (10). Despite these possible corrections, it is of interest to consider the constraints imposed by our model on the interaction length of neutrinos and compare them with a more conventional approach.

For this purpose, we considered the (water equivalent) interaction length, defined as

$$L_{\text{int}}(E_\nu )=\frac{1}{\sigma ^{\nu N}(E_\nu )N_A},$$

where $`N_A=6.022\times 10^{23}`$ cm-3 (water equivalent) is the Avogadro’s number, and evaluated the ratio of $`L_{\text{int}}(E_\nu )`$ to $`L_{\text{int}}(10^8\text{ GeV})`$ as obtained from our model, with $`b_q^+/a_q^+`$ ranging between 1/20 and 1/3 and setting $`Q_0^2=10`$ GeV2, $`\alpha _s(Q_0^2)=0.23`$ and $`x_M=0.01`$. In Fig. 2 our result is compared with the same ratio determined from the results of Ref. . At higher energies $`L_{\text{int}}`$ in our model could be larger than expected in more conventional approaches.

Experiments are planned to detect UHE neutrinos by observation of the nearly horizontal air showers (HAS) in Earth atmosphere resulting from $`\nu `$-air interactions. The expected rates are proportional to the neutrino-nucleon cross section. The resulting cross section at $`10^{20}`$ eV is rather high, $`10^{31}`$ cm2 according to the estimates of , and results from the extrapolation of parton distribution functions far beyond the reach of experimental data.

If the cross section is lower, then the event rate for neutrino-induced HAS is reduced by the same factor. This reduction would compromise the main detection signal that has been proposed for UHE neutrino experiments. A smaller cross section would however offer a double advantage for the planned experiments :

1. with a new search strategy, based on the observation of upgoing air showers (UAS) initiated by muon and tau leptons produced by neutrinos interacting just below the surface of Earth, the neutrino event rate with a small cross section is actually larger than the HAS rate with a large cross section;
2. the future detectors can also measure the neutrino cross section at energies far beyond those achievable in collider experiments thus providing important information for particle physics.

If the proposed model is correct, the above scenario may have better chances to be realized.

## 5 Conclusions

UHE neutrinos became an extremely rich and interesting testing field for various physical phenomena, such as cosmology, astrophysics (origin of UHE neutrinos), theory of gravitation (production of black holes), unified theories, testing the standard theory and beyond, as well as the strong interaction theory, responsible for the hadron structure. We predict that the UHE neutrino cross sections rise moderately with energy, in accord with the above-mentioned “soft” Pomeron dynamics.

In the present paper we concentrated on the calculation of the asymptotic behavior of the structure functions and resulting cross sections. We are aware, however, about the important role of the large and intermediate $`x`$ behavior of the structure functions to be fitted to the existing, relatively low-energy, data. The cross sections calculated in this way should match with their UHE asymptotic behavior. We intend to address this point in a forthcoming publication.

## A Appendix

In this Appendix we include formulas that are not strictly necessary to understand the logical lines in the main text, but can be useful in reproducing some relevant results.

In order to determine the asymptotic behavior of $`\sigma ^{\nu N}`$ we needed some cumbersome integrals and the first one is

$`J_{}(s,Q^2)={\displaystyle _{Q^2/s}^1}{\displaystyle \frac{dx}{x}}\left[1+\left(1{\displaystyle \frac{Q^2}{xs}}\right)^2\right]xf_q^{}(x,Q^2)`$ (A.1)
$`=`$ $`{\displaystyle \frac{1}{2}}e^{d_{}(1)t}\{3a_q+(3C_1(\nu )+{\displaystyle \frac{7}{2}})b_q+{\displaystyle \frac{4Q^2}{s}}[a_q(C_1(\nu )+1)b_q]`$
$`+`$ $`{\displaystyle \frac{Q^4}{s^2}}\left[a_q+\left(C_1(\nu )+{\displaystyle \frac{1}{2}}\right)b_q\right]+4\mathrm{ln}\left({\displaystyle \frac{s}{Q^2}}\right)\left[a_q\left(C_1(\nu )+{\displaystyle \frac{3}{4}}\right)b_q\right]`$
$`+`$ $`2\mathrm{ln}^2\left({\displaystyle \frac{s}{Q^2}}\right)b_q\}.`$
Next, we prove that the integral

$$\frac{Q^2}{s}_{Q^2/s}^{x_M}\frac{dx}{x^2}\rho I_1(\sigma )$$
(A.2)
is non-leading, $`O(\rho )`$, with respect to the integral evaluated in Eq. (15). The identity

$$\frac{d}{dx}\left[\rho ^nI_n(\sigma )\right]=\frac{1}{x}\rho ^{n+1}I_{n+1}(\sigma ),n=0,1,2,\mathrm{}$$
(A.3)
can be easily obtained by using the definitions (8). Integrating repeatedly by part, we get

$$\frac{dx}{x^2}\rho I_1(\sigma )=\frac{1}{x}\rho I_1(\sigma )\frac{dx}{x^2}\rho ^2I_2(\sigma )=$$

$$=\frac{1}{x}\underset{k=1}{\overset{\mathrm{}}{}}(1)^k\rho ^kI_k(\sigma ).$$
(A.4)
The dominant contribution to the integral (A.2) comes from the lower limit of integration and, since the $`z`$ asymptotics of $`I_k(z)`$ does not depend on $`k`$, we get

$$\frac{Q^2}{s}_{Q^2/s}^{x_M}\frac{dx}{x^2}\rho I_1(\sigma )\left(1+O(1/\sqrt{\mathrm{ln}s})\right)\left[\rho I_1(\sigma )\right]_{x=Q^2/s}$$

that is non-leading with respect to the integral in Eq. (15). By using the same arguments, it is easy to show that

$$\frac{Q^4}{s^2}_{Q^2/s}^{x_M}\frac{dx}{x^4}\rho I_1(\sigma )\frac{1}{2}\left(1+O(1/\sqrt{\mathrm{ln}s})\right)\left[\rho I_1(\sigma )\right]_{x=Q^2/s}.$$

Acknowledgments L.J. thanks the Dipartimento di Fisica Teorica dell’Università di Torino, the Dipartimento di Fisica dell’Università della Calabria and the Dipartimento di Fisica dell’Università di Padova, together with the Istituto Nazionale di Fisica Nucleare (INFN), Sezioni di Padova e Torino e Gruppo collegato di Cosenza, where this work was done, for their warm hospitality and financial support. A.K. thanks for the same the Dipartimento di Fisica dell’Università di Padova and the INFN Sezione di Padova; he is supported by the INFN-LThPh agreement program. L.J. thanks A. Kusenko for an inspiring correspondence.


text:
# Many-body aspects of coherent atom-molecule oscillations

(February 18, 2003)

## Abstract

We study the many-body effects on coherent atom-molecule oscillations by means of an effective quantum field theory that describes Feshbach-resonant interactions in Bose gases in terms of an atom-molecule hamiltonian. We determine numerically the many-body corrections to the oscillation frequency for various densities of the atomic condensate. We also derive an analytic expression that approximately describes both the density and magnetic-field dependence of this frequency near the resonance. We find excellent agreement with experiment.

Introduction — One of the most remarkable applications of Feshbach resonances in doubly spin-polarized alkali gases stwalley ; eite ; inouye  is the observation of coherent atom-molecule oscillations elisabeth2 . In this last experiment Donley et al. used the Feshbach resonance at $`B_0154.9`$ Gauss in the $`|f=2;m_f=2`$ state of 85Rb to perform a Ramsey-type experiment, consisting of two short pulses in the magnetic field towards resonance separated by a longer evolution time. As a function of this evolution time an oscillation in the number of condensate atoms was observed. Over the investigated range of magnetic field during the evolution time, the frequency of this oscillation agreed exactly with the molecular binding energy found from a two-body coupled-channels calculation servaas , indicating coherence between atoms and molecules.

Very recently, Claussen et al. have performed a similar series of measurements over a larger range of magnetic fields claussen2003 . It was found that close to resonance the frequency of the oscillation deviates from the two-body molecular binding energy, which indicates that many-body effects play an important role in this regime. It is the main purpose of this Letter to present the theory that explains this deviation quantitatively.

The first mean-field theory for Feshbach-resonant interactions in Bose-Einstein condensed gases is due to Drummond et al. drummond1998  and Timmermans et al. timmermans  and introduces the physical picture of an interacting atomic condensate coupled to a noninteracting molecular condensate. Although this theory contains the correct resonant scattering amplitude for the atoms, it does not contain the correct molecular binding energy. By studying the fluctuations around this mean-field theory within the Hartree-Fock-Bogoliubov approximation, it is possible to also incorporate the correct binding energy servaas ; mackie2002 ; kohler2002 . This comes about because the latter approach explicitly contains a so-called anomalous density, or pairing field, for the atoms which, after elimination, leads to a shift in the coupling constants and the molecular binding energy. Unfortunately, however, the elimination does not lead to a proper renormalization of all the coupling constants. In particular the interaction between condensate atoms and non-condensate atoms, which is ultimately responsible for the many-body corrections to the molecular binding energy that are of interest to us here, is not correctly described within the Hartree-Fock-Bogoliubov approximation.

Relying on the anomalous density for the description of the molecular properties also makes the theory inapplicable above the critical temperature for Bose-Einstein condensation, since the anomalous density is proportional to the atomic condensate density. To overcome all these problems, it is convenient to formulate an effective quantum field theory that incorporates the exact two-body physics not at the mean-field level but at the quantum level and is hence applicable both above and below the critical temperature. We have recently derived such an effective quantum field theory by starting from the microscopic hamiltonian for an atomic gas with a Feshbach resonance and explicitly summing all the ladder diagrams rembert4 . Here, we apply the mean-field theory for the Bose-Einstein condensed phase of the gas that results from this effective quantum field theory, to the study of the recent experiments by Claussen et al. claussen2003 .

Atom-molecule coherence — The mean-field equations consist of coupled equations for the macroscopic wave function $`\psi _\mathrm{a}(𝐱,t)`$, describing the atomic Bose-Einstein condensate, and the macroscopic wave function $`\psi _\mathrm{m}(𝐱,t)`$ that describes the condensate of bare molecules, i.e., with the appropriate bound-state wave function in the closed channel potential of the Feshbach problem. The mean-field equations for the coupled atom-molecule system are ultimately given by

$`i\mathrm{}{\displaystyle \frac{\psi _\mathrm{a}(𝐱,t)}{t}}`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2m}}+T_{\mathrm{bg}}^{2\mathrm{B}}|\psi _\mathrm{a}(𝐱,t)|^2\right]\psi _\mathrm{a}(𝐱,t)`$
$`+2g\psi _\mathrm{a}^{}(𝐱,t)\psi _\mathrm{m}(𝐱,t),`$
$`i\mathrm{}{\displaystyle \frac{\psi _\mathrm{m}(𝐱,t)}{t}}`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{4m}}+\delta (B(t))\right]\psi _\mathrm{m}(𝐱,t)+g\psi _\mathrm{a}^2(𝐱,t)`$ (1)
$`g^2{\displaystyle \frac{m^{3/2}}{2\pi \mathrm{}^3}}i\sqrt{i\mathrm{}{\displaystyle \frac{}{t}}+{\displaystyle \frac{\mathrm{}^2^2}{4m}}2\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}}\psi _\mathrm{m}(𝐱,t),`$
where $`T_{\mathrm{bg}}^{2\mathrm{B}}=4\pi a_{\mathrm{bg}}\mathrm{}^2/m`$ is the off-resonant two-body T(ransition) matrix with $`a_{\mathrm{bg}}`$ the off-resonant background scattering length and $`m`$ the mass of one atom. The atom-molecule coupling $`g`$ is found from experiment by adiabatically eliminating the molecular wave function and using the fact that the resulting magnetic-field dependent scattering length of the atoms must be equal to

$$a(B)=a_{\mathrm{bg}}\left(1\frac{\mathrm{\Delta }B}{BB_0}\right),$$
(2)
with $`\mathrm{\Delta }B`$ and $`B_0`$ the experimental width and position of the resonance, respectively. This procedure results in $`g=\mathrm{}\sqrt{2\pi a_{\mathrm{bg}}\mathrm{\Delta }B\mathrm{\Delta }\mu /m}`$, where we have made use of the fact that the detuning of the bare molecular state is given by $`\delta (B)=\mathrm{\Delta }\mu (BB_0)`$ with $`\mathrm{\Delta }\mu 2.2\mu _\mathrm{B}`$ for 85Rb servaas  and $`\mu _\mathrm{B}`$ the Bohr magneton.

The mean-field equation for the molecular condensate contains a fractional derivative, corresponding to the retarded self energy of the molecules. In momentum and frequency space, this self energy reads $`\mathrm{}\mathrm{\Sigma }_\mathrm{m}^{(+)}(𝐤,\omega )=(g^2m^{3/2}/2\pi \mathrm{}^3)i\sqrt{\mathrm{}\omega \mathrm{}^2𝐤^2/4m2\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}}`$ 111The full expression for the self energy of the molecules is given by $`\mathrm{}\mathrm{\Sigma }_\mathrm{m}^{(+)}(𝐤,\omega )=(g^2m^{3/2}/2\pi \mathrm{}^3)i\sqrt{z}(1+ia_{\mathrm{bg}}\sqrt{mz}/\mathrm{})^1`$, with $`z=\mathrm{}\omega \mathrm{}^2𝐤^2/4m2\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}`$. Close to resonance only the square-root in the numerator is important. We will use this approximation throughout the paper.. The square-root behavior is a result of the Wigner threshold law for the decay of a molecule with total energy $`\mathrm{}\omega `$ and center-of-mass momentum $`\mathrm{}𝐤`$ into the two-atom continuum. Due to the mean-field interaction of the noncondensed atoms with the condensate, the decaying molecule has to overcome a mean-field barrier given by $`2\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}`$. Here, $`\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}`$ denotes the Hartree-Fock self energy for the noncondensed atoms. Neglecting the momentum dependence of this self energy, it is given by

$`\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}=2n_\mathrm{a}(T_{bg}^{2\mathrm{B}}`$
$`+{\displaystyle \frac{2g^2}{\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}+\mu \delta (B)g^2\frac{m^{3/2}}{2\pi \mathrm{}^3}\sqrt{\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}\mu }}}),`$ (3)
where $`n_\mathrm{a}=|\psi _\mathrm{a}|^2`$ is the density of the atomic condensate, $`\mu `$ its chemical potential, and we have used that a collision between a condensed atom and a noncondensed atom has a mean-field shift of $`\mu +\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}`$. Far from resonance the energy-dependence of the interactions can be safely ignored and the Hartree-Fock self energy becomes equal to $`8\pi a(B)\mathrm{}^2n_\mathrm{a}/m`$, as expected.

The molecular binding energy is given by the pole of the molecular propagator at zero momentum, which from Eq. (Many-body aspects of coherent atom-molecule oscillations) is seen to be given by

$$G_\mathrm{m}^{(+)}(\mathrm{𝟎},\omega )=\frac{\mathrm{}}{\mathrm{}\omega +i0\delta (B)+\frac{g^2m^{3/2}}{2\pi \mathrm{}^3}i\sqrt{\mathrm{}\omega 2\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}}}.$$
(4)
In the limit of vanishing condensate density $`n_\mathrm{a}`$ and for negative detuning, it has a pole at

$$ϵ_\mathrm{m}(B)=\delta (B)+\frac{g^4m^3}{8\pi ^2\mathrm{}^6}\left[\sqrt{1\frac{16\pi ^2\mathrm{}^6}{g^4m^3}\delta (B)}1\right].$$
(5)
Close to the resonance it thus follows that $`ϵ_\mathrm{m}(B)=\mathrm{}^2/[ma(B)^2]`$, which is the correct molecular binding energy in vacuum servaas . Note that in the absence of the molecular self energy the binding energy would be equal to the detuning, which is incorrect close to resonance. Moreover, the residue of the pole is given by

$`Z(B)`$ $`=`$ $`\left[1{\displaystyle \frac{\mathrm{\Sigma }_\mathrm{m}^{(+)}(\mathrm{𝟎},\omega )}{\omega }}\right]^1`$ (6)
$`=`$ $`\left[1+{\displaystyle \frac{g^2m^{3/2}}{4\pi \mathrm{}^3\sqrt{|ϵ_\mathrm{m}(B)|}}}\right]^1`$
and always smaller than one. Physically, the latter can be understood from the fact that the dressed molecular bound state near the Feshbach resonance is given by

$`|\chi _\mathrm{m};\mathrm{dressed}`$ $`=`$ $`\sqrt{Z(B)}|\chi _\mathrm{m};\mathrm{bare}`$ (7)
$`+{\displaystyle \frac{d𝐤}{(2\pi )^3}C(𝐤)|𝐤,𝐤;\mathrm{open}},`$
where the coefficient $`C(𝐤)`$ denotes the amplitude of the dressed molecular state to be in the open channel of the Feshbach problem and the two atoms having momenta $`𝐤`$ and $`𝐤`$, respectively. They are normalized as $`𝑑𝐤|C(𝐤)|^2/(2\pi )^3=1Z(B)`$. The dressed molecular state therefore only contains with an amplitude $`\sqrt{Z(B)}`$ the bare molecular state $`|\chi _\mathrm{m};\mathrm{bare}`$. Close to resonance we have that $`Z(B)1`$, whereas it approaches one far off resonance. With respect to this remark it is important to note that the result of the Hartree-Fock-Bogoliubov theory for the density of the molecular condensate should be multiplied by a factor $`1/Z(B)1`$ to obtain the density of real dressed molecules, since in this theory always the density of bare molecules is calculated rembert5 ; braaten2003 .

Collective modes — To study the many-body effects on the frequency of the coherent atom-molecule oscillations it is important to realize that these oscillations are in fact a collective mode where the atomic condensate density oscillates out-of-phase with the molecular condensate density. It is thus worthwhile to study the collective modes of the mean-field equations in Eq. (Many-body aspects of coherent atom-molecule oscillations) and look for solutions of the form

$`\psi _\mathrm{m}(𝐱,t)`$ $`=`$ $`\left[\psi _\mathrm{m}+u_𝐤^{}e^{i\omega t+i𝐤𝐱}+v_{𝐤}^{}{}_{}{}^{}e^{+i\omega ti𝐤𝐱}\right]e^{i2\mu t/\mathrm{}},`$
$`\psi _\mathrm{a}(𝐱,t)`$ $`=`$ $`\left[\psi _\mathrm{a}+u_𝐤e^{i\omega t+i𝐤𝐱}+v_𝐤^{}e^{+i\omega ti𝐤𝐱}\right]e^{i\mu t/\mathrm{}}.`$ (8)
After substitution into the mean-field equations in Eq. (Many-body aspects of coherent atom-molecule oscillations) the eigenmodes are found by diagonalizing the resulting $`4\times 4`$ matrix. This yields a dispersion relation with two branches, one corresponding to the gapless Bogoliubov modes and one that corresponds to the atom-molecule oscillations. The zero-momentum part of the latter corresponds to the experimentally observed frequency of the coherent atom-molecule oscillations. Note that for this calculation the evaluation of the molecular self energy can be performed exactly because every term in the right-hand side of Eq. (Many-body aspects of coherent atom-molecule oscillations) is an eigenfunction of the operator under the square root.

In Fig. 1 we show the dispersion relations found by means of the above procedure, for a fixed atomic condensate density of $`n_\mathrm{a}=2\times 10^{12}`$ cm-3 at a fixed magnetic field of $`B=156`$ G. Physically, the upper branch corresponds to the phonon-like excitations in the atomic condensate and the lower branch to coherent atom-molecule oscillations. The dashed line denotes the Bogoliubov dispersion for the scattering length $`a(B)`$. At low momenta the phonon branch corresponds with the Bogoliubov dispersion, as expected. At higher momenta the dispersion starts to deviate from the Bogoliubov result, due to the energy and momentum dependence of the resonant interactions that reduce the scattering amplitude. The dispersion of the lower branch obeys $`\mathrm{}\omega _𝐤\mathrm{}\omega _\mathrm{J}+\mathrm{}^2𝐤^2/4m`$, where $`\mathrm{}\omega _\mathrm{J}`$ is the Josephson frequency that is observed by Claussen et al. in their Ramsey experiment. This dispersion is negative due to the fact that we are dealing with a metastable situation. For negative detuning the true ground state contains almost all atoms in the form of molecules.

In Fig. 2 the results are shown for the frequency of the atom-molecule oscillations as a function of the magnetic field, for three different densities. Clearly, for increasing density the frequency starts to deviate from the two-body result. In addition, we show in Fig. 3 the frequency of the atom-molecule oscillations relative to the two-body binding energy, as a function of the atomic condensate density. The calculation is performed in this case for several values of the magnetic field. Both the magnetic field and the atomic density dependence of the frequency can be understood as follows. We first observe that for the range of magnetic fields and atomic densities that are explored experimentally, the difference between the energy of the dressed molecular state and the threshold of the two-particle continuum is almost independent of density and equal to $`ϵ_\mathrm{m}(B)`$. It is interesting to note that only this difference is independent of density, whereas both quantities individually show a substantial mean-field shift of about $`2\mathrm{}\mathrm{\Sigma }^{\mathrm{HF}}`$ rembert4 . As a result of the cancellation, however, also the wave-function renormalization factor $`Z`$ is almost density independent and equal to $`Z(B)`$. Expressing our mean-field equations in Eq. (Many-body aspects of coherent atom-molecule oscillations) in terms of the condensate wave function for the dressed molecules by replacing $`\psi _\mathrm{m}(𝐱,t)`$ by $`\sqrt{Z(B)}\psi _\mathrm{m}(𝐱,t)`$, we see that the coupling between the atomic condensate and the dressed molecular condensate is reduced by a factor of $`\sqrt{Z(B)}1`$ rembert5 . We thus expect the frequency to obey

$$\mathrm{}\omega _\mathrm{J}\sqrt{16Z(B)g^2n_\mathrm{a}+(ϵ_\mathrm{m}(B))^2}.$$
(9)
This analytic expression indeed turns out to give a first approximation to the deviation of the frequency from the two-body result.

In order to confront our results with the experimental data we have to realize that the experiments are performed in a magnetic trap. Taking only the ground-states into account for both the atomic and the molecular condensates, this implies effectively that the atom-molecule coupling $`g`$ is reduced by an overlap integral. Hence we define the effective homogeneous condensate density by means of $`n_\mathrm{a}=N_\mathrm{a}\left[𝑑𝐱\varphi _\mathrm{a}^2(𝐱)\varphi _\mathrm{m}(𝐱)\right]^2`$, where $`N_\mathrm{a}`$ denotes the number of condensed atoms and $`\varphi _\mathrm{a}(𝐱)`$ and $`\varphi _m(𝐱)`$ denote the atomic and molecular ground state wave function, respectively. For the experiments of Claussen et al. this results in an effective density of $`n_\mathrm{a}2\times 10^{12}`$ cm-3 claussen2003 . Fig. 2 clearly shows an excellent agreement with the experimentally observed frequency for this density.

It is important to note that there are two hidden assumptions in the above comparison. First, we have used that the dressed molecules are trapped in the same external potential as the atoms. This is not obvious because the bare molecular state involved in the Feshbach resonance is high-field seeking and therefore not trapped. However, Eq. (6) shows that near resonance almost all the amplitude of the dressed molecule is in the low-field seeking open channel and its magnetic moment is therefore almost equal to twice the atomic magnetic moment. Second, we have determined the frequency of the coherent atom-molecule oscillations in equilibrium. In contrast, the observed oscillations in the number of condensate atoms is clearly a nonequilibrium phenomenon. This is, however, expected not to play an important role because the Ramsey-pulse sequence is performed on such a fast time scale that the response of the condensate wave function can be neglected.

Conclusions — With the linear-response calculation presented in this Letter we have obtained excellent agreement with the experimental results on the frequency of coherent atom-molecule oscillations. The next step is a more detailed understanding of other quantities that are of interest for the two-pulse experiments elisabeth2 ; claussen2003 . This includes a quantitative study of the number of condensed atoms as a function of time that goes beyond the linear approximation discussed here. Experimentally an overall decay of the number of condensed atoms is observed and also that the oscillations have a finite damping rate. Both effects increase as one approaches the resonance. We expect that an important contribution to these effects is due to the so-called rogue-dissociation process mackie2002 , which has its physical origin in the decay of a molecule into two noncondensed atoms that in the center-of-mass system have opposite momenta $`𝐤`$ and $`𝐤`$, respectively. In equilibrium this process is always forbidden due to energy conservation, since the molecular state lies below the two-atom continuum threshold. However, in the nonequilibrium setting of the experiments of interest, it may occur due to the fact that the detuning is strongly time-dependent. In the mean-field equations in Eq. (Many-body aspects of coherent atom-molecule oscillations) it is the fractional derivative term corresponding to the self energy of the molecules that automatically incorporates the rogue-dissociation effect. A full numerical simulation of these equation is challenging due to the fact that we are essentially dealing with a term that is nonlocal in time. Nevertheless, work in this direction is in progress and will be reported in a future publication.