text:
# Neutron production by cosmic-ray muons at shallow depth

## I Introduction

In the present work we report the results of a new measurement of the muon-induced neutron production at shallow depth, based on the Palo Verde neutrino detector . Even though the device was designed for a different purpose, namely the detection of reactor $`\overline{\nu }_e`$, it was possible to operate it in a mode suitable for the identification of muon-induced neutrons.

Neutrons and other hadrons produced by cosmic-ray muons in the earth are an important and unavoidable source of background for underground low rate experiments. Knowledge of the rates for these processes is an essential part of understanding backgrounds in neutrino detectors and other low counting rate experiments. For example, since muon-produced fast single or multiple neutrons can mimic the correlated signature of inverse neutron beta decay, searches for neutrino oscillations at nuclear reactors  must cope with this source of background. Other neutrino and proton decay experiments, as well as dark matter searches (even though often at greater depth), have to cope with this source of background as well. The CDMS experiment, for instance, is searching for cold dark matter , and is presently at shallow depth; muon-induced neutrons represent a major source of background. Low-energy accelerator neutrino oscillation searches  are usually performed near the Earth’s surface where muon-induced neutrons are a significant source of background. Another example is the proposed ORLaND neutrino detector at the Spallation Neutron Source at Oak Ridge National Laboratory  which is also planned to be at shallow depth.

Despite the importance of the subject, relatively little recent progress has been reported. Several measurements of the neutron production rates at various depths conducted in the past suggest a smooth dependence of the neutron yield on depth or, equivalently, on the average muon energy . However, a more recent measurement by the LVD collaboration , resulting in a smaller yield, is in disagreement with this simple dependence on depth. At shallow depth, where the average muon energy is $``$ 10 GeV, the smooth trend was confirmed by measurements with relatively small detectors without much shielding against the neutrons produced outside them . In the more recent experiment , the single and double neutron yields were determined separately, and pion production by muons was also observed. One drawback of that experiment, possibly present in the other ones as well, is the difficulty of distinguishing between the neutrons produced in the detector (the intended source of neutrons) and the neutrons produced outside by the muon-induced showers. In fact, in Ref.  it has been estimated that half of the detected single neutrons originated in the hadron showers outside the detector volume.

The theoretical description of this background process is usually based on the assumption that the electromagnetic interaction of high-energy muons with matter can be modeled by replacing the exchanged virtual photon by “equivalent” real photons, and using known photo-nuclear reaction cross sections . The analysis is complicated, since in order to relate the theoretical neutron production yield to measurement, the propagation and possible cascade multiplication of all reaction products must be understood. (For example, a $`\pi ^{}`$ produced by a muon will make more neutrons, as does the $`(n,2n)`$ reaction, etc.). While the smooth variation noted above for the neutron reaction rate versus depth is supported in some calculations , other approaches come to different conclusions. In particular, Ref.  was devoted, unlike the others, to the relatively shallow depth relevant for the present work. In that work, the calculated neutron yield was smaller than the measured one , while the calculated and measured $`\pi ^+`$ yields agreed with each other. One should also keep in mind that while the theoretical description quoted above deals with the muon-nucleus interaction involving the exchange of a virtual photon, neutrons can be also produced by nuclear interactions involving real bremsstrahlung photons, and electron/positron pairs created during the passage of muons through matter.(The present experiment cannot separate neutrons created by interactions involving virtual photons from those produced by interactions of real photons.)

## II Experiment

The measurement of the neutron production rate was performed at $`32\pm 3`$ meter-water-equivalent (mwe) depth using the Palo Verde neutrino detector. The detector and its operation is described in detail in Ref. . Briefly, the apparatus consists of 66 acrylic cells, each 900 x 12.7 x 25.4 cm3. These cells are filled with liquid scintillator loaded with Gd 0.1% by weight, for a total scintillator volume of 11.34 tons. The acrylic material has an aggregate mass of 3.48 tons. Muons can thus spall either in the scintillator volume or in the acrylic material, and the resulting neutrons cannot be distinguished by source. Neutrons, when moderated to thermal energies, are preferentially captured on Gd, resulting in a $`8`$ MeV gamma-ray cascade, the characteristic neutron capture signal. The central active detector volume is surrounded by a 1 m thick water shield, and the outermost layer of the detector is an active muon veto counter providing 4$`\pi `$ coverage.

Since the apparatus was intended for the detection of reactor neutrinos, it was not optimized for the neutron yield measurement. For the present purpose special runs were performed without the trigger rejection of the cosmic rays. The throughgoing muon sample was selected off-line such that at least two veto hits were recorded. Only muons which at the same time went through at least three cells of the central detector were included. The delayed neutron capture events, recognized by their energy deposit, are of two kinds: The “single bank” events, where only one neutron capture candidate event occurs following the veto hit; and the “two bank” events, which have two neutron capture candidate events following the muon. (The detector electronics is not capable of recording more than two correlated events.)

The measured quantities are the numbers of single ($`N_1`$) and double neutron-capture events ($`N_2`$) associated with $`N_\mu `$ muons traversing the central detector. The average path length of these muons in the central detector is $`X`$ (measured in units of g/cm2). Since some of the neutrons could have been created in the water shield or in other external detector parts, correction factors $`Q_k<1`$ are applied to $`N_k,k=1,2`$. Finally, the neutron detection efficiencies $`ϵ_{k,l}`$ are introduced. Here, e.g, $`ϵ_{1,1}`$ is the probability that one neutron was created and one detected, while $`ϵ_{1,2}`$ is the probability that two neutrons were created and only one detected. If all these quantities were known one could define neutron yields $`Y_l`$, etc. (per muon and g/cm2), where $`l`$ is the number of neutrons, that are independent of the detector properties and obey the relation

$$N_k=\frac{N_\mu Xϵ_k^{daq}}{Q_k}\underset{l=1}{\overset{3}{}}Y_lϵ_{k,l}.$$
(1)
Here $`ϵ_k^{daq}`$ is the detection livetime correction which in this case depends on the number of detected neutron-like events $`k`$.

Obviously, since only two quantities, $`N_1`$ and $`N_2`$, are measured, only two yields can be determined, and the system above cannot be solved without approximations. It is assumed further that the contributions from four or more produced neutrons can be neglected, as indicated in equation (1). Moreover, as will be shown below, the $`total`$ neutron yield

$$Y_{tot}=Y_1+2Y_2+3Y_3$$
(2)
is essentially independent of the ratio $`Y_3/Y_2`$, i.e., on the assumed value of $`Y_3`$, while the deduced single and double neutron yields $`Y_1`$ and $`Y_2`$ depend on that ratio significantly. Thus, the final results of the present experiment will be expressed as the measurement of $`Y_{tot}`$.

Note that the quantities $`Y_k`$, and thus also $`Y_{tot}`$, contain $`all`$ processes that lead to the production of $`k`$ neutrons by the muon. Thus, if the muon creates a $`\pi ^{}`$ (or any other particle except a neutron) which then, in turn, creates $`k`$ neutrons, all of the neutrons, regardless of source, contribute to $`Y_k`$. This somewhat awkward definition is necessary since the evaluation of the efficiencies $`ϵ_{k,l}`$ is based on a code that tracks just neutron propagation in the detector.

The neutron capture events were recorded in two runs, each about half of a day long. The raw muon veto rate was about 2 kHz and the rate of muons which went through at least three cells of the central detector and caused two or more detectable veto hits in the first (second) run was 270(275) Hz. These two runs are essentially equivalent, and the final total number of muons was determined to be $`N_\mu =1.42\times 10^7`$, very similar in both runs.

The average path length $`X`$ of the muons is estimated with a simple ray-tracing Monte Carlo simulation, starting from a $`\mathrm{cos}^2\theta `$ zenith angle distribution. The resulting $`X`$ was 125 g/cm2 for the central detector (i.e., the scintillator, the acrylic cells and the small amount of other materials (Cu, Fe) in the central detector) and 317 and 62 g/cm2 for the water shield and veto counter, respectively.

The neutron capture events $`N_k`$ were selected using cuts similar to the neutron part of the neutrino signal . For the single neutron events a time cut of 10 - 100 $`\mu `$s between the last muon event and the neutron capture event was used, while for the two neutron events only the lower limit of 10 $`\mu `$s was kept, the upper one being irrelevant. The time between the two neutron capture events was restricted to 5 - 100 $`\mu `$s. The energy distributions of all neutron events had identical shapes within statistics.

There are several contributions to the livetime correction factors $`ϵ_k^{daq}`$. The largest one, which affects primarily the single neutron events, arises since the initial information is overwritten when another muon strikes the veto after the single neutron event was recorded and before the 450 $`\mu `$s time interval of accepting new events. Combining the correction for this effect with other deadtime corrections, the resulting $`ϵ_1^{daq}`$ = 30.8(27.8)% for the one-neutron events in the first(second) runs, and $`ϵ_2^{daq}`$ = 72.4(70.3)% for the two-neutron events. The difference in $`ϵ_k^{daq}`$ between the two runs is caused primarily by the $`8\%`$ change in the raw muon rate.

The distribution of the time interval between the muon and the single-neutron event is shown in the top panel of Fig. 1, while the lower panel shows the same distribution for the first of the two neutron events. The fitted curves are of the form

$$F(t)=(a_1e^{t/\tau _1}+a_2e^{t/\tau _2}+a_3)e^{t/500\mu \mathrm{s}}.$$
(3)
Here $`\tau _1=28.8\pm 1.0\mu `$s for the single neutron case while $`\tau _1=11.8\pm 0.4\mu `$s (not very far from half of the previous $`\tau _1`$, the value one expects for two neutrons) for the double neutron case, are the characteristic neutron capture times. Both agree quite well with the Monte Carlo (MC) simulation, as seen in Fig. 1. The presence of the time constant $`\tau _2`$ is related to the inhomogeneous nature of the detector. Some neutrons enter the acrylic, where there is no Gd available for captures, thus prolonging their capture time. Finally, the last term, with the fixed time constant of 500 $`\mu `$s, represents the accidental background (500 $`\mu `$s is the average time interval between successive muons). Using such a fit one can subtract the accidental background in the relevant time window, and calculate the corresponding correction for the described choice of the time cuts.

In Fig. 2 the distribution of the elapsed time between the two neutron capture events is shown, and again fitted to the same functional dependence, Eq. (3). (The distinction between the top and bottom panels is explained below.) Here, the time constant $`\tau _1=18.1\pm 1.0\mu `$s in the top panel is substantially smaller than the $`\tau _1`$ obtained by the MC simulation ($`\tau _1^{MC}29\mu `$s, as expected if only two neutrons are involved). This feature suggests that the data sample contains a nonnegligible component with more than two neutrons. We will return to this point later.

## III Efficiencies

While the evaluation of $`N_\mu `$, $`N_1,N_2`$, and $`X`$ is rather straightforward, the determination of the efficiency matrix $`ϵ_{k,l}`$ is somewhat model dependent. This is because, in order to calculate $`ϵ_{k,l}`$, one needs the initial energy and angular distributions of the created neutrons, both of which are poorly known. In practice, a number of physically plausible assumptions are made about these distributions, and the spread among the resulting efficiencies yields a measure of the systematic uncertainty. Once the neutron initial distribution is chosen, the Monte Carlo simulation of the neutron transport developed for the neutrino experiment  is used. The detector geometry, materials and electromagnetic interactions are simulated using GEANT . Hadronic interactions are simulated by FLUKA , and the low energy neutron transport by GCALOR .

To evaluate the quantities $`ϵ_{1,1}`$ and $`ϵ_{2,1}`$ initial single neutrons were distributed randomly in position and initial angles through the detector volume. For the initial energy distribution, several possibilities were used: the exponential distribution $`\mathrm{exp}^{E/39\mathrm{M}\mathrm{e}\mathrm{V}}`$ as proposed by the Karmen collaboration , and functions with the power dependence $`E^x,0.5x2`$ (see ). These cover the shape following $`E^{1.86}`$ suggested by the experience with photo-nuclear processes . All of these distributions result in similar efficiencies. Using the spread of the calculated values as a measure of the systematic error, and taking a simple average as the most probable value, one obtains

