text:
# 1 Introduction

## 1 Introduction

Understanding the QCD phase diagram in the temperature – chemical potential ($`T,\mu `$) plane has many important implications in cosmology, in astrophysics and in the phenomenology of heavy-ion collisions (for a review, see ). In particular, it is extremely important to localize with high accuracy the critical lines in this plane and to determine if the transition across them is first order, second order or just a crossover .

The formulation of the theory on a space-time lattice offers a unique tool for such task. However, at non-zero chemical potential the determinant of the fermion matrix becomes complex, thus preventing to perform standard Monte Carlo simulations. One way out is to take advantage of physical fluctuations in the thermal ensemble generated at $`\mu =0`$ in order to extract from it information at (small) non-zero $`\mu `$, after suitable reweighting . However, it is not possible to know a priori to what extent the reweighted ensemble overlaps the true one at non-zero $`\mu `$.111In the (2+1)-dimensional U(1) Gross-Neveu model a clear disagreement was evidenced  in the $`\mu `$-dependence of the chiral condensate and of the number density between the determinations from reweighting and from standard Monte Carlo techniques. Of course, this does not imply the ineffectiveness of the method in other theories and in different temperature regimes.

An alternative approach is based on the idea to perform numerical simulations at imaginary chemical potential , for which the fermion determinant is again real, and to analytically continue the results to real $`\mu `$ .222Some interesting arguments about the QCD phase diagram for imaginary chemical potential have been discussed in Ref. . In practice, the method consists in fitting the numerical data obtained for imaginary chemical potential $`\mu =i\mu _I`$ by a polynomial, which is then prolongated to real $`\mu `$. The advantage of this method lies in the fact that its applicability is not restricted by large lattice volumes, since it makes no use of reweighting. Therefore, size scaling analysis can be made after the analytical continuation to clarify the nature of the transition across the critical lines. On the other side, a limitation comes from the presence of non-analyticities in $`\mu `$, thus implying, as we will see later on, that the method works practically only for $`\mu T`$.

In particular, if non-analyticities are singled out for imaginary $`\mu `$ and even if numerical data are obtained inside a region of analyticity in the $`(T,\mu _I)`$ plane, it is necessary that during the “rotation” from $`\mu =i\mu _I`$ to real $`\mu =\mu _R`$ no other non-analyticities are met. It is understood that non-analyticities show up only in the infinite volume limit and that all observables are analytic on finite lattices. However, the presence of non-analyticities in the thermodynamical limit is indicated on (large) finite volumes by sharp changes in the behavior of physical observables, which is practically difficult to interpolate by polynomials, if not by including a large number of terms. The papers where the method of analytical continuation has been used so far limited the attention to regions of temperature and chemical potential where this situation does not occur.

Of course, the compatibility between the results from two very different approaches, such as reweighting and analytical continuation, would give reasonable confidence on the validity of both. Such consistency has been shown in the case of 4-flavor QCD in Ref. .

Nevertheless, a test of both methods by comparison with direct determinations from standard Monte Carlo is interesting, especially because it could shed light also on the applicability outside the regions usually considered in literature. For this purpose, it is necessary to turn to a theory where direct Monte Carlo simulations in presence of real chemical potential are feasible.

It is well known that the determinant of the fermion matrix of 2-color QCD or SU(2) gauge theory is real even in presence of a real chemical potential (see, for instance ). Moreover, although being very different in the structure of the spectrum , 2-color QCD shares many similarities with true QCD, therefore represents the best candidate for the mentioned test.

In this paper we adopt 2-color QCD to compare determinations of the Polyakov loop and of the chiral condensate from the method of analytical continuation with results from direct Monte Carlo simulations. In particular, we characterize different regions of the phase diagram in the $`(T,\mu _I)`$ plane according to the applicability of this method. Our findings can be helpful to apply and to understand the results of the method of analytical continuation in the theory of physical interest, i.e. in true QCD.

The paper is organized as follows: in Section 2 we describe the structure of the phase diagram of 2-color QCD in the $`(T,\mu _I)`$ plane and discuss the implications on the method of analytical continuation; in Section 3 we present the numerical results; in Section 4 we draw some conclusions and anticipate some future developments.

## 2 The phase diagram of 2-color QCD

In order to draw the phase diagram of 2-color QCD we make use of known theoretical and numerical results on the behavior with temperature and chemical potential of the Polyakov loop $`L`$ and of the chiral condensate $`\overline{\psi }\psi `$. These objects represent true order parameters only in two limiting cases: the Polyakov loop is the order parameter of the confinement-deconfinement transition in the limit of infinite quark masses, while the chiral condensate is that of the chiral transition in the limit of vanishing quark masses. However, there is a broad numerical evidence that they both exhibit a rapid variation in correspondence of certain values of the parameters of the theory, which therefore are taken as critical ones.

We start the discussion from the case of imaginary chemical potential, $`\mu =i\mu _I`$. Roberge and Weiss have shown  that the partition function of any SU(N) theory is periodic in the parameter $`\theta =\mu _I/T`$ (the Boltzmann constant is set equal to one) as consequence of a remnant of the Z(N) symmetry valid in absence of fermions. Moreover, they have shown by weak coupling perturbation theory that above a certain temperature $`T_E`$ the free energy and the Polyakov loop show cusps or discontinuities for $`\theta =2\pi (k+1/2)/N`$, with $`k=0,1,\mathrm{}`$. This suggests that there are first order vertical critical lines in the $`(T,\theta )`$ plane located at $`\theta =2\pi (k+1/2)/N`$ and extending from $`T_E`$ to infinity (see Fig. 1). They argue also that below $`T_E`$, there is no phase transition for any value of $`\theta `$, this being true also in the physical case of $`\theta =\mu =0`$.

The latter conclusion cannot be true for all values of $`n_f`$ and of the fermion masses. In the chiral limit of QCD, there are arguments  according to which a chiral phase transition for $`\mu =0`$ does exist and is first order for $`n_f3`$ and second order for $`n_f=2`$. For SU(2), where a strong coupling analysis could be applied, according to Ref.  there is a second order chiral phase transition for $`n_f=1`$ and $`n_f=2`$. For non-vanishing masses, this phase transition probably becomes a crossover. Furthermore, a numerical study on a $`12^3\times 4`$ lattice in SU(2) with $`n_f=8`$ and $`am=0.07`$ has shown a peak in the susceptibility of the chiral condensate at $`\beta _c1.41`$, suggesting the possibility of a chiral phase transition near this point .

In view of all this, the worst case for the purposes of the analytical continuation (but the most interesting from the point of view of the present study) is that the Roberge-Weiss (RW) critical line at $`\theta =\pi /N`$ could be prolongated downwards and bent up to reaching the $`\theta =0`$ axis in correspondence of a critical temperature $`T_c`$ (see also ). For the reasons mentioned above, this prolongation of the RW critical line will be referred to in the following as the chiral critical line . The RW periodicity and the symmetry for $`\mu \mu `$, resulting from CP invariance, imply that knowing the expectation value of an observable for $`T`$ and $`\theta `$ inside the strip $`\{0\theta \pi /N,\mathrm{\hspace{0.33em}0}T<\mathrm{}\}`$ in the $`(T,\theta )`$ plane is enough to fix it on the remaining part of the plane. Applying this argument also to the critical lines, we end up with the qualitative structure of the phase diagram drawn in Fig. 2. The actual nature of the transition across the branches of the critical lines of Fig. 2 which lie below $`T_E`$ is to a large extent unknown. Other scenarios of phase diagrams are possible, depending on $`n_f`$ and on the fermion masses, such as, for instance, that where a chiral critical line emanating from $`T_c`$ ends somewhere before joining the vertical lines .

Inside the mentioned strip we can single out two regions where the (non-perturbative) determinations of an observable can be interpolated by a Taylor series in $`\theta ^2=(\mu _I/T)^2`$: the region above and the region below the chiral critical line. This Taylor series can be analytically continued to real $`\mu =\mu _R`$ by the “rotation” $`\mu _Ii\mu _R`$, if no non-analyticities are met during this rotation. Such a condition is not a priori guaranteed. It is evident, that the crucial point is how the critical line in the $`(T,\theta )`$ plane, implicitely defined by $`T=f(\mu _I/T)`$, evolves during the rotation from imaginary to real $`\mu `$. A reasonable assumption is that, at the end of the rotation, the critical line in the plane with real $`\mu `$ is given by $`T=f(i\mu _R/T)g(\mu _R/T)`$, with $`g`$ a monotonically decreasing function in $`\mu _R/T`$. Another assumption, which needs confirmation by numerical tests, is that the critical line evolves smoothly during the rotation, i.e. it spans in the $`(T,\mu _R/T,\mu _I/T)`$ space a 2-dimensional surface, connecting the critical line in the $`(T,\mu _I/T)`$ plane with the one in the $`(T,\mu _R/T)`$ plane, which does not possess bumps or ridges or valleys.

