format: text

text:
# Full quantum solutions to the resonant four-wave mixing of two single-photon wave packets

## I Introduction

The cancellation of resonant linear absorption and refraction via electromagnetically induced transparency (EIT) harris1997 ; Marangos  has led to a range of new possibilities in non-linear optics. One important application is optical frequency mixing close to atomic resonances which allows the use of the enhanced nonlinear interaction without suffering from linear absorption and refraction. It has been predicted that EIT could lead to a new regime of nonlinear optics on the level of few light quanta Harris1998 ; HauHarris99 ; Imamoglu96 ; Lukin-AdvAMO-2000 ; johnsson2002a . Several schemes for resonant nonlinear processes have been proposed and analyzed, both theoretically and experimentally Lukin-AdvAMO-2000 . A particularly interesting system is resonant four-wave mixing using atoms with a double-lambda configuration Hemmer1994 ; Babin1996 ; Popov1997 ; Lu1998 . Efficient frequency conversion, generation of squeezing LukinPRL1999  as well as the possibility of mirrorless oscillations Lukin1998  with extremely low thresholds and narrow linewidth have been predicted Fleischhauer-PRL-1999  and in part experimentally observed Zibrov-PRL-2000 .

Most theoretical and experimental studies of resonant nonlinear processes have been done for classical fields. For one-dimensional set-ups where common co-moving frames exist, full analytical solutions have been derived for the interaction of classical pulses in the adiabatic limit Korsunsky2002 . In view of the potential for an efficient nonlinear interaction on the level of few photons, however, a full quantum theoretical analysis of these systems is necessary. In addition, in order to take into account finite size effects which become increasingly important in the few-photon regime, and to analyze the potential for quantum information processing, such a quantum analysis has to go beyond the usual linearization approaches. Here very little work has been done, the few exceptions being the integrable models of resonantly enhanced Kerr interaction Imamoglu-Kerr  and photon blockade Harris1998 ; blockade .

In a previous paper johnsson2002a  we have shown numerically that if single-mode fields are considered, it is possible to use an atomic vapor in a four-wave mixing double-lambda configuration to obtain full conversion from two input fields into two generated fields within a few centimeters of interaction length, even if the input fields only consist of single light quanta. The treatment was fully quantum but being restricted to a single-mode analysis, propagation effects of wave packets were not considered. In the present paper we extend this study by undertaking a multimode analysis in one spatial dimension and for co-propagating fields. In order to keep the problem tractable we restrict ourselves to the important special case of two single-photon wave packets as inputs. For this case we obtain simple analytic solutions and compare them to numerical simulations. We show that soliton-like quantum solutions exist that display an oscillatory exchange of excitation between the two input and the two generated fields. Finally we discuss briefly possible applications of single-photon four-wave mixing to quantum information processing and entanglement generation.

## II System and effective field equations

The situation we consider is resonantly enhanced four-wave mixing in the modified double-lambda system shown in Fig. 1. Note that a five-level atomic system is used instead of the original four-level system put forward in Hemmer1994 ; Lukin1998 . This is due to the fact that in the four-level system the finite detuning $`\mathrm{\Delta }`$ is associated with an AC Stark effect, which leads to intensity dependent dynamical phase shifts of the fields. These phase shifts are of minor consequence in the case where the fields are counter-propagating Fleischhauer00b , but for co-propagating fields have a detrimental influence leading to impaired phase matching and inefficient frequency conversion. As shown in johnsson2002b  these phase shifts are eliminated in the five-level scheme when the relative sign of the dipole moments of the $`|4|2,|1`$ transitions is opposite to that of the $`|3|2,|1`$ transitions.

The two fields with carrier frequencies $`\nu _{\mathrm{\Omega }1}`$ and $`\nu _{\mathrm{\Omega }2}`$ and slowly varying amplitudes $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ are initially excited and form the pump fields. The other fields of carrier frequencies $`\nu _{E1}`$ and $`\nu _{E2}`$ and slowly varying amplitudes $`E_1`$ and $`E_2`$ are generated during the interaction process and are assumed to be initially zero. $`\mathrm{\Omega }_1`$ and $`E_1`$ are taken to be exactly on resonance, i.e. $`\nu _{\mathrm{\Omega }1}=\omega _{52}`$, $`\nu _{E_1}=\omega _{51}`$, while the other two fields are detuned by an amount $`\mathrm{\Delta }`$, i.e. $`\nu _{\mathrm{\Omega }2}=\omega _{41}\mathrm{\Delta }=\omega _{31}+\mathrm{\Delta }`$ and $`\nu _{E2}=\omega _{42}\mathrm{\Delta }=\omega _{32}+\mathrm{\Delta }`$. A finite detuning $`\mathrm{\Delta }`$, large compared to the Rabi frequencies, the Doppler broadening and the decay rates from the excited states, is necessary to maximize the ratio of nonlinear gain to linear absorption. Decay from the two lower levels is considered to be negligible and all fields have the same propagation direction.

Because of energy conservation there is overall four-photon resonance, i.e. $`\nu _{\mathrm{\Omega }1}+\nu _{\mathrm{\Omega }2}=\nu _{E1}+\nu _{E2}`$. It can be shown that the contributions of the resonant transitions to the linear refractive index vanish if the fields are pairwise in two-photon resonance. Phase matching will thus favor two-photon resonance and we assume that this condition is fulfilled for the carrier frequencies of the four pulses, i.e. $`\nu _{E1}\nu _{\mathrm{\Omega }1}=\nu _{\mathrm{\Omega }2}\nu _{E2}=\omega _{21}`$.

Extending the analysis of johnsson2002a  to a multimode description, the interaction can be described by the effective adiabatic Hamiltonian hamiltonian

$`H_{\mathrm{int}}={\displaystyle \frac{\mathrm{}gc}{\mathrm{\Delta }}}{\displaystyle dz\frac{\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_2+\widehat{E}_1^{}\widehat{E}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2}{\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1+\widehat{E}_1^{}\widehat{E}_1}}.`$ (1)
Here $`\widehat{\mathrm{\Omega }}_1(z),\widehat{\mathrm{\Omega }}_1^{}(z)`$ etc. denote dimensionless, slowly-varying (both in time and space) positive and negative-frequency components of the corresponding electric fields

$`\widehat{E}_j(z)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{L}}}{\displaystyle \underset{k}{}}\widehat{a}_{jk}\mathrm{e}^{ikz}\mathrm{e}^{i\nu _{Ej}(tz/c)},`$ (2)
$`\widehat{\mathrm{\Omega }}_j(z)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{L}}}{\displaystyle \underset{k}{}}\widehat{b}_{jk}\mathrm{e}^{ikz}\mathrm{e}^{i\nu _{\mathrm{\Omega }j}(tz/c)}.`$ (3)
$`L`$ is the quantization length and $`k\nu /c`$ is the wave-vector component in $`z`$-direction relative to $`\nu /c`$. The four fields are assumed to have either sufficiently different carrier frequencies or different polarizations. Thus operators corresponding to different fields commute. The commutator between positive and negative frequency components can be approximated by a spatial delta function

$`[\widehat{E}_i(z),\widehat{E}_j^{}(z^{})]`$ $`=`$ $`{\displaystyle \frac{\delta _{ij}}{L}}{\displaystyle \underset{k}{}}\mathrm{e}^{ik(zz^{})}\delta _{ij}\delta (zz^{}).`$ (4)
We furthermore assume in (1) that all four transitions give rise to the same coupling $`g=Nd_i^2\omega _i/(2\mathrm{}cϵ_0)=3N\lambda ^2\gamma /8\pi `$, where $`N`$ is the atomic number density, $`\lambda `$ the typical wavelength of the fields and $`\gamma `$ the typical radiative decay rate. It should be noted that numerator and denominator of Eq. (1) commute and thus there is no ambiguity with respect to the ordering of the terms.

The structure of the denominator results from the saturation of the two-photon transition $`|1|2`$, whose coherence lifetime is taken to be infinite. If a finite decay rate $`\gamma _0`$ of the $`|1|2`$ coherence is taken into account, a term proportional to $`\gamma \gamma _0`$ has to be added in the denominator.

The non-polynomial character of the interaction Hamiltonian causes the nonlinear coupling to behave unusually: As shown in johnsson2002a , the interaction increases with decreasing pump field strength, making effective nonlinear frequency conversion possible even for single photons. In the derivation of the effective Hamiltonian in johnsson2002a  adiabatic conditions were assumed. This limits the applicability of (1) in the multimode case to sufficiently long pulses. A discussion of non-adiabatic corrections and their effect on the propagation of the pulses is, however, outside of the scope of the present paper and will be discussed elsewhere.