$$ϵ_{1,1}=0.20\pm 0.03,ϵ_{2,1}=0.007\pm 0.004.$$
(4)
Note that the probability that a single initial neutron will result in two neutron captures, characterized by $`ϵ_{2,1}`$, is quite small.

The determination of two neutron efficiencies is even more difficult, since it depends on the energy distribution of both initial neutrons. To evaluate the quantities $`ϵ_{1,2}`$ and $`ϵ_{2,2}`$ two neutrons were created at the same random point and with random directions for each of them. Two extremes were considered. In the first of them, the total energy of the neutrons followed one of the previously described functions used above in the single neutron case. This initial energy was distributed randomly between them. (Essentially the same efficiency was obtained if the two neutrons shared the initial available energy equally.) The other extreme is obtained if both neutrons each have the energy distribution used in the single neutron case, i.e., the total neutron energy is on average twice as large. The spread (not very large) was again used in the averaging and in the estimate of the systematic error, resulting in

$$ϵ_{1,2}=0.22\pm 0.01,ϵ_{2,2}=0.06\pm 0.01.$$
(5)
Note that there is a sizable probability that only one neutron is detected when two were initially created. None of the previous analyses, in particular Ref. , took that into account. Next, one has to consider the possibility that three neutrons were spalled by the muon, but only one or two were detected. Again, using the average and the two extreme possibilities to divide the available energy among the three neutrons, one obtains

$$ϵ_{1,3}=0.19\pm 0.01,ϵ_{2,3}=0.10\pm 0.01.$$
(6)
It is important to realize that these efficiencies are not much smaller than those for initial one or two neutrons. Thus, the effect of $`Y_3`$ should be considered. (The effect of four and more neutron spallation will be neglected, however.)

The correction factors $`Q_k`$ that exclude neutrons created outside the central detector volume (i.e. in the water shield since the effect of the veto is negligible), but captured there, must be determined also. To do that, the efficiencies for each passive volume were determined using the same Monte Carlo code as for the case of the central detector. Then $`Q_k`$ were determined from

$$Q_k=\frac{(X_lY_lϵ_{k,l})_{centraldet.}}{(X_lY_lϵ_{k,l})_{centraldet.}+(X_lY_lϵ_{k,l})_{water}},$$
(7)
resulting in $`Q_1=0.80\pm 0.10`$ and $`Q_2=0.94\pm 0.07`$. In order to obtain the above values a crude assumption $`Y_1=2Y_2=2Y_3`$, was used. It is important to note that the uncertainties of $`Q_k`$ are strongly correlated with the error in the efficiencies.

The effect of the neutron component of hadronic cascades created outside the veto, and coincident with the muon which created them, remains to be determined. The ratio of muons making three hits in the veto, thus involving more than one particle, to the prevalent case of the two hits was employed as a measure of the frequency of showers. This ratio was $`6.7\pm 0.3\%`$ and $`7.3\pm 0.3\%`$ in the first and second runs, respectively, i.e., about 7% of muons entering the veto were accompanied by a shower.

To see that the ‘triple-veto-hit’ events are really different from the standard throughgoing muons with just two veto hits, one can form ratios $`R_k`$

$$R_k=\frac{N_k^{shower}/N_k^{total}}{N_\mu ^{shower}/N_\mu ^{total}},$$
(8)
which represent a quantitative measure of the neutron content of the showers. The measured values are $`R_1=2.0\pm 0.1`$ and $`R_2=4.4\pm 0.2`$, both significantly larger than unity. Thus the shower events indeed are richer in neutrons, particularly in the two neutron sample.

Moreover, the double neutron events with three veto hits have distinctly different time structure than the more common ones with just two veto hits. This is shown in the bottom part of Fig. 2 for the interevent time, i.e. the elapsed time between the first and second neutron captures. The corresponding time constant $`\tau _1`$ for the three-veto-hit events is significantly smaller than for the two-veto-hit events, showing that the shower events have a large multineutron component. A similar effect is present when the capture time of the first neutron is considered. It is therefore likely that the discrepancies between data and simulation in the time dependence of the interevent time interval, noted above, are the consequence of a multineutron ($`k>2`$) component in the two bank events.

Since the shower events contain an unknown number of neutrons created outside the detector volume, they are excluded from further consideration.

## IV Results

In the first run, the observed numbers of neutron captures, corrected for the random background and with the effect of showers subtracted, were $`N_1=3916\pm 66`$, $`N_2=828\pm 29`$, while in the second run $`N_1=3451\pm 62`$, $`N_2=829\pm 29`$. The two runs, which were separated in time by 10 months, give consistent results when the differences in $`ϵ^{daq}`$ is taken into account, proving that the experiment is stable.

As a first step, the yields are analysed as in Ref. , i.e., only the ‘diagonal’ efficiencies $`ϵ_{k,k}`$ are taken into account, thus

$$Y_k^{simple}=N_kQ_k/(N_\mu Xϵ_{k,k}ϵ_k^{daq}).$$
(9)
The resulting yields obtained by averaging the two runs are:

$`Y_1^{simple}`$ $`=`$ $`(2.94\pm 0.04(stat)\pm 0.50(syst))\times 10^5{\displaystyle \frac{\mathrm{neutrons}}{\mu \mathrm{g}\mathrm{cm}^2}}`$ (10)
$`Y_2^{simple}`$ $`=`$ $`(0.98\pm 0.03(stat)\pm 0.14(syst))\times 10^5{\displaystyle \frac{\mathrm{neutrons}}{\mu \mathrm{g}\mathrm{cm}^2}}.`$ (11)
The exclusion of external showers, characterized by the three-veto-hit events, has resulted in the reduction of $`Y_1^{simple}`$ by 14%, and of $`Y_2^{simple}`$ by 45%. Thus, presumably due to the presence of the sizable shield $`X380\mathrm{g}/\mathrm{cm}^2`$, the effect of the external showers, while still clearly present, was reduced compared to the findings of Ref. , where (in the same units) $`Y_1^{simple}=4.3(2.0)\times 10^5`$ and $`Y_2^{simple}=1.6(0.5)\times 10^5`$ without (with) the correction for external showers.

However, the proper analysis should include the full efficiency matrix, i.e., the possibility that two neutrons were initially produced but only one neutron capture was recorded and vice versa, as well as the possibility that three neutrons were originally produced and only two or one neutron captures were observed. As pointed out earlier, with only two measured capture rates, $`N_1`$ and $`N_2`$, it is impossible to deduce all the relevant information without further constraints. To avoid these difficulties, and to make the comparison with other experiments easier, the $`total`$ number of neutrons produced per muon was evaluated. This quantity, $`Y_{tot}`$, defined earlier in Eq. (2), has the further advantage that it is essentially independent of the ratio $`Y_3/Y_2`$ for three to two neutron production, as demonstrated in Table I.

Thus, the final result of the present measurement can be expressed as

$$Y_{tot}=(3.60\pm 0.09\pm 0.31)\times 10^5\frac{\mathrm{neutrons}}{\mu \mathrm{g}\mathrm{cm}^2},$$
(12)
where the systematic error is an estimate based on the spread of values in Table I added in quadrature to the spread of the evaluated efficiencies.

All available data on neutron yields are collected and compared in Fig. 3. (It is assumed that the other measurements () also are really $`Y_{tot}`$ measurements.) For the measurements of Ref.  both results, with and without the correction for external showers, are shown. In some of the other measurement the shower contribution was excluded, but it is not clear how well this has been done, since at least the results at 25 and 316 mwe were obtained with relatively small and unshielded detectors.

At shallow depth there are now three measurements, with essentially consistent results. Clearly, still better and more complete measurements are desirable, in which the full neutron multiplicity and energy spectra are determined and the neutrons produced externally are reliably identified.

### A Acknowledgment

We would like to thank the Arizona Public Service Company for the generous hospitality provided at the Palo Verde plant. The important contributions of M. Chen, R. Hertenberger, K. Lou, and N. Mascarenhas in the early stages of this project are gratefully acknowledged. We also acknowledge the generous financial help from the University of Alabama, Arizona State University, California Institute of Technology, and Stanford University. Finally, our gratitude goes to CERN, DESY, FNAL, LANL, LLNL, SLAC, and TJNAF who at different times provided us with parts and equipment needed for the experiment. This project was supported in part by the USDoE.


format: text

text:
# Plasma Oscillations and Expansion of an Ultracold Neutral Plasma

## Abstract

We report the observation of plasma oscillations in an ultracold neutral plasma. With this collective mode we probe the electron density distribution and study the expansion of the plasma as a function of time. For classical plasma conditions, i.e. weak Coulomb coupling, the expansion is dominated by the pressure of the electron gas and is described by a hydrodynamic model. Discrepancies between the model and observations at low temperature and high density may be due to strong coupling of the electrons.

One of the most interesting features of neutral plasmas is the rich assortment of collective modes that they support. The most common of these is the plasma oscillation , in which electrons oscillate around their equilibrium positions and ions are essentially stationary. This mode is a valuable probe of ionized gases because the oscillation frequency depends solely on the electron density.

In an ultracold neutral plasma as reported in , the density is nonuniform and changing in time. A diagnostic of the density is thus necessary for a variety of experiments, such as determination of the three-body recombination rate at ultralow temperature , and observation of the effects of strong Coulomb coupling  in a two-component system. A density probe would also aid in the study of the evolution of a dense gas of cold Rydberg atoms to a plasma , which may be an analog of the Mott insulator-conductor phase transition .

In this work we excite plasma oscillations in an ultracold neutral plasma by applying a radio frequency (rf) electric field. The oscillations are used to map the plasma density distribution and reveal the particle dynamics and energy flow during the expansion of the ionized gas.

The creation of an ultracold plasma has been described in . A few million metastable xenon atoms are laser cooled to approximately $`10\mu `$K. The peak density is about $`2\times 10^{10}`$ cm-3 and the spatial distribution of the cloud is Gaussian with an rms radius $`\sigma 220\mu `$m. These parameters are determined with resonant laser absorption imaging . To produce the plasma, up to $`25`$% of the atoms are photoionized in a two-photon excitation. Light for this process is provided by a Ti:sapphire laser at $`882`$ nm and a pulsed dye laser at $`514`$ nm ($`10`$ns pulse length). Because of the small electron-ion mass ratio, the resulting electrons have an initial kinetic energy ($`E_e`$) approximately equal to the difference between the photon energy and the ionization potential. In this study we vary $`E_e/k_B`$ between 1 and 1000 K. The initial kinetic energy of the ions varies between 10 $`\mu `$K and 4 mK.

For detection of charged particles, a small DC field (about $`1`$ mV/cm) directs electrons to a single channel electron multiplier and ions to a multichannel plate detector. The amplitude of the rf field that excites plasma oscillations, $`F`$, varies between $`0.220`$mV/cm rms. All electric fields are applied to the plasma with grids located above and below the laser-atom interaction region.

In the absence of a magnetic field, the frequency of plasma oscillations is given by $`f_e=(1/2\pi )\sqrt{e^2n_e/ϵ_0m_e}`$ . Here, $`e`$ is the elementary charge, $`n_e`$ is the electron density, $`ϵ_0`$ is the permittivity of vacuum, and $`m_e`$ is the electron mass. This relation is most often derived for an infinite homogeneous plasma, but it is also valid in our inhomogeneous system for modes which are localized in regions of near resonant density. Corrections to $`f_e`$ due to finite temperature  depend on the wavelength of the collective oscillation, which is difficult to accurately estimate. Such corrections are not expected to be large and will be neglected. We observe plasma oscillations with frequencies from $`1`$ to $`250`$ MHz. This corresponds to resonant electron densities, $`n_r`$, between $`1\times 10^4`$ cm-3 and $`8\times 10^8`$ cm-3. The oscillation frequency is sensitive only to $`n_e`$, but, as explained in , the core of the plasma is neutral. This implies that plasma oscillations measure electron and ion densities in this region ($`n_e=n_in`$).

