text:
# Measurement of the Electric Form Factor of the Neutron at 𝑄²=0.5 and 1.0 GeV2/c2

The Jefferson Lab E93-026 Collaboration

(January 30, 2004)

## Abstract

The electric form factor of the neutron was determined from measurements of the $`\stackrel{}{d}(\stackrel{}{e},e^{}n)p`$ reaction for quasielastic kinematics. Polarized electrons were scattered off a polarized deuterated ammonia (15ND3) target in which the deuteron polarization was perpendicular to the momentum transfer. The scattered electrons were detected in a magnetic spectrometer in coincidence with neutrons in a large solid angle detector. We find $`G_E^n=0.0526\pm 0.0033(stat)\pm 0.0026(sys)`$ and $`0.0454\pm 0.0054\pm 0.0037`$ at $`Q^2=0.5`$ and 1.0 (GeV/c)2, respectively.

preprint: HEP/123-qed

The electric form factor of the neutron $`G_E^n`$ is a fundamental quantity in nuclear physics. Knowledge of $`G_E^n`$ over a large range of momentum transfer leads to insight to physics beyond the simple SU(6) symmetric models, for which it must vanish. Beyond nucleon structure, our understanding of $`G_E^n`$ has an impact on a broad range of topics, which vary from the charge radii of nuclei Friar:pp  to extracting the strangeness content of the nucleon Kumar:eq .

Historically, measurements of $`G_E^n`$ have suffered from large uncertainties due to experimental technique and nuclear model dependence. Early attempts to measure $`G_E^n`$ from cross sections sometimes failed to determine even the sign. Until recently, the best determination of $`G_E^n`$ came from elastic electron-deuteron measurements, but the errors were large, $``$30%, due to their dependence on the nucleon-nucleon potential model platchkov90 . The first polarization measurements of $`G_E^n`$, conducted with 2H and 3He targets, differed significantly because final state interactions were not addressed becker99 ; ostrik99 . These experiences highlight the importance of measuring $`G_E^n`$ using different reactions.

In the last few years the experimental understanding of $`G_E^n`$ has improved considerably. The results from polarization experiments described in Ref. HG01 -bermuth03  are consistent and provide good accuracy in the kinematic region of the four-momentum transferred squared $`Q^20.7`$ (GeV/c)2 (henceforth units of $`Q^2`$ are assumed to be (GeV/c)2). Until now, values of $`G_E^n`$ determined from the deuteron quadrupole form factor SS02  provide the only information in the kinematic region above $`Q^20.7`$ and leave the $`Q^2>1.6`$ region undetermined.

This letter describes a first measurement of $`G_E^n`$ at $`Q^2=1.0`$ using a polarized target. In addition, the result for $`Q^2=0.5`$ has half the uncertainty as the previous measurement HG01  and is the most precise result near the peak of $`G_E^n`$. These results are largely insensitive to the model of the nucleon-nucleon (NN) potential so that, compared to those of Ref. SS02 , they are more reliable.

To determine $`G_E^n`$, the helicity dependent rate asymmetry in electron scattering was measured. In the ideal case of a polarized electron scattering elastically off a free polarized neutron, with the neutron polarization vector in the scattering plane and perpendicular to the momentum transfer $`\stackrel{}{q}`$, $`G_E^n`$ is related to the beam-target asymmetry term $`A_{en}^V`$ donnelly  by

$$A_{en}^V=\frac{2\sqrt{\tau (\tau +1)}\mathrm{tan}(\theta _e/2)G_E^nG_M^n}{(G_E^n)^2+\tau [1+2(1+\tau )\mathrm{tan}^2(\theta _e/2)](G_M^n)^2},$$
(1)
where $`\tau =Q^2/4M_n^2`$, $`M_n`$ is the mass of the neutron, $`G_M^n`$ is the magnetic form factor of the neutron, and $`\theta _e`$ is the electron scattering angle.

For lack of a free neutron target, the actual measurements were performed on a polarized deuterium target. The measured experimental asymmetries were due to a combination of several physics asymmetries scaled by the electron and target vector and tensor polarizations (that are described in detail in our previous work HG01 , not repeated here for brevity’s sake). The proper averaging of the asymmetry (symmetrically around $`\stackrel{}{q}`$) and the negligible contributions from the target tensor asymmetry simplify the relationship of the measured asymmetry $`ϵ`$ to the deuteron vector asymmetry $`A_{ed}^V`$ so that $`ϵ=fP_eP_1^dA_{ed}^V`$, where $`P_e`$ is the beam polarization, $`P_1^d`$ is the deuteron vector polarization and $`f`$ is the dilution factor due to scattering from nucleons other than polarized deuterons in the target. Calculations show that $`A_{ed}^V`$ depends linearly on $`G_E^n`$ for the kinematics of the experiment arenhoevel88 .

The measurements were conducted in Hall C of the Thomas Jefferson National Accelerator Facility in a setup similar to that of the previous measurement HG01 . The longitudinally polarized electron beam source  was scattered off a polarized frozen deuterated ammonia (15ND3) target. The scattered electrons were detected by the High Momentum Spectrometer (HMS), and the neutrons were detected by a dedicated neutron detector. The central kinematics for the two measurements as well as the average beam and target polarizations are listed in Table 1. The average deuteron luminosity was $`10^{35}/cm^2/s`$.

The polarized target target  consisted of ammonia granules submerged in liquid He and maintained at 1 K by a 4He evaporation refrigerator. The target spins were aligned by a 5 T magnetic field generated by a pair of superconducting coils. The polarization was enhanced via dynamic nuclear polarization Crabb94  and measured with a continuous-wave NMR system nmr . To minimize localized heating and depolarization of the target material, the electron beam was rastered uniformly in a circular pattern with a 1 cm radius.

A two-magnet chicane compensated for the deflection of the electron beam by the target field. After traversing the target, the beam passed through a helium bag to a special dump in the experimental hall. The polarization of the beam was measured at regular intervals throughout the experiment with a Møller polarimeter moller . The beam helicity was changed in a pseudorandom sequence 30 times per second to minimize charge fluctuations correlated with helicity.

The HMS was operated in its standard mode for the detection of electrons. The established reconstruction algorithms were augmented to account for the large beam rastering and the effects of the target magnetic field on the scattered electrons.

The neutron detector consisted of multiple vertical planes of segmented plastic scintillators. Two planes of thin scintillators served to distinguish charged particles. Behind these were six planes of thick scintillators to detect the neutrons. All scintillators were equipped with photomultipliers on both ends to provide spatial and time information for the detected particle. The 88 thick scintillators provided a neutron detection volume that was roughly 160 cm wide, 160 cm tall and 90 cm deep. The front of the detection volume was approximately 4.2 m and 6.2 m from the target during the $`Q^2=0.5`$ and 1.0 measurements, respectively. The detector was shielded from direct gamma rays from the target by a 2.5 cm lead curtain, and the entire assembly was housed in a thick-walled concrete hut, which was open to the target.

The trigger was set up so that the neutron detector was read out for every electron trigger in the HMS. Coincidences between the electrons and the knock-out nucleon were determined offline.

The experiment was simulated using Monte Carlo (MC) techniques. The simulation package, based on MCEEP MCEEP , included the charged particle transport through the target’s magnetic field and the optical and aperture model of the HMS from the Hall C simulation package SIMC SIMC . The MC served two principal functions: to determine the dilution factor and to average the theoretical asymmetries over the acceptance. For the dilution factor, all target materials were included in the simulation: the deuterium and the nitrogen in the ammonia, the liquid helium bath, the aluminum target cell windows and the NMR coil.

Contributions from electron-neutron events originating from quasielastic scattering and pion production were included in the MC for all target materials. It was found that the $`(e,e^{}n)\pi `$ contribution to the event sample was negligible for $`Q^2=0.5`$ and less than 0.5% for $`Q^2=1.0`$. Two-body knockout contributions were examined in the MC (following ryckebusch ) and were also negligible.

The simulated rates were normalized to the measured ammonia rates, and the variations in the ratios of MC rate to observed rate for various target materials were used to determine the uncertainty in the dilution factor. A comparison of event distributions of the data and MC for several kinematic variables for $`Q^2=1.0`$ is shown in Fig. 1; agreement is very good. The comparison for the $`Q^2=0.5`$ data is similar and slightly better than in the previous experiment HG01 .

Several corrections to the measured asymmetry were applied. The leading correction was for charge exchange of protons in the lead shielding: $`3.8\pm 1.1\%`$ ($`3.3\pm 3.0\%`$) for $`Q^2=0.5`$ ($`Q^2=1.0`$). The charge exchange in the target material itself was small: $`0.3\pm 0.3\%`$ ($`0.3\pm 0.3\%`$). The accidental background rate was $`1.9\pm 0.1\%`$ ($`0.5\pm 0.2\%`$) with no statistically significant asymmetry. Internal radiation effects on the asymmetry were also corrected: $`0.9\pm 0.5\%`$ ($`0.6\pm 0.5\%`$). The effect of external radiation was insignificant compared to the 0.5% MC statistical uncertainty. The contamination from neutral pions generated by protons in the lead shielding of the neutron detector was also found to be insignificant.

The physics model of $`A_{ed}^V`$ used in the MC was based on the calculations of Arenhövel, Leidemann and Tomusiak arenhoevel88 . It included a non-relativistic description of the $`np`$ system in the deuteron using the Bonn $`R`$-space $`NN`$ potential bonn  for both the bound state and the final state interactions. The full calculation included meson exchange, isobar configuration currents and relativistic corrections.

The model assumed a scaled Galster parameterization galster  for $`G_E^n`$ and the dipole parameterization for $`G_M^n`$. The Galster parameterization is $`G_E^n(Q^2)=\mu _n\tau /(1+5.6\tau )G_D(Q^2),`$ where $`G_D=1/(1+Q^2/0.71)^2`$ is the dipole form factor. Various scale factors of this parameterization, ranging from 0.5 to 1.5, were examined. The potential impact of the $`Q^2`$ dependence of the $`G_E^n`$ parameterization was found to be negligible because the $`Q^2`$ acceptance was not large and the $`Q^2`$ dependence of $`G_E^n`$ over the acceptance was mostly linear. The acceptance-averaged value of $`A_{ed}^V`$ has a linear dependence on $`G_M^n`$, so one can trivially correct for more accurate $`G_M^n`$ values.

The value of $`G_E^n`$ was determined by comparing the acceptance averaged $`A_{ed}^V`$ of the data to that of the MC. The theoretical asymmetries were determined for a range of scaling factors of the Galster parameterization to determine the corresponding $`G_E^n`$: $`G_E^n/\mathrm{Galster}=1.003\pm 0.064`$ and $`1.172\pm 0.140`$ for $`Q^2=0.5`$ and 1.0, respectively. The uncertainties are statistical only. To account for the slight deviations of $`G_M^n`$ from the dipole form factor, the recent fit to the $`G_M^n`$ data baselgmn  was used: $`G_M^n/\mu _nG_D=1.007\pm 0.005`$ and $`1.072\pm 0.014`$ for $`Q^2=0.5`$ and 1.0, respectively. Then the values for $`G_E^n`$ are:

$`G_E^n(Q^2=0.5)`$ $`=`$ $`0.0526\pm 0.0033\pm 0.0026,`$
$`G_E^n(Q^2=1.0)`$ $`=`$ $`0.0454\pm 0.0054\pm 0.0037,`$
where the first uncertainty is statistical and the second is systematic. The $`Q^2=0.5`$ result agrees well with the previous result reported in Ref. HG01 .

Other systematic uncertainties for $`Q^2=0.5`$ ($`Q^2=1.0`$) were: dilution factor, 3.4% (3.0%); target polarization, 2.9% (4.6%); central kinematic values, 1.2% (3.4%); beam polarization, 1.1% (3.3%).

The results as compared to recent measurements are shown in Fig. 2. The new data at $`Q^2=1.0`$ provides the important experimental confirmation of the decline of $`G_E^n`$ following the Galster form. The data shown in the figure were fit to the form $`G_E^n(Q^2)=\mu _na\tau /(1+p\tau )G_D(Q^2),`$ The parameter $`a=0.895\pm 0.039`$, the slope at $`Q^2=0`$, is fixed by atomic measurements of the neutron charge radius kopecky . The fit yielded $`p=3.69\pm 0.40`$. The errors obtained for the fit parameters are uncorrelated. The one-sigma error region of the fit is shown in Fig. 2 as the shaded band.

Many recent models lomon miller  have attempted to predict or fit the nucleon electromagnetic form factors. Fig. 2 compares the data with two recent calculations that use covariant formulations of the constituent quark model with quark-quark interactions fitted to spectroscopic data. The point-form spectator approximation (PFSA) of pfsa  (dot-dash) uses a Goldstone boson exchange interaction with pointlike constituent quarks while the light-front (LF) calculation of lightfront  (short dash) uses a one-gluon exchange interaction with constituent quark form factors fitted to data for $`Q^2<1`$. The use of constituent form factors improves the fit to the nucleon magnetic form factors at larger $`Q^2`$, but the PFSA seems to describe $`G_E^n`$ better at low $`Q^2`$ with fewer parameters.

Also shown in Fig. 2 are the results from a hybrid model that interpolates between vector-meson dominance at low $`Q^2`$ and perturbative QCD at high $`Q^2`$ lomon (solid line), from a light front model where the nucleon is considered a system of three bound quarks surrounded by a cloud of pionsmiller  (long dash) and from a soliton modelholzwarth  (dotted) whose basic features include an extended object, partial coupling to the vector mesons and relativistic recoil corrections. While all these models agree qualitatively with the data, none agree with the data for the entire range of $`Q^2`$.

It is remarkable that in the last five years, the experimental precision in $`G_E^n`$ measurements has improved to the 10% level. This significant improvement provides a rigorous challenge for models of the nucleon structure because this electromagnetic form factor is the most sensitive to physics beyond the simplistic SU(6) symmetric picture.

In conclusion, the electric form factor of the neutron at $`Q^2=0.5`$ and $`1.0`$ has been determined from measurements of the beam-target asymmetry. This experiment provides the highest $`Q^2`$ measurement to date using a polarized target and the most precise measurement near the maximum of $`G_E^n`$.

We wish to thank the Hall C technical and engineering staff at TJNAF as well as the injector, target and survey groups for their outstanding support. This work was supported by the Commonwealth of Virginia, the Schweizerische Nationalfonds, the U.S. Department of Energy, the U.S. National Science Foundation, the U.S.-Israel Binational Science Foundation and the Deutsche Forschungsgemeinshaft. The Southeastern University Research Association (SURA) operates the Thomas Jefferson National Accelerator Facility for the U.S. Department of Energy under contract DE-AC05-84ER40150.


format: text

text:
# Evidence of two-dimensional macroscopic quantum tunneling of a current-biased DC-SQUID

## Abstract

The escape probability out of the superconducting state of a hysteretic DC-SQUID has been measured at different values of the applied magnetic flux. At low temperature, the escape current and the width of the probability distribution are temperature independent but they depend on flux. Experimental results do not fit the usual one-dimensional (1D) Macroscopic Quantum Tunneling (MQT) law but are perfectly accounted for by the two-dimensional (2D) MQT behaviour as we propose here. Near zero flux, our data confirms the recent MQT observation in a DC-SQUID Li02 .

Quantum dynamics of a current-biased Josephson junction (JJ) have been extensively studied during the past two decades Leggett92 . The phase difference across the JJ obeys the same dynamics as the position of a fictitious particle moving in a one dimensional tilted washboard potential, whose average slope is proportional to the current bias. For bias below the critical current, the particle is trapped in a local minimum. It can escape out by Thermal Activation (TA) or by Macroscopic Quantum Tunneling (MQT). Many experiments have demonstrated these two escape regimes. Quantized energy levels and resonant activation have been observed by applying microwaves Leggett92 ; Esteve86 . Very recently Rabi oscillations and coherent temporal oscillations were observed in a current biased JJ Yu02 ; Martinis02  or in a Cooper pair transistor coupled to a biased JJVion02 . Experiments on RF-SQUID have also clearly demonstrated quantum dynamics in form of MQT Leggett92 , resonant tunneling or quantum superposition of states Friedman00 .

By contrast and surprisingly, very few studies have been performed in a current biased DC-SQUID in the quantum regime. The dynamics of a DC-SQUID are described by two degrees of freedom and new phenomena can arise. The fictitious particle moves in this case in a two-dimensional (2D) potential. At zero magnetic field it has been predicted that escape out of a local minimum can occur via two different trajectories Ivlev87 ; Morais94 . At non zero magnetic field, MQT was derived only in the limit of small inductance of the loop Chen86 , where the SQUID is nearly equivalent to a single junction with a flux-dependent critical current. In this case, the dynamics are one dimensional. To our knowledge only three previous experiments were reported in the past discussing escape from a current-biased DC-SQUID. In the TA regime good agreement with 2D transition state theory Brinkman57  has been obtained Seguin92 . Similar results were obtained in an RF-SQUID Han89 . In the quantum regime no experimental evidence of two-dimensional behaviour has been demonstrated up to now. Sharifi et al. Sharifi88  claimed MQT observations but neither the amplitude of quantum fluctuations nor the TA regime were understood. Recent measurements Li02  have clearly observed MQT and TA regimes in a DC-SQUID at zero magnetic field. 2D behaviour could, however, not be observed because at zero flux the escape rate of their SQUID behaves exactly as that of a single Josephson junction, ie. a 1D system. Understanding the dynamics of a hysteretic DC-SQUID is very important in interpreting some recent quantum experiments VanderWal00 ; Saito02 , or a recently proposed experiment of a one shot quantum measurement in a superconducting charge qubit using a DC-SQUID Buisson03 .

In this letter, we report on an escape measurement in a DC-SQUID both in the thermal and in the quantum regime. In contrast to till now reported results on MQT, 1D tunneling escape formula can not describe our experimental data. Since no theory has described MQT in a DC-SQUID in the 2D regime at non zero magnetic flux, we propose a method based on the zero magnetic field result derived in Ref. Ivlev87  but taking into account the flux dependence of the 2D potential. Our experimental data support quantitatively this approach.

The current-biased DC-SQUID consists of two Josephson junctions in parallel in an inductive superconducting loop. The whole system has two degrees of freedom described by $`\phi _1`$ and $`\phi _2`$, the phase difference across the two junctions. Assuming identical JJs, each having a critical current equal to $`I_0`$ and a capacitance equal to $`C_0`$, the dynamics of the system can be treated as that of a fictitious particle moving in a two dimensional potential given by Tesche77

$$U(x,y)=U_0\left(sx\mathrm{cos}x\mathrm{cos}y\eta sy+b(y\pi \mathrm{\Phi }/\mathrm{\Phi }_0)^2\right)$$
(1)
where $`U_0=\mathrm{\Phi }_0I_0/\pi `$ is the sum of the Josephson energy of the two JJs, $`s=I_{DC}/(2I_0)`$, $`x=(\phi _1+\phi _2)/2`$, $`y=(\phi _1\phi _2)/2`$, $`b=\mathrm{\Phi }_0/(2\pi LI_0)`$ and $`I_{DC}`$ is the bias current. The asymmetry of inductances is parametrized by $`\eta =(L_2L_1)/L`$ where $`L_1`$ and $`L_2`$ are the inductances of the two branches and $`L=L_1+L_2`$ is the inductance of the SQUID loop. $`\mathrm{\Phi }`$ is the applied flux and $`\mathrm{\Phi }_0=h/2e`$ is the flux quantum.

The potential presents valleys and mountains with local minima separated by saddle points along which the particle can escape. The barrier height between the minima and the saddle points depends on the bias current and magnetic flux and vanishes at the critical current $`I_c`$. $`I_c`$ follows the usual flux dependence with $`\mathrm{\Phi }_0`$ periodicity. In the following we consider $`b1`$ close to that of our DC-SQUID. In this limit, there exists only one optimal trajectory that passes through the saddle point. Therefore we will not discuss escape out of a local minimum through two different trajectories Ivlev87 ; Morais94 . Yet the 2D signature still exists because the optimal trajectory is not a straight line along the $`x`$-direction but follows a path in the landscape with the tilt in the potential.

In the limit $`ss_c`$ where $`s_c=I_c/(2I_0)`$ is the normalized effective critical current of the DC-SQUID, we can expand the potential (Eq. (1)) to the third order. The escape occurs in a current and flux dependent angle $`\theta `$ with respect to the $`x`$-axis. Along this direction, we obtain the potential barrier height, $`\mathrm{\Delta }U=\mathrm{\Delta }U_0\sqrt{s_c\mathrm{cos}^3\theta /u}(1+\eta \mathrm{tan}\theta )^{3/2}`$, and the plasma frequency, $`\omega _p=\omega _{p0}(u\mathrm{cos}\theta /s_c)^{1/4}(1+\eta \mathrm{tan}\theta )^{1/4}`$ where $`u`$ is the third derivative of the potential along the escape direction Seguin92 . The quantities $`\omega _{p0}=2^{1/4}(\pi I_c/\mathrm{\Phi }_0C_0)^{1/2}(1s/s_c)^{1/4}`$ and $`\mathrm{\Delta }U_0=(4\sqrt{2}/3)(\mathrm{\Phi }_0/2\pi )I_c(1s/s_c)^{3/2}`$ are the usual plasma frequency and potential height, respectively, of a single JJ with an effective critical current $`I_c`$ and a capacitance $`2C_0`$. The 2D character of the escape is exhibited in these expressions by the renormalization of the plasma frequency and the potential barrier height of a single JJ. Therefore the ratios $`\omega _p/\omega _{p0}`$ and $`\mathrm{\Delta }U/\mathrm{\Delta }U_0`$ show deviations from 1D behavior.

Using the previous results of the renormalized plasma frequency and potential height, the escape rate out of a 2D potential well, in the TA regime, $`\mathrm{\Gamma }^{TA}`$, is predicted to be Brinkman57

$$\mathrm{\Gamma }^{TA}=\frac{\omega _p}{2\pi }\frac{\omega _w}{\omega _s}\mathrm{exp}\left(\frac{\mathrm{\Delta }U}{k_BT}\right)$$
(2)
where $`\omega _w`$ and $`\omega _s`$ are respectively the transverse frequency modes in the well and at the saddle point Seguin92 . We numerically analyzed the modulation of the escape current as a function of flux in the TA regime taking into account the transverse frequency modes in the escape rate. For our SQUID parameters, the deviation of $`\omega _w/\omega _s`$ from unity has been estimated to be $`2\%`$ at zero flux and it quickly decreases to zero when increasing $`\mathrm{\Phi }/\mathrm{\Phi }_0`$. Since this small error appears only in the prefactor of Eq. (2), we will assume in the following that $`\omega _w/\omega _s=1`$.