As shown in the Appendix the slowly varying amplitudes of the electric field obey

$`_t\widehat{E}_j(z,t)=c_z\widehat{E}_j(z,t)+{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}},\widehat{E}_j(z,t)]`$ (5)
and similarly for $`\widehat{\mathrm{\Omega }}_j`$. Thus from (1) we arrive at the following equations of motion

$`\left(_t+c_z\right)\widehat{E}_1`$ $`=`$ $`i\kappa c\widehat{\mathrm{\Lambda }}\left(\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_1\widehat{E}_2\widehat{\mathrm{\Omega }}_1^{}\widehat{E}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\right)\widehat{\mathrm{\Lambda }}`$
$`\left(_t+c_z\right)\widehat{\mathrm{\Omega }}_1`$ $`=`$ $`i\kappa c\widehat{\mathrm{\Lambda }}\left(\widehat{E}_1^{}\widehat{E}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{E}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_1\widehat{E}_2\right)\widehat{\mathrm{\Lambda }}`$
$`\left(_t+c_z\right)\widehat{E}_2`$ $`=`$ $`i\kappa c\widehat{\mathrm{\Lambda }}\widehat{E}_1^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2`$
$`\left(_t+c_z\right)\widehat{\mathrm{\Omega }}_2`$ $`=`$ $`i\kappa c\widehat{\mathrm{\Lambda }}\widehat{\mathrm{\Omega }}_1^{}\widehat{E}_1\widehat{E}_2`$ (6)
where $`\kappa =g/\mathrm{\Delta }`$ and $`\widehat{\mathrm{\Lambda }}=(\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1+\widehat{E}_1^{}\widehat{E}_1)^1`$.

These equations admit four independent constants of motion:

$`(_t+c_z)(\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1+\widehat{E}_1^{}\widehat{E}_1)`$ $`=`$ $`0`$
$`(_t+c_z)(\widehat{\mathrm{\Omega }}_2^{}\widehat{\mathrm{\Omega }}_2+\widehat{E}_2^{}\widehat{E}_2)`$ $`=`$ $`0`$
$`(_t+c_z)(\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2^{}\widehat{\mathrm{\Omega }}_2)`$ $`=`$ $`0`$
$`(_t+c_z)(\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_2+\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{E}_1^{}\widehat{E}_2^{})`$ $`=`$ $`0,`$ (7)
which represent the quantum analogs of the Manley-Rowe relations and additionally the conservation of the relative phase between the fields.

If we assume that the input fields consist of two, single-photon wave packets in $`\widehat{\mathrm{\Omega }}_1`$ and $`\widehat{\mathrm{\Omega }}_2`$, then it is clear, due to the constants of motion (7), that the state of the system can be represented at all times by

$$|\phi (t)=\underset{k,k^{}}{}\xi _{kk^{}}(t)|1_k\mathrm{\hspace{0.17em}1}_k^{}\mathrm{\hspace{0.17em}0\hspace{0.17em}0}+\underset{k,k^{}}{}\eta _{kk^{}}(t)|\mathrm{0\hspace{0.17em}0\hspace{0.17em}1}_k\mathrm{\hspace{0.17em}1}_k^{},$$
(8)
where $`|n_km_k^{}p_{k^{\prime \prime }}q_{k^{\prime \prime \prime }}`$ denotes $`n`$ photons in the $`k`$th mode of $`\widehat{\mathrm{\Omega }}_1`$, $`m`$ photons in the $`k^{}`$th mode of $`\widehat{\mathrm{\Omega }}_2`$, and so on.

## III Field intensities

Since for the case of a single-photon input the expectation values of all fields vanish at all times, all relevant information about the state of the system is given by the mean intensities of the fields

$`\phi (t)|\widehat{\mathrm{\Omega }}_j^{}(z)\widehat{\mathrm{\Omega }}_j(z)|\phi (t),\phi (t)|\widehat{E}_j^{}(z)\widehat{E}_j(z)|\phi (t),`$ (9)
and the two-photon wave functions

$`\psi _\mathrm{\Omega }(z,z^{},t)`$ $`=`$ $`0|\widehat{\mathrm{\Omega }}_1(z)\widehat{\mathrm{\Omega }}_2(z^{})|\phi (t)`$ (10)
$`=`$ $`{\displaystyle \underset{k,k^{}=\mathrm{}}{\overset{\mathrm{}}{}}}e^{2\pi ikz/L}e^{2\pi ik^{}z^{}/L}\xi _{kk^{}}(t)`$
$`\psi _E(z,z^{},t)`$ $`=`$ $`0|\widehat{E}_1(z)\widehat{E}_2(z^{})|\phi (t)`$ (11)
$`=`$ $`{\displaystyle \underset{k,k^{}=\mathrm{}}{\overset{\mathrm{}}{}}}e^{2\pi ikz/L}e^{2\pi ik^{}z^{}/L}\eta _{kk^{}}(t).`$
$`\psi _\mathrm{\Omega }(z,z^{},t)`$ and $`\psi _E(z,z^{},t)`$ represent the amplitude of finding the $`\widehat{\mathrm{\Omega }}_1`$ ($`\widehat{E}_1`$) photon at position $`z`$ and simultaneously the $`\widehat{\mathrm{\Omega }}_2`$ ($`\widehat{E}_2`$) photon at position $`z^{}`$. The $`\psi `$’s are 2D Fourier transforms from the $`k`$-space representations $`\xi _{kk^{}}`$ and $`\eta _{kk^{}}`$ into a real-space representation.

We first discuss the dynamics of the mean intensities of the fields. Due to the constants of motion it is sufficient to calculate say $`\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1`$.

$`(_t+`$ $`c_z)\phi |\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1|\phi =\phi |{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H},\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1]|\phi `$ (12)
$`=i\kappa c\left(\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_2\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{E}_1^{}\widehat{E}_2^{}\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\right)\widehat{\mathrm{\Lambda }}`$
$`=i\kappa c\widehat{E}_1^{}\widehat{E}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_2`$
where we have dropped the common spatial coordinate $`z`$ and used $`\widehat{H}|0=0`$, $`\widehat{\mathrm{\Lambda }}|\phi =|\phi `$ as well as the fact that three annihilation operators acting on $`|\phi `$ give zero. Similarly we find

$`\left(_t+c_z\right)`$ $`\widehat{E}_1^{}\widehat{E}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_2`$ (13)
$`=2i\kappa c\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{E}_1^{}\widehat{E}_2^{}\widehat{E}_1\widehat{E}_2`$
$`\left(_t+c_z\right)`$ $`\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{E}_1^{}\widehat{E}_2^{}\widehat{E}_1\widehat{E}_2`$ (14)
$`=2i\kappa c\widehat{E}_1^{}\widehat{E}_2^{}\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_2.`$
Consequently the differential equation we must solve is

$$\left(_t+c_z\right)^3\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1=4\kappa ^2c^2\left(_t+c_z\right)\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1.$$
(15)
We take the input to be two independent single-photon wave packets in $`\widehat{\mathrm{\Omega }}_1`$ and $`\widehat{\mathrm{\Omega }}_2`$ with the same spatial envelope $`f_0(z)`$ and vacuum in the other two fields, which corresponds to a separable initial state of the form

$`|\phi _{\mathrm{in}}`$ $`=`$ $`\left[{\displaystyle \underset{k}{}}\xi _k^0|1_k_{\mathrm{\Omega }_1}\right]\left[{\displaystyle \underset{k}{}}\xi _k^0|1_k_{\mathrm{\Omega }_2}\right]`$ (16)
$`|0_{E1}|0_{E2}.`$
The $`\xi _k^0`$ are Fourier-transforms of $`f_0(z)`$

$`f_0(z)={\displaystyle \underset{k}{}}\xi _k^0\mathrm{e}^{2\pi ikz/L}.`$ (17)
Thus

$`\widehat{\mathrm{\Omega }}_1^{}(z,t)\widehat{\mathrm{\Omega }}_1(z,t)_{\mathrm{in}}`$ $`=`$ $`\widehat{\mathrm{\Omega }}_2^{}(z,t)\widehat{\mathrm{\Omega }}_2(z,t)_{\mathrm{in}}`$ (18)
$`=`$ $`\psi _0(zct)`$
where $`\psi _0(z)=f_0^2(z)`$. With these initial conditions one finds