Figure 1a shows electron signals from an ultracold neutral plasma created by photoionization at time $`t=0`$. Some electrons leave the sample and arrive at the detector at about $`1\mu `$s, producing the first peak in the signal. The resulting excess positive charge in the plasma creates a Coulomb potential well that traps the remaining electrons . In the work reported here, typically $`9099`$% of the electrons are trapped. Debye shielding maintains local neutrality inside a radius $`r_e`$ beyond which the electron density drops to zero on a length scale equal to $`\lambda _D`$. The value of $`r_e`$ depends on the fraction of electrons that has escaped, and $`\lambda _D`$ is the Debye screening length, $`\lambda _D=\sqrt{ϵ_0k_BT_e/e^2n_e}`$, where $`T_e`$ is the electron temperature. For our conditions $`r_e`$$`\stackrel{>}{}`$$`\mathrm{\hspace{0.17em}2}\sigma `$, and $`\lambda _D\sigma `$. As the plasma expands, the depth of the Coulomb well decreases, allowing the remaining electrons to leave the trap. This produces the broad peak at $`25\mu `$s.

In the presence of an rf field an additional peak appears in the electron signal (Fig. 1a). We understand the generation of this peak as follows: The applied rf field excites plasma oscillations only where the frequency is resonant. Energy is thus pumped into the plasma in the shell with the appropriate electron density ($`n=n_r`$). The amplitude of the collective electron motion is much less than $`\sigma `$, but the acquired energy is collisionally redistributed among all the electrons within $`101000`$ ns , raising the electron temperature. This increases the evaporation rate of electrons out of the Coulomb well, which produces the plasma oscillation response on the electron signal.

The resonant response at a given time, $`S(t)`$, is proportional to the number of electrons in the region where the density equals $`n_r`$. If we make a simple local density approximation and neglect decoherence of the oscillations, $`S(t)F^2d^3rn(𝐫,t)\delta [n(𝐫,t)n_r]`$. The width in time of the observed signal (Fig. 1a) reflects the density distribution of the sample . At early times when the density is higher than $`n_r`$ almost everywhere, $`S(t)`$ is negligibly small. As the cloud expands and the density decreases, the response grows because the fraction of the plasma which is in resonance increases. The peak of the response appears approximately when the average density, $`\overline{n}`$, becomes resonant with the rf field. $`S(t)`$ vanishes when the peak density is less than $`n_r`$.

The resonant response arrives later for lower frequency (Fig. 1b) as expected because $`\overline{n}`$ decreases in time. Assuming that the plasma density profile remains Gaussian during the expansion, $`S(t)`$ can be evaluated and its amplitude scales as $`F^2/n_r`$. In Fig. 1b the data have been normalized by this factor and the resulting amplitudes are similar for all conditions.

By equating $`\overline{n}`$ to $`n_r`$ when the response peak arrives, we can plot the average plasma density as a function of time (Fig. 2). The data are well described by a self similar expansion of a Gaussian cloud, $`\overline{n}=N/[4\pi (\sigma _0^2+v_0^2t^2)]^{3/2}`$, where $`\sigma _0`$ is the initial rms radius and $`v_0`$ is the rms radial velocity at long times. $`N`$ is determined independently by counting the number of neutral atoms with and without photoionization. The extracted values of $`\sigma _0`$ are equal to the size of the initial atom cloud. In such an expansion, the average kinetic energy per particle is $`3mv_0^2/2`$.

Figure 3 shows the dependence of $`v_0`$ on density and initial electron energy. We first discuss data with $`E_e70`$K, for which the expansion velocities approximately follow $`v_0=\sqrt{E_e/\alpha m_i}`$, where $`m_i`$ is the ion mass and $`\alpha =1.7`$ is a fit parameter. For the plasma to expand at this rate, the ions must acquire, on average, a velocity characteristic of the electron energy. This is much greater than the initial ion thermal velocity. Electron-ion equipartition of energy would yield $`v_0=\sqrt{E_e/3m_i}`$, close to the observed value. However, due to the large electron-ion mass difference, this thermalization requires milliseconds . The observed expansion, in contrast, occurs on a time scale of tens of microseconds. One might expect the expansion to be dominated by the Coulomb energy arising from the slight charge imbalance of the plasma, but this energy is about an order of magnitude less than the observed expansion energy. Also, by Gauss’ law, it would only be important in the expansion of the non-neutral outer shell of the plasma. The oscillation probe provides information only on the neutral core because it relies on the presence of electrons.

A hydrodynamic model , which describes the plasma on length scales larger than $`\lambda _D`$, shows that the expansion is driven by the pressure of the electron gas. The pressure is exerted on the ions by outward-moving electrons that are stopped and accelerated inward in the trap. For the hydrodynamic calculation, ions and electrons are treated as fluids with local densities $`n_a(𝐫)`$ and average velocities $`𝐮_a(𝐫)=𝐯_a(𝐫)`$. Here, $`a`$ refers to either electrons or ions, and $`\mathrm{}`$ denotes a local ensemble average. Particle and momentum conservation lead to the momentum balance equations

$$m_an_a\left[\frac{𝐮_a}{t}+(𝐮_a)𝐮_a\right]=(n_ak_BT_a)+𝐑_{ab}.$$

Here $`n_ak_BT_a`$ represents a scalar pressure . The ion and electron equations are coupled by $`𝐑_{ab}`$, which is the rate of momentum exchange between species $`a`$ and $`b`$. The exact form of this term is unimportant for this study, but $`𝐑_{ab}=𝐑_{ba}`$. Plasma hydrodynamic equations typically have electric and magnetic field terms, but applied and internally generated fields are negligible when describing the expansion.

We can make a few simplifying approximations that are valid before the system has significantly expanded. The directed motion is negligible, so we set $`𝐮_a0`$ everywhere. Because $`n_en_i=n`$, $`𝐮_e/t𝐮_i/t`$. Due to the small electron mass, the rate of increase of average electron momentum is negligible compared to that of the ions. The electron momentum balance equation then yields $`(nk_BT_e)𝐑_{ei}`$, which describes a balance between the pressure of the electron gas and collisional interactions. This is the hydrodynamic depiction of the trapping of electrons by the ions.

In the ion momentum balance equation, we eliminate $`𝐑_{ie}`$ using the electron equation, and we drop the pressure term because the ion thermal motion is negligible. Thus $`m_in𝐮_i/t(nk_BT_e)`$, which shows that the pressure of the electron gas drives the expansion. This result implies that the ions acquire a velocity of order $`\sqrt{k_BT_e/m_i}`$, which is in qualitative agreement with the high $`E_e`$ data of Fig. 3. To calculate the expansion velocity more quantitatively, one must consider that as electrons move in the expanding trap, they perform work on the ions and cool adiabatically. The thermodynamics of this process  is beyond the scope of this study.

The data in Fig. 3 indicate that about $`90`$ % of the initial kinetic energy of the electrons is transferred to the ions’ kinetic energy, $`3m_iv_0^2/2=3E_e/2\alpha `$. This does not imply that the temperature of the ions becomes comparable to $`E_e/k_B`$ in this process. For the ions, $`𝐮_i`$ increases, but $`m_i|𝐯_i𝐮_i|^2`$, which measures random thermal motion and thus temperature, is expected to remain small. This follows from slow ion-electron thermalization  and correlation between position and velocity during the expansion .

We now turn our attention to systems with $`E_e<70`$K (Fig. 3). They expand faster than expected from an extrapolation of $`v_0=\sqrt{E_e/\alpha m_i}`$, and thus do not even qualitatively follow the hydrodynamic model. A relative measure of the deviation is $`(m_iv_0^2E_e/\alpha )/(E_e/\alpha )`$. Figure 4 shows that the relative deviation increases with increasing electron Coulomb coupling parameter , $`\mathrm{\Gamma }_e=(e^2/4\pi \epsilon _0a)/k_BT_e.`$ Here, $`a=(4\pi n/3)^{1/3}`$ is the Wigner-Seitz radius, $`n`$ is the peak density at $`t=0`$, and the temperature is calculated by $`3k_BT_e/2=E_e`$.

The fact that the relative deviation depends only on $`\mathrm{\Gamma }_e`$, and that it becomes significant as $`\mathrm{\Gamma }_e`$ approaches 1, suggests that we are observing the effects of strong coupling of the electrons . The hydrodynamic model of the plasma is only valid when $`\mathrm{\Gamma }_e1`$. When $`\mathrm{\Gamma }_e`$$`\stackrel{>}{}`$$`\mathrm{\hspace{0.17em}1}`$, electron and ion spatial distributions show short range correlated fluctuations that are not accounted for in a smooth fluid description . Correlations between the ion and electron positions would provide the excess kinetic energy observed in the expansion by lowering the potential energy of the plasma. This satisfies overall energy conservation and it may also explain the systematically poor fits of the data for high $`\mathrm{\Gamma }_e`$ (See Fig. 2).

Strong coupling is also predicted to alter the relation for the frequency $`f_e`$ , with which we extract the plasma density, size, and expansion velocity. The trend of this effect agrees qualitatively with the observed deviation, but knowledge of the wavelength of the collective oscillation is needed for a quantitative comparison.

Other possible explanations for the deviation are related to how the ultracold plasma is created. The $`10`$ns duration of the photoionization pulse is long compared to the time required for electrons to move an interparticle spacing. Photoionization late in the pulse thus occurs in the presence of free charges, which will depress the atomic ionization threshold by $`\mathrm{\Delta }E_{IP}\frac{1}{2}k_BT_e\{[(3\mathrm{\Gamma }_e)^{3/2}+1]^{2/3}1\}`$ . This effect might increase the electron kinetic energy by $`\mathrm{\Delta }E_{IP}`$ above what has been assumed. However, as shown in Fig. 4, the calculated $`\mathrm{\Delta }E_{IP}`$ is about an order of magnitude smaller than the observed effect. The random potential energy of charged particles when they are created may also yield a greater electron energy than $`E_e`$ .

High $`\mathrm{\Gamma }_e`$ (high density and low temperature) conditions are desirable for studying the three-body recombination rate in an ultracold plasma. The theory  for this process was developed for high temperature, and is expected to break down in the ultracold regime . Measuring or setting an upper limit for the recombination rate is not possible until the dynamics of high $`\mathrm{\Gamma }_e`$ systems is understood. We are currently studying this problem with molecular dynamics calculations.

We have shown that plasma oscillations are a valuable probe of the ion and electron density in an ultracold neutral plasma. This tool will facilitate future experimental studies of this novel system, such as the search for other collective modes in the plasma and further investigation of the effects of correlations due to strong coupling.

We thank Lee Collins for helpful discussions and Michael Lim for assistance with data analysis. S. Kulin acknowledges funding from the Alexander-von-Humboldt foundation. This work was funded by the ONR.


text:
# References

Reply to the comment of Chudnowski and Garanin on Phonon Assited Tunneling in Mn12 (PRL 83, 416 (1999))
G. Bellessa1, N. Vernier1, B. Barbara2 and D. Gatteschi3

1 Laboratoire de Physique des Solides, Bâtiment 510, Université Paris-Sud, 91405 Orsay, France
2 Laboratoire de Magnétisme Louis Néel, CNRS, BP 166, 38042 Grenoble, France
3 Dipartimento di Chimica, Università di Firenze, Via Maragliano 77, 50144 Firenze, Italia

Chudnovsky and Garanin have attacked harshly our paper  in a comment in arXiv:cond-mat . They write in their comment that “all formulas of Ref.  are incorrect”. They question “ how can $`\mathrm{\Delta }`$ be the magnetic dipole matrix element and the tunneling splitting at the same time?”. These two $`\mathrm{\Delta }`$ are obviously different and we write it explicitlely in the text:

* The first $`\mathrm{\Delta }`$ appears in the expression of the fundamental doublet:

$$\sqrt{\mathrm{\Delta }^2+ϵ^2}$$
(1)
$`\mathrm{\Delta }`$ is called the tunneling rate and $`ϵ`$ is a Zeeman term. We seize the opportunity of restablishing the historical truth of Eq. 1. It is Korenblit and Shender who the first have pointed out that high spins with large anisotropy give tunneling states, have resolved analytically the Hamiltonian in the high-spin approximation and have given the explicite form of $`\mathrm{\Delta }`$ in Eq. 1 . This is often ignored by the MQT community.
* In our paper there is a second $`\mathrm{\Delta }`$ which is obviously different from the first one. We write (Eq. 3 of our paper):