At low temperature, the escape is dominated by quantum tunneling. This regime has been theoretically studied at zero flux in Ivlev87 ; Morais94  discussing the splitting of the two escape trajectories, or flux dependence in the low inductance limit ($`b>>1`$) Chen86 . There is so far, to our knowledge, no theoretical treatment of flux dependence of 2D MQT, ie. in the regime $`b1`$. Therefore we propose to use the MQT escape rate at zero magnetic flux Ivlev87 ; Li02  in the limit $`b1`$. We introduce in this escape rate, $`\mathrm{\Gamma }^{MQT}`$, the 2D behaviour by taking the plasma frequency $`\omega _p`$ and the potential barrier height $`\mathrm{\Delta }U`$ along the escape direction, analogously to what was done in the TA regime earlier Seguin92 :

$$\mathrm{\Gamma }^{MQT}=f_{2D}\frac{\omega _p}{2\pi }\sqrt{864\pi \mathrm{\Delta }U/\mathrm{}\omega _p}\mathrm{exp}(36\mathrm{\Delta }U/5\mathrm{}\omega _p).$$
(3)
This 2D MQT rate thus depends on the magnetic flux and on the bias-current. In the following, because $`b1`$, two escape trajectories are impossible and $`f_{2D}1`$. Note, however, that the 2D behaviour remains through $`\mathrm{\Delta }U`$ and $`\omega _p`$. At zero magnetic flux $`\mathrm{\Gamma }^{MQT}`$ equals that of a single JJ of critical current $`2I_0`$ and capacitance $`2C_0`$.

The escape probability $`P_{esc}`$ from the zero voltage state during time $`\mathrm{\Delta }t`$ is a function of the amplitude of the bias current, the flux and temperature. Measured escape probability using current DC-pulses of duration $`\mathrm{\Delta }t=50`$ $`\mu `$s is plotted in Fig. 3. Each point corresponds to 4000 current pulses. By increasing the bias current, $`P_{esc}`$ starts from zero far below the critical current where the barrier is high, and it becomes unity at currents close to and above the critical one. In the following the escape current $`I_{esc}`$ is defined as the bias current corresponding to $`P_{esc}=0.5`$ and the width of the escape probability distribution $`\mathrm{\Delta }I`$ as the difference between the currents at $`P_{esc}=0.9`$ and $`P_{esc}=0.1`$. The measured escape rate is extracted from $`\mathrm{\Gamma }(I_{DC})=\mathrm{ln}(1P_{esc}(I_{DC}))/\mathrm{\Delta }t`$ and it ranges from about 10 s-1 to 105 s-1 (inset of Fig. 2).

The sample is a DC-SQUID consisting of two $`Al/AlO_x/Al`$ tunnel junctions obtained by e-beam lithography and shadow evaporation (inset of Fig. 1). The area of each JJ is about 11 $`\mu `$m2. The tunnel resistance of the SQUID, $`R_T80`$ $`\mathrm{\Omega }`$, and the superconducting gap $`\mathrm{\Delta }=230`$ $`\mu `$eV are deduced from the dissipative branches of the IV characteristic and yield $`I_02.3`$ $`\mu `$A using the Ambegaokar-Baratoff relation assuming two identical JJs. From MQT of a single JJ with the same area, we deduced $`C_0=0.46`$ pF, which agrees with the typical specific capacitance of about 50 fF/($`\mu `$m)2 for an $`AlOx`$ tunnel junction. The rectangular SQUID loop is about 10 $`\mu `$m by 40 $`\mu `$m giving an estimated magnetic inductance of about 70 pH. The kinetic inductance of the SQUID, deduced from the thin film resistivity, $`\rho 5.3\mu \mathrm{\Omega }`$cm, is coarsely estimated to be about 120 pH. The total inductance $`L`$, the sum of the kinetic and magnetic inductances, is therefore about 190 pH. An inductance asymmetry exists because the two SQUID branches have different lengths. We estimate $`\eta 0.2`$ by assuming proportionality between length of a branch and inductance. The immediate electromagnetic environment of the DC-SQUID is determined by two 20 nm thin, 1 $`\mu `$m wide and 700 $`\mu `$m long superconducting aluminium wires terminated by two large pads for wire bonding. The estimated total inductance of the in-situ wires, $`L_e=4.6`$ nH, results again of the sum of 1.2 nH magnetic inductance of the coplanar strips and of 3.4 nH kinetic inductance of the wires. Although these SQUID parameters have been independently estimated, the accuracy in determining them this way is not sufficient. Therefore we extracted the inductance of the DC-SQUID $`L=244`$ pH and its asymmetry $`\eta =0.26`$, the zero-flux, the critical current $`I_0=2.33`$ $`\mu `$A and the capacitance $`C_0=0.46`$ pF from the flux modulation characteristics of the escape current in the low temperature regime (see Fig. 3). These parameters are in very good agreement with the rough estimates given above. Our $`b=0.6`$ is much smaller than $`b3`$ of the earlier works Sharifi88 ; Li02  on MQT in a DC-SQUID, indicating that the fictitious particle is much more affected by the two-dimentionality of the potential.

The sample is enclosed in a cavity with resonances above 15 GHz, and it is anchored to the mixing chamber of the dilution refrigerator. This cavity is enclosed in a low temperature filter Balestro03  consisting of thermocoax microwave filters Zorin95  and $`\mathrm{\Pi }`$-filters. The attenuation of the DC-lines is estimated to be at least 200 dB above 1 GHz at the mixing chamber temperature. Thermocoax filters were also introduced from 1.5 K to the low temperature filter, and low pass LC filters were used at room temperature. Current bias and voltage probe lines of the sample are separated from the computer using differential amplifiers. A superconducting shield was inserted inside the vacuum jacket of the cryostat and $`\mu `$-metal surrounds the dewar to protect the DC-SQUID from external flux noise. Using this measurement set-up, we separately measured a width of $`\mathrm{\Delta }I=11`$ nA on a 1 $`\mu `$A Josephson junction in its MQT regime.

Figure 2 shows escape probability and escape rate at $`\mathrm{\Phi }/\mathrm{\Phi }_0=0.07`$ and $`T=42`$ mK. At this particular flux value, the escape current is maximum (see Fig. 3) and the ratios $`\omega _p/\omega _{p0}`$ and $`\mathrm{\Delta }U/\mathrm{\Delta }U_0`$ are equal to unity. Therefore the SQUID behaves as a single JJ. We can perfectly fit our data using 1D-MQT predictions (continuous line in Fig. 2) without free parameters which demonstrates that 1D-MQT theory describes the behaviour of a DC-SQUID at the maximum escape current. This first result confirms the recent MQT observation at zero flux performed by Li $`et`$ $`al.`$ Li02 . It does not, however, demonstrate 2D escape behaviour.

The escape current and the width versus applied external flux are directly extracted from the measured switching current distribution and they are plotted in Fig. 3 in the MQT regime at T= 42mK . The width is periodic in flux with periodicity $`\mathrm{\Phi }_0`$. It is minimum near $`\mathrm{\Phi }/\mathrm{\Phi }_0=0.07`$ and increases to reach maximum at $`\mathrm{\Phi }/\mathrm{\Phi }_0=\pm 0.5`$. There is asymmetry between negative and positive flux values and discountinuities at $`\mathrm{\Phi }/\mathrm{\Phi }_0=\pm 0.5`$. The usual 1D tunneling model cannot account for the main feature of our results. Namely the 1D model predicts a decrease of the width, $`\mathrm{\Delta }I(2I_0s_c)^{3/5}`$ in the MQT regime (dashed line in Fig. 3) which is contrary to our observation. Our proposed 2D formula in the MQT regime suggests increase of the width and it agrees with our data perfectly also quantitatively. It has no free parameters since all of them were already fixed by the escape current characteristics. We claim that these results are the first experimental observation of 2D escape from the zero voltage state of a DC-SQUID in its MQT regime.

The escape current and width against temperature are shown in Fig. 4 for $`\mathrm{\Phi }/\mathrm{\Phi }_0=0.25`$. The escape current increases as temperature drops and saturates at low temperature. The width decreases with temperature and saturates at low temperature. These features are perfectly accounted for by Eq. (3) in the MQT regime and by Eq. (2) in the TA regime (dashed lines). Yet the model seems to underestimate the measured width in the thermal regime by about 10 %, indicating that influence of the environment must be taken into account to better fit with the data. Finally, the observed cross-over temperature ($`80`$ mK) between MQT and TA is consistent with that predicted by the well-known formula $`\mathrm{}\omega _p/(2\pi k_B)`$ Leggett92 .

In order to estimate the effect of the environment on the escape process, we take into account the quadratic current fluctuations $`\delta I^2`$. Since the SQUID is connected to the external circuit through inductive wires on the chip, we assume the electrical environment to be the inductance $`L_e`$ in series with an external resistor $`R_{ext}`$. For temperatures $`k_BT>\mathrm{}R_{ext}/L_e`$, the current fluctuations are given by $`\delta I^2=k_BT/L_e`$. For $`k_BT<\mathrm{}\omega _p`$ thermal fluctuations are slow as compared to the SQUID dynamics, and escape rate can be corrected by $`\mathrm{\Gamma }_{env}^{TA}=\mathrm{\Gamma }^{TA}\mathrm{exp}((\mathrm{\Delta }U/I_{DC})\delta I^2/2(k_BT)^2)`$ in the TA regime and $`\mathrm{\Gamma }_{env}^{MQT}=\mathrm{\Gamma }^{MQT}\mathrm{exp}(((\mathrm{\Delta }U/\mathrm{}\omega _p)/I_{DC})\delta I^2/2)`$ in the MQT regime following the idea in Martinis88 . Using the estimated wires inductance, our data perfectly are fitted both in the quantum and in the TA regime by the fluctuation corrected expressions above (continuous line in Fig. 4). We notice that like in the recent work in Li02  the measured $`\mathrm{\Delta }I`$ crosses over from MQT to TA much more gradually than what we predict.

In summary, we have shown for the first time evidence of two dimensional macroscopic quantum escape of a hysteretic DC-SQUID from its zero voltage state. Moreover, we have proposed a way to theoretically analyze the flux dependence of escape in the MQT regime, and we show that our data perfectly fit this model, also in the TA regime with the same values of the circuit parameters.

We thank K. Hasselbach, F. Hekking, Ph. Lafarge, L. Lévy and A. Niskanen for useful discussions. This work was supported by the French ACI program. JP acknowledges support from CNRS and Joseph Fourier University.


text:
# Non-Gaussianity of Large-Scale Cosmological Perturbations

## I Introduction

Inflationary models inflation  explain the origin of cosmological perturbations, and predict the primordial super-horizon non-Gaussianity which is small in the standard one-field slow-roll case maldacena . A simple and natural extension of the standard model, however, predicts a much larger non-Gaussianity dgz ; zal . To measure the primordial non-Gaussianity one needs to compute the current sub-horizon perturbations in terms of primordial super-horizon perturbations up to second order in the amplitude. This is done here for the universe filled with matter and cosmological constant. The pure matter case has been computed by gks ; in that paper we also discuss observations of non-Gaussian effects in galaxy surveys.