$$\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_1=\widehat{\mathrm{\Omega }}_2^{}\widehat{\mathrm{\Omega }}_2=\psi _0(zct)\mathrm{cos}^2(\kappa z),$$
(19)
i.e. a sinusoidal exchange of excitation between the two pump and the two generated fields. Complete conversion is achieved at $`z=\pi /\kappa `$. It is worthwhile noting that as shown in johnsson2002a , and contrary to the classical dynamics, complete conversion can only be achieved in the quantum case for initial Fock states with one or two photons in the two pump modes.

## IV Dynamics of the two-photon wave function; quantum solitons

We proceed by calculating the two-photon wave functions. The calculation of $`\psi _\mathrm{\Omega }(z,z^{},t)`$ can be split into two distinct cases, depending on whether $`z=z^{}`$ or $`zz^{}`$. Let us first assume $`z=z^{}`$.

$`(_t+`$ $`c_z)0|\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2|\phi =0|{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H},\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2]|\phi `$ (20)
$`=i\kappa c0|\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\left(\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2^{}\widehat{E}_1\widehat{E}_2+\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2\widehat{E}_1^{}\widehat{E}_2^{}\right)\widehat{\mathrm{\Lambda }}|\phi `$
$`=i\kappa c0|\widehat{E}_1\widehat{E}_2|\psi `$
where we have again used the results that $`\widehat{H}|0=0`$, $`\widehat{\mathrm{\Lambda }}|\phi =|\phi `$ and that three annihilation operators acting on $`|\phi `$ give zero. Operating on $`0|\widehat{E}_1\widehat{E}_2|\psi `$ we find

$$(_t+c_z)0|\widehat{E}_1\widehat{E}_2|\phi =i\kappa 0|\widehat{\mathrm{\Omega }}_1\widehat{\mathrm{\Omega }}_2|\phi .$$
(21)
Thus, to determine the wave function, we need to solve the differential equation

$$\left(_t+c_z\right)^2\psi _\mathrm{\Omega }(z,z,t)=\kappa ^2c^2\psi _\mathrm{\Omega }(z,z,t).$$
(22)
For an input consisting of two independent single-photon wave packets with the same spatial shape, as considered above, the initial two-photon wave function reads

$`\psi _\mathrm{\Omega }(z,z^{},t)_{\mathrm{in}}=f(zct)f(z^{}ct).`$ (23)
We also have at the entrance of the medium

$$\left(_t+c_z\right)\psi _\mathrm{\Omega }(z,z,t)|_{z=0}=i\kappa c0|\widehat{E}_1\widehat{E}_2|\phi |_{z=0}=0$$
(24)
so that the solution to (22) is given by

$$\psi _\mathrm{\Omega }(z,z,t)=\psi _0(zct)\mathrm{cos}(\kappa z),$$
(25)
with $`\psi _0(z)=f^2(z)`$. Thus the two-photon wave function at equal spatial points propagates through the medium modulated by a factor of $`\mathrm{cos}(\kappa z)`$. We see that after one full conversion cycle $`z=\pi /\kappa `$ the phase of the wave function has changed sign. This agrees with a numerical simulation of the quantum problem, the results of which are shown in Fig. 2.

We can use a similar procedure to find $`\psi _E(z,z,t)=0|\widehat{E}_1\widehat{E}_2|\phi `$, the wave function of the generated fields. We obtain

$$\psi _E(z,z,t)=i\psi _0(zct)\mathrm{sin}(\kappa z).$$
(26)
The evolution of this field is shown in Fig. 3. Note that it is $`\pi /2`$ out of phase with the drive field wave functions, as expected.

For full knowledge of the state of the system, we also need to find $`\psi (z,z^{},t)`$ where $`zz^{}`$. From (4) we see that operators at different points in space commute, so that

$`(_t+`$ $`c_z+c_z^{})0|\widehat{\mathrm{\Omega }}_1(z)\widehat{\mathrm{\Omega }}_2(z^{})|\phi `$
$`={\displaystyle \frac{i}{\mathrm{}}}0|[\widehat{H},\widehat{\mathrm{\Omega }}_1(z)]\widehat{\mathrm{\Omega }}_2(z^{})+\widehat{\mathrm{\Omega }}_1(z)[\widehat{H},\widehat{\mathrm{\Omega }}_2(z^{})]|\phi `$
Expanding the Hamiltonian, normal ordering the expression and again noting the form of (8) we eventually arrive at

$$\left(_t+c_z+c_z^{}\right)0|\widehat{\mathrm{\Omega }}_1(z)\widehat{\mathrm{\Omega }}_2(z^{})|\phi =0.$$
(28)
Thus the wave function of the system in the case where $`zz^{}`$ is simply given by the corresponding input expression

$$\psi (z,z^{},t)=f(zct)f(z^{}ct),zz^{}.$$
(29)
It is evident that the initial two-photon wave function at different coordinates propagates undisturbed throughout the medium, and that there is a discontinuous change in the behavior when moving away from the line $`z=z^{}`$. This is clearly seen in Fig. 4, which shows the two dimensional wave function $`\psi _\mathrm{\Omega }(z,z^{},\tau )`$ at time $`\tau =\pi /\kappa `$. This is the time required for one full conversion from the pump fields to the generated fields and back again. Essentially we have at this time

$`\psi _\mathrm{\Omega }(z,z^{},\tau )=\psi _\mathrm{\Omega }(z,z^{})_{\mathrm{in}}`$ $`z=z^{}`$ (30)
$`\psi _\mathrm{\Omega }(z,z^{},\tau )=\psi _\mathrm{\Omega }(z,z^{})_{\mathrm{in}}`$ $`zz^{}`$ (31)
This behavior can easily be understood. The two-photon wave function represents the joint probability density that two photons are located at positions $`z`$ and $`z^{}`$. As both photons propagate at the same speed they will never meet at the same point in space unless they did initially. Since the nonlinear interaction in (1) is local, two photons at different spatial points do not interact. Thus only if $`z=z^{}`$ is there a conversion from the pump fields into the generated fields due to the local nonlinearity.

The simple sign change in the wave function for $`z=z^{}`$ hints at the possibility of using the system as a phase gate for quantum computation, as was mentioned in the single-mode case considered in johnsson2002a . One chooses the length of the nonlinear medium such that exactly one full conversion cycle can occur. If the entire wave function changes sign, we would have a system that behaves as a true phase gate: if only one of the two inputs is populated, the pulse exits the medium unchanged, but if both are present, a phase shift of $`\pi `$ occurs.

However, due to the multimode nature of this problem, it is $`\psi (z,z^{},t)`$ for all $`z`$ and $`z^{}`$ that is a true reflection of the system rather than $`\psi (z,z,t)`$ and there is no sign change for $`zz^{}`$. To see whether we can still use the system as an approximate phase gate, we consider the behavior of the state vector coefficients in (8). Returning briefly to a picture with $`2N+1`$ discrete modes, by using Eq. (38) of the Appendix one can show that (30) and (31) imply

$$\xi _{kk^{}}(\tau )=\xi _{kk^{}}(0)\frac{2}{2N+1}\underset{mn}{}\delta _{m+n,k+k^{}}\xi _{mn}(0).$$
(32)
This gives a good picture of how the different modes have mixed among themselves, and shows that in general a simple sign change of all the coefficients cannot exist. Only if we can enforce that at least temporarily only a single effective mode of the two pump fields is excited, e.g. by using a resonator set-up is it possible to use the system as a phase gate.

We see however from (32) that after one full conversion cycle the initially factorized state (16) evolves into an entangled state between all modes of the fields $`\widehat{\mathrm{\Omega }}_1`$ and $`\widehat{\mathrm{\Omega }}_2`$. Thus the nonlinear interaction generates entanglement.

On the other hand if the initial state is a two-photon wave packet in an entangled state, such that the initial two-photon wavefunction has only a contribution for $`z=z^{}`$

$`\psi _\mathrm{\Omega }(z,z^{})|_{\mathrm{in}}=\varphi _0(z)\delta (zz^{})`$ (33)
only diagonal components of the two-photon wave function will ever be nonzero. According to (25) and (26) they undergo sinusoidal oscillations