$$\chi ^{\prime \prime }=CN\mathrm{\Delta }^2T_2tanh(\frac{\mathrm{}\omega }{2k_BT})$$
(2)
and then we write:“where $`\mathrm{\Delta }`$ is the magnetic dipole matrix element betwen the two states of the fundamental doublet”. We take the same definition as Abragam and Bleaney in their formula that has been called “ trivial” by Chudnowsky and Garanin . Hence, this $`\mathrm{\Delta }`$ is explicitly different from the first one which is the splitting (equal to $`\mathrm{}\omega _0`$, where $`\omega _0`$ is the resonance frequency and depends on the applied magnetic field).

Another way to show that the two $`\mathrm{\Delta }`$ are different is to treat the fundamental doublet as an effective spin $`S=1/2`$ and to write the Hamiltonian:

$$=_0+^{}$$
(3)
where $`_0`$ is the unperturbed hamiltonian in the representation of the unperturbated eigenstates $`|1>`$ and $`|2>`$. The eigenstates are coupled to the electromagnetic field, which gives $``$’. In our case the magnetic field is perpendicular to anisotropy axis. So, the Zeeman term is zero and the splitting $`E=\sqrt{\mathrm{\Delta }^2+ϵ^2}`$ is $`\mathrm{\Delta }`$. Then the Hamiltonian is:

$$=1/2\left(\begin{array}{cc}\mathrm{\Delta }& 0\\ 0& \mathrm{\Delta }\end{array}\right)+1/2\left(\begin{array}{cc}\delta & M\\ M& \delta \end{array}\right)$$
(4)
where $`M`$ is the matrix element $`<|^{}|+>`$ and $`\delta `$ is negligeable as compared with $`\mathrm{\Delta }`$. The second $`\mathrm{\Delta }`$ in our paper is now $`M`$. It is this matrix element which induces absorption (and/or emission) of electromagnetic quanta and transitions within the fundamental doublet, these transitions being obviously tunneling transitions.

Eq. 2 gives a susceptibility $`\chi ^{\prime \prime }`$ constant at low temperature ($`\mathrm{}\omega <k_BT`$). However, we observe a strong increase of $`\chi ^{\prime \prime }`$ with increasing temperature. The only parameter in Eq. 2 is $`\mathrm{\Delta }`$ (the second one of our paper or M now), which means that this one increases with temperature, i.e. when phonons are added. This is an experimental fact.

To conclude, the attack of Chudnovski and Garanin upon the tunneling rate in $`Mn_{12}`$ is absolutely irrelevant.


text:
# Diffusive mass transfer by non equilibrium fluctuations: Fick’s law revisited

## Abstract

Recent experimental and theoretical works have shown that giant fluctuations are present during diffusion in liquid systems. We use linearized fluctuating hydrodynamics to calculate the net mass transfer due to these non equilibrium fluctuations. Surprisingly the mass flow turns out to coincide with the usual Fick’s one. The renormalization of the hydrodynamic equations allows us to quantify the gravitational modifications of the diffusion coefficient induced by the gravitational stabilization of long wavelength fluctuations.

Recent experiments have shown that giant non equilibrium concentration fluctuations are a universal feature of diffusion in binary liquid mixtures and macromolecular solutions vailati97\_2 ; brogioli00 . The presence of long-range correlations in a diffusing liquid mixture is a surprising result, as diffusion is usually thought to give rise to a homogeneous mixing at the molecular level, other length scales playing no role in the process. The results of these experiments have been interpreted by means of an extension of Landau’s Fluctuating Hydrodynamics Landau  to time-dependent diffusion processes vailati98 .

Linearized hydrodynamics shows that the fluctuations are generated by thermal velocity fluctuations, due to the presence of a macroscopic non equilibrium concentration gradient. The mechanism is analogous to that first predicted law89\_2 ; segre93\_3  and then observed segre93 ; li94 ; vailati96 ; li98  for stationary thermal diffusion in a binary liquid mixture. Velocity fluctuations parallel to the macroscopic gradient displace parcels of fluid into regions with different concentration, thus giving rise to concentration fluctuations weitz97 . Since the fluctuations displace mass along the concentration gradient, one might wonder if they contribute significantly to the macroscopic diffusive mass flow.

We will show that, quite unexpectedly, the net mass transfer due to non equilibrium fluctuations corresponds to the macroscopic Fick’s mass flow. Traditionally the diffusive flow has been interpreted as due to the random hops of the molecules in the fluid. Here we show that Fick’s flow is due to random hops of parcels of fluid whose dimensions span all the length scales from the molecular up to the macroscopic one. This is a surprising result, as usually thermal fluctuations are thought of as small perturbations of a macroscopic state, while here we show that the macroscopic flow is determined by a second-order effect in the fluctuations. The introduction of new degrees of freedom, the hydrodynamic fluctuations, determines an ultraviolet divergence of the diffusion coefficient, which we regularize by imposing a cut-off at large wave vectors. The regularization procedure involves the introduction of an arbitrary length scale. By imposing that the cut-off length scale is of the order of the correlation length of the mixture we find that the diffusion coefficient displays the usual Stokes-Einstein form and is of the same order of magnitude of the phenomenological one. However the introduction of an arbitrary length scale into the theory can be avoided by a renormalization of the hydrodynamic equations. The renormalization procedure allows us to quantify the modifications of the phenomenological diffusion coefficient due to the gravitational stabilization of long wavelength fluctuations.

We will now derive the net contribution of the non equilibrium fluctuations to the mass transfer. The sample considered is a binary mixture of concentration $`c(𝐱,t)`$ under rest. A local fluctuation in the current of mass can be written as the second order product

$$\delta 𝐣_f(𝐱,t)=\rho \delta c(𝐱,t)\delta 𝐮(𝐱,t),$$
(1)
where $`\delta c(𝐱,t)`$ is the local concentration fluctuation and $`\delta 𝐮(𝐱,t)`$ the local velocity fluctuation.

The net contribution of the fluctuations to the diffusive current is obtained by averaging over the configurations of the system:

$$𝐣_f(𝐱,t)=\rho \delta c(𝐱,t)\delta 𝐮(𝐱,t).$$
(2)
Equation (1) shows that a current fluctuation takes place whenever a concentration fluctuation and a velocity fluctuation occur simultaneously at the same place. At thermodynamic equilibrium fluctuations in the velocity and concentration are not correlated. Therefore the probability of a current fluctuation is small and with random sign, and the net flow $`𝐣_f`$ averages to zero. Under non equilibrium conditions, however, velocity fluctuations are correlated to concentration fluctuations vailati98 ; law89\_2 ; segre93\_3 . The presence of a macroscopic concentration gradient determines the existence of velocity-induced concentration fluctuations: a velocity fluctuation determines a small vortex inside the fluid and this brings parcels of fluid with a certain concentration into layers of fluid with a different concentration. On the contrary, the gravity force determines the opposite mechanism: whenever a concentration fluctuation is produced, the buoyancy force tends to displace the fluctuation. Therefore terms in the average of Eq. (2) are correlated and this determines a non vanishing mass flow $`𝐣_f`$ .

We will shortly derive an explicit expression for $`𝐣_f`$ by determining the cross-correlation properties of concentration and velocity fluctuations in the presence of a macroscopic concentration gradient. This task can be accomplished by using Landau’s linearized hydrodynamics Landau . Basically one writes down the relevant hydrodynamic equations for a binary mixture. The equations are then linearized for small fluctuations of the hydrodynamic variables around their macroscopic values and random source terms are added to them.

The equations thus obtained for a binary mixture where a uniform concentration gradient is present are vailati98

$$\{\begin{array}{c}\delta c(𝐪,\omega )\left(i\omega +D_0q^2\right)=\delta 𝐮(𝐪,\omega )\mathbf{}c+i𝐪𝐅(𝐪,\omega )\hfill \\ \delta 𝐮(𝐪,\omega )\left(i\omega +\nu q^2\right)=\beta \delta c(𝐪,\omega )\left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)+\frac{i}{\rho }𝐪𝐒(𝐪,\omega )\hfill \end{array},$$
(3)
where $`D_0`$ is the bare diffusion coefficient, $`\nu `$ the kinematic viscosity, $`𝐠`$ the gravitational acceleration, $`\beta =\rho ^1\left(\rho /c\right)`$ and $`c(𝐪,\omega )`$ is the space-time Fourier transform for the concentration fluctuations defined by

$$\delta c(𝐪,\omega )=\frac{1}{\left(2\pi \right)^4}e^{i𝐪𝐱+i\omega t}\delta c(𝐱,t)d𝐱dt,$$
(4)
and an analogous expression holds for the velocity fluctuations.

In Eq. (3)$`𝐅(𝐪,\omega )`$ and $`𝐒(𝐪,\omega )`$ represent random source currents for concentration fluctuations and velocity fluctuations, respectively. As customary, their correlation properties are assumed to retain their equilibrium values cohen71 :

$$F_i(𝐪,\omega )F_j(𝐪^{},\omega ^{})=\delta _{ij}\delta (𝐪+𝐪^{})\delta (\omega +\omega ^{})\frac{K_BT}{8\pi ^4}\frac{D_0}{\rho }\left(\frac{c}{\mu }\right)_{p,T},$$
(5)
$$S_i(𝐪,\omega )S_j(𝐪^{},\omega ^{})=\delta _{ij}\delta (𝐪+𝐪^{})\delta (\omega +\omega ^{})\frac{K_BT}{8\pi ^4}\nu \rho ,$$
(6)
$$F_i(𝐪,\omega )S_j(𝐪^{},\omega ^{})=0.$$
(7)
Let us briefly comment Eq. (3). The left terms specify that concentration and velocity fluctuations decay with lifetimes $`1/D_0q^2`$ and $`1/\nu q^2`$, respectively. The right terms are source terms for the fluctuations. Beyond the random currents $`𝐅(𝐪,\omega )`$ and $`𝐒(𝐪,\omega )`$, they contain the cross-coupling between velocity and concentration fluctuations. The cross-coupling term in the concentration equations shows that a velocity fluctuation $`\delta 𝐮`$ couples with the macroscopic concentration gradient $`\mathbf{}c`$, thus giving rise to a concentration fluctuation. The cross coupling term in the velocity equation represents the acceleration acting on a concentration fluctuation due to the buoyancy force. We have assumed that the fluid is incompressible; therefore the acceleration acting on the fluctuation is the component of gravity perpendicular to the wave vector $`𝐪`$.

The explicitation of $`\delta c`$ and $`\delta 𝐮`$ in Eq. (3) yields

$$\delta c(𝐪,\omega )=\frac{i\left(i\omega +\nu q^2\right)𝐪𝐅(𝐪,\omega )\frac{i}{\rho }𝐪𝐒(𝐪,\omega )\mathbf{}c}{\left(i\omega +D_0q^2\right)\left(i\omega +\nu q^2\right)+\beta \left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)\mathbf{}c}$$
(8)
and

$`\delta 𝐮(𝐪,\omega )={\displaystyle \frac{i𝐪𝐅(𝐪,\omega )\beta \left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)+\frac{i}{\rho }\left(i\omega +D_0q^2\right)𝐪𝐒(𝐪,\omega )}{\left(i\omega +D_0q^2\right)\left(i\omega +\nu q^2\right)+\beta \left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)\mathbf{}c}}`$
$`+{\displaystyle \frac{\frac{i}{\rho }\frac{\beta }{i\omega +\nu q^2}\left[𝐪𝐒(𝐪,\omega )\mathbf{}c\left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)𝐪𝐒(𝐪,\omega )\mathbf{}c\left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)\right]}{\left(i\omega +D_0q^2\right)\left(i\omega +\nu q^2\right)+\beta \left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)\mathbf{}c}}.`$ (9)
We are now in the position to evaluate the mass flow $`𝐣_f`$. In the reciprocal space Eq. (2) becomes