Neglecting tensor (gravitational) perturbations, we assume that primordial perturbations are pure scalar. Following maldacena ; bst ; sb , we use a gauge invariant variable $`\zeta `$ that remains constant outside the horizon to characterize the primordial perturbations. On super-horizon scales $`\zeta `$ is defined as follows: $`e^\zeta `$ is proportional to the local scale factor measured on uniform local Hubble parameter hypersurfaces.

The degree of non-Gaussianity of primordial perturbations $`f_{nl}`$ can be defined by assuming that $`\zeta `$ is of the form komsper

$$\zeta =g(3/5)f_{nl}g^2,$$
(1)
where $`g`$ is Gaussian. The standard one-field slow-roll inflation predicts $`f_{nl}0.05`$ maldacena ; spergel . Other inflationary models dgz ; zal  predict $`|f_{nl}|`$ few. The WMAP results are $`58<f_{nl}<134`$ at 95% confidence komatsu . The Sloan Digital Sky Survey should provide similar accuracy, and future galaxy surves can be much better roman .

The non-Gaussianity of sub-horizon density perturbations is due to Newtonian evolution, relativistic effects, and primordial non-Gaussianity. As we will see, the ratio of relativistic to primordial contributions is $`1/f_{nl}`$; accurate determination of $`f_{nl}`$ must include relativistic effects.

In §2 we derive the second-order expression for the density perturbation in terms of $`\zeta `$. In §3 we generalize the standard Newtonian calculation of the density skewness.

## II Second-order density perturbation in terms of $`\zeta `$.

### II.1 Results

Consider flat matter plus cosmological constant universe, with no primordial gravitational waves. We use the comoving synchronous gauge. The first-order metric is

$$ds^2=dt^2a_0^2(t)(1+2\zeta )(\delta _{ij}+2\tau \zeta _{ij})dx_idx_j.$$
(2)
Here $`a_0(t)`$ is the unperturbed scale factor; $`\zeta (x)`$ is the primordial perturbation, which is a function of $`x`$ only; $`\zeta _{ij}`$ indicates spatial derivatives; $`\tau (t)`$ is the linear growth factor given in the Appendix. For zero cosmological constant, $`\mathrm{\Lambda }=0`$, we have $`\tau t^{2/3}`$ – the familiar linear growth of perturbations in the matter dominated universe.

The density perturbation to second order is

$$\delta =\tau \zeta _{ii}+\frac{\tau ^2}{2}\left((1+𝒞)\zeta _{ii}^2+(1𝒞)\zeta _{ij}^2\right)+\frac{\tau }{2}\left(4\zeta \zeta _{ii}\zeta _i^2\right).$$
(3)
Here $`𝒞(t)`$ is a dimensionless numerical factor given in the Appendix. For $`\mathrm{\Lambda }=0`$, $`𝒞=3/7`$ and we recover gks .

The first two terms in (3) represent the second-order Newtonian perturbation in the synchronous coordinate system peebles ; matarrese . The third term represents the relativistic effects. The first-order metric can be derived along standard lines mukh . The second order density perturbation is obtained in the next subsection. The applications of this expression are discussed in §3 and §4 .

### II.2 Derivation

We calculate the second order density perturbation under spherical symmetry. Then we simply generalize to the non-spherical case.

Consider the well known exact solution of the Einstein equations – the spherically symmetrical dust dingle ; peebles  with cosmological constant. The metric is

$$ds^2=dt^2\frac{(ax)^2}{1x^2/R^2}dx^2a^2x^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2).$$
(4)
The matter density is

$$\rho =\frac{3(x^3F)^{}}{8\pi G(x^3a^3)^{}}.$$
(5)
Here $`R(x)`$ and $`F(x)`$ are arbitrary functions of $`x`$ but not of $`t`$, while $`a(x,t)`$ and $`\rho (x,t)`$ are functions of both $`x`$ and $`t`$, and $`a`$ satisfies the Friedmann equation

$$a\dot{a}^2+a/R^2=F+\mathrm{\Lambda }a^3/3.$$
(6)
Dot and prime are time and $`x`$ derivatives.

We can choose

$$Fe^{3\zeta },$$
(7)
then at $`t0`$ equation (6) gives $`ae^\zeta t^{2/3}`$. If we also choose

$$R^2=\frac{2}{x}\zeta ^{}\zeta ^2,$$
(8)
equation (4) gives the early time (superhorizon) metric of the form

$$ds^2=dt^2a_0^2(t)e^{2\zeta }d𝐱^2,$$
(9)
where $`a_0(t)t^{2/3}`$, $`t0`$, is the unperturbed scale factor. We see that what we have called $`\zeta `$ is indeed the standard gauge invariant quantity that inflation theories like maldacena ; dgz  claim to predict to second or arbitrary order.

Having expressed the two free parameters of the spherical model, $`R(x)`$ and $`F(x)`$, in terms of $`\zeta (x)`$, we can now calculate the spherical density perturbations in terms of $`\zeta `$ to second order. The dimensionless density perturbation $`\delta `$ is defined by

$$\rho (1+\delta )\rho _0,$$
(10)
$`\rho _0`$ is the unperturbed matter density. It is convenient to define $`\delta _a`$ by

$$ae^\zeta a_0(t)(1+\delta _a).$$
(11)
Then equation (5) gives

$$\delta =3\delta _ax\delta _a^{}+6\delta _a^2+4x\delta _a\delta _a^{}+x^2\delta _a^2+x^2\zeta ^{}\delta _a^{},$$
(12)
and equation (6) gives

$$\delta _a=\delta _{a1}(1𝒞\delta _{a1}2\zeta +\frac{1}{2}x\zeta ^{}),$$
(13)
where

$$\delta _{a1}=\tau \frac{1}{x}\zeta ^{}.$$
(14)
Finally we put (13) into (12) and get the spherically symmetrical density perturbation to second order

$$\delta =\tau (\zeta ^{\prime \prime }+\frac{2}{x}\zeta ^{})+\tau ^2(\zeta ^{\prime \prime 2}+\frac{2+2𝒞}{x}\zeta ^{}\zeta ^{\prime \prime }+\frac{3+𝒞}{x^2}\zeta ^2)+\frac{\tau }{2}(4\zeta \zeta ^{\prime \prime }+\frac{8}{x}\zeta \zeta ^{}\zeta ^2).$$
(15)
To get the corresponding expression without the spherical symmetry assumption, we write down the most general expression for a scalar which contains the same powers of $`\zeta `$, $`\tau `$ and $``$ as the expression (15):

$$\delta =\tau C_1\zeta _{ii}+\tau ^2(C_2\zeta _{ii}^2+C_3\zeta _{ij}^2)+\tau (C_4\zeta \zeta _{ii}+C_5\zeta _i^2).$$
(16)
Now assume that $`\zeta `$ is spherically symmetrical. Then (16) should reproduce (15). We have 5 arbitrary coefficients in (16) in terms of which we should get the 8 known coefficients of (15). If we could do it, it would be a check of our spherically symmetrical expression. We can, and the result is given by (3).

## III The skewness of the density perturbation

The skewness of the density perturbation $`S<\delta ^3>/<\delta ^2>^2`$. The computation should include the Jacobian of (2), $`J=1+3\zeta +\tau \zeta _{ii}`$. To lowest order,

$$S=(4+2𝒞)\frac{1}{\tau k^2}\frac{36}{5}(f_{nl}+\frac{35}{24})$$
(17)
Here $`k^2<\zeta _{ii}^2>/<\zeta _i^2>`$. For $`\mathrm{\Lambda }=0`$, neglecting primordial and relativistic contributions of the second term, we obtain the Newtonian result peebles , $`4+2𝒞=34/7`$.

The actual perturbation spectrum is scale-free, and our formula is not directly applicable. However it may serve to compare primordial and relativistic contributions to the sub-horizon non-Gaussianity.

## IV Conclusions

We have expressed the large-scale cosmological perturbations in terms of the primordial scalar perturbations up to second order. Our expression may allow to measure primordial non-Gaussianities using galaxy surveys roman .

I thank Raisa Karasik and Roman Scoccimarro for useful discussions, and P.J.E. Peebles for telling me about the spherical perturbation case. This work was supported by the David and Lucile Packard Foundation.

## Appendix A  $`\tau (t)`$ and $`𝒞(t)`$

The linear growth rate in the synchronous comoving gauge coincides with the Newtonian expression for perturbations in a universe with matter and cosmological constant. We just need to re-derive it by perturbing eq. (6) according to (11). To first order we get

$$\dot{a}_0\dot{\delta }_{a1}+4\pi G\rho _0a_0\delta _{a1}=\frac{1}{a_0}\frac{\zeta ^{}}{x}.$$
(18)
We write the solution in the form (14). Then

$$\tau (t)D(t)_0^t\frac{dt_1}{D(t_1)a_0(t_1)\dot{a_0}(t_1)},$$
(19)
where $`D(t)`$ is another growth factor (the homogeneous solution of (18)):

$$\mathrm{ln}D(t)=4\pi G^t\frac{dt_1\rho _0(t_1)a_0(t_1)}{\dot{a_0}(t_1)}.$$
(20)
These equation coincide, after simple rearrangements, with the expressions given in peebles . On the other hand, the cosmological perturbation theory mukh  gives the growth factor of perturbations in the comoving synchronous gauge

$$\dot{\tau }=a_0^3𝑑ta_0,$$
(21)
which is equal to (19 ).

At early time the universe is matter dominated, $`\rho _0=1/(6\pi Gt^2)`$ and $`a_0t^{2/3}`$. Then eq.(20) gives $`D=t^1`$, and eq.(19) gives the familiar growth factor $`\tau =\eta ^2/10t^{2/3}`$, where $`\eta `$ is the conformal time. At a later time, when the cosmological constant accelerates the expansion of the universe, the growth slows down.

When equation (6) is solved to second order, a non-linear growth factor will appear. The second-order solution can be written in the form (13, 14), with a dimensionless coefficient

$$𝒞(t)=\frac{D(t)}{\tau ^2(t)}_0^t𝑑t_1\frac{2\tau +a_0^2\dot{\tau }^2/2+2a_0\dot{a_0}\dot{D}\tau ^2/Da_0\dot{a_0}\tau \dot{\tau }}{Da_0\dot{a_0}},$$
(22)
with all quantities under the integral being at time $`t_1`$. At early time, $`𝒞=3/7`$ and we recover the zero cosmological constant result gks .


text:
# Entropy Production at RHIC

## Abstract

For central heavy ion collisions at the RHIC energy, the entropy per unit rapidity $`dS/dy`$ at freeze-out is extracted with minimal model dependence from available experimental measurements of particle yields, spectra, and source sizes estimated from two-particle interferometry. The extracted entropy rapidity density is consistent with lattice gauge theory results for a thermalized quark-gluon plasma with an energy density estimated from transverse energy production at RHIC.

Lattice quantum chromodynamics (QCD) calculations indicate lattice  that at zero net baryon density at a critical temperature of $`T_c170`$ MeV and energy density $`ϵ_c1`$ GeV/fm3 a color deconfined and a chirally restored quark-gluon plasma (QGP) phase is formed. These energy densities are attained in relativistic heavy ion collisions, where it is believed that a dense system of quarks and gluons is created which undergoes rapid collective expansion before the partons hadronize and eventually decouple.