$`\psi _\mathrm{\Omega }(z,t)`$ $`=`$ $`\varphi _0(zct)\mathrm{cos}(\kappa z),`$ (34)
$`\psi _E(z,t)`$ $`=`$ $`i\varphi _0(zct)\mathrm{sin}(\kappa z).`$ (35)
The superposition of pump and generated fields

$`\mathrm{\Phi }(z,t)`$ $``$ $`\mathrm{cos}(\kappa z)\psi _\mathrm{\Omega }(z,t)+i\mathrm{sin}(\kappa z)\psi _E(z,t)`$ (36)
$`=`$ $`\varphi _0(zct)`$
propagates in a form-stable manner and represents a quantum soliton solution. The two-photon wavefunction $`\mathrm{\Phi }(z,t)`$ of the quantum soliton corresponds to a quasi-particle excitation

$`\widehat{\mathrm{\Psi }}^{}(z,t)`$ (37)
$`\mathrm{cos}(\kappa z)\widehat{\mathrm{\Omega }}_1^{}(z,t)\widehat{\mathrm{\Omega }}_2^{}(z,t)i\mathrm{sin}(\kappa z)\widehat{E}_1^{}(z,t)\widehat{E}_2^{}(z,t).`$
## V summary

In the present paper we have presented a full quantum theoretical treatment of resonant forward four-wave mixing for pulses. For this we have used an effective Hamiltonian derived from the interaction of the fields with an ensemble of atoms in a double-lambda configuration using an adiabatic approximation johnsson2002a . We were particularly interested in the few-photon regime since here quantum effects dominate the dynamics and because of potential applications in quantum information processing. Thus we have restricted ourselves to the important special case of an input consisting of two single-photon wave packets. For this case we were able to analytically solve the propagation equations for the field intensities and two-photon wave functions which contain all relevant information about the quantum state.

We found that there is an oscillatory energy exchange between the two pump and generated fields with 100% conversion at periodic intervals of interaction. This result is characteristic for a few-photon Fock-state input; for a coherent input complete conversion can only be achieved asymptotically for very large input power.

We have also shown that after even multiples of the conversion length the two-photon wave function $`\psi (z,z^{},t)`$ regains its initial form, while after odd multiples there is a sign flip for $`z=z^{}`$.

If the two input wave packets are not independent but in a highly entangled state, the two-photon wave function can be made zero outside of the diagonal. It was shown that such a pair of input wave packets form a formstable soliton-like quantum solution which is a superposition of pump and generated fields with oscillating coefficients.

The process of resonant four-wave mixing was shown to generate large entanglement between the modes forming the two, single-photon wave packets. Furthermore the nonlinear interaction strength is large enough to generate a controlled phase shift of a single photon by the presence of another one. Thus, if the number of relevant modes is at least temporarily restricted by some external means like a resonator, the system could have interesting applications as a photonic phase gate.

## Acknowledgement

The authors acknowledge financial support by the Deutsche Forschungsgemeinschaft through the special program on quantum information and the Graduiertenkolleg 792 at the University of Kaiserslautern.

## Appendix

To simplify the transition from a single-mode description to a multimode one, we will first consider one single-mode quantum field $`\widehat{a}`$, and then generalize the result to our four-field system. Suppose that the single-mode Hamiltonian governing the evolution of $`\widehat{a}`$ is given by $`\widehat{H}=\mathrm{}\omega _0\widehat{a}^{}\widehat{a}+\widehat{H}_{\mathrm{int}}(\widehat{a},\widehat{a}^{})`$. To go over to the multimode description we consider an interaction region of length $`L`$, divided into $`2N+1`$ cells, and consider a discrete set of modes around the carrier frequency of the field, i.e. $`k_n=k_0+2n\pi /L`$, $`NnN`$. We now define localized field operators (denoted by a tilde) via

$`\widehat{\stackrel{~}{a}}_l`$ $`=`$ $`{\displaystyle \underset{k=N}{\overset{N}{}}}\widehat{a}_k\mathrm{exp}\left[{\displaystyle \frac{2\pi ikl}{2N+1}}\right]`$ (38)
$`\widehat{a}_k`$ $`=`$ $`{\displaystyle \frac{1}{2N+1}}{\displaystyle \underset{l=N}{\overset{N}{}}}\widehat{\stackrel{~}{a}}_l\mathrm{exp}\left[{\displaystyle \frac{2\pi ikl}{2N+1}}\right]`$ (39)
where the $`\widehat{a}_k`$ are annihilation operators for mode $`k`$. $`\widehat{\stackrel{~}{a}}_l`$ is related to the field strength in cell $`l`$.

These localized field operators have the commutation relations

$$[\widehat{\stackrel{~}{a}}_l,\widehat{\stackrel{~}{a}}_l^{}^{}]=(2N+1)\delta _{ll^{}}.$$
(40)
The multimode Hamiltonian can now be written as

$`\widehat{H}`$ $`=`$ $`\mathrm{}{\displaystyle \underset{k}{}}\omega _k\widehat{a}_k^{}\widehat{a}_k+{\displaystyle \underset{l}{}}\widehat{H}_{\mathrm{int}}(\widehat{\stackrel{~}{a}}_l,\widehat{\stackrel{~}{a}}_l^{})`$
$`=`$ $`{\displaystyle \frac{\mathrm{}\omega _0}{2N+1}}{\displaystyle \underset{l}{}}\widehat{\stackrel{~}{a}}_l^{}\widehat{\stackrel{~}{a}}_l\mathrm{}{\displaystyle \underset{ll^{}}{}}\omega _{ll^{}}\widehat{\stackrel{~}{a}}_l^{}\widehat{\stackrel{~}{a}}_l^{}+{\displaystyle \underset{l}{}}\widehat{H}_{\mathrm{int}}(\widehat{\stackrel{~}{a}}_l,\widehat{\stackrel{~}{a}}_l^{})`$
where $`\omega _0`$ is the carrier frequency and

$$\omega _{ll^{}}=\underset{k=N}{\overset{N}{}}\frac{2\pi kc}{(2N+1)^2L}\mathrm{exp}\left[\frac{2\pi ik(ll^{})}{2N+1}\right].$$
(42)
Commutation relations yield the following Heisenberg equation of motion

$$\frac{}{t}\widehat{\stackrel{~}{a}}_l=i\omega _0\widehat{\stackrel{~}{a}}_li(2N+1)\underset{l^{}}{}\omega _{ll^{}}\widehat{\stackrel{~}{a}}_l^{}+\frac{i}{\mathrm{}}[\widehat{\stackrel{~}{H}}_{\mathrm{int}},\widehat{\stackrel{~}{a}}]$$
(43)
Now, as we let $`N\mathrm{}`$ we find

$`lL/(2N+1)`$ $``$ $`z`$ (44)
$`\widehat{\stackrel{~}{a}}_l`$ $``$ $`\widehat{a}(z,t)`$ (45)
$`i(2N+1){\displaystyle \underset{l^{}}{}}\omega _{ll^{}}\widehat{\stackrel{~}{a}}_l^{}`$ $``$ $`c_z\widehat{a}(z,t)`$ (46)
$`[\widehat{a}(z,t),\widehat{a}^{}(z^{},t)]`$ $``$ $`L\delta (zz^{}).`$ (47)
Introducing slowly time-varying amplitudes we obtain

$$\left(_t+c_z\right)\widehat{a}(z,t)=\frac{i}{\mathrm{}}[\widehat{H}_{\mathrm{int}}(t),\widehat{a}(z,t)]$$
(48)
Thus the multimode equations of motion look exactly like the single-mode equations of motion with the exception that a $`c_z`$ term has been added and the fields now have a spatial dependence.

Returning to the four-wave mixing situation described in this paper, we see that the interaction Hamiltonian given by (1) is of the form shown above, with the summation replaced by an integral in the continuum limit. The multimode annihilation and creation operators $`\widehat{E}(z)`$ and $`\widehat{\mathrm{\Omega }}(z)`$ defined in (2) and (3) are analogous to the localized field operators $`\widehat{a}(z,t)`$ defined above, except for the factor of $`\sqrt{L}`$ inserted to ensure correct commutation relations regardless of the quantization length. Thus the equation of motion (5) follows from (48) above.


text:
# Generalized Parton Distributions from Lattice QCD

DESY 03-051

LU-ITP 2003/006

Edinburgh-2003/04

LTH-573

MIT-CTP-3365

## Abstract