$$𝐣_f=\rho \delta c(𝐪,\omega )\delta 𝐮(𝐪^{},\omega ^{})d𝐪d𝐪^{}d\omega d\omega ^{}.$$
(10)
The integrand in Eq. (10) can be decomposed as the sum of terms containing correlations of the random currents $`𝐅`$ and $`𝐒`$. Due to (7) the only non vanishing terms are those originated from self-correlations:

$$\delta c(𝐪,\omega )\delta 𝐮(𝐪^{},\omega ^{})=\delta c(𝐪,\omega )\delta 𝐮(𝐪^{},\omega ^{})_F+\delta c(𝐪,\omega )\delta 𝐮(𝐪^{},\omega ^{})_S,$$
(11)
where

$$\delta c(𝐪,\omega )\delta 𝐮(𝐪^{},\omega ^{})_F=\delta \left(𝐪+𝐪^{}\right)\delta \left(\omega +\omega ^{}\right)\frac{K_BTD_0}{8\pi ^4\rho }\left(\frac{c}{\mu }\right)_{p,T}\beta \frac{q^2\left(i\omega +\nu q^2\right)\left(𝐠𝐪\frac{𝐠𝐪}{q^2}\right)}{\left|\left(i\omega +D_0q^2\right)\left(i\omega +\nu q^2\right)+\beta \left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)\mathbf{}c\right|^2}$$
(12)
and

$$\delta c(𝐪,\omega )\delta 𝐮(𝐪^{},\omega ^{})_S=\delta \left(𝐪+𝐪^{}\right)\delta \left(\omega +\omega ^{}\right)\frac{K_BT\nu }{8\pi ^4\rho }\frac{q^2\left(i\omega +D_0q^2\right)\left(\mathbf{}c𝐪\frac{\mathbf{}c𝐪}{q^2}\right)}{\left|\left(i\omega +D_0q^2\right)\left(i\omega +\nu q^2\right)+\beta \left(𝐠𝐪\frac{𝐪𝐠}{q^2}\right)\mathbf{}c\right|^2}.$$
(13)
The denominator of Eqs. (12-13) is an even function of $`\omega `$; therefore in the evaluation of the integral in Eq. (10) we can neglect odd functions of $`\omega `$ in the numerator.

By combining Eqs. (10-13) and by assuming that $`\widehat{𝐳}=𝐠/g=\mathbf{}c/\left|c\right|`$, the current $`𝐣_f`$ due to non equilibrium fluctuations is finally determined:

$$𝐣_f=\rho D_f\left[\mathbf{}c\beta 𝐠\left(\frac{c}{\mu }\right)_{p,T}\right],$$
(14)
where

$$D_f=\frac{K_BT\nu D_0}{8\pi ^4\rho }\frac{q^4\left(\widehat{𝐳}𝐪\frac{\widehat{𝐳}𝐪}{q^2}\right)\widehat{𝐳}}{\left|\left(i\omega +D_0q^2\right)\left(i\omega +\nu q^2\right)+\beta gc\left(\widehat{𝐳}𝐪\frac{\widehat{𝐳}𝐪}{q^2}\right)\widehat{𝐳}\right|^2}d𝐪d\omega .$$
(15)
Surprisingly the flow $`𝐣_f`$ displays the customary dependence of Fick’s flow on the driving force $`\mathbf{}c`$ and on the barodiffusion flow.

Equation (15) can be rewritten by factorizing its denominator as $`\left(\omega ^2+\omega _+^2\right)\left(\omega ^2+\omega _{}^2\right)`$ and by approximating the two roots as $`\omega _+^2=D_0^2q^4\left\{1+\left[q/q_{RO}\left(g\right)\right]^4\left(\widehat{𝐳}𝐪\frac{\widehat{𝐳}𝐪}{q^2}\right)\widehat{𝐳}\right\}^2`$ and $`\omega _{}^2=\nu ^2q^4`$. Basically this approximation holds true for $`\nu >>D`$, a condition fulfilled by most binary liquid mixtures; further details are discussed in Ref. vailati98 :

$$D_f=\frac{K_BT}{\left(2\pi \right)^3\rho \nu }\frac{q^2\left(\widehat{𝐳}𝐪\frac{\widehat{𝐳}𝐪}{q^2}\right)\widehat{𝐳}}{q^4+q_{RO}\left(g\right)^4\left(\widehat{𝐳}𝐪\frac{\widehat{𝐳}𝐪}{q^2}\right)\widehat{𝐳}}d𝐪,$$
(16)
where

$$q_{RO}\left(g\right)=\sqrt[4]{\frac{g\beta c}{\nu D_0}}$$
(17)
is a roll-off wave vector due to the gravity force segre93\_3 . As thoroughly discussed in Ref. vailati98 , the roll-off wave vector marks the transition from the diffusive relaxation of fluctuations at large wave vectors to the gravitational stabilization of fluctuations at smaller ones. At the roll-off wave vector the diffusional and gravitational fluctuation relaxation times coincide.

By introducing polar coordinates Eq. (16) becomes

$$D_f^Q=\frac{K_BT}{\left(2\pi \right)^3\rho \nu }_0^Qdq_0^\pi d\theta _0^{2\pi }d\phi \frac{q^4\mathrm{sin}^3\left(\theta \right)}{q^4+q_{RO}\left(g\right)^4\mathrm{sin}^2\left(\theta \right)}.$$
(18)
The integral on $`q`$ is divergent. Therefore we have regularized it by introducing a brute-force cut-off at the arbitrary wave vector $`Q`$. The evaluation of the integrals in $`\phi `$ and $`q`$ yields

$$D_f^Q=\frac{K_BT}{\left(2\pi \right)^2\rho \nu }_0^\pi d\theta Q\left[1\frac{\sqrt{2}\pi }{4}\frac{q_{RO}}{Q}\sqrt{\mathrm{sin}\theta }+O\left(\frac{q_{RO}\sqrt{\mathrm{sin}\theta }}{Q}\right)^4\right]\mathrm{sin}^3\theta .$$
(19)
By assuming that $`q_{RO}\left(g\right)/Q1`$ Eq. (19) becomes

$$D_f^Q\left(g\right)=\frac{4}{3}\frac{K_BT}{\left(2\pi \right)^2\rho \nu }Q\sqrt{2}I\frac{\pi }{4}\frac{K_BT}{\left(2\pi \right)^2\rho \nu }q_{RO}\left(g\right)+QO\left(q_{RO}\left(g\right)/Q\right)^4,$$
(20)
where

$$I=_0^\pi d\theta \mathrm{sin}^{7/2}\theta 1.248599.$$
(21)
In the absence of gravity $`q_{RO}=0`$, and Eq.(20) becomes

$$D_f^Q=\frac{2}{3}\frac{K_BT}{\pi \rho \nu \mathrm{\Lambda }},$$
(22)
where $`\mathrm{\Lambda }=2\pi /Q`$ can be roughly assumed to be of the order of correlation length of the fluid. It should be noticed that the diffusion coefficient $`D_f`$ of Eq. (22) does not depend on the bare mass diffusion coefficient $`D_0`$ which appeared in the linearized fluctuating hydrodynamic equations (3) and (5).

Quite remarkably, Eq. (22) exhibits the usual Stokes-Einstein dependence on temperature, viscosity and correlation length, the proportionality factor being of the right order of magnitude. This result, combined with Eq.(14), strongly suggests that the mass transfer due to fluctuations and Fick’s mass flow can be identified. This is a striking result, as it shows for the first time that the diffusive mass transfer is entirely determined by non equilibrium fluctuations occurring at length scales ranging from the microscopic up to the macroscopic one.

Although it is physically reasonable to assume that wave vectors larger than a typical molecular one have no meaning, the imposition of a cut-off involves the introduction of an arbitrary length scale, which in turn requires to introduce a microscopic model for the molecular interactions. This necessity can be eliminated by a renormalization of the linearized hydrodynamic equations. The normalization procedure we will discuss in the following closely mirrors the renormalization of mass of Quantum Field Theory lebellac  and of Classical Electrodynamics feynman . The overall diffusive mass flow is the superposition of the contribution due to fluctuations $`𝐣_f`$ and the bare fickean mass flow:

$$𝐣\left(g\right)=𝐣_f^Q\left(g\right)+𝐣_0^Q\left(g\right)=\rho \left(D_f^Q\left(g\right)+D_0^Q\right)\left[\mathbf{}c\beta 𝐠\left(\frac{c}{\mu }\right)_{p,T}\right].$$
(23)
The renormalization condition is

$$𝐣\left(g=0\right)=\rho D\mathbf{}c,$$
(24)
where $`D`$ is the phenomenological diffusion coefficient in the absence of gravity.

By combining Eq. (23) and Eq. (24) the renormalization condition can be rewritten as

$$D_0^Q=DD_f^Q\left(g=0\right),$$
(25)
and, by inserting Eqs. (20,25) into Eq. (23), the renormalized mass current is

$$𝐣_R\left(g\right)=\rho D_R\left(g\right)\left[\mathbf{}c\beta 𝐠\left(\frac{c}{\mu }\right)_{p,T}\right],$$
(26)
where $`D_R`$ is the renormalized diffusion coefficient

$$D_R\left(g\right)=D\sqrt{2}\frac{\pi }{4}I\frac{K_BT}{\left(2\pi \right)^2\rho \nu }q_{RO}\left(g\right).$$
(27)
In Eq. (27) the $`Q`$ dependence has been eliminated by letting $`Q\mathrm{}`$, because $`D_R`$ does not diverge. The microscopic details have been reabsorbed into $`D`$, the only experimental parameter needed by the theory in order to obtain the value of the diffusion coefficient in the presence of gravity.

Equation (27) shows that the diffusion coefficient is smaller due to the presence of gravity, the additional term having the Stokes-Einstein form, the only length scale involved being $`2\pi /q_{RO}`$.

This modification of the diffusion coefficient is due to the fact that, in the presence of gravity, fluctuations with wave vector smaller than $`q_{RO}`$ are stabilized and therefore do not contribute to the mass transfer. Therefore $`q_{RO}`$ acts as an infrared cut-off for the fluctuations. As the length scales stabilized by gravity are large, the relative correction introduced by gravity is very small, being of the order of $`q_{RO}\xi `$, where we have assumed a Stokes-Einstein dependence of the phenomenological diffusion coefficient $`D`$ on the correlation length $`\xi `$ of the fluid. For an ordinary liquid mixture $`q_{RO}`$ can be as high as $`10^3\mathrm{cm}^1`$ and $`\xi `$ is of the order of $`10^7\mathrm{cm}`$. Therefore the gravitational correction is estimated to be roughly one part in $`10^4`$. In the case of critical fluids and macromolecular solutions the correlation length can grow very large, thus suggesting that the gravitational correction to the diffusion coefficient can get large. However, in order for our hydrodynamic description to be valid, the correlation length has to be significantly smaller than $`1/q_{RO}`$. Therefore the largest correction to the diffusion coefficient reliably predicted by this description can be roughly estimated to be of the order of 10%. Such a variation in the diffusion coefficient should be observable by performing microgravity experiments on critical fluids.

The authors wish to thank Marzio Giglio for stimulating discussion. Work partially supported by the Italian Space Agency (ASI).


text:
# Many–particle entanglement with Bose–Einstein condensates

## Abstract

We propose a method to produce entangled states of several particles starting from a Bose-Einstein condensate. In the proposal, a single fast $`\pi /2`$ pulse is applied to the atoms and due to the collisional interaction, the subsequent free time evolution creates an entangled state involving all atoms in the condensate. The created entangled state is a spin-squeezed state which could be used to improve the sensitivity of atomic clocks.

The possibility of creating and manipulating entangled states of many–particle systems has recently boosted the field of quantum information since it may yield new applications which rely on the basic principles of Quantum Mechanics . So far, up to four atoms have been entangled in a controlled way . The experimental achievement of Bose–Einstein condensation  has also raised a lot of attention since it may lead to applications in several fields of Science. Some of these applications are based on the fact that condensates can be considered as pure states at the single particle level, which is a crucial requirement for the production of entangled state. Thus, a natural question is to investigate whether Bose–Einstein condensates can also be used in some applications of quantum information. In this paper we show that this is indeed the case. We propose a method to obtain substantial entanglement of a large number of atoms with present technology. In our proposal, a single resonant pulse is applied to all atoms in the condensate, and the collisional interaction entangles the atoms in the subsequent free evolution.