If these assumptions are true and the scenario of Fig. 2 holds, we can identify three regions in the strip $`\{0\theta \pi /N,\mathrm{\hspace{0.33em}0}T<\mathrm{}\}`$ of the $`(T,\mu _I/T)`$ plane according to the possibility to apply the analytical continuation:

* $`T>T_E`$: analytical continuation possible for $`0\mu _I/T<\pi /N`$;
* $`T_c<T<T_E`$: analytical continuation possible for $`0\mu _I/T<\mu _I^{}/T`$, where $`\mu _I^{}/T`$ is the solution of $`T=f(\mu _I^{}/T)`$, i.e. the value of $`\mu _I/T`$ on the chiral critical line at the given temperature;
* $`T<T_c`$: analytical continuation possible for $`0\mu _I/T<i\mu _R^{}/T`$, where $`\mu _R^{}/T`$ is the solution of $`T=g(\mu _R^{}/T)`$, i.e. the value of $`\mu _R/T`$ on the critical line of the $`(T,\mu _R/T)`$ plane at the given temperature.

So far the method of analytical continuation has been applied, implicitely making the assumptions we have mentioned, to the chiral condensate in the region $`T>T_E`$ for QCD with $`n_f=4`$  and to the critical line itself for QCD with $`n_f=4`$  and with $`n_f=2,3`$ . In the case of QCD with $`n_f=4`$, the analytically continued critical line has been found in agreement with the determination of Ref. .

In the next Section we will compare the scenario sketched in the above list with numerical simulations in all the three temperature regions. As a test-field theory we use 2-color QCD with $`n_f=8`$ and as observables we adopt the Polyakov loop $`L`$ and the chiral condensate $`\overline{\psi }\psi `$.

## 3 Numerical results

We consider the formulation of 2-color QCD or SU(2) gauge theory on a $`N_\sigma ^3\times N_\tau `$ lattice with staggered fermions for non-zero temperature and chemical potential of the baryonic density.

The finite temperature is realized by compactifying the time direction and by imposing (anti)periodic boundary conditions on (fermion) boson fields. The connection between the temperature $`T`$ and the lattice size in the time direction is $`T=1/(N_\tau a)`$, where $`a=a(\beta )`$ is the lattice spacing and $`\beta 4/g^2`$ is the bare coupling constant. The chemical potential is introduced by the replacements

$`U_\tau `$ $``$ $`e^{a\mu }U_\tau `$
$`U_\tau ^{}`$ $``$ $`e^{a\mu }U_\tau ^{},`$ (1)
where $`U_\tau `$ and $`U_\tau ^{}`$ are the gauge link variables in the forward and backward time direction, respectively. When $`\mu `$ is purely imaginary these replacements amount to add a constant background U(1) field to the original theory. For SU(2) the fermionic determinant, appearing in the partition function after integration of the fermion fields, is real for any complex value of the chemical potential $`\mu `$, thus allowing the Monte Carlo importance sampling without any limitation.

We adopted in this work the same hybrid simulation algorithm of Ref. . Since for $`\mu 0`$ it is not possible to use the even-odd partitioning, we were forced to simulate eight degenerate continuum flavors. We performed all the simulations on a $`8^3\times 4`$ lattice, setting the quark bare mass to $`am=0.07`$ and the microcanonical time step to $`dt=0.02`$. We made refreshments of the gauge fields every 5 steps of the molecular dynamics and of the pseudofermion fields every 3 steps. In order to reduce autocorrelation effects, we made “measurements” every 50 steps and analyzed data by the jackknife method combined with binning.

Before presenting the numerical results, it is convenient to redraw the phase diagram discussed in the previous Section in units more suitable for the lattice. In particular, on the vertical axis we put the inverse coupling $`\beta `$ which is connected to the temperature $`T`$ by a monotonic function in the region of interest (increasing $`\beta `$ means increasing temperature). On the horizontal axis we put the chemical potential in lattice units, $`\widehat{\mu }a\mu `$. The qualitative structure of the phase diagrams in the $`(\beta ,\widehat{\mu }_I)`$ and in the $`(\beta ,\widehat{\mu }_R)`$ planes is the same as in the $`(T,\mu _I/T)`$ and in the $`(T,\mu _R/T)`$ planes. In the new units, however, the RW transition lines in the $`(\beta ,\widehat{\mu }_I)`$ plane appear at $`\widehat{\mu }_I=2\pi (k+1/2)/(NN_\tau )`$. For $`N=2`$ the first RW transition line is located at $`\widehat{\mu }_I=\pi /(2N_\tau )`$. Of course, the scenario illustrated at the end of the previous Section must be reformulated in terms of $`\beta `$, $`\beta _E`$ and $`\beta _c`$.

An estimate of the value of $`\beta _c`$ can be deduced from the literature. Indeed, in Ref. , a study of the chiral condensate on a $`12^3\times 4`$ lattice for bare quark mass $`am=0.07`$ lead to the result $`\beta _c=1.41\pm 0.03`$. Since we have a different lattice, we expect a small difference between our $`\beta _c`$ and the determination of Ref. . However, for our purposes we can safely adopt this value of $`\beta _c`$.

As a first step we took two $`\beta `$ values, one (much) above $`\beta _c`$ and likely above $`\beta _E`$, the other below $`\beta _c`$ in order to possibly single out a different behavior in $`\widehat{\mu }_I`$ of lattice observables across $`\widehat{\mu }_I=\pi /(2N_\tau )`$. We took $`\beta `$=1.80 for the higher $`\beta `$ value and $`\beta =1.35`$ for the lower and determined the Polyakov loop and the chiral condensate for $`\widehat{\mu }_I`$ ranging between 0 and the value corresponding to the second RW critical line, i.e. $`3\pi /(2N_\tau )`$. The results are shown in Figs. 3 and 4. We can clearly see the RW periodicity in all cases; however, for $`\beta =1.80`$ we have discontinuity for the Polyakov loop and cusp for the chiral condensate across the first RW critical line, while for $`\beta =1.35`$ the behavior is smooth for both observables. We have verified that the overshooting of the first RW critical line seen in Fig. 3 is due to an hysteresis effect, which confirms the first order nature of the transition.

In order to proceed with the test of the method of analytical continuation we need also an estimation of $`\beta _E`$. There is a nice procedure to do it, suggested in Ref. . It relies on performing simulations at a critical value of $`\widehat{\mu }_I`$, say $`\widehat{\mu }_I=\pi /(2N_\tau )`$, and on changing $`\beta `$ from a very high value downwards. If the zero field configuration is chosen as the starting one, the system remains always in the phase on the left of the RW line until $`\beta >\beta _E`$, if the flip to the other phase is prevented by a large space volume. Due to our limited computational power, we are constrained to a space volume which is too small to prevent flips between the two phases. An alternative procedure to determine $`\beta _E`$ could be to study the behavior of the Polyakov loop (or of its phase) and of the chiral condensate across the RW critical lines for different $`\beta `$ values and to see when the discontinuities or the cusps are smoothed out. Another possibility could be to exploit the hysteresis effect in the Polyakov loop or in the chiral condensate determinations across the RW critical line. The procedure consists in building two sequences of expectation values of the Polyakov loop; in the first (second) sequence the starting configuration of any new run is taken to be a thermal equilibrium configuration being always in the phase on the left (right) side of the first RW transition. Across the RW critical line the comparison of these two sequences of data shows a hysteresis effect as long as $`\beta >\beta _E`$. This effect disappears below $`\beta _E`$. Of course, both the above procedures demand for a high statistics and, so far, we can only say that $`\beta _E`$ is bounded between about 1.53 and about 1.57. This level of accuracy is, however, enough for the purposes of the present study.