The ongoing experimental program at RHIC is devoted to the search for the signals of QGP. Unfortunately, measurements are confined to final-state particles which are mostly hadrons. One interesting quantity that may provide valuable insight into the state of the matter in the early stages of the collision is the net entropy which is roughly conserved between initial thermalization and freeze-out Siemens ; Letes . After freeze-out, when particles are free-streaming, the entropy is fixed. Entropy would also be conserved during the expansion stage if the mean free paths were very short. Viscous effects, shock waves, and the decoupling process can only increase the entropy. Thus, the entropy from the final state provides an upper bound for the entropy of the initial state.

Compared to a pion gas at the same energy density, the QGP is a high-entropy state as can be understood by considering the equation,

$$\mathrm{\Delta }S=\frac{Q}{T},$$
(1)
where $`\mathrm{\Delta }S`$ is the change in entropy resulting from the addition of an amount of heat $`Q`$ into a fixed volume at temperature $`T`$. Lattice calculations show that the transition is nearly a first order with a distinct peak in the specific heat near $`T_c`$. In this region, the temperature stays relatively fixed while the energy density changes substantially. Thus, the existence of a latent heat keeps the temperature lower than a pion gas at the same energy density, which results in a higher entropy as can be seen from Eq. (1). At high energy densities, the temperature of a QGP remains lower than that of the pion gas due to the disparity in the number of degrees of freedom. Asymptotically, the energy density of a QGP approaches that of a massless gas of quarks and hadrons,

$$ϵ=\left(N_B+\frac{7}{8}N_F\right)\frac{\pi ^2}{30}T^4,$$
(2)
where $`N_B`$ and $`N_F`$ are the number of Bosonic and Fermionic degrees of freedom. Whereas a pion gas has 3 bosonic degrees of freedom, a QGP has 47.5 effective degrees of freedom if one considers gluons and the three light quark flavors. Thus, relative to a pion gas, the QGP has a lower temperature and higher entropy. If the partonic degrees of freedom were not fully liberated at RHIC, one would thus expect a lower entropy for the same energy density. Furthermore, if one were to consider a full hadron gas, incorporating all the resonances from the particle data book at energy densities far above the point where the hadronic models are justified, one would expect to have a state of matter with even more entropy than a QGP for temperatures near or above 250 MeV.

The entropy density per unity rapidity $`dS/dy`$ can be calculated from the final phase space density, $`f(𝐩,𝐫,t_f)`$, evaluated at the freeze-out time $`t_f`$. Spectra provide information only about the final momenta, $`𝐩`$. The spatial information in $`f`$ must be determined by correlations which can provide the characteristic dimensions of $`f(𝐩,𝐫,t)`$ for given values of $`𝐩`$ Pratt ; Bertsch2 . In fact, it has been suggested that the softer equation of state associated with a first order phase transition to QGP with a large latent heat results in long lived and/or large source sizes Bertsch ; Pratt2 ; Rischke ; Soff ; teaneyhydro . The fact that the QGP models tend to accurately reproduce spectra while over-predicting source sizes Soff ; teaneyhydro ; heinzkolb ; kolbheinz ; kolbheinz2  suggests that the entropy at RHIC might fall below the expectations for a quark-gluon plasma.

In this letter we extract $`dS/dy`$ for hadrons at freeze-out from yields and spectra combined with the source dimensions in central Au+Au collisions at the RHIC energy of $`\sqrt{s}=130A`$ GeV. By comparing this estimated net hadron $`dS/dy`$ with the entropy at the early stage in the parton medium, we address several key issues concerning the initial stage of the collisions and the equilibration of partonic degrees of freedom.

In the Landau quasi-particle approximation the entropy is given by

$`S`$ $`=`$ $`(2J+1){\displaystyle \frac{d^3rd^3p}{(2\pi )^3}\left[f\mathrm{ln}f\pm (1\pm f)\mathrm{ln}(1\pm f)\right]}`$ (3)
$``$ $`(2J+1){\displaystyle \frac{d^3rd^3p}{(2\pi )^3}\left[f\mathrm{ln}f+f\pm f^2/2\right]},`$
where the $`\pm `$ sign corresponds to bosons/fermions with total spin $`J`$. In the classical limit $`f1`$, the $`f^2`$ and higher orders of $`f`$ can be neglected. The classical limit is sufficient for all particles except for pions. Even for pions, the $`f^2`$ term contributes only on the 3% level. Assuming a three-dimensional Gaussian form for the density profile in coordinate space, the single-particle phase space density at any time $`t`$ at or after freeze-out Bertsch ; Tomasik  is

$$f(𝐫,𝐩,t)=f_{\mathrm{max}}(𝐩)\mathrm{exp}\left(\frac{x^2}{2R_x^2}\frac{y^2}{2R_y^2}\frac{z^2}{2R_z^2}\right),$$
(4)
where $`R_x`$, $`R_y`$ and $`R_z`$ are implicit functions of $`𝐩`$. Since spectra are obtained from integrating the phase space density, $`f_{\mathrm{max}}`$, and therefore $`f(𝐩,𝐫,t)`$, can be determined from the radii and spectra.

$`{\displaystyle \frac{dN}{d^3p}}`$ $`=`$ $`(2J+1){\displaystyle \frac{d^3r}{(2\pi )^3}f(𝐫,𝐩)},`$
$`f_{\mathrm{max}}(𝐩)`$ $`=`$ $`{\displaystyle \frac{(2\pi )^{3/2}}{(2J+1)}}{\displaystyle \frac{dN}{d^3p}}{\displaystyle \frac{1}{R^3}},`$ (5)
where $`R^3`$ is the product of the three radii. Chemical equilibration has not been invoked as the expressions do not involve temperature or chemical potentials. Instead the phase space densities and therefore the entropies will be functions of the spectra and the source sizes which can be inferred from experiment.

In Hanbury-Brown Twiss (HBT) correlation measurements for identical particles, if the single-particle emission source is characterized by a Gaussian as in Eq. (4), the Gaussian widths are each inversely proportional to the source dimensions in the canonically conjugate spatial variables Pratt3 . One then obtains the two-particle correlation function in terms of the radii,

$`C(𝐊,𝐪)`$ $`=`$ $`{\displaystyle \frac{P(𝐩_1,𝐩_2)}{P(𝐩_1)P(𝐩_2)}}`$ (6)
$`=`$ $`1+\mathrm{exp}\left(q_o^2R_o^2q_s^2R_s^2q_l^2R_l^2\right),`$
where $`𝐊=(𝐩_1+𝐩_2)/2`$ and $`𝐪=𝐩_1𝐩_2`$. The momentum difference $`q`$ has been decomposed into “out-side-long” components defined by the beam direction $`q_l`$, parallel to the total momentum of the pair in the transverse plane $`q_o`$, and the direction orthogonal to the above two direction $`q_s`$ Bertsch ; Pratt2 . Data are generally analyzed in the longitudinal co-moving system (LCMS) frame in which the longitudinal component of the pair momentum vanishes.

The Gaussian widths of Eq. (6) in the LCMS frame are related to the radii in the two-particle rest frame via a contraction factor $`\gamma =E_T/m`$. Phase space densities depend on the product of the three radii. In the rest frame of the pair, this product is referred to as $`R_{\mathrm{inv}}^3=\gamma R_oR_sR_l`$. Substituting the HBT radii in Eq. (Entropy Production at RHIC), the entropy rapidity density from Eq. (3) in the classical limit is then

$`{\displaystyle \frac{dS}{dy}}`$ $`=`$ $`{\displaystyle d^2p_TE\frac{dN}{d^3p}\left[\frac{5}{2}\mathrm{ln}(𝐩)\pm \frac{1}{2^{5/2}}(𝐩)\right]},`$
$`(𝐩)`$ $`=`$ $`{\displaystyle \frac{1}{m}}{\displaystyle \frac{(2\pi )^{3/2}}{(2J+1)}}E{\displaystyle \frac{dN}{d^3p}}{\displaystyle \frac{1}{R_{\mathrm{inv}}^3}}.`$ (7)
Thus, $`dS/dy`$ of the final state can be estimated from the Gaussian source-size parameters and spectra.

The two-pion correlation function has been measured by STAR 2pis  and PHENIX 2pip  for $`11\%`$ central Au+Au collisions at $`\sqrt{s}=130A`$ GeV at midrapidity. In spite of about 70% increase in pion multiplicity from SPS to the RHIC energy the extracted radius, $`R_{\mathrm{inv}}6.3`$ fm at RHIC, was very similar to what was measured at the SPS na44 . Moreover, the prediction of $`R_o`$ larger than $`R_s`$ as a QGP signature was not observed. Instead $`R_o/R_s1`$ was measured over a broad range of transverse momenta $`K_T=(𝐩_{T_1}+𝐩_{T_2})/2`$. The smaller values of $`R_l`$ and $`R_o`$, as compared to hydrodynamic models, may be suggestive of a rather high final phase space density Bertsch2 ; Ray , and therefore, a low entropy. In our estimate of $`\pi ^{}`$ entropy we use the scaling relation $`R_l=A_l/\sqrt{m_T}`$ where $`m_T=\sqrt{K_T^2+m_\pi ^2}`$ with the fitted value $`A_l=3.32\pm 0.03`$ fm GeV1/2 2pip . While for the $`s`$-$`o`$ direction we employ the blast-wave inspired parameterization, $`R_s^2=R_o^2=R_{\mathrm{geom}}^2/(1+\beta ^2m_T/T)`$, where $`R_{\mathrm{geom}}=6.7\pm 0.2`$ fm, $`\beta =0.69`$ and the freeze-out temperature is $`T=125`$ MeV 2pip . The measured corresponding single particle $`\pi ^{}`$ spectra at RHIC spectrap  (after correcting for electromagnetic and weak decays) is fitted to a Bose-Einstein distribution $`EdN/d^3p=A^{}/[\mathrm{exp}(m_T/T_{\mathrm{eff}})1]`$ where $`A^{}=739`$ GeV-2 and effective (slope) temperature $`T_{\mathrm{eff}}=197`$ MeV. The entropy per unit rapidity estimated from Eq. (Entropy Production at RHIC) is $`dS/dy|_\pi ^{}=680`$, and is shown in Table I.

Pions resulting from electromagnetic decays, i.e., those from $`\eta `$ and $`\eta ^{}`$ decays, should be subtracted from the spectra. From a thermal model with a chemical equilibrium temperature of 170 MeV, we estimate that 12% of the pions come from such decays and correct the spectra accordingly. We add in the entropy of the $`\eta `$ and $`\eta ^{}`$ by assuming that the $`\eta `$ contributes the same amount of entropy as a single species of kaon which has a similar mass, and that the $`\eta ^{}`$ contributes half the amount of entropy as the average of protons and anti-protons. The factor of 1/2 is applied to account for the lower spin degeneracy of the $`\eta ^{}`$. The entropies of the $`\eta `$ and $`\eta ^{}`$ are displayed in Table I.

For a dilute gas in chemical equilibrium, it was demonstrated Siemens  that the entropy per nucleon carried by the nucleons is related to the ratio of the deuteron to the proton yield,