Consider a set of $`N`$ two–level atoms confined by some external trap. In order to describe the internal properties of these atoms, it is convenient to let the internal states $`|a_n`$ and $`|b_n`$ of the $`n`$-th atom represent the two states of a fictitious spin 1/2 particle with angular momentum operators $`j_z^{(n)}=1/2(|aa|_n|bb|_n)`$, $`j_x^{(n)}=1/2(|ba|_n+|ab|_n)`$, and $`j_y^{(n)}=i/2(|ba|_n|ab|_n)`$. We consider collective effects of the atoms which are described by total angular momentum operators, $`\stackrel{}{J}=_{n=1}^N\stackrel{}{j}^{(n)}`$. The entanglement properties of the atoms can be expressed in terms of the variances and expectation values of these operators. In the appendix we show that if

$$\xi ^2\frac{N(\mathrm{\Delta }J_{\stackrel{}{n}_1})^2}{J_{\stackrel{}{n}_2}^2+J_{\stackrel{}{n}_3}^2}<1,$$
(1)
where $`J_\stackrel{}{n}\stackrel{}{n}\stackrel{}{J}`$ and the $`\stackrel{}{n}`$’s are mutually orthogonal unit vectors, then the state of the atoms is non–separable (i.e. entangled). The parameter $`\xi ^2`$ thus characterizes the atomic entanglement, and states with $`\xi ^2<1`$ are often referred to as “spin squeezed states” . In the following we show how to reduce $`\xi ^2`$ by several orders of magnitude using the collisional interactions between atoms in a Bose–Einstein condensate.

We consider a two component weakly interacting Bose-Einstein condensate, which is present in several laboratories, and we assume that the interactions do not change the internal state of the atoms. This situation is described by the second quantized Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{j=a,b}{}}{\displaystyle d^3r\widehat{\mathrm{\Psi }}_j^{}(𝐫)H_{0,j}\widehat{\mathrm{\Psi }}_j(𝐫)}`$ (2)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=a,b}{}}U_{jj}{\displaystyle d^3r\widehat{\mathrm{\Psi }}_j^{}(𝐫)\widehat{\mathrm{\Psi }}_j^{}(𝐫)\widehat{\mathrm{\Psi }}_j(𝐫)\widehat{\mathrm{\Psi }}_j(𝐫)}`$ (3)
$`+`$ $`U_{ab}{\displaystyle d^3r\widehat{\mathrm{\Psi }}_a^{}(𝐫)\widehat{\mathrm{\Psi }}_b^{}(𝐫)\widehat{\mathrm{\Psi }}_a(𝐫)\widehat{\mathrm{\Psi }}_b(𝐫)},`$ (4)
where $`H_{0,j}`$ is the one particle Hamiltonian for atoms in state $`j`$ including the kinetic energy and the external trapping potential $`V_j(𝐫)`$, $`\widehat{\mathrm{\Psi }}_j(𝐫)`$ is the field operator for atoms in the state $`j`$, $`U_{jk}=4\pi \mathrm{}^2a_{jk}/m`$ is the strength of the interaction between particles of type $`j`$ and $`k`$, parameterized by the scattering length $`a_{jk}`$, and $`m`$ is the atomic mass.

Assume that we start with a Bose–Einstein condensate in state $`|a`$ at very low temperature ($`T0`$), so that all the atoms are in a single particle (motional) state $`|\varphi _0`$. A fast $`\pi /2`$ pulse between the states $`|a`$ and $`|b`$ prepares the atoms in the state $`|\varphi _0^N(|a+|b)^N/2^{N/2}`$ which is an eigenstate of the $`J_x`$ operator with eigenvalue $`N/2`$. If we choose in (1) $`\stackrel{}{n}_1=[0,\mathrm{cos}(\theta ),\mathrm{sin}(\theta )]`$ and $`\stackrel{}{n}_2`$ along the $`x`$ axis of the fictitious angular momentum, we have $`\xi _\theta ^2=1`$ at $`t=0`$. By using the equations of motion of the angular momentum operators in the Heisenberg picture we find the time derivative of $`\xi _\theta ^2`$ at $`t=0`$

$$\frac{d}{dt}\xi _\theta ^2=\mathrm{sin}(2\theta )\frac{(N1)\left(U_{aa}+U_{bb}2U_{ab}\right)}{2\mathrm{}}d^3r\left|\varphi _0\right|^4.$$
(5)
This equation immediately shows that spin squeezing will be produced for certain angles $`\theta `$ if $`U_{aa}+U_{bb}2U_{ab}`$ as it is in the Na experiments .

To quantify the amount of squeezing which may be obtained, we assume identical trapping potentials $`V_a(𝐫)=V_b(𝐫)`$ and identical coupling constants for interaction between atoms in the same internal state $`a_{aa}=a_{bb}>a_{ab}`$. Physically, this could correspond to the $`|F=1,M_F=\pm 1`$ hyperfine states of Na trapped in an optical dipole trap. Due to the symmetry of these states their scattering lengths and trapping potentials will be identical; moreover, due to angular momentum conservation there are no spin exchanging collisions between these states as required by our model. This is exactly the experimental setup used in Ref. . In this experiment it is shown that these states have an antiferromagnetic interaction $`a_{ab}<a_{aa}`$ which according to (5) enables the production of squeezed states. To avoid spin changing collisions that populate the state $`|M_F=0`$ one has to slightly modify such an experimental set–up. If the $`F=1`$ manifold is coupled to the $`F=2`$ with a far off-resonant blue detuned $`\pi `$-polarized microwave field, the $`|F=1,M_F=0`$ state is raised in energy with respect to the $`|F=1,M_F=\pm 1`$ states, since the Clebsch-Gordan coefficient for the $`|F=1,M_F=0|F=2,M_F=0`$ is larger than the coefficients for the $`|F=1,M_F=\pm 1|F=2,M_F=\pm 1`$ transitions, and spin exchanging collisions become energetically forbidden. If for instance one chooses a detuning $`\delta =(2\pi )25`$ MHz and a resonant Rabi frequency for the $`11`$ transition of $`\mathrm{\Omega }=(2\pi )2`$ MHz, the energy separation is $`\mathrm{\Delta }E=640`$nK. With a typical chemical potential $`\mu 220`$ nK, this energy separation is much higher than the available energy in the collisions and the $`|F=1,M=0`$ state is completely decoupled. In Ref.  a much smaller energy difference is shown to exclude spin exchanging collision and we therefore expect that much weaker fields will suffice.

The assumption $`a_{aa}=a_{bb}`$ has several advantages. Firstly, it reduces the effect of fluctuations in the total particle number. If $`a_{aa}a_{bb}`$ the mean spin performs a $`N`$ dependent rotation around the $`z`$-axis, and fluctuations in the number of particles introduces an uncertainty in the direction of the spin which effectively reduces the average value and introduces noise into the system. With $`a_{aa}=a_{bb}`$ the mean spin remains in the $`x`$-direction independent of the number of atoms in the trap. Secondly, this condition ensures a large spatial overlap of different components of the wavefunction. After the $`\pi /2`$ pulse the spatial wavefunction is no longer in the equilibrium state. That is, due to the atomic repulsions (which are now different than before since the atoms are in different internal states), the spatial distribution of the atomic cloud will start oscillating. Furthermore, since the state of the system is now distributed over a range of number of particles in the $`|a`$ state ($`N_a`$), and since this number enters into the time evolution, the $`N_a`$ dependent wavefunctions $`\varphi _a`$ and $`\varphi _b`$ are different for particles in the states $`|a`$ and $`|b`$. With $`a_{aa}=a_{bb}`$, $`\varphi _a`$ and $`\varphi _b`$ are identical if $`N_a`$ equals the average number $`N/2`$. In the limit of large $`N`$, the width of the distribution on different $`N_a`$’s is much smaller than $`N_a`$ and all the spatial wavefunctions are approximately identical $`\varphi _a(N_a,t)\varphi _b(N_a,t)\varphi _0(t)`$. This relation is only true if $`a_{aa}>a_{ab}`$ where small deviations from the average wavefunction perform small oscillations. In the opposite case the deviations grow exponentially  resulting in a reduction of the overlap of the $`a`$ and the $`b`$ wavefunctions and a reduced squeezing. The advantages mentioned above could also be achieved with $`a_{aa}a_{bb}`$ by using the breath–together solutions proposed in Ref. .

Before analyzing quantitatively the complete system, we estimate the amount of spin squeezing we can reach with our proposal by using a simple model. Assuming the same wavefunction $`\varphi _0`$ for both $`|a`$ and $`|b`$ atoms is constant and independent of the number $`N_a`$, the spin dependent part of the Hamiltonian (4) may be written as $`H_{\mathrm{spin}}=\mathrm{}\chi J_z^2`$, where $`\chi `$ depends on the scattering lengths and the wavefunction $`\varphi _0`$. The spin squeezing produced by this Hamiltonian can be calculated exactly . In the limit of large $`N`$, the minimum obtainable squeezing parameter is $`\xi _\theta ^2=\frac{1}{2}\left(\frac{3}{N}\right)^{2/3}`$, which indicates that our proposal might produce a reduction of $`\xi ^2`$ by a factor of $`N^{2/3}`$ which would be more than three orders of magnitude with $`10^5`$ atoms in the condensate.

In contrast to the simplified Hamiltonian $`H_{\mathrm{spin}}`$ the real Hamiltonian (4) will also entangle the internal and motional states of the atoms which is a source of decoherence for the spin squeezing. To quantify this effect we have performed a direct numerical integration following the procedure developed in Ref. . We split the whole Hilbert space into orthogonal subspaces containing a fixed number of particles $`N_a`$ and $`N_b=NN_a`$ in each of the internal states, respectively. In each subspace we make a Hartree-Fock variational ansatz in terms of three-dimensional spatial wavefunctions $`\varphi _a(N_a,t)`$ and $`\varphi _b(N_b,t)`$, which are evolved according to the time dependent coupled Gross-Pitaevskii equations. This is an approximation to the full problem which is valid in the limit of low temperatures and short times, where the population of the Bogoliubov modes is small. Particularly it is a good approximation in our case with $`a_{ab}<a_{aa}`$, $`a_{bb}`$ where there are no demixing instabilities. With this procedure, the decoherence induced by the entanglement with the motional state is effectively taken into account. Together with the prediction from the simple Hamiltonian $`H_{\mathrm{spin}}`$ the result of the simulation is shown in Fig. 1. The two curves are roughly in agreement confirming that the system is able to approximate the results of the Hamiltonian $`H_{\mathrm{spin}}`$. The numerical solution shows fluctuations due to the oscillations of the spatial wavefunction. The large dips at $`\omega t`$ 4, 9, 13, and 18 are the points where the atomic cloud reach the initial width. At these instants the overlap of the wavefunctions is maximal and the two curves are very close (up to a factor of two). The small dips at $`\omega t=`$2, 7, 11, and 16 corresponds to the points of maximum compression. With the realistic parameters used in the figure, our simulation suggests that three orders of magnitude squeezing is possible. Also, note the time scale in the figure. The maximally squeezed state is reached after approximately two oscillation periods in the trap. For a fixed ratio of the scattering lengths $`a_{ab}/a_{aa}`$, the optimal time scales as $`(a_{aa}/d_0)^{2/5}N^{1/15}`$, where $`d_0=\sqrt{\mathrm{}/(m\omega )}`$ is the width of the ground state of the harmonic potential.

The analysis so far has left out a number of possible imperfections. Specifically, we have assumed that all scattering length are real so that no atoms are lost from the trap and we have not considered the role of thermal particles. To estimate the effect particle losses we have performed a Monte Carlo simulation of the evolution of squeezing from the Hamiltonian $`H_{\mathrm{spin}}`$. The particle loss is phenomenologically taken into account by introducing a loss rate $`\mathrm{\Gamma }`$ which is identical for atoms in the $`|a`$ and the $`|b`$ state. In Fig. 2 we show the obtainable squeezing in the presence of loss. Approximately 10 % of the atoms are lost at the time $`\chi t610^4`$ where the squeezing is maximally without loss. With the parameters of Fig. 1 this time corresponds to roughly two trapping periods. Such a large loss is an exaggeration of the loss compared to current experiments and the simulation indicate that even under these conditions, squeezing of nearly two orders of magnitude may be obtained. On the other hand, the effects of thermal particles can be suppressed at sufficiently low temperatures but due to the robustness with respect to particle loses shown in Fig. 2, we expect to obtain high squeezing even at some finite temperatures.