We turn now to the main part of this work, i.e. the test if the scenario presented at the end of the previous Section holds. We selected three values of $`\beta `$ (1.90, 1.45, 0.90), each representative of one of the regions $`\beta >\beta _E`$, $`\beta _c<\beta <\beta _E`$ and $`\beta <\beta _c`$. For each $`\beta `$ we have performed simulations for both real and imaginary chemical potential, varying $`\widehat{\mu }_I`$ and $`\widehat{\mu }_R`$ between 0 and $`\pi /(2N_\tau )`$ and have determined expectation values of the Polyakov loop and of the chiral condensate. Then, we have interpolated the data obtained for imaginary $`\mu `$ with a truncated Taylor series of the form $`a+b\widehat{\mu }_I^2+c\widehat{\mu }_I^4`$. After “rotating” this polynomial to real chemical potential, thus leading to $`ab\widehat{\mu }_R^2+c\widehat{\mu }_R^4`$, we have compared it with the determinations obtained directly for real $`\mu =\mu _R`$.

For $`\beta =1.90`$ we have interpolated with the fourth order polynomial the data obtained for $`0\widehat{\mu }_I\pi /(2N_\tau )`$ and have found that the rotated polynomial perfectly interpolates data obtained directly for real $`\mu =\mu _R`$ in the same range (see Fig. 5). The uncertainty in the rotated polynomial has been obtained by propagating in standard way the uncertainty on the fitted coefficients $`a`$, $`b`$ and $`c`$. We have verified that the same situation occurs for other values of $`\beta `$ larger than $`\beta _E`$. From this outcome we can conclude that the analytical continuation works in the region $`\beta >\beta _E`$ for all the $`\widehat{\mu }_I`$’s before the first RW critical line.

For $`\beta =1.45`$ we are in the region $`\beta _c<\beta <\beta _E`$ for which, if the scenario presented at the end of the previous Section holds, we expect that by varying $`\widehat{\mu }_I`$ we meet the chiral critical line, while no transitions should be met at the same $`\beta `$ by varying $`\widehat{\mu }_R`$. Therefore this time we have interpolated with the fourth order polynomial the data obtained for the real chemical potential $`0\widehat{\mu }_R\pi /(2N_\tau )`$ and have rotated this polynomial to its counterpart in $`\widehat{\mu }_I`$. The comparison with data obtained directly for imaginary chemical potential shows agreement for $`\widehat{\mu }_I`$ below $`0.18`$, which can be taken as an estimate of the chiral critical value in the $`(\beta ,\widehat{\mu }_I)`$ plane at the given $`\beta `$ (see Fig. 6).

For $`\beta =0.90<\beta _c`$, we expect analyticity for all possible values of $`\widehat{\mu }_I`$, while the chiral critical line in the $`(\beta ,\widehat{\mu }_R)`$ plane is met by varying $`\widehat{\mu }_R`$ at the given value of $`\beta `$. Therefore we have interpolated with the fourth order polynomial the data obtained for $`0\widehat{\mu }_I\pi /(2N_\tau )`$ and have rotated this polynomial to its counterpart in $`\widehat{\mu }_R`$. The comparison with data obtained directly for real chemical potential shows agreement for $`\widehat{\mu }_R`$ below $`0.12`$, which can be taken as an estimate of the critical value in the $`(\beta ,\widehat{\mu }_R)`$ plane at the given $`\beta `$ (see Fig. 7).

We can conclude from these results that the scenario of the end of the previous Section is consistent, thus giving support to the assumptions underlying it.

## 4 Conclusions and outlook

In this paper we have presented a possible structure of the phase diagram in the temperature – imaginary chemical potential in 2-color QCD and have argued how this structure can be continued to real chemical potential. Then, we have shown that the scenario we have described is consistent with Monte Carlo numerical determinations of the Polyakov loop and of the chiral condensate for real and imaginary chemical potential.

In particular, we have found that the method of analytical continuation, largely exploited now in the literature for the physically interesting case of QCD, works fine within the restrictions posed by the presence of non-analyticities.

High precision determinations were not the aim of this paper. We limited ourselves to give convincing numerical support to the proposed scenario, postponing massive calculations to forthcoming works.

We plan, moreover, to further develop the present ideas by (a) determining the location of the critical lines for real and imaginary chemical potential by the standard method of susceptibilities, (b) verifying that the critical lines for real chemical potential can be obtained from analytical continuation of those for imaginary chemical potential, (c) following the “evolution” of the critical lines during the rotation from imaginary to real chemical potential by simulations with complex chemical potential.

Acknowledgments We have the pleasure to thank M.-P. Lombardo for many stimulating discussions and for providing us with the Monte Carlo code. We acknowledge also several interesting conversations with R. Fiore.


format: text

text:
# Observation of Large Atomic-Recoil Induced Asymmetries in Cold Atom Spectroscopy

## Abstract

The atomic recoil effect leads to large (25 %) asymmetries in simple spectroscopic investigations of Ca atoms that have been laser-cooled to $`10\mu `$K. Starting with spectra from the more familiar Doppler-broadened domain, we show how the fundamental asymmetry between absorption and stimulated emission of light manifests itself when shorter spectroscopic pulses lead to the Fourier transform regime. These effects occur on frequency scales much larger than the size of the recoil shift itself, and have not been observed before in saturation spectroscopy. These results are relevant to state-of-the-art optical atomic clocks based on freely expanding neutral atoms.

When we recently started using much colder Ca atoms for our ultra-high resolution optical clock spectroscopy, we discovered unexpectedly large asymmetries in the resulting spectra, which we found to be the direct result of atomic recoil.Curtis et al. (2003) Most surprising was that these effects appeared on a frequency scale more than ten times larger than the size of the recoil effect itself. In this letter we use a simple saturation spectroscopic configuration to show how the small effect of atomic recoil can cause these large asymmetries in the spectra of samples of freely expanding cold atoms. Not only are these effects of interest from the point of view of basic physics, they also have important implications for future optical clocks based on laser-cooled neutral atoms.Curtis et al. (2003); Wilpers et al. (2002)

The effects of atomic recoil have been of considerable interest to laser spectroscopists for more than 30 years Kol’chenko et al. (1969) and are the foundation of laser cooling of atoms. The first experimental evidence that saturation absorption spectroscopy produced a recoil splitting in the Doppler-free spectra of atomic and molecular lines was demonstrated by Hall, Bordé, and Uehara using a high-resolution laser spectrometer.Hall et al. (1976) Since then there have been several beautiful demonstrations of recoil splitting via various forms of saturation spectroscopy (see for example refs. Barger et al. (1979); Sterr et al. (1992); Riehle et al. (1988); Bagayev et al. (1989)). In all of these studies the closely spaced recoil components were superimposed on a broad Doppler background, thereby obscuring the recoil-induced asymmetries inherent in light-atom interactions. In the work presented here, the width of the Doppler background is only seven times that of the recoil splitting, thereby clearly exposing the asymmetric aspect of the recoil components in the Doppler-broadened regime (see Fig 1).

Moreover, by gradually changing the resolution, we show how this leads to large spectral asymmetries in the Fourier-transform limit. These effects differ from asymmetries described in earlier papers, which resulted from spontaneous emission or collisional quenching.Kol’chenko et al. (1969); Hall et al. (1976)

As described elsewhere, the recoil doublet can be readily understood from conservation of momentum and energy.Kol’chenko et al. (1969); Hall et al. (1976) When an atom initially at rest absorbs a photon of frequency $`\nu `$, it recoils with a velocity, $`v_\mathrm{r}=h\nu /mc`$, where $`m`$ is the atomic mass. For a two-level atom with an energy level spacing of $`E_0=h\nu _0`$, this absorption resonance occurs at a frequency,

$$\nu _\mathrm{A}=\nu _0+\frac{1}{2}\frac{h\nu _0^2}{mc^2}=\nu _0+\frac{1}{2}\nu _\mathrm{r}.$$
(1)
Eq. 1 shows that the incident photon needs to supply energy for both the atomic excitation and recoil. For an atom at rest but starting in the excited state, a similar analysis shows that the resonance for stimulated emission occurs at a frequency $`\nu _{\mathrm{SE}}=\nu _01/2\nu _\mathrm{r}`$. This resonance is red-shifted relative to that for absorption by $`\nu _\mathrm{r}`$, the splitting of the recoil doublet, which is typically tens of kilohertz for optical transitions. This frequency splitting ensures that the alternating absorption and stimulated emission cycles in Rabi flopping are simultaneously resonant even when the Doppler shift associated with the atomic recoil is included. Thus, one cannot observe the recoil splitting with a singly-passed laser beam. Instead, one must use two counter-propagating laser beams, the usual configuration for saturation spectroscopy. Then the atomic recoil due to photon absorption from one beam pushes an atom toward the counter-propagating beam, effectively reversing the sign of the recoil shift. In this way it is possible to see two distinct sub-Doppler features split by $`\nu _\mathrm{r}`$, although this small splitting can only be resolved in ultra-high-resolution experiments.