$$S=5/2\mathrm{ln}R_{dp}/3\sqrt{2}.$$
(8)
If, within a given phase space cell, the ratio of a species $`a`$ to the product of two species $`b`$ and $`c`$, which can undergo a reaction $`a+bc`$, is determined by the ratio of available phase space for bound and unbound states, one can determine the Gaussian source parameters from coalescence arguments prattllope ,

$`R_{\mathrm{inv}}^3`$ $`=`$ $`\pi ^{3/2}{\displaystyle \frac{(2J_c+1)}{(2J_a+1)(2J_b+1)}}`$ (9)
$`\times ({\displaystyle \frac{E_a}{m_a}}{\displaystyle \frac{dN_a}{d^3p_a}})({\displaystyle \frac{E_b}{m_b}}{\displaystyle \frac{dN_b}{d^3p_b}})/\left({\displaystyle \frac{E_c}{m_c}}{\displaystyle \frac{dN_c}{d^3p_c}}\right).`$
This expression is based on the assumption that the relative populations within a given phase space cell are determined by the ratios of the available phase space for bound and unbound states; whereas Eq. (8) requires the assumption that all phase space cells are equilibrated according to a single chemical potential. Once $`R_{\mathrm{inv}}`$ is obtained as a function of the momentum, one can extract the entropy from Eq. (Entropy Production at RHIC).

Assuming that $`\overline{p}`$ and $`\overline{n}`$ possess identical phase space densities, which in turn is identical to that for $`\overline{d}`$ at the same relative velocity, one can determine a source size for $`\overline{p}`$ using measured anti-proton and anti-deuteron spectra. Figure 1 shows the $`R_{\mathrm{inv}}`$ for antiprotons (open triangles) estimated (using Eq. (9)) from the measured dbars  spectra of $`\overline{d}`$ and weak-decay-corrected $`\overline{p}`$s at three $`K_T`$ bins for the $`18\%`$ most central Au+Au collisions at $`\sqrt{s}=130A`$ GeV. As expected, the radii for the massive $`\overline{p}`$ exhibits a smaller decrease with $`K_T`$ compared to that of pions. The $`\overline{p}`$ transverse momentum spectra for this centrality from the STAR pbars  and the PHENIX spectrap  data are corrected for the weak hyperon decay lambda  and matched in the overlap region. We fit this spectra in the range $`0.25<p_T<3.2`$ GeV/c with a Gaussian distribution $`EdN/d^3p=dN/dy\mathrm{exp}(p_T^2/2\sigma ^2)/(2\pi \sigma ^2)`$ resulting in a yield of $`dN/dy=11.4`$ and width $`\sigma =872`$ MeV. Given the limited range and poor statistics for anti-deuterons, we choose a constant average value, $`R_{\mathrm{inv}}=5.34`$ fm, independent of $`K_T`$. Applying Eq. (Entropy Production at RHIC), the entropy carried by the anti-protons is then $`dS/dy|_{\overline{p}}=118`$.

The procedure to estimate the $`\overline{p}`$ source radii can be identically applied to extract the kaon radii using the $`\varphi `$ meson, $`\varphi K^+K^{}`$, in place of the deuteron. The coalescence arguments used in Eq. (9) ignore the binding energy of the composite particle, which is justified if the binding energies are much less than the breakup temperatures. This is an excellent assumption for the deuteron and a reasonable assumption for the $`\varphi `$ which is a resonance with a binding energy of $`B=30`$ MeV. Although we neglect the binding energy in this analysis, one could adjust the expression for $`R_{\mathrm{inv}}^3`$ by a factor of $`e^{B/T}`$ to account for the binding energy prattllope .

Figure 1 displays the $`K^{}`$ invariant radii extracted from measured midrapidity $`\varphi `$ phis  and $`K^{}`$ kaons  yields for the 11% most central Au+Au collisions at $`\sqrt{s}=130A`$ GeV. As expected the kaon radii fall between those for the proton and pion. The $`K_T`$ dependence for $`K^{}(K^+)`$ radii is parameterized as $`R_{\mathrm{inv}}=A/\sqrt{m_T}`$, with $`A=4.45(4.74)`$ fm GeV1/2. The estimated kaon radii are consistent with that of the preliminary three-dimensional measurements Eno  at the highest RHIC energy of $`\sqrt{s}=200A`$ GeV. We note that the $`m_T`$ dependence of the HBT radii $`R_{\mathrm{inv}}=A/\sqrt{m_T}`$ for these hadrons at RHIC are not consistent with one scaling coefficient as observed at the SPS energy Murray ; the fitted values here are $`A=3.81,4.45,5.32`$ fm GeV1/2 for $`\pi ^{},K^{},\overline{p}`$. The measured single particle $`m_T`$ spectra for $`K^{}`$ and $`K^+`$ kaons  are fitted with an exponential function $`EdN/d^3p=dN/dy\mathrm{exp}[(m_Tm_K)/T]/(2\pi T(m_K+T))`$ that results in a rapidity density of $`dN/dy=38.2(42.1)`$ for $`K^{}(K^+)`$ using a slope parameter of $`T=280`$ MeV. Inserting the spectra along with the the parameterization of $`R_{\mathrm{inv}}`$ into Eq. (Entropy Production at RHIC) yields $`dS/dy|_K^{}=256`$ and $`dS/dy|_{K^+}=286`$. We assume that the neutral (anti)kaons carry the same entropy as the corresponding charged (anti)kaons.

The measured lambda  spectra for $`\overline{\mathrm{\Lambda }}`$ for the 11% most central collisions at midrapidity has feed-down from the electromagnetic decay of the $`\overline{\mathrm{\Sigma }}^0`$ and the weak decays of the $`\overline{\mathrm{\Xi }}^+`$, $`\overline{\mathrm{\Xi }}^0`$, and $`\overline{\mathrm{\Omega }}^+`$. The feed-down from the charged $`\overline{\mathrm{\Xi }}`$ and $`\overline{\mathrm{\Omega }}`$ states is corrected by using the exponential parametrization of the the measured spectra for $`\overline{\mathrm{\Xi }}^+`$ and $`\overline{\mathrm{\Omega }}^+`$ of Ref. casom . The corrected $`\overline{\mathrm{\Lambda }}`$ spectra is then fitted with an exponential function resulting in $`dN/dy=7.4`$ and slope parameter $`T=390`$ MeV. The entropy rapidity density for the antihyperons estimated by assuming that they have the same source radii as the antiprotons are presented in Table I. To estimate $`dS/dy`$ for $`\overline{\mathrm{\Sigma }}^\pm `$, we note that the particle density in the statistical model in the Boltzmann approximation is $`N/V=gm^2TK_2(m/T)\mathrm{exp}(\mu /T)/(2\pi ^2)`$, where the notations are as usual. Assuming isospin symmetry $`\overline{\mathrm{\Sigma }}^+\overline{\mathrm{\Sigma }}^0\overline{\mathrm{\Sigma }}^{}`$, and same temperature $`T_{\overline{\mathrm{\Lambda }}}T_{\overline{\mathrm{\Sigma }}}=173`$ MeV, the thermodynamical model gives the ratio $`\overline{\mathrm{\Lambda }}/\overline{\mathrm{\Sigma }}^{\pm ,0}=1.47`$ and thus $`\overline{\mathrm{\Sigma }}^\pm =(\overline{\mathrm{\Lambda }}+\overline{\mathrm{\Sigma }}^0)/2.94`$. The $`dS/dy`$ extracted under this assumption for $`\overline{\mathrm{\Sigma }}^\pm `$ is listed in Table I.

To extract the entropy density for baryons, we employ the measured ratio for the yields $`\overline{p}/p=0.72`$, $`\overline{\mathrm{\Lambda }}/\mathrm{\Lambda }=0.71`$, $`\overline{\mathrm{\Xi }}^+/\mathrm{\Xi }^{}=0.83`$, $`\overline{\mathrm{\Omega }}^+/\mathrm{\Omega }^{}=0.95`$ ratio , and assume that all (anti)baryons have the same source radii as antiprotons. Furthermore, the ratios are found to be independent over a large range a $`p_T`$ suggesting that the transverse momentum distribution and thus the slope temperatures of the baryons and their antibaryons can be taken to be identical. The $`dS/dy`$ for the baryons thus estimated are shown in Table I.

The values of $`R_{\mathrm{inv}}`$ shown in Fig. 1 (top panel) may be used in conjunction with the corresponding spectra of the hadrons to calculate at a given $`p_T`$ the average phase space density at freeze-out $`f(p_T)=/2\sqrt{2}`$, where $``$ is given by Eq. (Entropy Production at RHIC). This is displayed in Fig. 1 (bottom panel) where the pion $`f`$ at RHIC is remarkably higher than the values extracted at SPS na44 ; ferenc  and AGS e877  energies.

Under the assumption of isospin symmetry the total entropy rapidity density for the hadron-gas at freeze-out for central collisions at RHIC is 4451. Despite the large uncertainty in source sizes, this number is surprisingly robust since the radii only enter logarithmically. A 30% change of $`R_{\mathrm{inv}}^3`$ changes the entropy by 0.26 units per particle. For baryons, which have approximately 10 units of entropy per particle, the overall entropy would change by only a few percent. Another source of error comes from estimating the contribution of the $`\eta `$ and $`\eta ^{}`$. However, since an $`\eta `$ carries two to three times more entropy than a pion, and since the $`\eta `$ decays into two or three pions, the net entropy is rather insensitive to the fraction of pions assigned to $`\eta `$s. The same would be true if one were to explicitly include the entropy of the $`\varphi `$ meson while correspondingly reducing the kaon yields. Since the yields are not measured to better than 10% accuracy, the systematic error in estimating the final-state entropy is estimated $`10\%`$.

We now estimate the entropy per rapidity density generated in the QCD partonic matter in the early stage at RHIC. Assuming boost-invariant longitudinal expansion, the energy density averaged over transverse coordinate is Bjorken

$$\overline{ϵ_{\mathrm{Bj}}(\tau )}=\frac{dE_T}{dy}\frac{1}{\tau \pi R^2}.$$
(10)
Using $`dE_T/dy542`$ GeV etp  and $`\pi R^2=120`$ fm2 for the $`11\%`$ most central Au+Au collisions, one obtains $`\overline{ϵ(\tau )}=4.5`$ GeV/fm3 at a proper time $`\tau =1`$ fm/c. Present data from lattice QCD simulation at $`\mu _B=0`$ shows a rapid increase of $`ϵ/T^4`$ by about an order of magnitude within a narrow temperature range $`\mathrm{\Delta }T<0.1T_c16`$ MeV near the critical temperature $`T_c`$ lattice . This may indicate a first order transition, or nearly first-order transition. For three degenerate quark flavors, $`n_f=3`$, with a mass $`m/T=0.4`$, transition to a plasma phase occurs at $`T_c=154`$ MeV at a critical energy density $`ϵ/T_c^46`$. Assuming the initial matter at RHIC is already thermalized at 1 fm/c, lattice data indicates that $`\overline{ϵ(\tau )}=4.5`$ GeV/fm3 would correspond to an (average) initial temperature of $`T230`$ MeV and pressure density $`p3T^4=1.1`$ GeV/fm3. The resulting entropy density for the 3 flavor system is $`s=(ϵ+p)/T=24.3`$ fm-3, and the average entropy rapidity density at $`\tau =1`$ fm/c is