In conclusion we believe that we have presented a simple and robust method to produce entangled states of a large number of atoms with present technology. The produced entangled state are interesting in fundamental physics and they also have possible technological applications in atomic clocks , where the projection noise $`(\mathrm{\Delta }J_\theta )^2`$ is currently the main source of noise . In particular Boyer and Kasevich  have shown how to use an entangled state with half of the particles $`N/2`$ in one internal state $`a`$, and the other half in the other state $`b`$ to obtain the Heisenberg limit in atomic interferometric measurements. On the other hand, it has been shown that if the atoms are prepared in a state with $`\xi ^2<1`$ before they are injected into an atomic clock, one can reduce the frequency noise (variance in the frequency measurements) or the measuring time to obtain a desired precision by a factor $`\xi ^2`$ as compared to the case in which one uses atoms in an uncorrelated state . If the squeezing produced by our proposal is transferred into a suitable clock transition, and the atomic cloud is allowed to expand so that the role of collisions is reduced during the frequency measurement, the states could be directly applied in the current set–up of atomic clocks. Other theoretical proposals for noise reduction in neutral atom systems have been made , and a weak squeezing of the spin has recently been produced experimentally . However, our proposal has the considerable advantage that it offers a very strong noise reduction, and it is directly applicable to existing experimental set–ups. In future experiments with negligible particle loss even for very long interaction times, the Hamiltonian $`H_{\mathrm{spin}}`$ could also be used to produce maximally entangled state of any number of atoms , which could reduce the frequency noise to the fundamental limit of quantum mechanics .

Appendix
Here we present the derivation of Eq. (1) as a criterion for entanglement. An $`N`$-particle density matrix $`\rho `$ is defined to be separable (non-entangled) if it can be decomposed into

$$\rho =\underset{k}{}P_k\rho _k^{(1)}\rho _k^{(2)}\mathrm{}\rho _k^{(N)},$$
(6)
where the coefficients $`P_k`$ are positive real numbers fulfilling $`_kP_k=1`$, and $`\rho _k^{(i)}`$ is a density matrix for the $`i`$’th particle. The variance of $`J_z`$ may be described as $`(\mathrm{\Delta }J_z)^2=\frac{N}{4}_kP_k_ij_z^{(i)}_k^2+_kP_kJ_z_k^2J_z^2`$, and using Cauchy-Schwarz’s inequality and $`j_x^{(i)}_k^2+j_y^{(i)}_k^2+j_z^{(i)}_k^21/4`$ we find three inequalities for separable states $`_kP_kJ_z_k^2J_z^2`$, $`_kP_k_ij_z^{(i)}_k^2\frac{N}{4}+_kP_k_ij_x^{(i)}_k^2+j_y^{(i)}_k^2`$, and $`J_x^2N_kP_k_ij_x^{(i)}_k^2`$. From these inequalities we immediately find that any separable state obeys $`\xi ^21`$ and hence any state with $`\xi ^2<1`$ is non-separable.

Acknowledgments

This work was supported by the Austrian Science Foundation, the European Union project EQUIP, the TMR european network, the ESF under the PESC program ”Quantum Information”, the Institute for Quantum Information GmbH, and Thomas B. Thriges Center for Kvanteinformatik. A. S. acknowledges the hospitality of the University of Innsbruck.


text:
# Quantum mechanics can be led by metrics in Minkowski’s time-space world

## Abstract

A body of theory is completely different between relativity theory and quantum mechanics. Most targeting physical phenomena are different between them as well. Despite that both theories describe our time-space world, no one has succeeded in arriving at one of these theories by using the other one, deductively. Fusion of these two theories produced quantum field theory. However this theory impresses us as a theory that was made by merging the two theories by force. So this theory cannot indicate the reason why both relativity theory and quantum mechanics exist in this world at the same time. Here I show you that quantum mechanics can be led by the equation of Minkowski’s metrics, which stands up in the relativistic world of space and time. This is the first report proving mathematically that quantum mechanics can be born reductively in our relativistic world.
PACS numbers: 03.65.Bz, 03.75.-b, 11.90.+t
Keywords: Quantum mechanics, Relativity theory, Minkowski’s time-space world

The Institute for Enzyme Research, The University of Tokushima 3-18-15 Kuramoto-cho, Tokushima, 770-8503, Japan
Present temporary address: Amgen Institute/ Ontario Cancer Institute, 620 University Avenue, Suite 706, Toronto, Ontario, M5G 2C1, Canada
Correspondence should be addressed to the author (e-mail: sakai@ier.tokushima-u.ac.jp, Tel.:+81-88-633-7430, Fax:+81-88-633-7431)
The strong success of quantum mechanics is not surprising considering its much greater degree of efficacy, on an atomic scale, compared to any other theory. However, it has also been difficult to understand the fundamental implication of quantum mechanics because of its specificity from the structure of the theory including wave-particle duality. We are able to see such specificity from the aspect that the theory has borne lots of paradoxical questions. On the other hand, the theory includes another important aspect: that it seems to be incompatible with aspects of relativity theory that exert great validity against our macroscopic scale. Here I would like to propose some novel techniques and ideas in order to understand the relationship between these two theories and to solve the paradoxical questions of quantum mechanics, such as wave-particle duality or the measurement problem. Here it is shown that quantum mechanics can be led from the equation in Minkowski’s space-time world.

The differential distance $`ds`$ in Minkowski’s space-time world (Minkowski’s metrics) is given by

$$(ds)^2=(cdt)^2+(dx)^2+(dy)^2+(dz)^2$$
(1)
For our simplified spacetime, with just one space dimension, this reduces to

$$(ds)^2=(cdt)^2+(dx)^2$$
(2)
This can also be rewritten by

$`(ds)^2`$ $`=`$ $`(cdt)^2(\mathrm{cos}^2\theta _t(i\mathrm{sin}\theta _t)^2)+(dx)^2(\mathrm{cos}^2\theta _x(i\mathrm{sin}\theta _x)^2)`$ (3)
$`=`$ $`(cdt\mathrm{cos}\theta _t,icdt\mathrm{sin}\theta _t)(cdt\mathrm{cos}\theta _t,icdt\mathrm{sin}\theta _t)`$
$`+(dx\mathrm{cos}\theta _x,idx\mathrm{sin}\theta _x)(dx\mathrm{cos}\theta _x,idx\mathrm{sin}\theta _x)`$
$`=`$ $`d\psi _td\psi _t^{}+d\psi _xd\psi _x^{}`$
$`i=\sqrt{1}`$
$`d\psi (\theta )`$ $`=`$ $`dx(\mathrm{cos}\theta ,i\mathrm{sin}\theta )`$
$`d\psi ^{}(\theta )`$ $`=`$ $`dx(\mathrm{cos}\theta ,i\mathrm{sin}\theta )=d\psi (\theta )`$
These $`d\psi `$ and $`d\psi ^{}`$ are expressed as the coordinates on a circumference in a complex space of real and imaginary axis (Figure 1). I will discuss later concerning this complex space. First of all, you consider a complex world consisting with an imaginary space lying at right angles to a real space for either its time dimension or its space dimension. Then you think about a circle consisting of a minute segment ($`dx`$) as a radius (Figure 1). The “$`d\psi (\theta )=dx(\mathrm{cos}\theta ,i\mathrm{sin}\theta )`$ ” or the “$`d\psi ^{}(\theta )=dx(\mathrm{cos}\theta ,i\mathrm{sin}\theta )`$ ” shown above represents points on the circumference. If this condition $`\psi `$ is involved in a quantum behavior, $`\theta `$ must be a function of the coordinates of the real space, $`x`$. When you consider the existence of a particle in this minute region $`\mathrm{\Delta }x`$, the condition must have a feature of a wave in the imaginary space. If you use $`\lambda `$ for its wavelength, you can substitute $`{\displaystyle \frac{2\pi }{\lambda }}x`$ for $`\theta `$. So if you write the particle’s condition as the function of $`x`$ like $`\psi (x)`$, the condition in the minute region $`d\psi (x)`$ is written as $`d\psi (x)=dx(i\mathrm{sin}{\displaystyle \frac{2\pi }{\lambda }}x,\mathrm{cos}{\displaystyle \frac{2\pi }{\lambda }}x)`$. This can be rewritten as $`{\displaystyle \frac{d\psi }{dx}}=(i\mathrm{sin}{\displaystyle \frac{2\pi }{\lambda }}x,\mathrm{cos}{\displaystyle \frac{2\pi }{\lambda }}x)`$. And this leads to $`\psi =i{\displaystyle \frac{\lambda }{2\pi }}(i\mathrm{sin}{\displaystyle \frac{2\pi }{\lambda }}x,\mathrm{cos}{\displaystyle \frac{2\pi }{\lambda }}x)=i{\displaystyle \frac{\lambda }{2\pi }}{\displaystyle \frac{d\psi }{dx}}`$. Assuming that the $`\lambda `$ is equivalent to the de Broglie wavelength, the relationship between the $`\lambda `$ and the momentum of the particle ($`P`$) can be described as $`{\displaystyle \frac{\lambda }{2\pi }}={\displaystyle \frac{\mathrm{}}{P}}`$. If you substitute this, you get the equation shown below.

$$P\psi =i\mathrm{}\frac{\psi }{x}$$
(4)
Next, you think about time dimension as well. You combine the equation $`d\psi (t)=dt(i\mathrm{sin}\theta ,\mathrm{cos}\theta )`$ and $`\theta =2\pi \nu t`$.
Then you get $`{\displaystyle \frac{d\psi }{dt}}=c(i\mathrm{sin}(2\pi \nu t),\mathrm{cos}(2\pi \nu t))`$. One of the simplest solutions of the equation is $`\psi (t)=i{\displaystyle \frac{c}{2\pi \nu }}(i\mathrm{sin}(2\pi \nu t),\mathrm{cos}(2\pi \nu t))=i{\displaystyle \frac{1}{2\pi \nu }}{\displaystyle \frac{d\psi }{dt}}`$. By substituting $`\nu ={\displaystyle \frac{E}{h}}`$ into this equation, you get

$$E\psi (t)=i\mathrm{}\frac{d\psi }{dt}$$
(5)
Equivalent results can be obtained concerning $`\psi ^{}(x)`$ and $`\psi ^{}(t)`$ as well. So, all the results can be put together with the simple equations shown below.

$$(ds)^2=g_{ij}dx_idx_j=g_{ij}d\psi _id\psi _j^{}$$
(6)
($`g_{ij}`$: metric tensor)

$$\psi ^{}(x)=\psi (x)=\frac{\mathrm{}}{P}(i\mathrm{cos}\frac{P}{\mathrm{}}x,\mathrm{sin}\frac{P}{\mathrm{}}x)$$
(7)
Here I succeeded in demonstrating that quantum mechanics can be led by Minkowski’s space-time world by assuming the existence of an imaginary space lying at right angles to a real space for either its time dimension or its space dimension. You can also see from the equation of $`\psi `$ or $`\psi ^{}`$ that the quantum mechanical area size of a particle is restricted in $`{\displaystyle \frac{\mathrm{}}{P}}`$, in which the particle exists as two waves written by complex numbers ($`\psi `$ and $`\psi ^{}`$). This means that any matter produces a quantum field by its momentum, and that the size of the field is restricted by the momentum. In addition, both of the quantum mechanical functions led by two wave functions $`\psi `$ and $`\psi ^{}`$ are equivalent, so these quantum mechanical features are undeniably with each other. The conclusion is that any particle moves as two quantum mechanical waves of $`\psi `$ and $`\psi ^{}`$ in its specific small areas. On the other hand, if you write the condition of a particle without the idea of an imaginary space, you cannot see such quantum mechanical features. Here you get a solution against the measurement problem of quantum theories. The measurement itself is equivalent to denying the existence of an imaginary space because we cannot measure any movement in imaginary spaces. So if you describe the movement according to the measurement, you never reach the functions of $`\psi `$ and $`\psi ^{}`$. So you can only measure the space by $`(dx)^2=g_{ij}dx_idx_j`$, which is an incomplete equation of the equation (6).
Next, I discuss the meaning of the complex world of real and imaginary axes. The imaginary axis lying at right angles to a real axis (Figure 1) can be interpreted as an axis of the fifth dimension, which is different from our four-dimensional world. However, more simple and easy alternative idea may be to interpret it as a time axis. On the contrary, an imaginary axis against time axis can be interpreted as an axis of the fifth or sixth dimensions. However, the more simple and easy way is to interpret it as a space axis. For instance, you put a space axis and an imaginary axis lying at right angles to the space axis as $`x`$ and $`iy`$, respectively ($`i=\sqrt{1}`$), metric tensor $`(ds)^2`$ is given by