In this investigation, we take advantage of the capabilities of our optical clock apparatus Curtis et al. (2003) to perform the simplest form of saturation spectroscopy. We excite a laser-cooled ($`T=10\mu `$K) sample of neutral Ca atoms using the closed 657 nm transition between the 1S0 ground state and the meta-stable 3P1 excited state, which has a lifetime of $`340\mu `$s. With transit-time broadening and spontaneous emission negligible, we can change the spectroscopic resolution just by changing the duration of the square probe pulses. Finally, we can excite the atoms sequentially with laser pulses of equal intensity and frequency detuning. Sequential excitation greatly simplifies the analysis, since it removes the possibility of events containing more than two photons from different directions that can distort the lineshapes Hall et al. (1976); Barger et al. (1979); Ishikawa et al. (1994).

We realize these experimental conditions with the following measurement cycle.Curtis et al. (2003) We first load Ca into a magneto-optic trap using the strongly-allowed 423 nm cooling transition. We then turn off the 423 nm light and use a 3-dimensional quenched narrow-line cooling scheme based on the clock transition to reduce the temperature of the atomic sample ($`10^6`$ atoms) to less than 10 $`\mu `$K.Curtis et al. (2003); Binnewies et al. (2001) When the atoms are cold, we switch off the trap and turn on a greater than 2 G magnetic bias field to perform the spectroscopy on the narrow $`m=0m=0`$ clock transition. We then probe the 657 nm clock transition with pulses derived from a cw diode laser, which is locked tightly to a narrow Fabry-Perot cavity fringe (producing a laser linewidth less than 70 Hz). Some light from this stabilized master laser is used to injection-lock a slave laser whose output is sent through two acousto-optic modulators to generate the pulses for the opposing directions. The deflected light is steered into optical fibers to spatially filter the beams. The beams coupled out of the fibers expand and are collimated to about 6 mm, so that the atoms see flat wavefronts (radius of curvature greater than 50 m). As much as 13 mW can be coupled into these beams, although we adjust the power to yield unit excitation on resonance (commonly called a $`\pi `$-pulse in Rabi flopping parlance). For these measurements we illuminate the atoms with a pulse from one direction, wait 6 $`\mu `$s (to make sure the first beam is completely turned off), and then illuminate the atoms with a pulse from the opposite direction. Finally, the fraction of atoms in the excited state is measured using a normalized shelving fluorescence detection technique.Curtis et al. (2003) We scan the probe frequency ($`\nu _\mathrm{L}`$) slowly (4 s sweep time) while continuously repeating the measurement cycle (duration 35 ms) to generate our spectra as a function of the laser detuning ($`\nu _\mathrm{L}\nu _0`$).

To make a connection with previous experiments, we first consider the Doppler-broadened regime, for which we choose a pulse duration ($`40\mu `$s) such that the spectroscopic resolution is about 21 kHz. This is less than the 150 kHz (FWHM) Doppler width of the 10 $`\mu `$K atoms and slightly less than the 23.1 kHz recoil splitting of the 657 nm clock transition. In Fig. 1 we show the resulting excitation spectrum, whose envelope is primarily determined by the Doppler background. In fact this envelope is simply twice the height of the curve seen with a single probe pulse, because for most laser frequencies the two counter-propagating laser beams excite non-overlapping velocity classes. However, at the absorption resonance frequency, $`\nu _\mathrm{L}=\nu _\mathrm{A}`$ (‘A’ in Fig. 1), there is not this doubling of the Doppler background, since the first laser pulse has already excited the majority of these atoms to the long-lived excited state. Note that the second pulse cannot de-excite many of these atoms either, since due to atomic recoil they have been shifted out of resonance. The second resonance (‘SE’ in Fig. 1), at $`\nu _\mathrm{L}=\nu _{\mathrm{SE}}`$, results from atoms that started with a velocity $`v=v_\mathrm{r}`$, and were thus excited by the first pulse (in the atoms’ rest frame, $`\nu _\mathrm{L}=\nu _\mathrm{A}`$); now at rest in the lab frame due to atomic recoil, they are resonant with the second pulse for stimulated emission, thereby reducing the net fraction of atoms excited.

The two dips resulting from these resonances are analogous to those seen in earlier experiments but with one important distinction: due to the small width of the Doppler background, we can readily see that the dips are asymmetrically located about the center of the Doppler background. Since the Doppler background results from absorption, it is naturally centered around the resonance at $`\nu _\mathrm{L}=\nu _\mathrm{A}`$, coincident with the absorption recoil dip. The dip associated with stimulated emission is located one recoil frequency ($`\nu _\mathrm{r}`$) below the absorption resonance, on the red side of the Doppler curve. It is important to emphasize that this asymmetry is not a result of the order of the laser pulses. If we reverse the temporal order, so that the pulse directions are reversed, we observe the identical lineshape, not its mirror image. Rather, this asymmetry is a fundamental feature of saturation spectroscopy, though one that is easily overlooked in experiments with broad Doppler backgrounds. As has been noted by other observers,Ishikawa et al. (1994); Dingler et al. (1994) this asymmetry can lead to undesired offsets in realizing optical frequency standards based on saturation absorption since the unperturbed line center of the transition (midway between the recoil components) is not centered on the background.

We now consider what happens to this spectrum as we broaden the Fourier spectrum of the probe pulse, not just beyond the recoil splitting but well beyond the width of the Doppler distribution itself. This regime has not previously been investigated experimentally, but is becoming important in state-of-the-art optical atomic clocks based on neutral atoms.Curtis et al. (2003); Wilpers et al. (2002) We access this regime by reducing the duration of the probe pulses (while maintaining the pulse area to keep the excitation probability constant). In Fig. 2, we show a set of spectroscopic lineshapes taken over probe durations ranging from 40 $`\mu `$s down to 2.6 $`\mu `$s.

Note that the fraction of atoms excited grows with decreasing resolution as our probe spectrum covers a larger fraction of the velocity distribution. For the shortest pulse length, the Fourier transform of the probe pulse has a spectral width approximately 2.5 times that of the Doppler distribution. As we see, changing the resolution from 21 kHz to 42 kHz (20 $`\mu `$s pulse duration) begins to obscure the recoil splitting, but the large asymmetry persists. As we move to 85 kHz resolution (10 $`\mu `$s pulse duration), the recoil splitting is no longer visible as the two dips merge into a single dip centered at $`\nu _0`$, but the spectrum now appears to consist of two peaks whose separation is determined by the spectroscopic resolution and whose amplitudes differ by more than 25 % ! The asymmetry persists for even the lowest resolution, where the maxima are separated by 280 kHz, more than 10 times that of the recoil splitting itself.

Examination of Fig. 2 provides an intuitive picture of how such a small effect can cause such large asymmetries. In the Fourier-transform regime we find the unusual spectroscopic condition where the widths of both the envelope and the dip are determined almost solely by the probe time, but are slightly offset from one another. This small offset leads to the envelope asymmetry by causing the dip to intersect the Doppler envelope at different heights on the two sides. Thus, the size of the asymmetry is related to the slope of the envelope multiplied by the size of the recoil effect.