$$\frac{dS}{dy}|_{\mathrm{QGP};n_f=3}=\tau d^2x_Ts2922.$$
(11)
The Bjorken estimate of the energy density only provides a lower bound since it ignores the energy of the local longitudinal (along the beam axis) motion and neglects the longitudinal work caused by the expansion of the plasma. The inclusion of local longitudinal motion would increase the estimate by a factor of $`4/\pi `$ if the emitted particles were massless and if there was no transverse collective flow. Accounting for transverse flow and masses should result into the increase being in the 10-15% range. For massless particles undergoing a purely one-dimensional Bjorken expansion, the energy density would fall as $`\tau ^{4/3}`$, and $`dE_T/dy`$ would change by a factor of $`\tau ^{1/3}`$. This would be a factor of 2 for times between 1 and 8 fm/c. The loss of $`E_T`$ from longitudinal work is proportional to the $`\pi R^2P𝑑\tau `$ in hydrodynamic models, or $`\pi R^2T_{zz}d\tau `$ if the stress tensor is non-isotropic, i.e., viscous effects are included. The effects of viscosity have been studied within the context of cascade models. Whereas, one-dimensional hydrodynamic models might show that $`dE_T/dy`$ is higher than that at breakup, the inclusion of viscosity reduces that increase to the 50% range pang ; cheng . After including the effects of longitudinal work and longitudinal motion, estimates of the initial average energy density might be in the range of $`7\pm 2`$ GeV/fm3, depending on the equation of state and the viscosity.

Of course, the initial entropy is less than the final state entropy. One expects the entropy to increase during the reaction by the order of 10% from viscous effects or shocks seibertshocks . If one were to lower the final-state estimate of $`dS/dy`$ by 10% it would come close to the value of $`dS/dy`$ predicted by the lattice gauge calculation. A lack of complete thermalization of the initial state would also point to a lack of entropy at early times. If the initial energy was stored in classical gluon fields cgc , the entropy would be negligible. An underpopulation of quark degrees of freedom would also lower the entropy. Since flavor and color charges cannot sample the entire volume, finite-size effects can also reduce the initial entropy, but only by a few percent joerg ; elze . Since the amount of entropy produced during the expansion is unknown, it is impossible to rule out that the initial state might have been of lower entropy.

Thus, our estimate of the final-state entropy, $`dS/dy=4451`$, does not seem inconsistent with expectations from lattice calculations from analyzing Fig. 2. It therefore seems puzzling that hydrodynamic calculations based on lattice-inspired equations of state matched spectra and failed to fit HBT sizes Soff ; teaneyhydro ; heinzkolb ; kolbheinz . An overestimate of a factor of two in $`R_{\mathrm{inv}}^3`$ would suggest an overestimate of the entropy of the order of $`\mathrm{ln}(2)`$ per pion, or $`400`$ units per rapidity. However, it appears from the data that other degrees of freedom, especially the baryonic degrees of freedom, compensate for the lack of entropy in the pions. Although the pions comprise two thirds of the final-state particles, they provide less than half the final-state entropy. The shift of entropy from pions to the (anti)baryon sector was discussed by Rapp and Shuryak rappshuryak  as a natural consequence of an increased phase space density of pions. This follows from the reaction $`p\overline{p}N\pi `$. Since pion number population is difficult to maintain, the overpopulation of pions was expected in an isentropic expansion greiner . An additional factor in increasing the overpopulation of pions would be a lowering of hadronic masses associated with chiral symmetry restoration brownrho . A lowering of hadron masses would lead to an increase in the population of heavier hadrons, many of which would decay back to multiple pions after the normal QCD vacuum is restored haglinpratt .

It should be stressed that the small pionic sources extracted from HBT remain a puzzle. Since the overall entropy of the final state is consistent with the lattice gauge theory, the possibility remains open that a quark-gluon plasma is indeed created at RHIC but that the dynamics of the expansion and decoupling somehow differ from current hydrodynamic descriptions Hirano ; KolbRapp .

To illustrate the sensitivity of the entropy to the number of degrees of freedom, Fig. 2 also displays $`dS/dy`$ for a pion gas. Since the pion gas has only three effective degrees of freedom, it has a much higher temperature for a given energy density than does the lattice equation of state. From Eq. 1, it therefore has a lower entropy when compared at the same energy density. At high temperatures, a hadron gas would effectively incorporate more degrees of freedom, and for energy densities near or above 10 GeV/fm3 would have more effective degrees of freedom than a plasma due to the extremely large number of baryonic resonances with masses between 1.0 and 2.0 GeV. However, since the volume of a baryon is in the neighborhood of 2 fm3, energy densities above 0.5 GeV/fm3 preclude a serious description in terms of a non-interacting hadron gas.

In conclusion, we have calculated the final-state entropy from $`130A`$ GeV Au+Au collisions at RHIC to be $`dS/dy=4451`$. We estimate that this number, which is based solely on measured spectra and extracted source-size radii of the principal final-state constituents, to have a systematic uncertainty of 10%. This number is not out of line with expectations from lattice gauge theory. Improved experimental measurements of spectra and source size might reduce the systematic uncertainties to the 5% level, but it is the uncertainty in the reaction dynamics between the initial and final stages of the collision that precludes our making a more definitive conclusion about the initial state entropy. Although our number represents an important upper bound on the initial state entropy, the amount of entropy generated throughout the reaction and the initial energy density remain objects of conjecture. These issues will be settled largely by other observables. For instance, elliptic flow measurements are sensitive to the viscosity of the initial stages molnar  and HBT measurements provide important three-dimensional constraints about the dynamics.

###### Acknowledgements.

This work was supported by the National Science Foundation under Grant No. PHY-0070818, and by the U.S. Dept. of Energy, Grant DE-FG02-03ER41259.


text:
# 1 Introduction

## 1 Introduction

Optical trapping, which is the trapping and manipulation of microscopic particles by a focussed laser beam or beams, is a widely used and powerful tool. The most common optical trap, the single-beam gradient trap, commonly called optical tweezers, consists of a laser beam strongly focussed by a lens, typically a high-numerical aperture microscope objective, with the same microscope being used to view the trapped particles (see fig. 1) . The trapped particle is usually in a liquid medium, on a microscope slide. Commonly used laser sources employed for trapping range from He-Ne lasers, through Ar ion and semiconductor lasers to TiS and NdYAG lasers. Varying laser powers are used in a broad range of applications of optical tweezers - from just a few milliwatts to hundreds of mW. For most of the lasers used, when the beam is passed through the objective lens, the focal spot of the trapping beam is of the order of a micron. The trapped objects can vary in size from hundreds of nanometres to hundreds of microns.

Although simple trapping and manipulation are sufficient for many applications, the use of optical trapping for quantitative research into physical, chemical, and biological processes, typically using a laser-trapped particle as a probe, requires accurate calculation of the optical forces and torques acting within the trap. The approximate methods commonly used for such calculations may prove inadequate for many applications . An accurate quantitative theory of optical trapping not only allows such calculations to be performed, but also greatly extends the usefulness of techniques such as optical force and torque measurement and optical particle characterisation.

The concept of optical trapping is based on a gradient force causing small particles to be attracted to regions of high intensity in a tightly focussed laser beam . Other optical forces, due to absorption, reflection, and scattering are termed scattering forces. Both the gradient and scattering forces result from the transfer of momentum from the trapping beam to the particle. Optical torques can also be produced by the transfer of angular momentum from the beam, which can result from birefringence, or particle shape, or absorption of a beam carrying non-zero angular momentum .

Although the optical forces and torques within a trap result in a straightforward manner from the change in momentum and angular momentum due to scattering processes, exact calculation is difficult, and a number of approximations are usually made. Approximate calculations use geometric optics for large particles (radius $`a>5\lambda `$), or assume that a small particle ($`a<\lambda /2`$) acts as a Rayleigh scatterer or a point-like polarisable particle. This leaves a large range of particle sizes without adequate results; an unfortunate gap since particles of size comparable to the trapping wavelength are of great practical interest because they can be readily and strongly trapped, and can be used to probe interesting microscopic and mesoscopic phenomena. Reecnt theoretical efforts have individually eliminated some of the deficiencies due to the various approximations used, but there still exists no general correct theory .

The lack of suitable theory is even more acute when non-spherical or anisoropic particles are considered. Non-spherical paricles are of particular interest due to their suitability for use as microscopic probes, and their frequent occurence in nature, for example, biological structures and crystals are usually non-spherical, and are often anisotropic. The possibility of rotating or controlling the orientation of such particles greatly extends the range of manipulation possible within an optical trap, introducing new applications, such as, for example, the investigation of microscopic rotational dynamics .

These theoretical deficiencies can be overcome by considering the scattering processes responsible for optical trapping. There exists a well-developed body of work on electromagentic scattering which can be applied to laser trapping in order to determine the scattered fields, from which, in turn, the optical forces and torques can be found.

There are some examples in the literature of attempts to develop a general theory of optical trapping. These are typically restricted to a limited range of particle types and sizes. The systematic application of scattering theory can eliminate these limitation, giving a general theory correct for all particles compositions, including transparent, conductive, absorbing, etc., all particle sizes and shapes, and for arbitrary trapping fields.

## 2 Light scattering in an optical trap

An optical trap in most ways presents a simple electromagnetic scattering problem, with usually only a single particle in a single orientation in the trap at any one time. The major problems are the representation of the beam, and the possibility that the particle is such that calculations will be difficult even for a single particle. A number of the available techniques require that the trapping beam (i.e. the incident field) be represented as a plane wave spectrum or in terms of vector spherical wavefunctions (VSWFs). The trapping beam is usually a strongly focussed (i.e. non-paraxial) Gaussian beam. An immediate problem is that the standard representations of Gaussian beams do not actually satisfy the Maxwell equations, leading to some difficulty in finding a plane wave or VSWF spectrum, though satisfactory methods exist . Additionally, the trapping beam cannot actually be Gaussian, but will have been truncated at some point in the optical system. Non-Gaussian beams are also used for trapping, for example, Laguerre-Gaussian “donut” beams .

We can note that while the existence of a wide variety of techniques for the calculation of scattering is indicative of the lack of a universally superior method, each technique has its own particular advantages and disadvantages, and we can find a usable or even ideal method for any particular case. A brief survey of commonly used methods follows, focussing on application to trapping problems. Computer codes implementing many of these methods are available .

In general, the solution of an electromagnetic scattering problem requires the solution of the Maxwell equations. Some geometries yield relatively simple solutions, more general cases require direct numerical solution of the Maxwell equations. The best known analytical method is Mie theory, restricted to scattering by a homogeneous isotropic sphere illuminated by a plane wave. Extensions of Mie theory, including the use of spheroidal expansions instead of spherical coordinates allow a broader, but still very limited, range of applicability. In general, numerical methods must still be used to obtain the final solutions in these cases .

Closely related to Mie theory is a family of numerical techniques where the incident and scattered fields are expressed in terms of VSWFs, and the expansion coefficients of the scattered field are found by the boundary conditions at the surface of the scatterer. These methods include the point matching method  and the T-matrix method . The T-matrix method is widely used, computationally efficient for axisymmetric particles, and a number of computer codes implementing this technique are freely available . The T-matrix method is of particular interest, since, for a given scatterer and illuminating wavelength, the T-matrix only needs to be computed once, and can then be used for repeated calculations. A surface integral over the particle must be computed, but in the case of a rotationally symmetric particle, this integral reduces to one dimension .