$$(ds)^2=(dx)^2+(idy)^2=(dx)^2(dy)^2$$
(8)
This indicates that the $`iy`$ axis is in a time dimension. In addition, the equation (3) indicates that both the space and time axes can be treated as the complex of real and imaginary axes in a minute region. And this operation can lead quantum mechanics as I showed. These facts lead us to imagine that the essence of quantum mechanics lies in the transformation between space and time regions in the world of space and time.


text:
# 1 Clebsch Decomposition

## 1 Clebsch Decomposition

The decomposition of vectors, as well as higher rank tensors, into irreducible parts is an ancient and extremely useful tool in fluid mechanics, electrodynamics, and gravity. The longitudinal/transverse split separates an arbitrary $`d=3`$ Euclidean vector into a scalar plus a transverse vector,

$$𝐀𝐀^L+𝐀^T=\mathbf{}\lambda +\mathbf{}\times 𝐖^T,$$
(1)
the vector $`𝐖^T`$ being defined up to a gradient. This linear orthogonal (upon spatial integration) parameterization naturally decomposes the 3-vector A into two, with 1 and 2 components, respectively. The completeness of $`(\mathbf{}\lambda ,𝐖^T)`$ in representing A is evidenced by their uniqueness and invertibility (up to zero modes); in particular, any variational principle yields the same Euler-Lagrange system whether we vary A or first decompose it and then vary $`(\lambda ,𝐖^T)`$ separately. It is easy to check that since the field strength is

$$𝐁\mathbf{}\mathbf{}\times 𝐀=^2𝐖^T,$$
(2)
the Chern-Simons (first used physically in ) and Maxwell actions

$$I_{CS}=\frac{1}{2}\mu 𝐀𝐁=\frac{1}{2}\mu (\times 𝐖^T)(^2𝐖^T)I_M=\frac{1}{2}𝐁^2=\frac{1}{2}(^2𝐖^T)^2,$$
(3)
vary with $`A`$ into

$$\delta I_{CS}/\delta 𝐀=\mu 𝐁,\delta I_M/\delta 𝐀=\mathbf{}\times 𝐁$$
(4)
which agrees completely with the variations of (3) with respect to $`(\lambda ,𝐖^T)`$. The extra curl that appears upon varying $`𝐖^T`$ does not lose information: the curl is not a projector in the space of transverse functions.

Another ancient parameterization, due to Clebsch , has long been used in fluid mechanics , but is less familiar elsewhere. This decomposition is strikingly different – it involves only scalars rather than vectors and is nonlinear:

$$𝐀=\mathbf{}\omega +ϵ^{ab}\varphi _a\mathbf{}\varphi _b=\mathbf{}\omega ^{}+2\varphi _1\mathbf{}\varphi _2,𝐁=\mathbf{}\times 𝐀=ϵ^{ab}\mathbf{}\varphi _a\times \mathbf{}\varphi _b.$$
(5)
The representations (1) and (5) share the pure gauge longitudinal term, but the role of $`𝐖^T`$ has been assumed by the pair of scalars $`\varphi _a`$, $`a`$=1,2; $`ϵ^{ab}=ϵ^{ba}`$ is the alternating symbol. The construction of these scalars is described in the texts ; for the moment we note only that the $`\varphi _a`$ are not uniquely determined by A or B, nor can all (A,B) be put in this form in terms of well-defined scalar functions.

Our main interest is the particularly strange phenomenon (not hitherto noted to our knowledge), that the Clebsch decomposition (5) is incomplete in the very concrete sense that it fails (differently for each), when inserted into the two actions of (3) and varied with respect to $`(\omega ,\varphi _a)`$, to give the field equations (4) implied by the generic A-variations. The formal source of this incompleteness may be traced by comparing the decompositions (1) and (5),

$$^2\lambda =\mathbf{}𝐀=^2\omega +ϵ^{ab}\varphi _a^2\varphi _b,^2𝐖^T=\mathbf{}\times 𝐀=ϵ^{ab}\mathbf{}\varphi _a\times \mathbf{}\varphi _b.$$
(6)
Note the demand of (6) that the curl of A be expressible as the cross-product of two gradients: the curl must then lie on the line orthogonal to the plane they define. It follows that the Chern-Simons action becomes a pure surface term,

$$𝐀\mathbf{}\times 𝐀=\mathbf{}\omega \mathbf{}\times 𝐀=\mathbf{}(\omega \mathbf{}\times 𝐀).$$
(7)
In $`\omega =0`$ gauge, A is locally orthogonal to B, and in any gauge the bulk integral of (7) vanishes111When the volume integral of $`𝐀\mathbf{}\times 𝐀`$ is nonzero, its entire contribution in the Clebsch parameterization for $`𝐀`$ must come from surface terms. Consequently, the Clebsch potentials $`\omega ,\varphi _a`$ must be ill behaved, either in the finite volume or at infinity, so that the surface integrals do not vanish. For an example and further details, see . A fortiori, the variation of $`I_{CS}`$ with respect to $`(\omega ,\varphi _a)`$ is trivial. Pure Chern-Simons dynamics here reduces entirely to the boundary222In this respect, when CS is used as the “helicity” in magnetohydrodynamics, it bears a resemblance to gravitational energy; there is no sensible local (gauge invariant) density, but a perfectly well defined integral over the boundary., rather than defining the (trivially $`𝐁`$=0) behavior in the bulk.

Variation of the Maxwell action is also incomplete:

$$\delta \frac{1}{2}𝐁^2𝐁ϵ^{cd}(\mathbf{}\varphi _c\times \mathbf{}\delta \varphi _d),$$
(8)
resulting in the pair of Euler-Lagrange equations

$$(\mathbf{}\times 𝐁)\mathbf{}\varphi _a\mathbf{}(\varphi _a\mathbf{}\times 𝐁)=0,a=1,2.$$
(9)
That is, $`\mathbf{}\times 𝐁`$, like B itself, is orthogonal to the plane defined by $`(\mathbf{}\varphi _a)`$: the field equation is then

$$\mathbf{}\times 𝐁+\mu (𝐫)𝐁=0.$$
(10)
The function of integration, $`\mu (𝐫)`$, is constrained by the transversality of $`(𝐁,\mathbf{}\times 𝐁)`$ in (10) to have its gradient orthogonal to them both: $`\mathbf{}\mu `$ must lie in the ($`\mathbf{}\varphi _a`$) plane.

The indeterminacies in the two variational principles are rather different: from $`I_{CS}`$, there is no field equation at all; this is in fact traceable to the Darboux incompleteness of the Clebsch decomposition , which rests on the triviality of the bulk contribution to $`I_{CS}`$ that it imposes1. From the Maxwell action, we obtain only the direction of $`\mathbf{}\times 𝐁`$. Indeed, any gauge invariant action $`I[𝐁]`$ will provide the same form

$$[\mathbf{}\times \delta I/\delta 𝐁(x)]\mathbf{}\varphi _a=0,$$
(11)
and hence allow the same indeterminate $`\mu (r)𝐁`$ addition, subject to $`\mathbf{}\mu 𝐁=0`$.

The topologically massive-like  “Maxwell” equation (10) is quite unusual: Consider (9) in “covariant” notation,

$$(_\nu F^{\mu \nu })_\mu \varphi ^a_\mu (_\nu F^{\mu \nu }\varphi ^a)=0$$
(12)
whose solution is

$$_\nu F^{\mu \nu }\varphi ^a=ϵ^{\nu \alpha \beta }_\alpha Z_\beta ^a$$
(13)
where $`Z_\alpha ^a`$ are two arbitrary transverse 3-vectors. This means in turn that there are two (conserved) “currents”, for example, $`j^\mu =(\varphi _b\varphi ^b)^1\varphi _aϵ^{\mu \alpha \beta }_\alpha Z_\beta ^a.`$ Equations (9,10) and (12,13) present the field equations in two equivalent ways. If we instead introduce the explicit form (6) into (9), we get two (cubic) equations in $`\varphi `$ alone,

$$_\nu [ϵ^{bc}\varphi ^b(_{\mu \nu }^2\delta _{\mu \nu }^2)\varphi ^c]_\mu \varphi ^a=0.$$
(14)
A nonabelian generalization of the Clebsch decomposition has recently been suggested . It shares the abelian property that the CS term is a pure divergence and a similar degeneracy in the field equations occurs. The basic mechanism in both cases is that the variation of an action with respect to the initial variable, here $`A_\mu `$ or its color generalization, is to be multiplied by $`\delta A_\mu /\delta \omega ,\delta A_\mu /\delta \varphi ^a`$. The degenerate contributions noted above again occur as a result of the $`\varphi ^a`$ variations.

It is amusing to compare the Clebsch vector decomposition to the “string-embedding” partition of a 2-tensor , available in any dimension $`d`$,

$$g_{\mu \nu }(x)=y_\mu ^A(x)y_\nu ^B(x)\eta _{AB},y_\mu ^A_\mu y^A.$$
(15)
The scalar string variables $`y^A`$, with the index $`A`$ ranging over $`1,\mathrm{},\frac{1}{2}d(d+1)`$, are now the independent ones. Inserting this decomposition into the Einstein action (or indeed any other covariant one) and varying the $`y^A`$,

$$\delta I=G^{\mu \nu }\delta g_{\mu \nu }/\delta y_A\delta y_AD_\mu D_\nu (G^{\mu \nu }y^A)\delta y^A$$
$`(16\mathrm{a})`$
(where $`G^{\mu \nu }\delta I/\delta g_{\mu \nu }`$ is the necessarily identically conserved “Einstein” tensor density of the action $`I[g]`$) gives Euler-Lagrange equations very similar to the “weak Maxwell” ones, in terms of the covariant derivative $`D_\mu (g)`$

$$G^{\mu \nu }(g(y))D_\mu D_\nu y^A=0.$$
$`(16\mathrm{b})`$
These are much weaker  than333One might conjecture that a quantity similar to the Chern-Simons form becomes degenerate in the metric case, for example, that yielding the Cotton tensor (see ) in $`d=3`$. The absence of CS generalizations in even $`d`$ may account for the lack of a higher $`d`$ Clebsch representation for vectors. For odd $`d`$, the extension of (5), $`A=\mathbf{}\omega +_{i=1}^{d2}\varphi _i\mathbf{}\varphi _{i+1}`$, trivializes the corresponding CS $`AF_1\mathrm{}F_{d2}`$. $`G^{\mu \nu }=0`$. In both vector and tensor cases, the problem is with the derivative nature of a parameterization by scalars: the equally nonlinear *vielbein* decomposition of the metric $`g_{\mu \nu }=e_\mu ^ae_\nu ^b\eta _{ab}`$ into the set of vectors $`(e_\mu )^a`$ is of course perfectly acceptable (if the *vielbeins* are invertible). Also true of both the Clebsch and the tensor decompositions is the inaccessibility of weak field, linearized, theory from these nonlinear parameterizations, unlike the (invertible) *vielbein* choice, where there is a vacuum possibility $`e_\mu ^a=\delta _\mu ^a`$ about which to expand.

This work was supported in part by National Science Foundation grant PHY99-73935 and US Department of Energy (DOE) contract \# DE-FC02-94ER40818.