Alternatively, we can think of this spectrum as resulting from two contributions: one symmetric (due solely to absorption) and one asymmetric (due to stimulated emission by the second pulse). This framework allows a simple model to describe our spectra well (an exact theoretical treatment of coherent saturation spectroscopy including multi-photon effects has been developed by Bordé and co-workersIshikawa et al. (1994); Bordé et al. (1984). We start by considering the excitation spectrum resulting from the first (idealized) square $`\pi `$-pulse of duration $`T`$ and Rabi frequency $`\mathrm{\Omega }`$ illuminating a sample of ground-state atoms initially at rest. This yields the well-known Rabi spectrum Rabi et al. (1939) for the excitation probability:

$$P(\mathrm{\Delta })=\frac{\mathrm{\Omega }^2T^2}{4}\left(\frac{\mathrm{sin}\left[\sqrt{\mathrm{\Omega }^2+\mathrm{\Delta }^2}T/2\right]}{\left[\sqrt{\mathrm{\Omega }^2+\mathrm{\Delta }^2}T/2\right]}\right)^2,$$
(2)
where $`\mathrm{\Delta }=\nu _\mathrm{L}\nu _\mathrm{A}`$ is the laser detuning from the absorption resonance. Illuminating the atoms with a square $`\pi `$-pulse from the opposite direction gives two sets of atoms to consider. First, there are atoms that remained in the ground state after the first pulse (a fraction equal to $`1P(\mathrm{\Delta })`$); they will be resonant with a second pulse at frequency $`\nu _\mathrm{L}=\nu _\mathrm{A}`$, so the above equation applies to these atoms as well. Second, the atoms that were excited by the first pulse (a fraction equal to $`P(\mathrm{\Delta })`$) now have a velocity $`v_\mathrm{r}`$ toward the second counter-propagating laser beam. The associated Doppler effect will shift the stimulated emission resonance down by one recoil, so these atoms will be resonant with light at frequency $`\nu _{\mathrm{SE}}^{}=\nu _{\mathrm{SE}}\nu _\mathrm{r}=\nu _0\frac{3}{2}\nu _\mathrm{r}`$, or $`\nu _{\mathrm{SE}}^{}=\nu _\mathrm{A}2\nu _\mathrm{r}`$. The fraction of these atoms transferred back to the ground state can then be derived from the same probability function $`P`$, but with the argument $`\mathrm{\Delta }^{}=\nu _\mathrm{L}\nu _{\mathrm{SE}}^{}=\mathrm{\Delta }+2\nu _\mathrm{r}`$.

To find the total fraction in the excited state after two pulses, we simply add the two contributions:

$$P_{\mathrm{final}}=P(\mathrm{\Delta })\left[1P(\mathrm{\Delta })\right]+\left[1P(\mathrm{\Delta }+2\nu _\mathrm{r})\right]P(\mathrm{\Delta }).$$
(3)
The first product in this expression gives the ground-state contribution and is symmetric (both pulses are resonant at the same frequency). The second term gives the excited-state contribution, but it is asymmetric due to the offset of $`2\nu _\mathrm{r}`$ between arguments of the multiplicands, thus yielding a net asymmetry for $`P_{\mathrm{final}}`$. We can easily connect this model with experiment by including the initial velocity distribution (via the detuning) and the unequal laser intensities seen by the atoms due to the spatial distribution of the atoms in the laser mode (via the Rabi frequency). The free parameters in our simulation were overall signal amplitude (imperfect normalization required a 10 - 15 % reduction in signal size) and atomic cloud size (which we fixed at the same value for all simulations). We measured the velocity distributions and the laser mode size separately and used the resulting values for the simulations. We see good agreement over a variety of probe resolutions (simulations are shown as solid lines in Figs. 1 and 3), although we see small differences in the wings at the lowest resolutions,

most likely resulting from the non-ideal square pulses used in the experiment.

Interestingly, we can isolate the second product in the expression for $`P_{\mathrm{final}}`$ experimentally by using recoil suppression.Riehle et al. (1988); Dingler et al. (1994) In this case we suppress the ground-state (absorption) contribution by illuminating the atoms with a resonant 423 nm pulse (duration 20 $`\mu `$s) after the first red probe pulse but before the second red pulse. This heats the atoms in the ground state to the point where excitation by the second pulse produces a nearly flat Doppler background, upon which sits the contribution from the atoms shelved in the excited state. Fig. 3 shows that the feature associated with stimulated emission is fully responsible for the lineshape asymmetry, as the absolute peak height differences are virtually identical. We emphasize that this asymmetric envelope persists in four-pulse Bordé-Ramsey saturation spectroscopy Bordé et al. (1984) in the Fourier-transform regime, which is used in ultra-high resolution studies and for cold atom optical clocks.Curtis et al. (2003); Wilpers et al. (2002) However, the good agreement between theory and experiments give us good confidence that these effects will not limit the accuracy of these clocks.

In summary, we have used ultra-cold two-level atoms and ultra-high resolution spectroscopy to probe the atomic recoil structure unique to saturated absorption spectroscopy at an unprecedented level. The resulting spectra clearly reveal the fundamental asymmetry in the location of recoil components on the Doppler background. In addition we have shown how this asymmetry in the recoil frequencies leads to large amplitude asymmetries in the Fourier transform limit, a regime that has not previously been investigated but is important for state-of-the-art precision metrology. Extending these studies to a regime with Ca atoms laser cooled to sub-recoil temperatures, as demonstrated by Curtis et al. Curtis et al. (2003), would allow us to individually address the recoil components via laser detuning. This could enhance the sensitivity of precision measurements such as that of the photon recoil or other atom interferometry experiments Weiss et al. (1993); Snadden et al. (1998); Gupta et al. (2002).

###### Acknowledgements.

We thank J. Bergquist and J. Kitching for their insightful comments and J. Hall for illuminating discussions. This work was supported in part by NASA and ONR-MURI. G. W. acknowledges support from the Alexander von Humboldt Foundation. Work of a US government agency; not subject to copyright.


text:
# 1 Reach in Λ_LL and Λ_RR at 95% C.L. vs. integrated luminosity ℒᵢₙₜ obtained from the model-independent analysis for 𝑒⁺𝑒⁻→{𝑒⁺𝑒⁻,𝜇⁺𝜇⁻} and 𝑒⁻𝑒⁻→𝑒⁻𝑒⁻ at 𝐸_{c.m.}=0.5 TeV, |𝑃⁻|=0.8 and |𝑃⁺|=0.6.

Contact interaction Lagrangians (CI) provide an effective framework to account for the phenomenological effects of new dynamics characterized by extremely high intrinsic mass scales $`\mathrm{\Lambda }`$, at the ‘low’ energies $`\sqrt{s}\mathrm{\Lambda }`$ attainable at current particle accelerators. For the Bhabha scattering process

$$e^++e^{}e^++e^{},$$
(1)
as well as for Møller scattering

$$e^{}+e^{}e^{}+e^{},$$
(2)
we consider the flavor-diagonal, helicity conserving $`eeff`$ contact-interaction effective Lagrangian :

$$_{\mathrm{CI}}=\frac{1}{1+\delta _{ef}}\underset{i,j}{}g_{\mathrm{eff}}^2ϵ_{ij}\left(\overline{e}_i\gamma _\mu e_i\right)\left(\overline{f}_j\gamma ^\mu f_j\right).$$
(3)
In Eq. (3): $`i,j=\mathrm{L},\mathrm{R}`$ denote left- or right-handed fermion helicities, $`\delta _{ef}=1`$ for processes (1) and (2) and, if we assumed lepton universality, the same Lagrangian, with $`\delta _{ef}=0`$, is relevant to the annihilation processes

$$e^++e^{}\mu ^++\mu ^{}.$$
(4)
The CI coupling constants in Eq. (3) are parameterized in terms of corresponding mass scales as $`ϵ_{ij}=\eta _{ij}/\mathrm{\Lambda }_{ij}^2`$ and, according to the previous remarks concerning compositeness, one assumes $`g_{\mathrm{eff}}^2=4\pi `$. Also, by convention, one takes $`|\eta _{ij}|=1`$ or $`\eta _{ij}=0`$, leaving the energy scales $`\mathrm{\Lambda }_{ij}`$ as free, a priori independent, parameters.

We notice that for the case of the Bhabha process (1), Eq. (3) envisages the existence of three independent CI models, each one contributing to individual helicity amplitudes or combinations of them, with a priori free, and nonvanishing, coefficients (basically, $`ϵ_{\mathrm{LL}},ϵ_{\mathrm{RR}}\mathrm{and}ϵ_{\mathrm{LR}}=ϵ_{\mathrm{RL}}`$ combined with the $`\pm `$ signs). The same is true for the Møller process (2). In general, apart from the $`\pm `$ possibility, for $`e^+e^{}\overline{f}f`$ with $`fe`$ there are four independent CI couplings, so that in the present case of processes (1) and (2) there is one free parameter less. Correspondingly, in principle, a model-independent analysis of the data should account for the situation where the full Eq. (3) is included in the expression for the cross section. Potentially, in this case, the different CI couplings may interfere and such interference could substantially weaken the bounds. To this aim, in the case of the processes (1), (2) and (4) at the Linear Collider (LC) considered here, a possibility is offered by initial beam polarization, that enables us to extract from the data the individual helicity cross sections (or their combinations) through the definition of particular, polarized integrated cross sections and, consequently, to disentangle the constraints on the corresponding CI constants . In this note, we wish to present a model-independent analysis of the CI that complements that of Refs. , and is based instead on the measurements of more ‘conventional’ observables (but still assuming polarized electron and positron beams) such as the differential distributions of the final leptons. We also make a comparison of the results from these three processes.

With $`P^{}`$ and $`P^+`$ the longitudinal polarization of the electron and positron beams, respectively, and $`\theta `$ the angle between the incoming and the outgoing electrons in the c.m. frame, the differential cross section of process (1), including $`\gamma `$ and $`Z`$ exchanges both in the $`s`$ and $`t`$ channels and the contact interaction (3), can be written in the following form :

$`{\displaystyle \frac{\text{d}\sigma (P^{},P^+)}{\text{d}\mathrm{cos}\theta }}`$ $`=`$ $`{\displaystyle \frac{(1+P^{})(1P^+)}{4}}{\displaystyle \frac{\text{d}\sigma _\mathrm{R}}{\text{d}\mathrm{cos}\theta }}+{\displaystyle \frac{(1P^{})(1+P^+)}{4}}{\displaystyle \frac{\text{d}\sigma _\mathrm{L}}{\text{d}\mathrm{cos}\theta }}`$ (5)
$`+`$ $`{\displaystyle \frac{1+P^{}P^+}{2}}{\displaystyle \frac{\text{d}\sigma _{\mathrm{LR},t}}{\text{d}\mathrm{cos}\theta }}.`$
In Eq. (5):

$`{\displaystyle \frac{\text{d}\sigma _\mathrm{L}}{\text{d}\mathrm{cos}\theta }}`$ $`=`$ $`{\displaystyle \frac{\text{d}\sigma _{\mathrm{LL}}}{\text{d}\mathrm{cos}\theta }}+{\displaystyle \frac{\text{d}\sigma _{\mathrm{LR},s}}{\text{d}\mathrm{cos}\theta }},`$
$`{\displaystyle \frac{\text{d}\sigma _\mathrm{R}}{\text{d}\mathrm{cos}\theta }}`$ $`=`$ $`{\displaystyle \frac{\text{d}\sigma _{\mathrm{RR}}}{\text{d}\mathrm{cos}\theta }}+{\displaystyle \frac{\text{d}\sigma _{\mathrm{RL},s}}{\text{d}\mathrm{cos}\theta }},`$ (6)
with

$`{\displaystyle \frac{\text{d}\sigma _{\mathrm{LL}}}{\text{d}\mathrm{cos}\theta }}`$ $`=`$ $`{\displaystyle \frac{2\pi \alpha ^2}{s}}\left|A_{\mathrm{LL}}\right|^2,{\displaystyle \frac{\text{d}\sigma _{\mathrm{RR}}}{\text{d}\mathrm{cos}\theta }}={\displaystyle \frac{2\pi \alpha ^2}{s}}\left|A_{\mathrm{RR}}\right|^2,`$
$`{\displaystyle \frac{\text{d}\sigma _{\mathrm{LR},t}}{\text{d}\mathrm{cos}\theta }}`$ $`=`$ $`{\displaystyle \frac{2\pi \alpha ^2}{s}}\left|A_{\mathrm{LR},t}\right|^2,{\displaystyle \frac{\text{d}\sigma _{\mathrm{LR},s}}{\text{d}\mathrm{cos}\theta }}={\displaystyle \frac{\text{d}\sigma _{\mathrm{RL},s}}{\text{d}\mathrm{cos}\theta }}={\displaystyle \frac{2\pi \alpha ^2}{s}}\left|A_{\mathrm{LR},s}\right|^2,`$ (7)
and

$`A_{\mathrm{RR}}`$ $`=`$ $`{\displaystyle \frac{u}{s}}\left[1+{\displaystyle \frac{s}{t}}+g_\mathrm{R}^2\left(\chi _Z(s)+{\displaystyle \frac{s}{t}}\chi _Z(t)\right)+2{\displaystyle \frac{s}{\alpha }}ϵ_{\mathrm{RR}}\right],`$
$`A_{\mathrm{LL}}`$ $`=`$ $`{\displaystyle \frac{u}{s}}\left[1+{\displaystyle \frac{s}{t}}+g_\mathrm{L}^2\left(\chi _Z(s)+{\displaystyle \frac{s}{t}}\chi _Z(t)\right)+2{\displaystyle \frac{s}{\alpha }}ϵ_{\mathrm{LL}}\right],`$
$`A_{\mathrm{LR},s}`$ $`=`$ $`{\displaystyle \frac{t}{s}}\left[1+g_\mathrm{R}g_\mathrm{L}\chi _Z(s)+{\displaystyle \frac{s}{\alpha }}ϵ_{\mathrm{LR}}\right],`$
$`A_{\mathrm{LR},t}`$ $`=`$ $`{\displaystyle \frac{s}{t}}\left[1+g_\mathrm{R}g_\mathrm{L}\chi _Z(t)+{\displaystyle \frac{t}{\alpha }}ϵ_{\mathrm{LR}}\right].`$ (8)
Here: $`\alpha `$ is the fine structure constant; $`t=s(1\mathrm{cos}\theta )/2`$, $`u=s(1+\mathrm{cos}\theta )/2`$ and $`\chi _Z(s)=s/(sM_Z^2+iM_Z\mathrm{\Gamma }_Z)`$ and $`\chi _Z(t)=t/(tM_Z^2)`$ represent the $`Z`$ propagator in the $`s`$ and $`t`$ channels, respectively, with $`M_Z`$ and $`\mathrm{\Gamma }_Z`$ the mass and width of the $`Z`$; $`g_\mathrm{R}=\mathrm{tan}\theta _W`$, $`g_\mathrm{L}=\mathrm{cot}2\theta _W`$ are the SM right- and left-handed electron couplings of the $`Z`$, with $`\theta _W`$ the electroweak mixing angle.

With both beams polarized, the polarization of each beam can be changed on a pulse by pulse basis. This would allow the separate measurement of the polarized cross sections for each of the four polarization configurations RR, LL, RL and LR, corresponding to the four sets of beam polarizations $`(P^{},P^+)=(P_1,P_2)`$, $`(P_1,P_2)`$, $`(P_1,P_2)`$ and $`(P_1,P_2)`$, respectively, with $`P_{1,2}>0`$. To make contact to the experiment we take $`P_1=0.8`$ and $`P_2=0.6`$, and impose a cut in the forward and backward directions. Specifically, we consider the cut angular range $`|\mathrm{cos}\theta |<0.9`$ and divide it into nine equal-size bins of width $`\mathrm{\Delta }z=0.2`$ ($`z\mathrm{cos}\theta `$). We also introduce the experimental efficiency, $`ϵ`$, for detecting the final $`e^+e^{}`$ pair, and according to the LEP2 experience $`ϵ=0.9`$ is assumed. The reach on the CI couplings, and the corresponding constraints on their allowed values in the case of no effect observed, can be estimated by performing $`\chi ^2`$ analysis, assuming the data to be well described by the SM ($`ϵ_{\alpha \beta }=0`$) predictions, i.e., that no deviation is observed within the foreseen experimental uncertainty.

The procedure, and the criteria, to derive numerical constraints from the Møller process and muon pair-production process (4) are quite similar. One should notice only that in the case of Møller scattering one can find for the cross section results similar to Bhabha scattering, that can be obtained by crossing symmetry except for the overall normalization factor 1/2 related to identical particles . Also, from Eqs. (5)-(8) one can obtain the cross section for muon pair-production process accounting that it proceeds solely via $`s`$-channel exchange .

As for the systematic uncertainty, we take $`\delta _{\mathrm{int}}/_{\mathrm{int}}=0.5\%`$, $`\delta ϵ/ϵ=0.5\%`$ and, regarding the electron and positron degrees of polarization, $`\delta P_1/P_1=\delta P_2/P_2=0.5\%`$. As a criterion to constrain the values of the contact interaction parameters allowed by the non-observation of the corresponding deviations, we impose $`\chi ^2<\chi _{\mathrm{CL}}^2`$, where the actual value of $`\chi _{\mathrm{CL}}^2`$ specifies the desired ‘confidence’ level. We take the values $`\chi _{\mathrm{CL}}^2=`$7.82 and 9.49 for 95% C.L. for a three- (Bhabha and Møller processes) and a four-parameter ($`\mu ^+\mu ^{}`$ pair production) fit, respectively.

In Figs. 1-2 we show the derived limits on the electron contact interactions at a LC with longitudinally polarized beams and using a model-independent analysis that allows to simultaneously account for all independent couplings as non-vanishing free parameters. From these figures one can conclude that the two processes, (1) and (2), are complementary as far as the sensitivity to the individual couplings in a model-independent data analysis is concerned: the sensitivity of Bhabha scattering to $`\mathrm{\Lambda }_{\mathrm{LR}}`$ is dramatically higher, while Møller scattering is the most sensitive to $`\mathrm{\Lambda }_{\mathrm{LL}}`$ and $`\mathrm{\Lambda }_{\mathrm{RR}}`$.


text:
# Charge Transport in Non-Diluted Conjugated Polymers

## Abstract

We present a unified calculation method for variable range hopping transport with a varying charge concentration and a varying applied electrical field. We demonstrate that the major differences between the transport properties measured at high concentration and low concentration can be explained within this framework. In particular the difference between the measured mobility, and mobility activation energy in polymer field effect transistors and polymer light emitting diodes is explained. A theoretical method to extract the charge carrier density of state from the transport measurements of non-diluted materials is proposed.

Charge transport in conjugated polymers is often described in terms of thermal assisted tunneling (hopping) r1 ; r2  or polaronic transfer mechanism r3  in a disordered array of localized states r4 . As the density of states of the charge carrier has not been directly measured for most of the conjugated polymers, a Gaussian density of states is often assumed. The Gaussian disorder model (GDM) r5  and related models with energetic spatial correlations between the states r6  have been solved previously for non-interacting carriers, or explicitly for low charge density. Such calculations are typically successful in describing the field and temperature dependence of the measured mobility in the steady state regime. The dynamic transport properties at the short time scales, where the charge carriers are far from thermodynamic equilibrium, have been described using different methods r7 ; r8 . However, all of the mentioned calculations are based on the assumption that the time average occupation of the sites is small and hence, the effect of the charge concentration is negligible.

Two types of experiments demonstrated that the transport depends on the charge density: the super linear increase of current with increasing molecular doping, and the trans-conductance measurements in polymeric field effect transistors (FETs) r9 . Recently, we demonstrated a mobility enhancement in poly-\[2-methoxy-5-(2’-ethyl-hexiloxy)-p-phenylenevinylene\] (MEH-PPV) FET by a factor of $``$10 due to charge concentration build up r10  reaching values of up to $`310^5`$ cm2/V sec. In light emitting diodes (LEDs) made from the same polymer the mobility depends strongly on the applied field, with a zero field mobility value of $`110^6`$ cm2/V sec r11  (see also r12 ). Typically, FETs operate at much lower electrical fields compared to LEDs and the charge residence time in the device is much longer. Therefore, the increase of the mobility in FETs cannot be explained in terms of electric field dependent mobility or by the dynamical transport behavior at short time scales.

Recently, the mobility r13  and the diffusion coefficient r14  dependence on charge concentration has been calculated using a low field approximation of the GDM, successfully describing the molecular doping effect. Using the energy correlation model with traps, the mobility dependence on electrical field at different charge concentration was calculated r15 , resulting in a similar field dependence of the mobility for different carrier densities. In this letter we suggest a unified calculation of the concentration and electrical field dependence, using a semi-classical mean medium approximation (MMA) r16  and demonstrate that the fundamental behavior of the different devices is recovered. Moreover, we propose a technique to extract the charge carrier density of state.

A specific spatial and energetic configuration of sites is described by its density of states (DOS) and by its spatial distribution function. Here we assume a homogeneous spatial distribution and a Gaussian DOS r17 ,

$$g\left(\epsilon \right)=\frac{1}{2\pi \sqrt{\sigma }}\mathrm{exp}\left(\frac{\left(\epsilon \epsilon _0\right)^2}{\sqrt{2\sigma ^2}}\right)$$
(1)
where $`\sigma `$ is defined as the width of the DOS, $`\epsilon `$ is the energy, and $`\epsilon _0`$ is the center of the DOS. We use the Gaussian functional form to enable comparison with other methods. However, we do not expect a change in the principle results discussed here, while replacing the Gaussian DOS with a more accurate (experimentally determined) DOS. The second assumption we use is that the charge transfer process drives the charge carrier population towards equilibrium energy-distribution. To clearly define an equilibrium energy-distribution we assume that each state cannot contain more than one charge carrier (due to columbic repulsion). Therefore the equilibrium distribution function is given by the following Fermi-Dirac form:

$$f(\epsilon ,\eta )=1\mathrm{/}\left[1+\mathrm{exp}\left(\beta \left(\epsilon _i\eta \right)\right)\right]$$
(2)
where $`\eta `$ is the chemical potential. A transfer rate which will drive the charge carriers toward a detailed equilibrium described in Eq. (2), has the general form of:

$$\upsilon _{ij}=\upsilon (\left|\epsilon _j\epsilon _i\right|,𝐑_{\mathrm{𝐢𝐣}})\mathrm{exp}\left[\frac{\beta }{2}\left[\left(\epsilon _j\epsilon _i\right)+\left|\epsilon _j\epsilon _i\right|\right]\right]$$
(3)
where $`\epsilon _i,\epsilon _j,𝐑_{\mathrm{𝐢𝐣}}`$ are the initial energy, final energy, and the initial-final states vector, respectively; $`\upsilon `$ is an envelope function that depends only on the absolute energy difference and the spatial coordinates, and $`\beta `$ is the inverse temperature ($`1/kT`$). In the absence of an additional energy dissipation process, and when the form of the transfer rate deviates from Eq. (3) (as in the adiabatic polaronic rate), the charge carriers will never reach thermal equilibrium (Eq. (2)). On the other hand, when the envelope function is decaying exponentially in space, with an inverse localization radii $`\gamma `$ ($`\upsilon =\upsilon _0\mathrm{exp}\left(\gamma \left|𝐑_{\mathrm{𝐢𝐣}}\right|\right)`$), the transfer rate becomes the known Miller-Abrahams rate r2 , and the charge carriers reach equilibrium after a sufficient time has lapsed.

In order to calculate the transport properties one needs to solve the Master Equation determined by the transfer rate and the energetic and spatial configuration of states, and to average over all possible configurations. Instead, we propose to use the mean medium approximation (MMA) where averaging on all possible discrete configurations is replaced by solving the transport properties of a homogeneous and continuous media, characterized by the energetic DOS function. Such an approximation is valid only for large enough space and time scales, and where the excess energy gained by the accelerating field is dissipated fast. Under the conditions for which the MMA is valid, one can deduce: 1) the charge carrier population is near detailed equilibrium, and the occupation probability of each site is determined by the local quasi-chemical potential and the equilibrium distribution function. 2) Since the medium is homogeneous all points in space are characterized by the same charge concentration and transfer rates and hence, the relaxation of a charge excess must follow the Poissonian Debye pattern. Namely, the transport has to be Gaussian (Markovian) and can be characterized by the first two spatial moments of the transport Green function which are related to the mobility ($`\mu `$) and the diffusion coefficient ($`D`$). Based on detailed equilibrium and given a DOS and a charge concentration, the Einstein relation ($`D/\mu `$) can be uniquely derived. As the system is Markovian, one needs only to derive a value for the mobility as a function of charge density and electric field to complete the transport description. We perform this calculation below using the MMA framework.