These surface-based techniques are, in their simple forms, restricted to homogeneous particles, though extensions to layered particles exist . Other surface methods include the generalised multipole technique  and the method of moments .

If the scattering particle is such that techniques such as those above cannot be used, there are a number of general techniques, in principle usable for any scattering problem. In general, these methods are computationally intensive .

Since the Maxwell equations are a set of differential equations, finite difference or finite element methods can be used. The finite difference time domain method (FDTD) is widely used in computational electromagnetics, and can be applied to scattering problems . A discrete grid of points in space is set up, and the fields at successive time steps are calculated. Since the discretisation in space must be much smaller than the wavelength ($`\lambda /20`$), and a correspondingly small time step must be used, only a relatively small volume can be used for the calculation, so the boundary conditions at the edge of the computational volume must be carefully chosen, and if the far-field is desired, it must be found from the near-field via a suitable transformation. Finite element methods (FEM) can also be used, again using spatial discretisation to obtain a numerical solution to the system of differential equations . Both of these methods a conceptually simple, and can represent a particle of any shape or composition, but are computationally intensive and require special care with the boundary conditions.

An alternative method is the disrete dipole approximation (DDA) where the particle itself is divided into small volumes, each of which can be treated as a simple dipole with a polarisability depending on the composition of the particle. An initial guess of the final field is iteratively improved until convergence is obtained . Computer codes implementing DDA are publicly available .

We can note that the characterisation of a laser trap will require repeated calculations with the same particle in different positions and orientations within the trap. Thus, the T-matrix method is attractive, since the T-matrix need only be found once. The T-matrix method requires that the incident beam be represented in terms of VSWFs; this can be done directly, or via an intermediate plane wave expansion. The need for repeated calculations makes the general methods (such as FDTD, FEM and DDA) less attractive since the entire calculation must be repeated if the incident beam is changed. In these cases, the total calculation required can be minimised by representing the trapping beam by a plane wave spectrum, and calculating the scattering for the different angles of plane wave illumination. This represents no greater computational effort than calculating the scattering for all orientations of the particle at a single location within the trap, and gives results that can be used to find the sccattering at all points and all orientations.

## 3 Calculation of forces

Once the scattering has been calculated, the resulting force and torque can be found from the field, by integrating around the particle. The momentum of an electromagnetic field is given by

$$𝐏_{\mathrm{field}}=ϵ_0_V𝐄\times 𝐁d^3x$$
(1)
and the angular momentum by

$$𝐋_{\mathrm{field}}=ϵ_0_V𝐱\times (𝐄\times 𝐁)d^3x$$
(2)
If we have the far-field scattered field, which is an outgoing radiation field, we can choose a spherical surface for the integration, and the momentum and angular momentum fluxes are normal to the surface. Then, the field can be represented by the electric field alone, with two orthogonal components, with complex amplitudes $`E_\theta `$ and $`E_\varphi `$, $`𝐄_\theta =E_\theta \widehat{\theta }`$ and $`𝐄_\varphi =E_\varphi \widehat{\varphi }`$, tangential to the surface of integration. The rate of transfer of linear and angular momentum per unit area can be found using

$$𝐏_{\mathrm{flux}}=\frac{ϵ_0}{2}(E_\theta E_\theta ^{}+E_\varphi E_\varphi ^{})\widehat{𝐧}$$
(3)
and

$$𝐋_{\mathrm{flux}}=\frac{icϵ_0}{2\omega }(E_\theta E_\varphi ^{}E_\theta ^{}E_\varphi )\widehat{𝐧}$$
(4)
where $`\widehat{𝐧}`$ is the unit vector normal to the surface of integration. Equation (3) can be readily recognised as the familiar result of the Poynting vector divided by $`c`$.

From the conservation of momentum and angular momentum, the integrated momentum and angular momentum fluxes give the force and torque acting on the trapped particle.

If the discrete dipole approximation (DDA) is used to calculate the scattering, the particle is represented by a number of dipoles, within a known field (after the scattering calculations have been performed). Therefore, the force and torque on each dipole can be found, giving the force and torque on the entire particle, and incidentally, the stresses within the particle. This allows the force and torque to be found without needing to integrate around the particle .

The basic methods can be summarised as:

* If the T-matrix method can be used:

1. calculate T-matrix
    2. find VSWF representation of beam at the desired position within the trap, for the desired orientation
    3. find scattered field (in terms of VSWFs)
    4. integrate around particle
* If a general method must be used

1. calculate plane wave scattering by the particle for all angles of incidence
    2. find plane wave spectrum of trapping beam at the desired position
    3. find the scattering for each spectral component, and combine to find the total scattered field
    4. integrate around particle
* If DDA is used

1. find scattering at desired point
    2. find forces and torques on all dipoles, and combine to find the total force

Some sample force calculations are presented in figure 2. The trapped particle is a glass prolate spheroid ($`n=1.5+0.02i`$), with $`a=b=0.5\mu \mathrm{m}`$ and $`c=1.0\mu \mathrm{m}`$, trapped in water ($`n=1.33`$) by a Gaussian beam of free space wavelength $`1064\mathrm{n}\mathrm{m}`$, focussed to a spot width of $`1\mu \mathrm{m}`$. The T-matrix code developed by Mishchenko  was used to find the T-matrix an amplitude matrix. The trapping beam was decomposed into a plane wave spectrum consisting of 97 components. The size parameter of the particle is $`4.95`$. The wavelength and beam profile and particle size, shape and composition were chosen to best model the most typical situation occurring in the experimental realisation of the technique.

Figure 2 shows the forces exerted by the trapping beam, which is propagating from the left to the right, focussed to a small spot in the centre. The magnitude of the force is proportional to the length of the arrows, and the direction is the same as the direction of the arrows. The contours in the background are equal-intensity contours for the trapping beam. The trapping beam intensity falls of with increasing radial distance due to the narrow width of the beam. The intensity falls off axially due to the spreading of the beam. The strong radial trapping forces are clearly visible, and it can be seen that the forces are approximately normal to the intensity contours. Since the radial intensity gradient is so large, the radial forces are the largest optical forces acting on the particle. In the regions of low intensity (at the top and bottom of the middle of the graph), the forces are very small.

## 4 Optical measurement of forces and torques

The forces and torques acting on the trapped particle result from the scattering within the trap, and can be calculated if the scattered fields are known. So far, we have considered techniques for the calculation of the the scattered field. In many cases, this will not be possible, for example, if our particle is a biological structure within a living cell, with unknown optical properties, or a particle with known optical properties at an unknown position within the trap.

If the scattered radiation field can be experimentally measured, equations (3) and (4) can be used to find the force and torque by integrating over a surface around the particle. In principle, the required measurement of the scattered field is possible, but the application of this will be limited by a number of practical considerations, namely:

1. The normal method for measuring scattered light from trapped particles is to measure either the back-scattered light , or the forward scattered light . Back-scattering measurements usually record the total power received by the detector and do not give the directional resolution required for measuring the scattered field (which, in any case, would be virtually impossible to do with any accuracy after the return passage of the back-scattered light through the microscope optics). Forward scattering measurements usually use low resolution detectors, such as quadrant detectors, and do not make use of all the possible measurable information. It is generally impossible to collect all of the scattered light due to the spatial limitations imposed by a typical optical tweezers setup. To collect more than the forward scattered light emergent from the bottom of the trap will be very difficult. The placement of the detector(s) is limited by the design of the laser trap. Since the trapping cell, contained a microscope slide, will restrict collection of side-scattered light, and the focussing lens (the objective) will restrict the collection of back-scattered light, the detector must be placed below the trapping cell (see fig. 3). Since the detector must be capable of spatial resolution, and must be large enough to measure the forward-scattered light that passes through the trapping cell, a CCD array is ideal.

A CCD array, by itself, is not sufficient to measure the field; only the intensity will be measured, and it will still be necessary to measure the state of polarisation. This can be done by incorporating a polariser acting as an analyser (see fig. 4).
2. The bottom of the trap will reflect and refract the scattered light. Since the scattered light is initially in a higher refractive index medium (the fluid within the trapping cell and the glass microscope slide), some of this light will be totally internally reflected, and will be unmeasurable (see fig. 5). For typical values of refractive index, the maximum scattering angle for which the light can be measured is approximately $`45^{}`$. Reflection can also occur at interfaces internal to the trapping cell, although these will be smaller since the refractive index differences will be relatively small. The reflections that do occur will depend on the polarisation of the light. If too much light is scattered at large angles (and therefore unmeasurable), it will not be possible to accurately determine the optical force. Being restricted to measuring forward-scattered light will also make it difficult to measure torques other than that acting about the beam axis. The partial reflection of light that exits the trapping cell and is measured can, and in the interests of accuracy, should be compensated for.
3. The resolution and size of the detector will limit the measurements.

The restriction of the measurable scattering to the forward-scattered light limits the applicability of the technique. It remains useful, however, in the two cases of most interest, one case being organisms and biological structures, the second probe particles.

Organisms and biological structures in a trap will typically have unknown optical properties, and often have complex shapes. Thus, it will be necessary to measure the scattering since calculation will be impossible. Since such particles are optically soft (with a relative refractive index $`m1`$), very little light will be scattered at large angles, and the majority of the scattered light can be measured, and the force and torque acting on the particle determined with reasonable accuracy.

If the particle in the trap is not optically soft, with a large refractive index difference compared with the surrounding medium, it will not be possible to collect enough of the scattered light to accurately determine the force and torque in all cases. If, however, the particle in question is a probe particle of known size, shape, and optical properties, but at an unknown position within the trap, the portion of the scattered light that is measured can be used to determine what the total scattering pattern must be. Thus, the position of the probe particle within the trap, and the optical force and torque acting on it, can be determined. In this way, the external non-optical forces acting on the probe particle can be determined .

A special case of the use of a trapped particle as a probe is the rotating probe particle, which can be used to measure the viscosity of a fluid, colloid, or suspension on a microscopic scale. The rotating particle will typically be birefringent, and will remain in the centre of the trap. The change in polarisation of the light on scattering will cause the particle to rotate. Both the optical torque acting on the particle and its rotation speed can be measured .

## 5 Conclusion

Approaching laser trapping as a scattering problem allows the calculation of forces and torques using electromagnetic scattering theory. Such calculations can be performed for all types of particles: transparent, abosorbing, conductive, reflective, anisotropic, complex shapes, etc. For particles for which efficient computational methods (such as the T-matrix method) can be used, calculations are fast and relatively simple, and can be performed on readily available PCs.

This means that the gap in previous calculations, where no adequate results were available for particles comparable in size to the trapping wavelength, can be closed.

Optical measurement of forces and torques acting on particles within the trap can also be performed, by measuring the scattered light. Apart from being free of the usual calibration difficulties for force measurement in optical traps, it can be used where traditional force measurements are impossible, such as measurement of the forces acting on structures within living cells.

Measurement of the scattered light from a known probe particle, coupled with calculation of the scattering in different positions of the trap, allows measurement of the position of, and force acting on, the probe particle. Thus, external non-optical forces can be determined. A special case of this, where the probe is a rotating trapped particle, is particularly simple, and the rotation speed and optical torque can be simultaneously measured.