We perform a quenched lattice calculation of the first moment of twist-two generalized parton distribution functions of the proton, and assess the total quark (spin and orbital angular momentum) contribution to the spin of the proton.

Generalized parton distributions GPD  (GPDs) provide a deeper understanding of the internal structure of hadrons in terms of quarks and gluons. While ordinary parton distributions measure the probability $`|\psi (x)|^2`$ of finding a parton with fractional momentum $`x`$ in the hadron, GPDs describe the coherence of two different hadron wave functions $`\psi ^{}(x+\xi /2)\psi (x\xi /2)`$, one where the parton carries fractional momentum $`x+\xi /2`$ and one where this fraction is $`x\xi /2`$, from which information about parton-parton correlation functions can be deduced. As a consequence, GPDs depend on the momentum transfer $`\mathrm{\Delta }^2`$ between the initial and final hadron, which provides further information on the transverse location of quarks and gluons Diehl . Spatial images of hadrons can thus be obtained, where the resolution is determined by the virtuality $`Q^2`$ of the incoming photon. Last, but not least, GPDs allow us to isolate the contribution of the quark orbital angular momentum to the spin of hadrons. Lattice QCD is the only known method that is able to compute moments of GPDs from first principles.

We will restrict ourselves to the GPDs $`H_q`$ and $`E_q`$ of the nucleon, where $`q=u,d,\mathrm{}`$ denotes the flavor of the struck quark. We will not consider the gluon sector here. The lowest, zeroth moments of $`H_q`$ and $`E_q`$ are given by the Dirac and Pauli form factors:

$`{\displaystyle _1^1}\text{d}xH_q(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`F_1^q(\mathrm{\Delta }^2),`$ (1)
$`{\displaystyle _1^1}\text{d}xE_q(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`F_2^q(\mathrm{\Delta }^2).`$ (2)
Both form factors have been computed on the lattice in a similar calculation QCDSF  to the present one and found to be well described by a dipole ansatz

$$F_{1,2}^q(\mathrm{\Delta }^2)=F_{1,2}^q(0)/(1\mathrm{\Delta }^2/M_{1,2}^2)^2$$
(3)
for sufficiently small (and accessible) momenta, with dipole masses $`M_{1,2}`$ of the order of the $`\rho `$, $`\omega `$ mass, when extrapolated to the physical pion mass.

The first moments of $`H_q`$ and $`E_q`$ are of the form GPD

$`{\displaystyle _1^1}\text{d}xxH_q(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`A_2^q(\mathrm{\Delta }^2)+\xi ^2C_2^q(\mathrm{\Delta }^2),`$ (4)
$`{\displaystyle _1^1}\text{d}xxE_q(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`B_2^q(\mathrm{\Delta }^2)\xi ^2C_2^q(\mathrm{\Delta }^2),`$ (5)
where $`A_2^q(\mathrm{\Delta }^2)`$, $`B_2^q(\mathrm{\Delta }^2)`$ and $`C_2^q(\mathrm{\Delta }^2)`$ are generalized form factors, which are given by the nucleon matrix elements of the energy-momentum tensor (EMT):

$$\begin{array}{cc}\hfill p^{}|𝒪_{\{\mu \nu \}}^q|p& \frac{\text{i}}{2}p^{}|\overline{q}\gamma _{\{\mu }\stackrel{}{D}_{\nu \}}q|p\hfill \\ & =A_2^q(\mathrm{\Delta }^2)\overline{u}(p^{})\gamma _{\{\mu }\overline{p}_{\nu \}}u(p)\hfill \\ & B_2^q(\mathrm{\Delta }^2)\frac{\text{i}}{2m_N}\overline{u}(p^{})\mathrm{\Delta }^\alpha \sigma _{\alpha \{\mu }\overline{p}_{\nu \}}u(p)\hfill \\ & +C_2^q(\mathrm{\Delta }^2)\frac{1}{m_N}\overline{u}(p^{})u(p)\mathrm{\Delta }_{\{\mu }\mathrm{\Delta }_{\nu \}}.\hfill \end{array}$$
(6)
Here $`m_N`$ denotes the nucleon mass, $`\overline{p}=\frac{1}{2}\left(p^{}+p\right)`$, $`\mathrm{\Delta }=p^{}p`$, and curly brackets refer to symmetrization of indices and subtraction of traces. The EMT has twist two and spin two. It is assumed to be renormalized at the scale $`\mu `$, which makes $`A_2^q(\mathrm{\Delta }^2)`$, $`B_2^q(\mathrm{\Delta }^2)`$ and $`C_2^q(\mathrm{\Delta }^2)`$ scale and scheme dependent. For the classification of states of definite $`J^{PC}`$ contributing to (6) in the $`t`$-channel see Ji&Lebed . The so-called skewedness parameter $`\xi `$ is defined by $`\xi =n\mathrm{\Delta }`$, where $`n`$ is a light-like vector with $`n\overline{p}=1`$, and bounded by $`|\xi |2\sqrt{\mathrm{\Delta }^2/(\mathrm{\Delta }^24m_N^2)}`$. In the forward limit, $`\mathrm{\Delta }^20`$, we have

$$A_2^q(0)=x_q_0^1\text{d}xx\left(q_{}(x)+q_{}(x)\right),$$
(7)
where $`q_{()}(x)`$ are the usual quark distributions with spin parallel (antiparallel) to the spin of the nucleon. Furthermore, one derives Ji

$$\frac{1}{2}\left(A_2^q(0)+B_2^q(0)\right)=J_q,$$
(8)
where $`J_q`$ is the angular momentum of the $`q`$ quark, and $`J=_qJ_q`$ is the total angular momentum of the nucleon carried by the quarks. The angular momentum decomposes, in a gauge invariant way, into two pieces:

$$J_q=L_q+S_q,$$
(9)
where $`L_q`$ is the orbital angular momentum and

$$S_q=\frac{1}{2}\mathrm{\Delta }q\frac{1}{2}_0^1\text{d}x\left(q_{}(x)q_{}(x)\right)$$
(10)
is the spin of the quark. We know $`\mathrm{\Delta }q`$ from separate calculations QCDSF3 ; Deltaq , so that $`L_q`$ can be computed from (8).

In this Letter we perform a quenched lattice calculation of the generalized form factors $`A_2^q(\mathrm{\Delta }^2)`$, $`B_2^q(\mathrm{\Delta }^2)`$ and $`C_2^q(\mathrm{\Delta }^2)`$. The quenched approximation neglects fluctuations of virtual quark-antiquark pairs from the Dirac sea. The non-forward matrix elements (6) are computed from ratios of three- and two-point functions following QCDSF . Further details are given in MIT . To keep cut-off effects small, we use non-perturbatively $`𝒪(a)`$ improved Wilson fermions. The calculation is done on $`16^3\mathrm{\hspace{0.17em}32}`$ lattices at $`\beta =6.0`$ and for three different hopping parameters, $`\kappa =0.1324`$, $`0.1333`$ and $`0.1342`$, which allows us to extrapolate our results to the chiral limit. Using $`r_0=0.5\text{fm}`$ to set the scale, which results in the inverse lattice spacing $`1/a=2.12\text{GeV}`$, the corresponding pion masses are $`1070`$, $`870`$ and $`640\text{MeV}`$. If we use the nucleon mass extrapolated to the chiral limit to set the scale, the pion masses are $`930`$, $`760`$ and $`550\text{MeV}`$, and $`1/a=1.84\text{GeV}`$. The corresponding nucleon masses and the choice of nucleon momenta $`p,p^{}`$ can be inferred from QCDSF . For the EMT we consider two sets of (euclidean) operators:

$$\frac{1}{\sqrt{2}}\left(𝒪_{\mu \nu }+𝒪_{\nu \mu }\right),\mathrm{\hspace{0.33em}1}\mu <\nu 4$$
(11)
and

| $`{\displaystyle \frac{1}{2}}\left(𝒪_{11}+𝒪_{22}𝒪_{33}𝒪_{44}\right)`$ , |
| --- |
| $`{\displaystyle \frac{1}{\sqrt{2}}}\left(𝒪_{33}𝒪_{44}\right),{\displaystyle \frac{1}{\sqrt{2}}}\left(𝒪_{11}𝒪_{22}\right)`$ . |

(12)
Each set transforms irreducibly under the hypercubic group. The operators (11) and (12) are renormalized multiplicatively, $`𝒪(\mu )=Z(a\mu )𝒪(a)`$, with renormalization constants QCDSF3  $`Z_{v_{2a}}`$ and $`Z_{v_{2b}}`$, respectively. The renormalization constants are computed non-perturbatively QCDSF4  following Z . We obtain $`Z_{v_{2a}}^{\overline{MS}}(2\text{GeV})=1.10`$ and $`Z_{v_{2b}}^{\overline{MS}}(2\text{GeV})=1.09`$. The following results refer to the $`\overline{MS}`$ scheme at the renormalization scale $`\mu =2\text{GeV}`$.

In Figs. 1 and 2 we show the generalized form factors $`A_2^u(\mathrm{\Delta }^2)`$, $`B_2^u(\mathrm{\Delta }^2)`$, $`C_2^u(\mathrm{\Delta }^2)`$ and $`A_2^d(\mathrm{\Delta }^2)`$, $`B_2^d(\mathrm{\Delta }^2)`$, $`C_2^d(\mathrm{\Delta }^2)`$ of the proton for $`\kappa =0.1333`$. Data points with larger errors are not shown here but are included in the fit. The corresponding form factors of the neutron are obtained by interchanging $`u`$ and $`d`$. Similarly good results are found for $`\kappa =0.1324`$ and $`0.1342`$. The generalized form factors can be well described by the dipole ansatz

$$A_2^q(\mathrm{\Delta }^2)=A_2^q(0)/(1\mathrm{\Delta }^2/M^2)^2,$$
(13)
and similarly for $`B_2^q`$ and $`C_2^q`$. Fits of $`A_2^u(\mathrm{\Delta }^2)`$ and $`A_2^d(\mathrm{\Delta }^2)`$ give the same dipole mass $`M`$ within errors. The dipole masses obtained from separate fits of $`B_2^u(\mathrm{\Delta }^2)`$, $`B_2^d(\mathrm{\Delta }^2)`$, $`C_2^u(\mathrm{\Delta }^2)`$ and $`C_2^d(\mathrm{\Delta }^2)`$ are found to be consistent with that value. We therefore have decided to fit our data by a common dipole mass $`M`$. Our data do not favor a monopole behavior. The results of the fits are shown in Table I. For a reliable extrapolation to $`\mathrm{\Delta }^2=0`$ we find it important to cover a wide enough range of $`\mathrm{\Delta }^2`$ values. This may be the reason why our dipole masses turn out to be systematically larger than those found in a previous calculation Liu .

In Fig. 3 we show the dipole mass $`M`$ as a function of the pion mass. The mass values appear to lie on a straight line, as was observed already in the case of the nucleon form factors QCDSF . A linear extrapolation in $`m_\pi `$ to the physical pion mass gives $`M=1.1(2)\text{GeV}`$. This value is close to the physical masses of the $`f_2`$, $`a_2`$ mesons, which supports the hypothesis of tensor meson dominance. A quadratic extrapolation in $`m_\pi `$ leads to $`M=1.3(1)\text{GeV}`$. The form factor data $`A_2^q(0)`$, $`B_2^q(0)`$ and $`C_2^q(0)`$ show little variation with the quark mass and are extrapolated quadratically in $`m_\pi `$ to the physical pion mass. The results are shown in the bottom row of Table I. It should be stressed that all quantities refer (at best) to valence quark distributions, because sea quark effects have been neglected. In unquenched simulations there are also quark-line disconnected contributions. For an estimate see Liu .

If the dipole behavior (3), (13) continues to hold for the higher moments as well, and if we assume that the dipole masses continue to grow in a Regge-like fashion, we would obtain

$$_1^1\text{d}xx^nH_q(x,0,\mathrm{\Delta }^2)=x_q^n/(1\mathrm{\Delta }^2/M_{n+1}^2)^2,$$
(14)
with $`M_l^2=\text{const.}+l/\alpha ^{}`$, $`\alpha ^{}`$ being the slope of the Regge trajectory. This would mean that with increasing momentum transfer $`|\mathrm{\Delta }^2|`$ the lower moments of $`H_q(x,0,\mathrm{\Delta }^2)`$ are suppressed more than the higher ones, so that the observed peak in $`H_q(x,0,0)=q_{}(x)+q_{}(x)`$ around $`x0.2`$ is shifted towards the higher values of $`x`$. As a result, the $`\mathrm{\Delta }^2`$ dependence cannot be factorized in a simple way, as is sometimes assumed. Knowing $`x_q^n`$, we can reconstruct $`H_q(x,0,\mathrm{\Delta }^2)`$ from (14) by inverse Mellin transform. The $`\xi `$ dependence of both $`H_q`$ and $`E_q`$ appears to be rather weak, based on our knowledge of the first two moments, and in the isovector channel (corresponding to proton–neutron or $`u`$$`d`$ matrix elements) it largely cancels out.

In Fig. 4 we show the total angular momentum $`J=J_u+J_d`$. The dependence on the pion mass is rather flat, as expected Chen . The errors are due to the relatively large statistical errors of $`B_2^u`$ and $`B_2^d`$ and the fact that $`B_2^u`$ and $`B_2^d`$ cancel to a large extent. In Table II we give our results for $`J`$, and separately for $`J_q`$ and $`S_q`$, extrapolated quadratically (linearly in $`m_\pi ^2`$) to the physical pion mass. The numbers for $`S_q`$ refer to our latest results QCDSF4 , computed from the non-perturbatively improved axial vector current with non-perturbative renormalization factors. It turns out that the total angular momentum $`J`$ carried by the quarks amounts to $`70\%`$ of the spin of the (quenched) proton, leaving a contribution of $`30\%`$ for the gluons. The major contribution is given by the $`u`$ quark, while the contribution of the $`d`$ quark is found to be negligible, which hints at strong pairing effects. Our result for $`J`$ is somewhat smaller than that of Liu ; Gadiyak . We are able to compute $`L_q`$ now. The total orbital angular momentum of the quarks turns out to be consistent with zero:

$$LL_u+L_d=0.03(7).$$
(15)
This indicates that (at virtuality $`\mu =2\text{GeV}`$) the parton’s transverse momentum in the (quenched) proton is small. A similar conclusion can be drawn from our earlier finding d2  of a small twist-three contribution $`d_2`$ to the second moment of the polarized structure function $`g_2`$.

The generalized form factors $`C_2^q(\mathrm{\Delta }^2)`$ contribute to the beam charge asymmetry of deeply virtual Compton scattering. We obtain a rather small value: $`C_2^u(0)+C_2^d(0)=0.2(1)`$. This result is to be compared with the value $`0.8`$ obtained in the chiral quark-soliton model at $`\mu 0.6\text{GeV}`$ Kivel . For a discussion see also pheno .

As far as one can compare, quenched and unquenched results agree surprisingly well, and we do not expect to find significant differences here either. For a recent study of quenching artifacts, as well as cut-off effects, see QCDSF5 .

###### Acknowledgements.

This work is supported by DFG and under the Feodor Lynen program. Discussions with P. Hägler, J. Negele, D. Renner and C. Weiss are acknowledged. The numerical calculations have been performed at NIC (Jülich).


text:
# Measurements of Proton, Helium and Muon Spectra at Small Atmospheric Depths with the BESS Spectrometer

## 1 Introduction

Primary cosmic rays interact with nuclei in the atmosphere, and produce atmospheric secondary particles, such as muons, gamma rays and neutrinos. It is important to understand these interactions to investigate cosmic-ray phenomena inside the atmosphere. For precise study of the atmospheric neutrino oscillation , it is crucial to reduce uncertainties in hadronic interactions, which are main sources of systematic errors in the prediction of the energy spectra of atmospheric neutrinos. At small atmospheric depths below a few ten g/cm2, production process of muons is predominant over decay process, thus we can clearly observe a feature of the hadronic interactions. In spite of their importance, only a few measurements have been performed with modest statistics because of strong constraints of short observation time of a few hours during balloon ascending periods from the ground to the balloon floating altitude.

In 2001, using the BESS spectrometer, precise measurements of the cosmic-ray fluxes and their dependence on atmospheric depth were carried out during slow descending from 4.5 g/cm2 to 28 g/cm2 for 12.4 hours. The growth curves of the cosmic-ray fluxes were precisely measured. The results were compared with the predictions based on the hadronic interaction models currently used in the atmospheric neutrino flux calculations.

## 2 BESS spectrometer

The BESS ($`\underset{¯}{\mathrm{B}}`$alloon-borne $`\underset{¯}{\mathrm{E}}`$xperiment with a $`\underset{¯}{\mathrm{S}}`$uperconducting $`\underset{¯}{\mathrm{S}}`$pectrometer) detector  is a high-resolution spectrometer with a large acceptance to perform precise measurement of absolute fluxes of various cosmic rays , as well as highly sensitive searches for rare cosmic-ray components. Fig. 1 shows a schematic cross-sectional view of the BESS instrument. In the central region, a uniform magnetic field of 1 Tesla is provided by a thin superconducting solenoidal coil. The magnetic-rigidity ($`RPc/Ze`$) of an incoming charged particle is measured by a tracking system, which consists of a jet-type drift chamber (JET) and two inner-drift-chambers (IDC’s) inside the magnetic field. The deflection ($`R^1`$) is calculated for each event by applying a circular fitting using up-to 28 hit-points, each with a spatial resolution of 200 $`\mu `$m. The maximum detectable rigidity (MDR) was estimated to be 200 GV. Time-of-flight (TOF) hodoscopes provide the velocity ($`\beta `$) and energy loss ($`E/x`$) measurements. A 1/$`\beta `$ resolution of 1.4 % was achieved in this experiment. For particle identification, the BESS spectrometer was equipped with a threshold-type aerogel Cherenkov counter and an electromagnetic shower counter. The refractive index of silica aerogel radiator was 1.022, and the threshold kinetic energy for proton was 3.6 GeV. The shower counter consists of a plate of lead with two radiation lengths covering a quarter area of the lower TOF counters, whose output signal was utilized for $`e`$/$`\mu `$ identification.

The data acquisition sequence is initiated by a first-level TOF trigger, which is a simple coincidence of signals in the upper and lower TOF counter. In order to build a sample of unbiased triggers, one of every four events were recorded. The TOF trigger efficiency was evaluated to be 99.4 % $`\pm `$ 0.2 % by using secondary proton beam at the KEK 12 GeV proton synchrotron. In addition to the TOF trigger, an auxiliary trigger is generated by a signal from the Cherenkov counter to record particles above threshold energy without bias or sampling. An efficiency of Cherenkov trigger were evaluated as a ratio of the Cherenkov-triggered events among the unbiased trigger sample. It was 95.1 % for relativistic particles ($`\beta 1`$). For flux determinations, the Cherenkov-triggered events were used only above 9.5 GV for protons, 11.1 GV for helium nuclei, and 0.90 GV for muons. Below these rigidities, the TOF-triggered events were used.

## 3 Balloon flight

The BESS-2001 balloon flight was carried out at Ft. Sumner, New Mexico, USA (3449N, 10422W) on 24th September 2001. Throughout the flight, the vertical geomagnetic cut-off rigidity was about 4.2 GV. The balloon reached at a normal floating altitude of 36 km at an atmospheric depth of 4.5 g/cm2. After a few hours, the balloon started to lose the floating altitude and continued descending for more than 13 hours until termination of the flight. During the descending period, data were collected at atmospheric depths between 4.5 g/cm2 and 28 g/cm2. The atmospheric depth was measured with accuracy of $`\pm `$ 1 g/cm2, which comes mainly from an error in absolute calibration of an environmental monitor system. Fig. 2 shows a balloon flight profile during the experiment.

## 4 Data analysis

### 4.1 Data reduction

We selected “non-interacted” events passing through the detector without any interactions. The non-interacted event was defined as an event, which has only one isolated track, one or two hit-counters in each layer of the TOF hodoscopes, and proper d$`E`$/d$`x`$ inside the upper TOF counters. There is a slight probability that particles interact with nuclei in the detector material and only one secondary particle goes into the tracking volume to be identified as a non-interacted event. These events were rejected by requiring proper d$`E`$/d$`x`$ inside the upper TOF counters, because the interaction in the detector is expected to give a large energy deposit in the TOF counter. According to the Monte Carlo simulation with GEANT , these events are to be considered only below a few GeV. In order to estimate an efficiency of the non-interacted event selection, Monte Carlo simulations were performed. The probability that each particle can be identified as a non-interacted event was evaluated by applying the same selection criterion to the Monte Carlo events as was applied to the observed data. The systematic error was estimated by comparing the hit number distribution of the TOF counters. For muons, the simulated data were compared with muon data sample measured on the ground. The systematic error for protons below 1 GeV was directly determined by the detector beam test using accelerator proton beam . The resultant efficiency and its error of non-interacted event selection for protons was 83.3 $`\pm `$ 2.0 % at 1.0 GeV and 77.4 $`\pm `$ 2.5 % at 10 GeV, and that for helium nuclei was 71.6 $`\pm `$ 2.6 % at 1.0 GeV/n and 66.0 $`\pm `$ 2.9 % at 10 GeV/n. The efficiency for muons was 94.0 $`\pm `$ 0.9 % and 92.8 $`\pm `$ 0.8 % at 0.5 GeV and 10 GeV, respectively.

The selected non-interacted events were required to pass through the fiducial volume defined in this analysis. The fiducial volume of the detector was limited to the central region of the JET chamber for a better rigidity measurement. The zenith angle ($`\theta _z`$) was limited within $`\mathrm{cos}\theta _z0.90`$ to obtain nearly vertical fluxes. For the muon analysis, we used only the particles which passed through the lead plate to estimate electron contaminations. For the proton and helium analysis, particles were required not to pass through the lead plate so as to keep the interaction probability inside the detector as low as possible.

In order to check the track reconstruction efficiency inside the tracking system, the recorded events were scanned randomly. It was confirmed that 996 out of 1,000 visually identified tracks which passed through the fiducial volume were successfully reconstructed, thus the track reconstruction efficiency was evaluated to be 99.6 $`\pm `$ 0.2 %. It was also confirmed that rare interacted events are fully eliminated by the non-interacted event selection.

### 4.2 Particle identification

In order to select singly and doubly charged particles, particles were required to have proper d$`E`$/d$`x`$ as a function of rigidity inside both the upper and lower TOF hodoscopes. The upper TOF d$`E`$/d$`x`$ was already examined in the non-interacted event selection. The distribution of d$`E`$/d$`x`$ inside the lower TOF counter and the selection boundaries are shown in Fig. 3.

In order to estimate efficiencies of the d$`E`$/d$`x`$ selections for protons and helium nuclei, we used another data sample selected by independent information of energy loss inside the JET chamber. The estimated efficiency in the d$`E`$/d$`x`$ selection at 1 GeV was 98.3 $`\pm `$ 0.4 % and 97.2 $`\pm `$ 0.5 % for protons and helium nuclei, respectively. The accuracy of the efficiency for protons and helium nuclei was limited by statistics of the sample events. Since muons could not be distinguished from electrons by the JET chamber, the d$`E`$/d$`x`$ selection efficiency for muons was estimated by the Monte Carlo simulation to be 99.3 $`\pm `$ 1.0 %. The error for muons comes from the discrepancy between the observed and simulated d$`E`$/d$`x`$ distribution inside the TOF counters. The d$`E`$/d$`x`$ selection efficiencies were almost constant in the whole energy region discussed here.

Particle mass was reconstructed by using the relation of 1/$`\beta `$, rigidity and charge, and was required to be consistent with protons, helium nuclei or muons. An appropriate relation between 1/$`\beta `$ and rigidity for each particle was required as shown in Fig. 4. Since the 1/$`\beta `$ distribution is well described by Gaussian and a half-width of the 1/$`\beta `$ selection band was set at 3.89 $`\sigma `$, the efficiency is very close to unity (99.99 % for pure Gaussian).

### 4.3 Contamination estimation

#### 4.3.1 Protons

Protons were clearly identified without any contamination below 1.7 GV by the mass selection, as shown in Fig. 4. Above 1.7 GV, however, light particles such as positrons and muons contaminate proton’s $`1/\beta `$-band, and above 4 GV deuterons ($`D`$’s) start to contaminate it.

To distinguish protons from muons and positrons above 1.7 GV, we required that light output of the aerogel Cherenkov counter should be smaller than the threshold. This requirement rejected 96.5 % of muons and positrons while keeping the efficiency for protons as high as 99.5 %. The aerogel Cherenkov cut was applied below 3.7 GV, above which Cherenkov output for protons begins to increase rapidly.

Contamination of muons in proton candidates after the aerogel Cherenkov cut was estimated and subtracted by using calculated muon fluxes , which was normalized to the observed fluxes below 1.0 GeV where positive muons were clearly separated from protons. The positron contamination was calculated from the normalized positive muon fluxes and observed $`e`$/$`\mu `$ ratios. The resultant positive muon and positron contamination in the proton candidates was less than 0.8 % and 4.1 % below and above 3.7 GV, respectively, at 26.4 g/cm2 where observed $`(\mu `$+$`e)`$/$`p`$ ratio is largest during the experiment. The error in this subtraction was estimated from the ambiguity in the normalization of the muon spectra to be less than 0.2 % and 0.8 % below and above 3.7 GV, respectively.

Since the geomagnetic cut-off rigidity is 4.2 GV, most of deuterons contaminating proton candidates above 4 GV are considered to be primary cosmic-ray particles. The primary $`D`$/$`p`$ ratio was estimated by our previous measurement carried out in 1998 – 2000 at Lynn Lake, Manitoba, Canada, where the cut-off rigidity is as low as 0.5 GV. The $`D`$/$`p`$ ratio was found to be 2 % at 3 GV. No subtraction was made for deuteron contamination because there is no reliable measurement of the deuteron flux above 4 GV. Therefore, above 4 GV hydrogen nuclei selected which included a small amount of deuterons. The $`D`$/$`p`$ ratio at higher energy is expected to decrease  due to the decrease in escape path lengths of primary cosmic-ray nuclei  and the deuteron component is as small as the statistical error of the proton flux.

#### 4.3.2 Helium nuclei

Helium nuclei were clearly identified by using both upper and lower TOF d$`E`$/d$`x`$ as shown in Fig. 3. Landau tail in proton’s d$`E`$/d$`x`$ might contaminate the helium d$`E`$/d$`x`$ band, however this contamination from protons was as small as 3$`\times `$10-4. It was estimated by using another sample of proton data selected by d$`E`$/d$`x`$ in the JET chamber. No background subtraction was made for helium. Obtained helium fluxes include both 3He and 4He.

#### 4.3.3 Muons

Electrons and pions could contaminate muon candidates. To estimate electron contamination, we used the d$`E`$/d$`x`$ information inside the lower TOF counters covered with the lead plate. We calculated lower TOF d$`E`$/d$`x`$ distribution for muons and electrons using the Monte Carlo simulation. The most adequate e/$`\mu `$ ratio was estimated by changing weights both for muons and electrons so as to reproduce the observed d$`E`$/d$`x`$ distribution. The simulated distributions well-agreed with the observed data as shown in Fig. 5. The electron contamination was about 10 % of muon candidates at 0.5 GV, and less than 1% above 1 GV. Since in the lower energy region, the difference between the d$`E`$/d$`x`$ distributions for electrons and muons are small, rejection power of d$`E`$/d$`x`$ selection are lower. The error was estimated during the fitting procedure to be about 5 % at 0.5 GV and 0.1 % at 8.0 GV. Accuracy of the electron subtraction was limited by poor statistics of electron events. For the pion contamination we made no subtraction, thus the observed muon fluxes include pions. According to a theoretical calculation , the $`\pi `$/$`\mu `$ ratio at a residual atmosphere of 3 g/cm2 through 10 g/cm2 was less than 3 % at 1 GV and less than 10 % at 10 GV.

Between 1.6 and 2.6 GV, protons contamination in the positive muon candidates and its error were estimated by fitting the 1/$`\beta `$ distribution with a double-Gaussian function as shown in Fig 6. The estimated ratio of proton contamination in muon candidates was less than 3.3 $`\pm `$ 1.3 %.

### 4.4 Corrections

In order to determine the proton, helium and muon spectra, ionization energy loss inside the detector material, live-time and geometrical acceptance need to be estimated.

The energy of each particle at the top of the instrument was calculated by summing up the ionization energy losses inside the instrument with tracing back the event trajectory. The total live data-taking time was measured exactly to be 40,601 seconds by counting 1 MHz clock pulses with a scaler system gated by a “ready” status that control the first level trigger. The geometrical acceptance defined for this analysis was calculated as a function of rigidity by using simulation technique . In the high rigidity region where a track of the particle is nearly straight, the geometrical acceptance is 0.097 $`\mathrm{m}^2\mathrm{sr}`$ for protons and helium nuclei, and 0.030 $`\mathrm{m}^2\mathrm{sr}`$ for muons. The acceptance for muons is about 1/3 of that for protons and helium nuclei, because we required muons to pass through the lead plate while protons and helium nuclei were required not to pass through the lead plate. The simple cylindrical shape and the uniform magnetic field make it simple and reliable to determine the precise geometrical acceptance. The error which arose from uncertainty of the detector alignment was estimated to be 1 %.

## 5 Results and discussions

The proton and helium fluxes in energy ranges of 0.5–10 GeV/n and muon flux in 0.5 GeV/$`c`$–10 GeV/$`c`$, at small atmospheric depths of 4.5 g/cm2 through 28 g/cm2, have been obtained from the BESS-2001 balloon flight. The results are summarized in Table 1. The statistical errors were calculated as 68.7 % confidence interval based on Feldman and Cousins’s “unified approach” . The overall errors including both statistic and systematic errors are less than 8 %, 10% and 20 % for protons, helium nuclei and muons, respectively. The obtained proton and helium spectra are shown in Fig. 7. Around at 3.4 GeV for protons and 1.4 GeV/n for helium nuclei, a geomagnetic cut-off effect is clearly observed in their spectra. The proton spectrum measured by the AMS experiment in 1998  at the similar geomagnetic latitude ($`0.7<\mathrm{\Theta }_M<0.8`$, where $`\mathrm{\Theta }_M`$ is the corrected geomagnetic latitude ) is also shown in Fig. 7. The AMS measured proton spectra in space (at an altitude of 380 km), which are free from atmospheric secondary particles111The AMS observed substantial “second” spectra below the geomagnetic cut-off. Most of them follow a complicated trajectory in the Earth’s magnetic field, and could not be observed at balloon altitude.. In the BESS results the atmospheric secondary spectra for protons below 2.5 GeV are observed. Fig. 8 shows the observed proton and helium fluxes as a function of the atmospheric depth. Below 2.5 GeV the proton fluxes clearly increase as the atmospheric depth increases. It is because the secondary protons are produced in the atmosphere. In the primary fluxes above the geomagnetic cut-off, the fluxes attenuate as the atmospheric depth increases. In this energy region, the production of the secondary protons is much smaller than interaction loss of the primary protons. This is because the flux of parent particles of secondary protons is much smaller due to the steep spectrum of primary cosmic rays.

Figs. 10 and 10 show the observed muon spectra together with theoretical predictions. The predictions were made with the hadronic interaction model, DPMJET-III , which was used for the evaluation of atmospheric neutrino fluxes . The obtained proton fluxes were used to reproduce the primary cosmic-ray fluxes in the calculation. Fig. 11 shows the observed muon fluxes as a function of atmospheric depth together with the calculated fluxes. The calculated fluxes show good agreement with the observed data. Further detailed study of the hadronic interaction models will be discussed elsewhere.

## 6 Conclusion

We made precise measurements of cosmic-ray spectra of protons, helium nuclei and muons at small atmospheric depths of 4.5 through 28 g/cm2, during a slow descending period of 12.4 hours, in the BESS-2001 balloon flight at Ft. Sumner, New Mexico, USA. We obtained the proton and helium fluxes with overall errors of 8 % and 10 %, respectively, in an energy region of 0.5 – 10 GeV/n. The muon fluxes were obtained with an overall error of 20 % in a momentum region of 0.5 – 10 GeV/$`c`$. The results provide fundamental information to investigate hadronic interactions of cosmic rays with atmospheric nuclei. The measured muon spectra showed good agreement with the calculations by using the DPMJET-III hadronic interaction model. The understanding of the interactions will improve the accuracy of calculation of atmospheric neutrino fluxes.

We would like to thank NASA and the National Scientific Balloon Facility (NSBF) for their professional and skillful work in carrying out the BESS flights. We are indebted to M. Honda and T. Kajita of ICRR, the University of Tokyo for their kindest cooperation for Monte Carlo simulations and theoretical interpretations. We would like to thank ISAS and KEK for their continuous support and encouragement for the BESS experiment. The analysis was performed with the computing facilities at ICEPP, the University of Tokyo. This experiment was supported by Grants-in-Aid (12047206 and 12047227) from the Ministry of Education, Culture, Sport, Science and Technology, MEXT.