Before embarking on the mobility calculation we emphasize the importance of charge concentration effects in real devices. We calculated the charge carrier distribution in energy, $`p\left(\epsilon \right)`$, at equilibrium (Fig. 1) for several chemical potentials ($`\eta `$). As long as the chemical potential is below $`\epsilon _0\sigma \left(\beta \sigma +2\right)`$ the charge concentration has a Gaussian distribution centered at $`\epsilon _0\beta \sigma ^2`$, and the Boltzmann approximation is valid. For chemical potentials above this value the charge carrier distribution deviates from the Gaussian shape and the average energy increases. At this point the Boltzmann approximation is no longer valid, and one can expect the transport properties of the charge carriers to change. The relevant charge concentration can be then calculated for a given DOS. The total DOS can be estimated as the maximum concentration of electronic units (sites) or one over the molecular volume of such unit (5-10 monomers or 300 atoms volume equivalent to a total DOS of $``$1020 \[cm-3\]) r19 . For a DOS width of 5$`kT`$ (130 meV at room temperature), and a chemical potential of $`\epsilon _0\sigma \left(\beta \sigma +2\right)`$, the corresponding charge concentration is approximately 1011 \[cm-3\]. As polymer based electronic devices (as LEDs and FETs) operate at charge concentrations beyond 1015 \[cm-3\], the charge concentration effects should be taken into account. It is evident that for any practical concentration the charge energy distribution in a Gaussian DOS does not follow the Boltzmann distribution function. Previously we have shown that the $`D/\mu `$ ratio diverges from the classical $`kT/q`$ value r20 . In the following we will consider the effect of the charge concentration on the mobility.

The total current is calculated by integrating the current between each two sites while applying the MMA assumption that each point contains all possible energy states with a relative weight determined by the DOS function ($`g\left(\epsilon \right)`$):

$$\begin{array}{c}𝐉=\underset{\mathrm{}}{}𝑑𝐑_{\mathrm{𝐢𝐣}}\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\epsilon _i\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\epsilon _j\upsilon _{ij}(𝐑_{\mathrm{𝐢𝐣}},\epsilon _i,\epsilon _j)g\left(\epsilon _i\right)f(\epsilon _i,\eta )g\left(\epsilon _j𝐑_{\mathrm{𝐢𝐣}}𝐄\right)\hfill \\ \left[1f(\epsilon _j𝐑_{\mathrm{𝐢𝐣}}𝐄,\eta )\right]𝐑_{\mathrm{𝐢𝐣}}\widehat{𝐄}\hfill \end{array}$$
(4)
where $`𝐑_{\mathrm{𝐢𝐣}}\widehat{𝐄}`$ is the potential drop between sites induced by the electric field E. Inserting the Miller Abrahams hopping rate and selecting the $`(\left|𝐑_{\mathrm{𝐢𝐣}}\right|,z,\phi )`$ coordinate system, where $`z=𝐑_{\mathrm{𝐢𝐣}}\widehat{𝐄}`$ and $`\phi `$ is the angle in the plane perpendicular to z, the current integral becomes dependent on one spatial dimension, and reads:

$$\begin{array}{c}𝐉=2\pi \underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑z\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\epsilon _i\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\epsilon _jg\left(\epsilon _i\right)f(\epsilon _i,\eta )g\left(\epsilon _j\left(zE\right)\right)\left[1f(\epsilon _j\left(zE\right),\eta )\right]\hfill \\ \mathrm{exp}\left(\frac{\beta }{2}\left(\left(\epsilon _j\epsilon _i\right)+\left|\epsilon _j\epsilon _i\right|\right)\right)z\mathrm{exp}\left(\gamma \left|z\right|\right)\left(\gamma \left|z\right|+1\right)\hfill \end{array}$$
(5)
The mobility is calculated by the definition $`\mu 𝐉\mathrm{/}p𝐄`$, where p is the total charge carrier concentration: $`p\left(\eta \right)=g\left(\epsilon \right)f(\epsilon ,\eta )𝑑\epsilon `$ . The charge density dependence of the mobility for a range of electrical fields is shown in Fig. 2. We note that the effect of the charge density becomes pronounced for $`\eta >\epsilon _0\sigma \left(\beta \sigma +2\right)`$ where the system is degenerate, namely non-diluted (see also Fig. 1). At very high charge concentration the mobility value has a maximum due to ”over filling” of the states to the point where there are less states left to transfer into. The figure also shows that at higher electric fields the density dependence is less pronounced and the effect of the DOS form diminishes. At low electric fields, the curve describing the mobility is determined by the shape of the DOS. This can be shown by linearizing Eq. (5) in the limit $`\beta |zE|<<1`$:

$$\mu _{E0}=\frac{32\pi \beta \upsilon _0}{q\gamma ^5p(\eta ,g\left(\epsilon \right))}\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\epsilon _i\underset{\epsilon _i}{\overset{\mathrm{}}{}}𝑑\epsilon _jg\left(\epsilon _i\right)f(\epsilon _i,\eta )g\left(\epsilon _j\right)\left[1f(\epsilon _j,\eta )\right]\mathrm{exp}\left(\beta \left(\epsilon _j\epsilon _i\right)\right)$$
(6)
The low field mobility, given in Eq. (6), depends only on the exact shape of the DOS. Therefore, the DOS shape can be extracted from measurements of the low field mobility dependence on charge concentration, using this equation.

In Fig. 3 the mobility dependence on the electrical field is illustrated, for several charge concentrations. It shows that at the low concentration region all the mobility curves merge to one curve which depends strongly on the field, as expected from the Boltzmann approximation. However, as the charge concentration increases (to a practical value) and the chemical potential crosses the value of $`\epsilon _0\beta \sigma \left(\sigma +2\right)`$, the low field mobility increases, and the mobility field dependence weakens. At high electrical fields the curves merge again, as the potential drop due to the applied field shifts most of the final sites below that of the initial site (”state saturation”). In any case we expect that the theory used here would not be strictly valid at the extreme cases of high charge concentration and/or high electric fields. For example, it was suggested that at the high field regime the polaron dissociates into a free electron and a structural conformation which is left behind r22 , hence a different hopping mechanism should be applied. This can also be viewed as a violation of our assumption that any excess energy gained by the field is dissipated fast enough.

The inconsistency between the measured mobility in LEDs and FETs that was pointed out at the beginning of this paper can now be explained using the MMA calculation results. We propose that it steams from the vastly different operating conditions of each device expressed as the typical charge concentration and electric field. The ovals in Fig. 3 denote the typical operation field and concentration ranges in LEDs and FETs. The low concentration, high field condition in LEDs leads to a relatively low, strongly field dependant mobility. On the other hand, in FETs the mobility is typically higher and its dependence on the field is weak, as the gate bias induces a high charge concentration (as was demonstrated in r10 ). Similarly, another difference between FETs and LEDs is observed while examining the influence of charge carrier concentration on the activation energy of the mobility (Fig. 4). Since the transfer rate of the model is in the Miller-Abrahams form the only possible origin of the activation energy is the energetic disorder (contrary to other possible contribution as the polaronic effect, e.g. Ref. r23 ). While the activation energy at high applied fields is low, at low applied fields it is strongly dependent on the charge carrier concentration (insert in Fig. 4). The difference of 200 meV between ”low concentration” (LED) mobility activation energy and ”high concentration” (FET) activation energy was measured by us (215$`\pm `$10 meV r24  and 400$`\pm `$40 meV r11  in MEH-PPV based FET and LED, respectively) and in Ref. r12 . However, the activation energy measured at high concentrations (FET) is significantly higher than the GDM, Miller-Abrahams model prediction, indicating an additional source for the activation energy besides the energetic disorder. Explicitly, the activation energy in the polaronic model is constituted from the energetic disorder contribution and half of the polaronic energy. Therefore, at high concentration measurement, where the energetic disorder contribution is negligible, the polaronic binding energy can be estimated as twice the measured activation energy ($``$400 meV for MEH-PPV at the previous example).

In conclusion, we present a unified calculation method for variable range hopping of diluted and non-diluted materials (low and high charge carrier concentration), and a varying applied electric field. We demonstrate that the major differences between the transport properties measured at high concentration (in FETs) and low concentration (in LEDs) can be explained by this model, in particular the difference in the mobility and the activation energy of the mobility between the two. Moreover, the low value of the predicted high concentration activation energy, where compared to measurements, enables us to estimate the polaron binding energy in MEH-PPV as 400$`\pm `$40 meV. A theoretical method to extract the charge carrier DOS from the charge transport measurements is proposed, as well. Finally, we note that charge concentration effects play a role also in other experimental situations such as time of flight experiments r25 .

This research (No. 56/00-11.6) was supported by the Israel science foundation. Y.R. thanks Israel science Ministry and Israel science foundation for the generous scholarships.