text:
# New Experimental Limits on Macroscopic Forces Below 100 Microns

## I Introduction

Experimental tests of Newtonian gravity and searches for new weak forces in addition to gravity have been conducted over length scales ranging from light–years down to laboratory distances. Based on recent results adelberger  and reviews fischbach , new forces with strength weaker than or comparable to gravity have been excluded over distances ranging between 200 $`\mu `$m and a light–year. To date, only the single experiment in Ref. adelberger  has attained gravitational sensitivity below 1 mm, and limits on the strength of new interactions increase very rapidly below 100 $`\mu `$bordag . The sub–millimeter range attracts continuing experimental interest because many recent theoretical attempts at the unification of fundamental forces predict specific new phenomena in this regime. We present the results of an experiment which has attained a maximum sensitivity of about four times gravitational strength above 200 $`\mu `$m, and which provides the most sensitive limits on new forces between 10 and 100 microns.

### I.1 Current Limits

Results from experimental searches for new macroscopic forces are most commonly parameterized by a Yukawa interaction. The potential due to gravity and an additional Yukawa–type force between two mass densities $`\rho _1(\stackrel{}{r}_1)`$ and $`\rho _2(\stackrel{}{r}_2)`$ can be written

$$V=𝑑\stackrel{}{r}_1𝑑\stackrel{}{r}_2\frac{G\rho _1(\stackrel{}{r}_1)\rho _2(\stackrel{}{r}_2)}{r_{12}}[1+\alpha \mathrm{exp}(r_{12}/\lambda )],$$
(1)
where $`G`$ is the Gravitational constant, $`r_{12}`$ is the distance between $`\stackrel{}{r}_1`$ and $`\stackrel{}{r}_2`$, $`\alpha `$ is the strength of the Yukawa interaction relative to gravity, and $`\lambda `$ is the range. Current experimental limits on $`\alpha `$ between 1 $`\mu `$m and 1 cm are shown in Fig. 7, together with several recent theoretical predictions of new phenomena.

All previously published limits in Fig. 7 (bold, dashed curves) are obtained from torsion balance experiments. These curves represent 95% confidence level limits with the exception of the result from the Irvine 2–5 cm null experiment hoskins , which is a 1 $`\sigma `$ limit. The torsion balance experiment at the University of Washington has attained gravitational sensitivity or better for ranges above 200 $`\mu `$adelberger . Below 20 $`\mu `$m the best previous limit is derived from the Casimir force measurement by Lamoreaux lamoreaux ; the curve shown in Fig. 7 is based on the analysis by the authors long ; mostepanenko . From the figure, the previous experimental limits allow for forces in nature several million times stronger than gravity at ranges as large as 20 $`\mu `$m.

### I.2 Motivation

In addition to the unexplored parameter space, there are many specific motivations to search for new effects below 100 $`\mu `$disclaimer . Several arise from the decades–long effort to describe gravity and the other fundamental interactions in a unified theoretical framework. The leading candidates for such a theory are string or M theories, which must be formulated in more than three spatial dimensions. In the traditional view, the extra dimensions are taken to be compactified on a scale on the order of the Planck length ($`L_P10^{35}`$ m or $`L_P^110^{19}`$ GeV), a fundamental scale at which the unification is expected to occur. Recently, a class of models has been discovered in which the unification can occur at much lower energies, near the weak scale ($`L_W^11`$ TeV) hewett . In this picture, perhaps the most frequently cited development in particle physics in recent years, the discrepancy between the Planck scale and the weak scale (the so–called hierarchy problem) is removed as a consequence of some of the extra dimensions remaining large and accessible only to gravity, with the Standard Model fields confined to the usual three dimensions arkani . The size $`R`$ of the extra dimensions is given by $`R^n=L_W^{2+n}/L_P^2`$, where $`n`$ is the number of compact dimensions. For $`n=1`$, $`R10^{13}`$ m, clearly ruled out by astrophysics. However, the choice of $`n=2`$ implies $`R1`$ mm, with the consequence that the gravitational potential will behave according to a $`1/r^3`$ law below this scale. As the scale $`R`$ is approached from above, Yukawa corrections are predicted; the model illustrated in Fig. 7 predicts $`\alpha =4`$ floratos .

Another class of predictions arises in superstring theories in which supersymmetry (SUSY) is broken at low energies. These models are of interest as unification scenarios that account for the absence of flavor–changing neutral currents. Superstring theories generally contain gravitationally–coupled scalar fields called moduli, which are massless at the string scale but which can acquire mass from the same process which breaks SUSY. For models in which SUSY is broken by a gauge–mediated process between 10 and 100 TeV, the moduli acquire masses at the sub-eV level corresponding to interaction ranges between 100 $`\mu `$m and 10 cm, with values of $`\alpha `$ as high as $`10^6`$ dimopoulos . Predictions for the largest such effects (for the moduli which couple to the strange quark and gluon) with the uncertainties are shown in Fig. 7. In a related model (radius modulus) in which SUSY is broken near 1 TeV via weak–scale compactification antoniadis , the moduli can acquire Compton wavelengths in the range from 10 $`\mu `$m to 1 mm. The associated strength is predicted to be $`\alpha =1/3`$, as shown in Fig. 7.

The dilaton is another scalar field predicted by string theories, which can also acquire mass from SUSY breaking. Due to its more universal couplings, it is expected to acquire a mass too large to be observable in the scenario in Ref. dimopoulos , but other authors take the dilaton mass to be an unknown parameter. The estimated coupling for the dilaton from Ref. ellis  is shown in Fig. 7, where previous experimental limits exclude a mass less than about $`3.2\times 10^3`$ eV; more recent calculations can be found in Ref. kaplan .

The axion, a light pseudoscalar boson motivated by the strong CP problem of the Standard Model, is also predicted to mediate macroscopic forces. Astrophysical and laboratory bounds have left an allowed window for the axion mass corresponding to 200 $`\mu `$m $`<\lambda <`$ 20 cm rosenberg . The coupling of the associated long–range force between unpolarized test masses is constrained by measurements of the neutron electric dipole moment to the region indicated in Fig. 7. For pseudoscalars which are not derivative coupled (unlike the axion), double exchange leads to a $`1/r^3`$ potential between unpolarized test masses; short range experiments may set better limits on pseudoscalars via this effect krause .

Additional predictions are motivated by the cosmological constant ($`\mathrm{\Lambda }`$) problem, the discrepancy between the observed flatness of the universe versus the extreme curvature expected from the vacuum energy contributions of the Standard Model fields. Recent models sundrum ; beane  assert that if a sufficient fraction of the energy density in the universe is in the form of vacuum energy, consistency with local field theory implies the existence of new interacting quanta with mass on the order of $`\mathrm{\Lambda }^{1/4}`$. The predictions in Fig. 7 are derived from the conservative assumption of cosmological vacuum energy density $`\rho _\mathrm{\Lambda }`$ in the range $`0.1\rho _c<\rho _\mathrm{\Lambda }<\rho _c`$, where $`\rho _c`$ is the critical density of the universe.

## II Experiment

### II.1 Apparatus

The present experiment is illustrated in Fig. 1. The test masses consist of a planar resonant detector mass and a planar source mass driven at the detector resonant frequency. This geometry is chosen to concentrate as much test mass density as possible at the length scale of interest. It is also largely null with respect to $`1/r^2`$ forces and so is effective in suppressing the Newtonian background in the context of a new force search. The source mass consists of a $`35\text{mm}\times 7\text{mm}\times 0.305\text{mm}`$ node–mounted tungsten reed. It is driven at its second cantilever mode by a PZT bimorph at a frequency carefully matched to the 1 kHz resonant frequency of a normal mode of the detector mass. The detector mass consists of a 0.195 mm thick tungsten torsional oscillator of a double–rectangle design. The small rectangle measures $`11.455\text{mm}\times 5.080\text{mm}`$. For the resonant mode of interest (the 5th normal mode) the two rectangular sections counter–rotate about the torsional axis, with most of the amplitude confined to the motion of the small rectangle. This double–torsional arrangement provides significant isolation of the small rectangle from the detector mount, reducing mode damping kleiman .

When data are collected, the front edge of the source mass is aligned with the back edge of the small detector rectangle, and a long edge of the source is aligned above the detector torsion axis. This geometry (called “fiducial” below) maximizes the on–resonance torque on the detector if a mass–coupled force is present between the source and detector.

A stiff conducting shield suspended between the test masses suppresses electrostatic and acoustic backgrounds. The shield consists of a 0.060 mm thick sapphire plate with a 100 nm gold plating. The shield is mounted on a brass support and secured on opposite ends with epoxy.

The test masses and shield support are each suspended from modified optical mirror mounts, which are used as a tilt stages to level each element in the horizontal plane. Each tilt stage is in turn mounted to the bottom segment of a vibration isolation stack. Each stack consists of five solid brass disks connected by fine steel wires, as described in detail elsewhere chan . At the 1 kHz operational frequency of the experiment, each stack provides an attenuation of approximately 200 dB. The vibration isolation stacks of the source and detector masses are suspended from inverted 3-axis micrometer stages, which provide full translation control. The relative positions and coplanarity of the test masses and shield are ascertained by touching the elements against each other, or against a series of 0.5 mm diameter sapphire hemispheres attached to small rods at the end of the shield support.

To further reduce acoustic backgrounds, the apparatus is placed in a 75 liter vacuum bell jar, which is pumped to a pressure of approximately $`2\times 10^7`$ torr with a liquid nitrogen–trapped diffusion pump. In order to isolate the (non-magnetic) central apparatus from stray fields generated in the steel pump components, the bell jar is connected to the diffusion pump via a 1 m long, 10 cm diameter aluminum riser. The position–control micrometer stages are manually operated with torque rods which exit the bell jar via rotary feedthroughs.

In the absence of electromagnetic, acoustic, and vibrational backgrounds, the experiment is limited by thermal noise due to dissipation in the detector mass. To reduce this dissipation, the detector mass is annealed at $`1300^{}`$C for several hours in an induction furnace under a helium atmosphere before installation in the bell jar. This is observed to increase the detector mechanical quality factor ($`Q`$) typically by factors of 5 or more.

The detector temperature is stabilized with an electronic temperature controller using a silicon diode sensor and resistive heating element mounted to the bottom stage of the detector vibration isolation stack. The temperature is maintained at 305 K with fluctuations of about 0.1 K.

### II.2 Readout

The detector readout electronics and source mass PZT drive are illustrated in Fig. 2. Oscillations of the detector mass are read out with a capacitive transducer. The capacitive transducer probe consists of a 2.5 mm diameter brass cylinder supported with its flat end approximately 0.1 mm above a rear corner of the large rectangle of the detector mass. The probe is biased at 200 V through a 100 G$`\mathrm{\Omega }`$ resistor, whose value must be large to reduce current noise. The front end of the preamplifier consists of an SK152 JFET located in a small box immediately above the detector oscillator. The 100 mK noise temperature of the SK152 preamplifier ensures that this circuit is more than sufficient for detecting the 300 K thermal oscillations (amplitude $``$ 100 fm) of the detector mass. The JFET preamplifier is followed by a second preamplifier (Stanford Research SR560), filters, and finally a two phase lock–in amplifier. The total voltage gain from the capacitive probe to the lock–in input is about 1600.

A crystal–controlled oscillator provides a reference signal for the lock–in amplifier and drives the source mass PZT through a 1:10 step–up transformer.

## III Data Sample

### III.1 Force Measurement Data

The data were recorded in 108 “cycles” over five days. Each cycle consists of 7 runs containing both force measurement data and diagnostic data. The diagnostic data are used to continually monitor the performance of the instrument. For all runs the bandwidth of the lock-in amplifier was set to 250 mHz (roughly 5 times the width of the detector resonance) to include the noise power of the detector thermal oscillations, which is used for calibration. A data sampling rate of 1 Hz was used, corresponding to an over-sampling of about a factor of 10.

Within each cycle, the runs were taken in the following order: First, a stable dc bias (usually 5–10 V) was applied to the shield to induce a large resonant electrostatic signal. A run of 120 samples was taken with the PZT drive frequency set 30 mHz below the approximate value of the detector resonant frequency. Four subsequent runs were taken, each with the PZT drive frequency incremented by 15 mHz relative to the previous run in order to display the detector resonance. The precise value of the detector resonant frequency was determined from the five biased runs, and was then used as the drive frequency for the next run of the cycle. The shield was then grounded and a run of 720 samples (12 minutes) was recorded. This was followed by a shorter run of 288 samples (4.8 minutes) taken with the PZT drive set to 1171.000 Hz (at least 2 Hz below the detector resonance), in order to monitor the system for non–resonant systematic offsets. The entire cycle was repeated indefinitely until an intervention to service the experiment was necessary. A total of 108 such cycles were acquired yielding a total of 77760 on–resonance samples.

A plot of the data from the biased diagnostic runs is shown in Fig. 3 (only data from 10 consecutive cycles are shown for clarity). The signals from the two phase–quadrature channels of the lock-in amplifier are plotted against each other to show the phase behavior of the signal with drive frequency. In all cycles, the signal maximizes at very nearly the same phase and magnitude, indicating good stability of the detector resonant frequency, source mass amplitude, and system gain.

Fig. 4 shows histograms of the on– and off–resonance unbiased data. The plots combine data from all 108 on–resonance and off–resonance runs, each plot displaying data from a single channel of the lock-in amplifier. The data exhibit smooth gaussian behavior centered about common means, as expected in the exclusive presence of detector thermal noise and amplifier noise. The widths of the distributions of the on–resonance data (left-hand plots in Fig. 4) are roughly twice those of the off–resonance data (right-hand plots) due to the contribution of detector thermal noise.

The means of the distributions in Fig. 4 are shown in Fig. 5, in which the data from the separate lock-in amplifier channels are plotted against each other as in Fig. 3. The on–and off–resonance means agree within the 1 $`\sigma `$ standard deviations shown, indicating the absence of any resonant force signal.

The offset from the origin of the mean of the off–resonance distribution is due to leakage of the reference signal internal to the lock–in.

### III.2 Consistency Checks

After the acquisition of the session data, several checks are made to verify the stability of the results and their consistency with known backgrounds. First, several additional cycles are acquired with variations in the test mass geometry, including larger vertical gaps and different overlap configurations. No resonant signal is observed in any of these sessions. This strongly disfavors the possibility that the observed null result is due to a fortuitous cancellation of surface potential, magnetic, and/or acoustic effects, all of which are expected to have different dependencies on the geometry.

Next, several on–resonance runs are acquired with different transducer probe bias voltage settings and with the source mass drive turned off. The observed linear dependence of the rms fluctuation in these data on the probe bias voltage is consistent with detector motion due only to thermal noise and rules out additional motion from transducer back-action noise. This check is important because the magnitude of the detector thermal motion is used for calibration.

The data from the diagnostic runs with shield bias voltage applied can be used to estimate the minimum size of the residual potential difference between the shield and the (grounded) test masses needed to produce a resonant signal. From Fig. 3, the maximum signal with 9.6 V shield bias is about 1.3 mV. Further measurements show that this signal scales as the fourth power of the shield bias voltage, as expected for an electrostatic force between source and detector mediated by a deflected shield. Scaling the ratio of the on–resonance error bar in Fig. 5 to the 1.3 mV diagnostic implies that a residual shield bias of at least 1.5 V would be needed to generate a systematic effect above detector thermal noise. This is about an order of magnitude larger than the measured residual potential difference between the shield and test masses.

Magnetic effects can generate background signals though several mechanisms. The most important effect involves generation of eddy currents when the source mass moves in an external magnetic field. Fields produced by the source eddy currents create eddy currents in the detector, which then interact with the applied field. Studies of this effect with large applied fields show that the induced force varies as the square of the applied field (as expected), and extrapolation to the ambient field actually present indicates that this effect should be about 5 times smaller than the thermal noise limited sensitivity.

### III.3 Auxiliary Measurements

Before acquiring the session data samples, several measurements are made to support the analysis of the data for evidence of a net force. These include the relative phases of the motion of the test masses, the detector mechanical quality factor, and a precise survey of the test mass geometry.

The Yukawa interaction for the test mass geometry of the experiment is used to model any potential signal. Such an interaction can result in either a purely attractive or repulsive force between the source and detector masses. In order to ascertain the phase of the signal for such a force, the source and detector are placed in their fiducial positions with an additional vertical offset so that both are situated above the electrostatic shield. An electrostatic bias of 1.5 V is applied to the detector mass with a dry cell in order to induce a large, stable, purely attractive electrostatic coupling between the test masses. A set of runs is taken with small increments of the PZT drive frequency between successive runs. The position of the data point in the complex plane corresponding to the maximum signal as a function of drive frequency is recorded. The phase of this signal corresponds to that of a purely attractive force. The phase of a purely attractive signal is $`189^{}`$, just below the negative horizontal axis in the plots in Figs. 3 and 5.

The detector mechanical quality factor $`Q`$ is measured by applying a 100 mV resonant ac signal directly to the detector, switching the detector to ground, and observing the ring-down signal on the lock-in amplifier with the reference frequency tuned 1 Hz off the detector resonance. A value for $`Q`$ is obtained from a least-squares fit to the ring-down waveform. Several such measurements are made before and after the data cycles reported here, yielding a $`Q`$ of $`25522\pm 29`$.

In the initial step of the test mass geometry survey, the source and detector masses are positioned over the electrostatic shield. The source is tilted on its horizontal axes to ensure that the front left corner is the lowest point on the bottom surface. This corner is then used as a probe to make a map of the top surface of the detector mass, using a square grid of points with a point density of 1.0 mm<sup>-2</sup>. The tilt stage of the detector mass is then adjusted to optimize the level of the detector, and the surface is mapped again. A similar procedure is then carried out for the electrostatic shield. Finally, a map of the bottom surface of the source mass is obtained by touching this surface off against the sapphire hemisphere probe fixed to the end of the shield mount. The level is optimized with the source mass tilt stage, and the map is repeated.

In order to obtain the modeshape and amplitude of the driven source mass, the PZT is driven at the resonance of the detector and another grid of points is obtained by touching the bottom surface of the source mass against the sapphire probe. These measurements are done at atmospheric pressure and the touch–offs are determined acoustically.

The system is then closed and brought down to diffusion pump pressure. The relative horizontal positions of the source and detector masses are determined by touching an edge of the source against an edge of the detector, and then the horizontal positions are set to the fiducial locations. Next, the source is lowered to bring the test masses into momentary contact, and their vertical and horizontal positions are recorded. The detector mass is then positioned under the electrostatic shield. First it is brought as far forward under the shield as possible (without the detector and shield stacks making physical contact), then centered between the sides of the shield mount as determined by touch-off signals. The same net horizontal translations required to position the detector are applied to the source, bringing the test masses into their fiducial positions.

Finally, the minimum vertical gap is established. The top surface of the detector is brought into momentary contact with the bottom surface of the shield, then backed off by 20 $`\mu `$m. The bottom surface of the source is brought into momentary contact with the top surface of the shield, and the source is backed off by 20 $`\mu `$m (the amplitude at the end of the source mass) plus an additional 10 $`\mu `$m for safety.

### III.4 Calibration

Conversion of the voltage data to observed force is achieved by direct comparison of the mean signal to thermal noise. Comparison of the root mean square thermal oscillations of the detector mass to the mean displacement of the detector when driven by the source mass leads to the expression

$$\overline{V^D}=\sqrt{\overline{\left|V^T\right|^2}}\frac{Q}{\omega _0\sqrt{k_BT\rho _d}}\frac{d^3\stackrel{}{r}^{}\stackrel{}{z}^F(\stackrel{}{r}^{})\stackrel{}{f}(\stackrel{}{r}^{})}{\sqrt{d^3\stackrel{}{r}^{}\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2}},$$
(2)
where $`\overline{V^D}`$ is the mean of the distribution of voltages in the data sample, $`\sqrt{\overline{\left|V^T\right|^2}}`$ is the component of the standard deviation of the same data sample due to detector thermal noise, $`k_B`$ is Boltzmann’s constant, and $`Q,\omega _0,T`$, and $`\rho _d`$ are the detector mechanical quality factor, resonant frequency, temperature, and density. Both $`\overline{V^D}`$ and $`\sqrt{\overline{\left|V^T\right|^2}}`$ are evaluated at a phase of $`189^{}`$, corresponding to a purely attractive force. The expression $`\stackrel{}{z}^F(\stackrel{}{r}^{})\stackrel{}{f}(\stackrel{}{r}^{})`$ is the projection of the driving force density $`\stackrel{}{f}(\stackrel{}{r}^{})`$ times the free modal displacement $`\stackrel{}{z}^F(\stackrel{}{r}^{})`$ of the detector at arbitrary point $`\stackrel{}{r}^{}`$ on the detector, and $`\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2`$ the mean square free modal displacement; both are integrated over the detector volume. A derivation of Eq. 2 is given in Appendix A. This expression can be used to convert $`\overline{V^D}`$ directly to observed force. Alternately, the instrument may be calibrated using a calculable electrostatic force or by using the reciprocity of the transducer. Both methods agree with calibration based on the thermal motion, but the uncertainties involved are substantially larger.

## IV Analysis

Generally, the constraints on the Yukawa parameter $`\alpha `$ for a given range $`\lambda `$ are obtained by ascertaining how large of a Yukawa interaction between the test masses could be present and still be consistent with the data in Fig. 5. We present first a simplified analysis based on an idealized geometry and neglecting all systematic effects. As the data are consistent with a null result limited by detector thermal noise, this preliminary estimate can made by comparing the hypothetical Yukawa force between the source and detector with the detector thermal noise.

Following the previous analysis of this experiment in Ref. long , the detector modeshape is taken to be a pure rotation and the source mass motion a pure translation normal to its surface. The source mass is driven at the detector resonant frequency, and only the Fourier amplitude of the Yukawa torque at this frequency is effective in driving the detector. From Ref. long , the torque amplitude is

$$|N_Y(\omega _0)|=2\pi \alpha G\rho _d\rho _sA_dR\lambda ^2I_1(dz_s/\lambda )\mathrm{exp}(g_{sd}/\lambda )[1\mathrm{exp}(t_d/\lambda )][1\mathrm{exp}(t_s/\lambda )].$$
(3)
Here, $`\rho _s`$ is the source mass density, $`A_d`$ is the area of the detector under the source mass (half the area of the small rectangle), $`R`$ is the distance from the edge of the detector to the torsion axis (equal to $`w_d/2`$; see Table 2 below), $`dz_s`$ is the source mass amplitude, $`g_{sd}`$ is the average source–detector gap, $`I_1(dz_s/\lambda )`$ is the modified Bessel function, and $`t_d`$ and $`t_s`$ are the detector and source thicknesses. This expression neglects edge effects and does not include the true cantilever modeshape of the source mass.

The thermal noise torque in the experimental bandwidth on the small detector rectangle is found from the mechanical Nyquist expression:

$$N_T=\sqrt{\frac{4k_BT}{\tau }\left(\frac{mR^2\omega _0}{3Q}\right)}.$$
(4)
Here, $`m`$ is the mass of the small rectangle, $`\tau `$ is the inverse of the bandwidth, and $`Q`$ is the detector quality factor.

The signal–to–noise ratio is the ratio of Eq. 3 to Eq. 4. Setting this equal to unity and solving for $`\alpha `$ yields an approximation to the experimental limit curve shown in Fig. 7. Evaluating this expression using the mean values of the parameters listed in Tables 1 and 2 below, assuming an integration time of 77760 s and an average source amplitude equal to 1/2 the measured tip amplitude, yields a limit in the ($`\alpha ,\lambda `$) space stronger than the result in Fig. 7 by a small numerical factor for ranges $`\lambda >50\mu `$m. The discrepancy increases to about two orders of magnitude at $`\lambda =5\mu `$m due to the high sensitivity of the Yukawa exponential to systematic errors in the source amplitude and the average gap.

To take the precise geometry and the systematic error into account, the Yukawa force between the test masses is computed numerically and constraints on the Yukawa strength $`\alpha `$ in Fig. 7 are calculated using a maximum likelihood technique hagiwara . In this approach, the interval $`[\alpha _{\text{lo}},\alpha _{\text{up}}]`$ which contains the true value of $`\alpha `$ with probability $`CL`$ is given by

$$CL=_{\alpha _{\text{lo}}}^{\alpha _{\text{up}}}p(\alpha |𝒙)𝑑\alpha ,$$
(5)
where $`p(\alpha |𝒙)`$ is the posterior probability density function (p.d.f.) for $`\alpha `$ given the experimental data $`𝒙`$. The posterior p.d.f. is in turn calculated from Bayes’ theorem,

$$p(\alpha |𝒙)=\frac{L(𝒙|\alpha )\pi (\alpha )}{_{\mathrm{}}^{\mathrm{}}L(𝒙|\alpha ^{})\pi (\alpha ^{})𝑑\alpha ^{}},$$
(6)
where $`L(𝒙|\alpha )`$ is the likelihood function and $`\pi (\alpha )`$ is the prior p.d.f. for $`\alpha `$.

In order to account for the effects of various systematics such as test mass geometry, density, and other mechanical properties, the likelihood function in Eq. 6 is replaced by the expression

$$L^{}(𝒙|\alpha )=L(𝒙|\alpha ,𝝂)\pi (𝝂)𝑑𝝂,$$
(7)
where $`𝝂`$ represents the set of systematic variables and $`\pi (𝝂)`$ is their prior p.d.f.. In Eq. 7, $`L(𝒙|\alpha ,𝝂)`$ is given by

$$L(𝒙|\alpha ,𝝂)=\left(\frac{1}{\sigma \sqrt{2\pi e}}\right)^N\mathrm{exp}\left[\frac{(\overline{x}\mu (\alpha ,𝝂))^2}{2(\sigma /\sqrt{N})^2}\right],$$
(8)
where $`\overline{x}`$ is the average of the voltages $`𝒙`$ in the data session, $`\sigma `$ is their standard deviation, and $`N`$ is the effective number of uncorrelated samples. The term $`\mu (\alpha ,𝝂)`$ is the predicted mean voltage for a given $`\alpha `$ and set of systematics $`𝝂`$. It is equivalent to $`\overline{V^D}`$ in Eq. 2, where $`\stackrel{}{f}(\stackrel{}{r}^{})`$ in that equation is understood to be the theoretical value of the Yukawa force density for the experimental test mass geometry. A more detailed expression for the likelihood function is derived in Appendix B.

### IV.1 Evaluation of the Likelihood Function

Eq. 7 is evaluated by Monte Carlo integration over the systematics $`𝝂`$. For a fixed value of the Yukawa range $`\lambda `$, 400 points are thrown in the sample space defined by the volume $`d𝝂`$; the integrand is calculated for each point and added to a running total.

The parameter $`\overline{x}`$ in Eq. 8 represents the mean of the on–resonance data with respect to the mean of the off–resonance data (which measures the effective zero), at the phase corresponding to the predicted signal. To find $`\overline{x}`$, all data are first projected onto the phase of the purely attractive signal. The mean of the off–resonance data is then subtracted from that of the on–resonance data. The value of $`\overline{x}`$ obtained is $`0.44\mu `$V.

The quantity $`\sigma /\sqrt{N}`$, equal to the standard deviation of the mean of the data, is calculated in the following way: First, the projected data sets are partitioned into $`N_{\text{part}}`$ sets of equal samples, with $`N_{\text{part}}`$ large enough so that there is very little correlation between partitions. The data in each of the partitions are averaged separately and the standard deviation of the resulting averages $`\sigma _{\text{av}}`$ is computed. The standard deviation of the mean is then calculated from $`\sigma _{\text{av}}/\sqrt{N_{\text{part}}}`$. The effective number of uncorrelated samples is the value of $`N`$ required to make $`\sigma /\sqrt{N}=\sigma _{\text{av}}/\sqrt{N_{\text{part}}}`$. Finally, values of $`\sigma /\sqrt{N}`$ obtained separately for the on– and off–resonance data sets are added in quadrature. The resulting value of $`\sigma /\sqrt{N}`$ is $`0.82\mu `$V.

A list the parameters $`𝝂`$ needed to evaluate Eq. 2 is given in Table 1. Table 2 lists additional systematics which specify the test mass geometry used in the computation of the integral $`d^3\stackrel{}{r}^{}\stackrel{}{z}^F(\stackrel{}{r}^{})\stackrel{}{f}(\stackrel{}{r}^{})`$ in Eq. 2.

All parameters in Tables 1 and 2 are listed with their mean value and error. As the information about most of the $`𝝂`$ is limited, a uniform p.d.f. centered about the mean with a width equal to twice the error is assumed for each of the priors $`\pi (𝝂)`$.

The quantities $`G`$, $`k_B`$, $`\rho _d`$, and $`\rho _s`$ are taken from tabulated values hagiwara2 . To obtain the thermal noise $`\sqrt{\overline{\left|V^T\right|^2}}`$, the standard deviation of the off–resonance data (that is, the noise due to the JFET amplifier) is subtracted in quadrature from the standard deviation of the on–resonance data. The error on $`\sqrt{\overline{\left|V^T\right|^2}}`$ is estimated in the same way as the standard deviation of the mean $`\sigma /\sqrt{N}`$ of the projected data sets.

Several measurements of the detector mechanical quality factor $`Q`$ were made immediately before and after the data sessions, with the test masses in their fiducial positions. The values in Table 1 represent the resulting average and standard deviation of these measurements.

The quantity $`\sqrt{d^3\stackrel{}{r}^{}\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2}`$ in Eq. 2, the root mean square free modeshape of the detector integrated over its volume, is computed numerically from a complete finite element model of the detector mass and is treated as a single systematic with a uniform error of 10%. The free modeshape is normalized to a maximum amplitude of $`1\mu m`$.

For a particular set of the $`𝝂`$ in Table 2, a value of the expression $`d^3\stackrel{}{r}^{}\stackrel{}{z}^F(\stackrel{}{r}^{})\stackrel{}{f}(\stackrel{}{r}^{})`$ in Eq. 2 is obtained in a separate Monte Carlo integration. In this calculation, the detector curvature, source mass curvature and source mass modeshape are described by second order polynomials obtained from fits to the survey data. The detector modeshape is described by a similar function obtained from the finite element model. The curvature and modeshape functions and the $`𝝂`$ are used to calculate bounding boxes for the source mass and the small rectangle of the detector mass. For a particular position of the source mass, $`1\times 10^5`$ point pairs are thrown. For those pairs landing within the bounding boxes a force density $`\stackrel{}{f}(\stackrel{}{r}^{})`$ is calculated and added to a running total. A total of 30 different phase positions of the source mass between 0 and 2$`\pi `$ radians are sampled, with the range defined by the amplitude $`dz_s`$.

The vertical separation of the test masses is computed from the touch gap, the test mass thicknesses ($`t_d,t_s`$), and a correction calculated from the functions describing the test mass surface curvature. The touch gap ($`g_{sd}`$) is the vertical difference in test mass positions between where the opposing test mass surfaces touch and where they are in their fiducial positions. The error on these quantities is taken to be twice the resolution of the translation stages, or 6.4 $`\mu `$m, as two measurements are needed to determine each. This is also the case for the errors in the test mass widths ($`w_d,w_s`$) and detector length ($`l_d`$). An exception is the source mass tip amplitude ($`dz_s`$), which is the average of five measurements. Also, the accuracy with which the tip of the source mass is aligned with the back edge of the small detector rectangle is 320 microns, as it is aligned by eye and taken to have an error corresponding to a complete turn of the translation stage screw. The accuracy with which the long edge of the source mass is aligned to the torsion axis is the usual 6.4 microns, as this position is determined from touch–offs. All geometry measurements are repeatable over several weeks to within the resolution of the translation stages.

Fig. 6 shows the likelihood function (Eq. 7) as a function of $`\alpha `$ for the case $`\lambda =20\mu `$m, a range at which the experiment has good sensitivity relative to previous experiments. The factor $`\left(\sigma \sqrt{2\pi e}\right)^N`$ has been suppressed as it has no bearing on the posterior p.d.f. for $`\alpha `$ (Eq. 6). The maximum is shifted slightly toward negative (repulsive) values of $`\alpha `$, in accordance with the offset of the mean of the on–resonance data with respect to the off–resonance data (Fig. 5).

### IV.2 Limits on Yukawa Parameters

The function in Fig. 6 is substituted into Eq. 6 and integrated numerically (Eq. 5) over an interval $`[\alpha _{\text{lo}},\alpha _{\text{up}}]`$ to obtain the limit on $`\alpha `$. The interval is centered around the maximum of the likelihood function and the length adjusted by trial and error to obtain a confidence level of 95%. In Eq. 6, a uniform prior distribution for $`\alpha `$ is used. For the case $`\lambda =20\mu `$m, the resulting 95% confidence interval is \[$`5.60\times 10^3<\alpha <3.56\times 10^3`$ \].

The likelihood function is recomputed and the integration repeated for several values of $`\lambda `$ between 5 $`\mu `$m and 500 $`\mu `$m. The resulting 95% confidence level limit curve in the ($`\alpha ,\lambda `$) parameter space is shown in Fig. 7. For each value of $`\lambda `$, the slightly weaker limit corresponding to a repulsive interaction is shown. The new limit is close to three orders of magnitude more sensitive than the previous experimental limits in the range between 10 and 100 $`\mu `$m.

## V Conclusion

An experimental search for new macroscopic forces at short distances has been conducted using 1 kHz mechanical oscillators as test masses with a stiff conducting shield between them. No evidence for a resonant signal above detector thermal noise has been observed. Based on these measurements, roughly three new square decades in the ($`\alpha ,\lambda `$) parameter space are ruled out in the range $`10\mu \text{m}<\lambda <100\mu \text{m}`$. About half of the remaining parameter space for the gluon and strange modulus forces is now excluded. Our limit crosses the predicted line for the dilaton at 36 microns, corresponding to a lower limit on the dilaton mass of $`5.5\times 10^3`$ eV.

From Fig. 7, the sensitivity of the experiment maximizes at about 4 times gravitational strength above $`\lambda =200\mu `$m. For the planar test mass geometry of the current experiment, the Newtonian signal is roughly an order of magnitude below the thermal noise at this range. Detection of this signal should be possible with more statistics and an optimized geometry. The small test masses in this experiment allow for relatively easy control of the minimum separation, currently limited to 100 $`\mu `$m by the thickness of the electrostatic shield. Preliminary tests with stretched metal membranes suggest that shields of sufficiently low compliance can be made with thicknesses as small as 10 $`\mu `$m. Work is underway toward an experiment with a stretched membrane shield in place of the sapphire plate. If the backgrounds can be controlled, this experiment could improve limits between 10 and 50 $`\mu `$m by at least another order of magnitude.

###### Acknowledgements.

The authors would like to thank Elizabeth Lagae for her continuing work in the laboratory, and Cole Briggs, Tracy Buxkemper, Leslie Czaia, Hans Green, Sid Gustafson, and Hans Rohner of the University of Colorado and JILA instrument shops for technical assistance. This work is supported by NSF grant PHY00-71029.

## Appendix A Calibration

The mean signal from the lock-in amplifier data can be converted to observed force by direct comparison to the detector thermal noise. The mean square thermal displacement of the detector mass $`\overline{\left|\stackrel{}{z}^T(\stackrel{}{r})\right|^2}`$ at the position $`\stackrel{}{r}`$ of the capacitive transducer probe is related to the temperature $`T`$ via the equipartition theorem:

$$\frac{1}{2}m(\stackrel{}{r})\overline{\left|\stackrel{}{z}^T(\stackrel{}{r})\right|^2}\omega _0^2=\frac{1}{2}k_BT.$$
(9)
Here, $`\omega _0`$ is the resonant frequency of the 5th detector mode and $`m(\stackrel{}{r})`$ is the modal mass of the detector oscillator at the position of the probe.

For a distributed mass with free modal displacement $`\stackrel{}{z}^F(\stackrel{}{r}^{})`$ at an arbitrary point $`\stackrel{}{r}^{}`$ driven on resonance, the driven displacement amplitude $`d\stackrel{}{z}^D(\stackrel{}{r}^{})`$ due to a differential force amplitude $`d\stackrel{}{F}(\stackrel{}{r}^{})`$ is given by

$$d\stackrel{}{z}^D(\stackrel{}{r}^{})=j\frac{Q}{\omega _0^2m(\stackrel{}{r}^{})}\stackrel{}{z}^F(\stackrel{}{r}^{})\frac{\stackrel{}{z}^F(\stackrel{}{r}^{})d\stackrel{}{F}(\stackrel{}{r}^{})}{\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2},$$
(10)
where $`Q`$ is the detector mechanical quality factor and the dot product picks out the component of the force parallel to the modal displacement. (Note the use of $`e^{+j\omega t}`$ time dependence.) The differential amplitude at the probe is

$$d\stackrel{}{z}^D(\stackrel{}{r})=\left|d\stackrel{}{z}^D(\stackrel{}{r}^{})\right|\frac{\stackrel{}{z}^F(\stackrel{}{r})}{\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|}.$$
(11)
Substituting Eq. 11 into Eq. 10 and integrating over $`\stackrel{}{r}^{}`$ yields the amplitude of the detector at the position of the transducer due the total force on the detector:

$$\stackrel{}{z}^D(\stackrel{}{r})=j\frac{Q}{\omega _0^2\rho _d}\stackrel{}{z}^F(\stackrel{}{r})\frac{d^3\stackrel{}{r}^{}\stackrel{}{z}^F(\stackrel{}{r}^{})\stackrel{}{f}(\stackrel{}{r}^{})}{d^3\stackrel{}{r}^{}\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2}.$$
(12)
Here the force amplitude density $`\stackrel{}{f}(\stackrel{}{r}^{})`$ is defined by $`\stackrel{}{f}(\stackrel{}{r}^{})d\stackrel{}{F}(\stackrel{}{r}^{})/d^3\stackrel{}{r}^{}`$, $`\rho _d`$ is the detector mass density, and the definition of the modal mass

$$m(\stackrel{}{r})=\frac{\rho _d\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2d^3\stackrel{}{r}^{}}{\left|\stackrel{}{z}^F(\stackrel{}{r})\right|^2}$$
(13)
has been used.

Combining Eqs. 912, and 13 yields

$$\frac{\stackrel{}{z}^D(\stackrel{}{r})}{\sqrt{\overline{\left|\stackrel{}{z}^T(\stackrel{}{r})\right|^2}}}=\frac{Q}{j\omega _0\sqrt{k_BT\rho _d}}\frac{\stackrel{}{z}^F(\stackrel{}{r})}{\left|\stackrel{}{z}^F(\stackrel{}{r})\right|}\frac{d^3\stackrel{}{r}^{}\stackrel{}{z}^F(\stackrel{}{r}^{})\stackrel{}{f}(\stackrel{}{r}^{})}{\sqrt{d^3\stackrel{}{r}^{}\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2}}.$$
(14)
The linear response of the capacitive transducer insures that the ratio $`\stackrel{}{z}^D(\stackrel{}{r})/\sqrt{\overline{\left|\stackrel{}{z}^T(\stackrel{}{r})\right|^2}}`$ is related to the measured voltages on the lock-in amplifier by

$$\frac{\left|\stackrel{}{z}^D(\stackrel{}{r})\right|}{\sqrt{\overline{\left|\stackrel{}{z}^T(\stackrel{}{r})\right|^2}}}=\frac{\overline{V^D}}{\sqrt{\overline{\left|V^T\right|^2}}},$$
(15)
where $`\overline{V^D}`$ is the mean of the distribution of voltages in the data sample and $`\sqrt{\overline{\left|V^T\right|^2}}`$ is the component of the standard deviation of the same data sample due to detector thermal noise. Both quantities are evaluated at the phase of interest. Substitution of Eq. 15 into Eq. 14 yields Eq. 2.

## Appendix B Likelihood Function

In addition to the unknown parameter $`\alpha `$ corresponding to the Yukawa strength, the signal depends on a set of systematics $`𝝂`$ which includes test mass geometry, density, detector temperature and mechanical properties. The effects of uncertainties in these quantities can be accounted for in the posterior p.d.f. (Eq. 6) by integrating over the systematics:

$$p(\alpha |𝒙)=p(\alpha ,𝝂|𝒙)𝑑𝝂.$$
(16)
Assuming the prior joint p.d.f. for $`\alpha `$ and $`𝝂`$ factorizes, this is equivalent to replacing the likelihood function $`L(𝒙|\alpha )`$ with

$$L^{}(𝒙|\alpha )=L(𝒙|\alpha ,𝝂)\pi (𝝂)𝑑𝝂.$$
(17)
Here, $`\pi (𝝂)`$ is the prior p.d.f. of the $`𝝂`$. The function $`L(𝒙|\alpha ,𝝂)`$ is the joint p.d.f. for the data, regarded as a function of both the unknown Yukawa strength $`\alpha `$ and the systematics $`𝝂`$:

$$L(𝒙|\alpha ,𝝂)=\underset{i}{}F(x_i|\alpha ,𝝂).$$
(18)
For the present experiment, the data $`x_i`$ are the sampled voltages from the lock-in amplifier, which are Gaussian distributed (Fig. 4). Therefore the function $`F`$ is taken to be the Gaussian distribution $`G(x_i|\mu (\alpha ,𝝂),\sigma )`$, where $`\mu (\alpha ,𝝂)`$ is the predicted mean voltage for a given $`\alpha `$ and set of systematics $`𝝂`$, and $`\sigma `$ is the observed standard deviation:

$$F(x_i|\alpha ,𝝂)=\left(\frac{1}{\sigma \sqrt{2\pi }}\right)\mathrm{exp}\left[(x_i\mu (\alpha ,𝝂))^2/2\sigma ^2\right].$$
(19)
The product is

$`{\displaystyle \underset{i}{}}F(x_i|\alpha ,𝝂)`$ $`=`$ $`\left({\displaystyle \frac{1}{\sigma \sqrt{2\pi }}}\right)^N\mathrm{exp}\left[{\displaystyle \frac{1}{2\sigma ^2}}{\displaystyle \underset{i=1}{\overset{N}{}}}(x_i\mu (\alpha ,𝝂))^2\right]`$ (20)
$`=`$ $`\left({\displaystyle \frac{1}{\sigma \sqrt{2\pi e}}}\right)^N\mathrm{exp}\left[{\displaystyle \frac{(\overline{x}\mu (\alpha ,𝝂))^2}{2(\sigma /\sqrt{N})^2}}\right],`$
where $`\overline{x}`$ is the average of the voltages $`x_i`$ and $`N`$ is the effective number of uncorrelated samples (the value of $`N`$ such that $`\sigma /\sqrt{N}`$ is the standard deviation of the mean of the actual data).

The predicted mean $`\mu (\alpha ,𝝂)`$ can be replaced by Eq. 2 for $`\overline{V^D}`$ if $`\stackrel{}{f}(\stackrel{}{r}^{})`$ is identified with the theoretical value of the Yukawa force density for a test mass geometry corresponding to a particular set of parameters $`𝝂`$. As a Yukawa force density, $`\stackrel{}{f}(\stackrel{}{r}^{})`$ is linear in the parameter $`\alpha `$, the gravitational constant $`G`$, and the source mass and detector mass densities $`\rho _s`$ and $`\rho _d`$: $`\stackrel{}{f}(\stackrel{}{r}^{})=\alpha G\rho _s\rho _d\stackrel{}{\text{f}}(\stackrel{}{r}^{})`$. Substituting into Eq. 2 and factoring these terms out of the integral yields

$$\mu (\alpha ,𝝂)=\left(\alpha G\rho _s\sqrt{\rho _d}\right)\sqrt{\overline{\left|V^T\right|^2}}\left(\frac{Q}{\omega _0\sqrt{k_BT}}\right)\frac{d^3\stackrel{}{r}^{}\stackrel{}{z}^F(\stackrel{}{r}^{})\stackrel{}{\text{f}}(\stackrel{}{r}^{})}{\sqrt{d^3\stackrel{}{r}^{}\left|\stackrel{}{z}^F(\stackrel{}{r}^{})\right|^2}}.$$
(21)
Substituting Eq. 21 into Eq. 20 and the result into Eq. 17 yields the likelihood function $`L^{}(𝒙|\alpha )`$.


format: text

text:
# Classical, non-linear, internal dynamics of large, isolated, vibrationally excited molecules

## 1 Introduction

The present research was triggered by the question whether infrared chemiluminescence was possible with large molecules. This question, in turn, arose from the need to interpret the flurry of accurate IR (infrared) spectra delivered by dedicated astrophysics satellites in the past two decades. It is indeed established that the ISM (InterStellar Medium) is interspersed with large molecules or small grains (dust), with a wide size distribution (up to 1000 nm), and whose composition is essentially limited to hydrocarbons and silicates. We are interested, here, in that major family of amorphous hydrocarbons which lie so far from any illuminating star that they can only reach a low steady state temperature ($``$200 K) when the UV radiative power they absorb equals the IR radiative power they emit. Under such circumstances, the relatively high IR spectral intensities observed are incompatible with Planck’s law at that temperature. A possible alternative excitation mechanism is chemiluminescence, based on the fact that, in such environments, the dominant material, hydrogen, remains in atomic form, which notably favours interaction with hydrocarbon dust. Such interaction is believed to be the main cause of formation of molecular hydrogen, which is the main constituent of the ”clouds” that are detected within the ISM (Duley and Williams1). When a free H atom hits an H atom that is attached to a C atom at the periphery of a dust particle, they may form such a strong bond that the nascent H2 molecule can fly away, leaving a dangling bond behind. Upon a later encounter with an H atom, the latter can be captured by the dangling bond to form a long-lived C-H bond. In such an event, an energy of order 4.5 eV (the bond dissociation energy) is deposited *locally* in the molecule. This work addresses the fate of that energy and of the resulting atomic vibrations: what vibrational modes are excited and can they survive energy redistribution long enough for them to finally relax through the relatively slow process of IR emission?

Of course, IR chemiluminescence is known to occur upon reactions between small molecules (e.g. see Lee and Pimentel2). Is this possible with larger molecules as well? and how much larger? More precisely, does statistical equilibrium set in, with equal probability for all states of equal energy (also referred to as complete randomization)? If it does, a temperature can be defined, which scales like the inverse of the number of atoms; then, for a dust particle larger than, say, 100 atoms, the temperature would be so low as to preclude detectable emission of all but the very lowest vibrational frequencies, contrary to observations. If, on the contrary, complete randomization does not obtain, then the motions of the various atoms retain at least some degree of coherence and the concept of temperature becomes irrelevant. The distribution of energy among the various modes is then totally different, being determined by classical ( or quantum, for more accuracy) mechanics: even the highest vibrational frequencies can be excited to a notable extent and IR chemiluminescence becomes detectable.

The RRKM theory of unimolecular reactions (see Gilbert and Smith3) indeed assumes thermodynamical equilibrium, and is vindicated by innumerable experimental data. However, these data are mostly about yields of particular unimolar reactions (dissociations) averaged over assemblies of molecules under finite pressure, so that vibrational excitation is effected collisionally and randomly, according to equilibrium thermodynamics. This, according to Freed4 , is enough to ensure the validity of all the mathematical predictions of the RRKM theory, even assuming no randomization at all occurs within the molecule! Besides, more recent experiments have shown deviation from statistical behaviour, especially when the energy content is low (see McDonald and Marcus5). Other types of experiments, especially those designed in pursuit of ”selective chemistry”, suggest that, under suitable conditions, the molecule can be prepared in an excited *and coherent* state, which can survive for quite long times (see Zewail6). More recently, Williams and Leone7 followed the IR-*radiative* relaxation of gaseous naphtalene at very low pressures, after laser excitation of a C-H stretch vibration; they showed that their spectra are not compatible with calculations based on thermodynamic equilibrium. Also, Deak et al.8 excited the C-H stretch of nitromethane molecules in a tiny segment of a jet, and observed the prompt ($``$2 ps) establishment of a nearly uniform distribution among other vibrations, which prevailed for a time about equal to the collisional relaxation time of the bath (30 to 50 ps). Afterwards, a new, Boltzmannian, distribution clearly set in, again suggesting that collisions where responsible for this to occur. Similar conclusions were drawn by Stromberg et al.9 , working on W(CO)6. Time and frequency resolved experiments reveal in fact three cases of IVR: dissipative (energy sharing), restricted, and no IVR at all (see, for instance, Zewail10). Correspondingly, molecular motions range from statistical, ergodic, chaotic, to quasi-periodic and periodic, or regular.

Turning to theory, it is known that the normal modes of vibration of a system are orthogonal and, therefore, in principle do not interact. However, they do so in real life, because anharmonicity always provides some coupling between them, albeit weak. This coupling increases in strength as the molecule becomes larger, and the distance between modes correspondingly decreases in frequency space. There is therefore no doubt that any local disturbance is bound to decay in some way, except, in principle, in the ideal limit of the weakest excitation of the purest normal mode. This type of reasoning led Fermi, in the 30’s, to conjecture that a statistical (Boltzmannian) equilibrium must ultimately set in (see Ford11). He therefore initiated numerical simulations of the dynamics of long chains of oscillators, the outcome of which did not, however, confirm his conjecture: hence the Fermi-Ulam-Pasta paradox (Fermi12). Long afterwards it became clear that statistical intramolecular energy transfer is not ubiquitous (see Brumer13) even with many degrees of freedom; conversely, under suitable circumstances, ”chaos” can occur even with a few degrees of freedom. Thus, the size of the system and the corresponding number of normal modes are not good indicators of fast randomization 14.

Finally, a hierarchy of statisticality (stochasticity) and randomness was established, with ergodic, mixing, K and Bernouilli systems in order of increasing randomness15. Conditions favourable to chaos and randomness were recognized (see Ford11,16): existence and overlapping of non-linear resonances ($`\mathrm{\Sigma }`$nk$`\omega `$k$``$$`ϵ`$), high energy content, brief random excitations (as in a perfect gas or with billiard balls), etc.

In terms of phase space, Boltzmann’s ergodic hypothesis implies that all available phase space is uniformly visited by the system’s trajectory during the observation time. This leads to the usual *canonical prescription* for deriving time averages of physical quantities. It is found that energy is distributed among the degrees of freedom (here, the vibrational modes) according to Boltzmann’s law with a single temperature common to all. While this delivers the least biased estimates of observables *in the absence of enough information on the system*, very many well documented physical systems are not properly described in this way17. In these cases, the phase space is full of structures which induce ”intramolecular bottlenecks” to the phase flow (see Uzer18 for a review). Such systems are found to obey a *restricted* form of ergodicity, *broken ergodicity*17, meaning that, in these cases, *statisticality (randomness) obtains only in limited regions of phase space*.

For all these developments, theory can hardly predict which type of internal motion will prevail in systems with more than a few degrees of freedom. Surely, experiment remains the ultimate recourse. However, numerical models of coupled oscillators have been of great help in highlighting the role of different parameters (particularly excitation strength and bond anharmonicity), in assessing the degree of randomness, in interpreting experimental outcomes and in evaluating relaxation times (see Noid et al.19 and Thompson and co-workers20). With the advent of highly performing chemical software, numerical modeling can now more closely simulate the motion of large, real-life molecules. This is the course followed in the present work to address the problem mentioned above.

The main result of this research is as follows: in the type of structure considered here, and following the recombination of an impinging H atom with a dangling C bond, the energy deposited is ultimately, and on average, roughly equipartitioned among all the vibrational modes. However, chaos does not ensue and the molecule remains in a nearly coherent state, i.e. a coherent superposition of normal vibration modes . A Fourier analysis of the time variations of the total electric dipole moment then provides the spectral distribution of energy, which shows that chemiluminescence is indeed detectable.

The originality of this work does not lie in the method used nor in the conclusion that statistical equilibrium does not always obtain after relaxation. Rather, it is in the combination of very large molecules ($``$ 500 atoms), amorphous structure, chemical excitation, no collision at all before complete internal relaxation and very long computation runs (up to 500 ps). Such a combination, suggested by astronomical conditions, is difficult to achieve in the laboratory and, to my knowledge, was never considered before. Moreover, the outcome of this study may seem counterintuitive to many who are not familiar with the FUP paradox, and therefore may fuel interesting debate.

## 2 The chemical software

The calculations were made using the Hyperchem package released by Hypercube, Inc. in 1999, and implemented on a desk-top PC equipped with a Pentium (r)II processor (450 MHz) and MMX (TM) technology. First, the molecule is built on the screen by picking atom after atom in a table of elements. The CC bondings are specified graphically: -, =, $``$. A computational code and options are then selected. The code is then asked to optimize the molecular geometry, i.e. determine the architecture which minimizes the total potential energy, E0. A perturbation is then applied by displacing one or more atoms, or imparting a velocity to them, or creating a new bond (like in chemisorption). This brings it up to an energy state E\*. The subsequent motion is then followed using a Molecular Dynamics routine based on the Born-Oppenheimer approximation: at each calculation step, the nuclei are assumed to be clamped at their initial positions and the code determines the electronic charge distribution, either by solving Schroedinger’s equation ((a) below) or empirically ((b) below). Hence the forces acting on the nuclei and the new velocities of the latter, which determine the displacement for the next step. The step length is set at 10-16 s, much less than the shortest characteristic period of the system. To save time, the electric dipole moment is computed only once every cycle of 30 dynamics steps, still ensuring correct frequency analysis up to 3000 cm-1 .

For a frequency resolution of 30 cm-1, a Molecular Dynamics run must span 1 ps, i.e. 10000 steps or 300 cycles. However, much longer times are needed for transients to decay. The coordinates and velocity of each nucleus, and hence the potential and kinetic energies are recomputed at every step. The total energy is conserved all along, as it should be, even for the longest calculation. The energy content, or perturbation energy is Ep= E\*-E0 . The accuracy and speed of this computation depend on the choice of the code; while most state-of-the-art chemical codes are available with this package, only two of them were used in the present work:

a) the semi-empirical PM3 code, developed by J. Stewart21 along the same line as AM1. The semi-empirical methods use a rigorous quantum-mechanical formalism combined with empirical parameters obtained from comparison with experimental results. PM3 incorporates a much larger number and wider variety of experimental data than AM1. These codes compute approximate solutions of Schroedinger’s equation, and use experimental data only when the quantum mechanical calculations are too difficult or too lengthy. This sometimes makes them more accurate than poor ab initio methods, and they are always faster. PM3 also provides for a direct calculation of the normal mode spectrum, which does not require running the Molecular Dynamics. This is based on the quadratic term of the interatomic potentials and only involves matrix manipulation. The result applies, more generally, to weak excursions from the equilibrium geometry. The spectrum is discrete, i.e. composed of separate vertical bars of length equal to the ”IR intensity”, A, of the structure at the corresponding wavenumber. In real life, this must be corrected for broadening effects: thermal fluctuations, rotations, and slight differences between particles of the same family. For large particles, the first two effects are negligible in the astrophysical context.

The accuracy of the mode wavenumbers has been discussed at length22: although it depends on the type of vibration and on the structure, it is considered to be of order 10 cm-1 for the C-H stretch and 100-200 for the rest of the spectrum.

An interesting feature of this code is the possibility of animating the atoms on the screen in any selected normal mode. This is very useful in the process of tailoring the composition and structure of the particle to fit the astronomical spectra as best as possible.

This code is very much faster than ab initio codes (which are out of the question here). However, for a molecule of 100 atoms, each Molecular Dynamics step requires $``$30 s computation time, which amounts to about 100 cycles a day or 1 month for a run of minimum length. This code was used in the present work to determine, by cut-and-try, a structure capable of displaying the main IR features observed in the sky; this structure is shown in fig. 1 (Papoular23,24). Note that it is 3-D, heterogeneous, disordered, amorphous and closely resembles models of coal25 and kerogen26, a material also found in meteorites. In this study, we have used larger molecules built up by assembling randomly part or totality of fig. 1 (101 atoms) so as to retain the same chemical characters. Most of the results below refer to the largest of these molecules, designated by LM and consisting of 493 atoms. This is too large to be conveniently handled with semi-empirical codes.

b) The Molecular Mechanics code MM+. This is based on the Allinger27 MM2 code. Molecular Mechanics codes, like MM2, Amber, Charmm,… approximate the required generalized potential energy surface (force fields) by analytical functions which are deduced from, and subsequently parametrized via, accurate calculations and experimental results. The force field contains terms associated each with bond stretching, angle (in-plane) bending, torsional (out-of-plane) bending, stretch-bend cross terms, twisting, van der Waals forces and other terms of chemical significance. Anharmonicity is arbitrarily built in various terms; it is cubic for bond stretches.

MM+ uses the latest MM2 parameters and atom types (designed primarily for small organic molecules), complemented, where necessary, by specific calculations. When the ”atomic charge” option is selected, the charge on every atom is invariable and determined empirically by the type of atom and bonding. As a result, the vibration frequencies derived from MM+ calculations may be in error by up to 200 cm-1. Moreover, the absolute value of the electric dipole moment is not reliable. However, this code is very much faster than the semi-empirical codes: about 20000 cycles a day for a particle of 100 atoms. It is, therefore, highly useful in clarifying the general trends of the spectra as a function of parameters such as time, energy content, structure, etc. Such trends, which are our main interest here, are not likely to depend upon the exact numerical values of vibration frequencies in large molecules like our’s. This code was therefore used extensively for our present purposes; since it is not suited to simulating the hydrogen recombination with the molecule, the latter is approximately simulated by starting with the model particle at rest in its ground state, E0, as above . Then, instead of breaking a C-H bond, the bond is simply altered in length and direction on the screen, by cut and try, until the new potential energy of the particle is E\*$``$E0+4.5 eV. The Molecular Dynamics is launched thereafter. Under the perturbation just described, the real molecule is expected to remain in its electronic ground state. Energy is conserved to at least 3 significant digits and much better on average over a few steps.

The relevance and efficiency of force fields and semi-empirical methods are now well established in spite of their neglect of quantum effects. The rationale behind their extensive use lies in the Ehrenfest theorem (see 34 for a discussion).

## 3 The approach to steady state

The quantitative criterion adopted here for stationarity is

$$\text{R=M}^{0.5}\sigma \text{/}\text{E}_{kin}$$
(1)
where $`M=3N6`$ is the number of degrees of freedom for *N* atoms, $``$Ekin$``$ is the average total kinetic energy, and $`\sigma `$, its standard deviation. As a function of time, this tends towards a horizontal asymptote. This is shown in fig. 2 for our largest molecule (LM, 493 atoms). The relaxation is initially very fast ($`\tau `$$``$0.5 ps) but later looks more like an extended exponential (dashed curve). Ultimately, R$``$1 and all the degrees of freedom appear to have been activated.

The code allows us to follow the ”diffusion” of energy in time and space, from the perturbed bond into the molecule ; this is done by defining a circle/sphere centered on the captured H atom, ”selecting” the atoms within and asking for their total kinetic energy; this procedure is repeated for any desired radius and instant along the calculation. Concurrently, autocorrelation and FFT (Fast Fourier Transform) can be obtained from the time evolution of the position/velocity of any atom and of the length or orientation of any bond, revealing the nature of the dynamics. In this way, it is found that

1) the deposited chemical (or potential) energy is promptly transformed into stretch ($``$3000 cm-1)and bend ($``$1500 cm-1) vibrations of the central C-H bond; these motions are initially highly non-linear and carry strong harmonics;

2) spatially extended low-frequency vibrations quickly invade the whole molecule while high-frequency vibrations lag behind because they are more localized in space and preferentially exchange energy between neighbouring atoms;

3) overall, low frequencies grow at the expense of high frequencies, until all modes are activated and rough energy equipartition obtains between all atoms;

4) steady state is essentially established at about 100 ps;

One would also wish to know the energy in each normal mode. However, a normal mode analysis is hopeless because of the system size. A less rigorous, but simpler and quite helpful, approach consists in Fourier analyzing the *total* kinetic energy, Ekin, over an arbitrary period (here from 0.4 to 16 ps). Because of spectral congestion, this ”energy” spectrum is inconvenient for illustrative purposes. It is better to use the *cumulative* spectral kinetic energy distribution, $``$Ekin($`\nu `$), which is obtained by summing up the spectrum from 0 to $`\nu `$, the abscissa. The result is drawn in fig.3 for several successive times during the transient; the curves are all normalized to 1 at the highest frequency. Again because of spectral congestion, they seem to be roughly monotonous on this scale, except for the isolated C-H stretch band near 6000cm-1. Since Ekin is a quadratic form in the coordinates or momenta, the frequencies are, in fact, double the vibration frequencies. In case of a purely harmonic motion, one also expects sum and difference frequencies. These are washed out, however, in the presence of strong dephasing between all frequency components, which seems to be the case here, making for an easier interpretation.

Points 3) and 4) made above are synthetically illustrated in fig. 3. It appears that the spectrum is essentially invariable after a few picoseconds, although some ”reddening” goes on for a much longer time. It is remarkable, at this point, that the energy initially concentrated in the C-H stretch flows roughly from high to low frequencies: the main relaxation path is not direct to the ”bath” frequencies. Neither is it a strict ”vibrational cascade” down a ladder of levels: significant amounts of intermediate and low frequency energy are observed at the earliest stage; moreover, the relaxation time roughly increases with decreasing mode frequency. All this was already deduced from experiments on small molecules in liquid environment (see Deak et al. 19998, Sibert and Rey8). In the present study, the system is isolated so the energy decay of the C-H stretch modes (5000 to 6000 cm-1 in fig. 3) is only partial; as expected from anharmonic oscillators, it is accompanied by line narrowing and blue shift .

## 4 The steady state

We now take up the detail character of the motion in steady state. There are a number of different features which distinguish the various types of molecular motions (Noid et al.19). Most are associated with phase space coordinates; these are well suited to theoretical studies but can hardly be accessed in the present case. More manageable are the methods which involve looking at the frequency spectra or correlation functions of dynamical variables. As noted by Noid et al., the frequency spectrum of a dynamical variable, for quasi-periodic (regular) motion, consists of a series of sharp lines, while the chaotic frequency spectrum has broad lines forming bands near the location of the sharp lines of the corresponding regular spectrum. The transition from regular to chaotic is also characterized by the continuum between bands becoming less and less smooth and more and more rugged (grassy). In principle, for a given perturbation, every normal mode shows up, to a smaller or larger extent, in the fluctuations of position and velocity of each atom, as well as of the length of each bond, depending on the type of atom or bond and immediate environment. C-C bonds were found to carry more and stronger spectral lines than C-H bonds. Figure 4 shows the amplitude-FFT of the length of one C-C bond in the LM molecule, more than 220 ps after excitation. The time length of the data sample was 65 ps, so that the frequency resolution is 1 cm-1 (with Hanning-filtering). Most lines are indeed about this wide, and roughly retain their (simple) shape but vary in strength to some extent because they exchange energy with distant modes. Some are so close that they partly blend together. A case in point is the C-H stretch band, which extends from 2920 to 2970 cm-1 (somewhat below experimental value, because of imperfections in the model force fields). Since there are about 250 H atoms in this molecule, we expect about that number of C-H str modes, all in that range, but slightly displaced from one another because of differences of environments. Clearly, these modes are close enough to strongly interact and share energy, like *nearly* resonant coupled oscillators. This is indeed what transpires from the comparison of successive FFT’s over, say, a 16 ps interval each: the band profile changes but the integrated ”intensity” remains essentially constant. This is an example of selectivity in IVR (Felker and Zewail28). In this molecule, however, we found no case of strong resonance entailing ”thermalization” between the corresponding oscillators, as shown by Ford et al.16. By contrast, we did simulate such a resonance by exciting some of the stretches of the highly symmetrical benzene molecule, and observing the subsequent fluctuating energy distribution among the six C-H bonds (the resonant oscillators) in steady state.

Overall, spectra like the one in fig. 4 are not grassy but consist of distinct, narrow bands whose positions do not vary notably in time; the Fourier transforms of the corresponding autocorrelation functions deliver identical spectra, frequencywise. The motion can therefore be considered regular overall. This conclusion is confirmed by the autocorrelation functions: e.g., in fig. 5, for the same C-C bond length as in fig. 4; note the very strong and persistent correlations; the superimposed modulation indicates fluctuations of line intensities due to energy exchanges between oscillators. Similar results are obtained for C-H bonds. Thus, despite heavy spectral congestion, the molecule as a whole does not qualify as a thermalizing bath to any individual mode. Reasons for this are probably the amorphicity of the structure, which hinders close resonances, and the lack of a ”phonon continuum” at the lowest frequencies of the spectrum (see fig. 6 of24) as opposed to the case of molecular crystals29, for instance.

For our present purposes, it is important to note that this steady state vibrational energy distribution is not very sensitive to the details of the H atom capture or, in terms of the procedure adopted here, to the details of the initial C-H bond perturbation which is meant to simulate that capture, nor to the particular C-H bond that is perturbed, as long as the excess energy thus deposited remains of order 4.5 eV. From a more fundamental point of view, however, it is interesting to observe that, because of the non-linearity of the system, the steady state depends on that excess energy as well as on the type of bond which is perturbed (C-C, O-H, stretch, bend, torsion, etc.). None the less, even if the perturbation consists in altering the torsion angle of a C-O-H group, thus exciting one of the lowest-frequency modes, vibrational energy finds its way up the frequency ladder and most of the normal modes are ultimately activated. Figure 6 illustrates the relaxation of LM to equilibrium after the torsion angles of six such groups were altered by about 100 degrees. Note the differences between the steady state spectrum reached here (curve 6) and the one that was reached after altering the length of a C-H bond (curve 7, identical with curve 7 of fig. 3): the higher frequencies are distinctly weaker here, and the opposite is true below 2000 cm-1. The same trend is observed, but to a lesser degree, if energy is instead initially deposited in a C-C bond. This sensitivity of the final state to the initial perturbation is an indication that statistical (thermodynamic) relaxation did not take place.

It can be shown24 that the expected spectral density of IR radiation from the excited molecule is proportional to the product of the square of the frequency with the power spectral density of the dipole moment. The latter is the FFT of the sum of the autocorrelation functions of the fluctuations of the dipole moment components about their respective means. Finally, these fluctuations can be recorded during a run of the code well after steady state has settled. This procedure was applied to the steady state corresponding to fig. 3 (curve 7) and fig. 4, and the result is displayed in fig. 7. This shows that only a fraction of the vibrational modes is ”infrared active”, but one can still recognize the C-H stretch bands (aromatic and aliphatic, around 3000 cm-1), the in-plane C-H bends ($``$ 1000-1200 cm-1) and the C-H out-of-plane bends ($``$ 750-1000 cm-1), and far infrared bands mainly associated with oxygen atoms in various environments23. It should be stressed again that the present computations used the MM+ molecular dynamics software, which is known to yield only approximate frequencies and intensities. A more accurate spectrum can be obtained using the PM3 semi-empirical software, albeit on a smaller molecule (100 atoms), but still representative of the present amorphous structure24.

## 5 Chaotic motions

The conclusions drawn above are strengthened by comparison with conspicuously chaotic motions. Such motions are expected when the molecule is immersed in a thermal bath. The Hyperchem package allows the simulation of this situation by means of Langevin dynamics, in which the usual equations of motion are supplemented by terms representing both the friction against the bath molecules and the random changes in position and velocity upon collisions with these molecules (Allen and Tildesley30). The user may specify the temperature, T, and the friction coefficient, $`\gamma `$; the stronger the latter the tighter the coupling of the molecule to the bath. For purposes of comparison with fig. 5, we set $`\gamma `$=100 ps-1 and T=35 K, a temperature at which the energy content of the molecule is about the same as above. Figure 8a shows, under these conditions, the autocorrelation function of the length of the same C-C bond as in fig. 5. The difference is striking: very weak correlation, with no regular pattern.

Theoretical studies also show that chaos is bound to set in beyond a total energy threshold, even in an isolated system. To explore this route to chaos, we first ran Langevin dynamics at 3000 K for 0.3 ps. This dumped an energy of 12.5 eV uniformly in the molecule.The latter was then isolated from the bath and molecular dynamics was run in vacuum for 16.5 ps. The corresponding correlation function is shown as curve b in fig. 8, again for the same C-C bond length as above. The weak correlations and highly disordered modulation indicate the partial onset of chaos, but this cannot be carried much further because the molecules dissociate at higher excitation energies. Chaos is much easier to induce that way in a small and/or symmetric structure like benzene; I confirmed this using the same code.

The frequency spectra of these chaotic motions are also distinctly different from the regular spectra: they are very grassy and individual modes are no longer recognizable because they blend together, even at the highest resolution. In the extreme case of Langevin dynamics, all features are drowned in the grass which extends to much higher frequencies. This is apparent from the cumulative spectral kinetic energy distributions in fig. 9, where curves 9a,b correspond, respectively, to 8a,b while 9c is taken from fig. 3, the steady state regular motion. In such graphs, of course, the grass is not apparent because it is integrated out, so that the progression to chaos from c to a is made clearer.

## 6 Conclusion

Given the size and structure of our system, and the type and strength of the activation process, there is no theoretical way to determine, beforehand, the nature of the subsequent relaxation. The computations performed in this study help drawing a realistic classical model of the dynamics to be expected.

Many markers point to a regular (quasi-periodic) motion in steady state: the steady state depends on the initial perturbation (although not wildly); the spectral lines are narrow; correlations are strong and long lived; the spectral energy distribution (energy spectrum) is very different than for a molecule in a (model) thermal bath; in particular, there is no continuum here. However, this steady state is not *perfectly* coherent: the spectral lines vary *erratically* in intensity (although their short term averages are constant), i.e. there are *random*, reversible, energy exchanges between modes; and the combination frequencies are absent from the energy spectrum, i.e. there are no apparent phase relations between modes.

These results are all understandable in terms of more or less strongly and non-linearly coupled anharmonic oscillators (see 11,16 and references therein). They fit neatly in the picture of partial chaos or broken ergodicity17: the phase point can wander randomly, but only in restricted regions of the energetically available phase space. Such situations occur when the vibrational spectrum is too sparse and/or the energy content is too low for multiple overlaps to occur between non-linear resonances. In the present instance, this is likely due to the limited size, and to the amorphicity and inhomogeneity of the structure.

From comparisons between corresponding classical and quantal motions (see Noid et al.19), it appears that quantum chaos does not set in before classical chaos is well established. It seems, therefore, that, even under the more realistic, but ”corresponding”, laws of quantum mechanics, the vibrational motions of our activated molecule will not reach statistical (thermodynamic) equilibrium. If that is the case, then the classical spectral intensity (fig. 6) may be interpreted as the ensemble average quantal energy content at the same frequency. One can therefore expect that such molecules, isolated and chemically excited, will ultimately deliver a significant part of their initial energy content in the form of near and mid-IR radiation.

The question of the upper size limit beyond which these conclusions no longer apply remains open to discussion. It is known that the vibrational energy deposited in a molecular *crystals* will ultimately flow into the *phonon bath*29,31. Phonons are very low-frequency ($`<`$100 cm-1) plane propagating waves which can transport energy through the structure. They can interact upon hitting a structural or chemical defect, and exchange energy like perfect gas particles in a vessel, so that thermal equilibrium can set in. But, for the localized vibrations (vibrons), which we have been considering above, to be efficiently converted into phonons, the spectrum of the latter must be dense enough as to form a quasi-continuum to which the Fermi Golden Rule of transitions can be applied. This is only possible if the structure is regular and infinite, like ideal crystals; very irregular structures obviously can hardly support phonons, a fact which translates into very low thermal conductivities and long thermal relaxation times32. Indeed, the relaxation times of some glasses at low temperatures have been found to range in the hours33. It is likely, therefore, that a grain of our model material will have to be quite large before it can give rise to a phonon spectrum of sufficient density that it constitutes an efficient thermal bath.

1. W. Duley and D. Williams, *Interstellar Chemistry* (Academic Press, London, 1984).

2. Y.-P. Lee and G. Pimentel, J. Chem. Phys. $`\mathrm{𝟔𝟗}`$, 3063 (1978).

3. R. Gilbert and S. Smith, *Theory of Unimolecular Reactions* (Blackwell Scientific Publ., London, 1990).

4. K. Freed, Faraday Disc. Roy. Chem. Soc. $`\mathrm{𝟔𝟕}`$, 231 (1979).

5. J. McDonald and R. Marcus, J. Chem. Phys. $`\mathrm{𝟔𝟓}`$, 2180 (1976).

6. A. Zewail, Acc. Chem. Res. $`\mathrm{𝟏𝟑}`$, 360 (1980).

7. R. Williams and S. Leone, Astrophys. J. $`\mathrm{𝟒𝟒𝟑}`$, 675 (1995).

8. J. Deak, L. Iwaki and D. Dlott, J. Phys. Chem. A $`\mathrm{𝟏𝟎𝟑}`$, 971 (1999); E. Sibert III and R. Rey, J. Chem. Phys. $`\mathrm{𝟏𝟏𝟔}`$, 237 (2002).

9. C. Stromberg D. Myers and M. Fayer, J. Chem. Phys. $`\mathrm{𝟏𝟏𝟔}`$,3540 (2002).

10. A. Zewail, in *The Nobel Prizes* (Almqvist and Wiksell, Stockholm, 1999).

11. J. Ford, Phys. Rep. $`\mathrm{𝟐𝟏𝟑}`$, 271 (1992).

12. E. Fermi, *Collected Papers*, Vol. II (U. Chicago Press, 1965), p. 978.

13. P. Brumer, Adv. Chem. Phys. $`\mathrm{𝟒𝟕}`$, 201 (1981).

14. D. Nesbitt and R. Field, J. Chem. Phys. $`\mathrm{𝟏𝟎𝟎}`$, 12735 (1996).

15. J. Lebowitz and O. Penrose, Phys. Today $`\mathrm{𝟐𝟔}`$, 23 (1973).

16. J. Ford and J. Waters, J. Math. Phys. $`\mathrm{𝟒}`$, 1293 (1963); J. Ford and G. Lunsford, Phys. Rev. $`\mathrm{𝐀𝟏}`$, 59 (1970).

17. R. Palmer, Adv. Phys. $`\mathrm{𝟑𝟏}`$, 669 (1982).

18. T. Uzer, Phys. Reports 199, $`\mathrm{𝟕𝟑}`$ (1990).

19. D. Noid, M. Koszykowski and R. Marcus, Ann. Rev. Phys. Chem. $`\mathrm{𝟑𝟐}`$, 267 (1981).

20. T. Sewell, D. Thompson and R. Levine, J. Phys. Chem. $`\mathrm{𝟗𝟔}`$, 8006 (1992); B. Sumpter and D. Thompson, J. Chem. Phys. $`\mathrm{𝟖𝟔}`$, 2805 (1987) and $`\mathrm{𝟖𝟐}`$, 4557 (1985).

21. J. Stewart, J. Am. Chem. Soc. $`\mathrm{𝟏𝟎𝟕}`$, 3902 (1985).

22. D. Seeger, D. Korzeniewski, W. Kowalchyk, J. Phys. Chem. $`\mathrm{𝟗𝟓}`$, 6871 (1991).

23. R. Papoular, Astron. Astrophys. $`\mathrm{𝟑𝟔𝟐}`$, L9 (2000).

24. R. Papoular, Spectrochem. Acta, Part $`\mathrm{𝐀𝟓𝟕}`$, 947 (2001).

25. J. Speight, Appl. Spectr. Rev. $`\mathrm{𝟐𝟗}`$, 117 (1994).

26. F. Behar and M. Vandenbroucke, Rev. Inst. Fr. Petrole $`\mathrm{𝟒𝟏}`$, 173 (1986).

27. N. Allinger, J. Am. Chem. Soc. $`\mathrm{𝟗𝟗}`$, 827 (1977).

28. P. Felker and A. Zewail, Phys. Rev. Lett. $`\mathrm{𝟓𝟑}`$, 501 (1984).

29. D. Dlott, in *Laser Spectroscopy of Solids II*, ed. W. Yen, Topics in Applied Physics vol.65 (Springer-Verlag, 1989), p.167.

30. M. Allen and D. Tildesley, *Computer Simulation of Liquids* (Clarendon Press, Oxford, 1987), p. 261.

31. V. Kenkre, A. Tomakoff and M. Fayer, J. Chem.Phys. $`\mathrm{𝟏𝟎𝟏}`$, 10618 (1994).

32. G. Srivastava, *The physics of phonons* (Adam & Hilger,1990).

33. E. Donth, The Glass Transition (Springer, 2001).

34. M. Ben-Nun and R. Levine, J. Chem. Phys. $`\mathrm{𝟏𝟎𝟓}`$, 8136 (1996).


text:
# Identified particles in Au+Au collisions at √𝑠__{𝑁𝑁} = 200 GeV

## 1 INTRODUCTION

In the PHOBOS experiment, particles are measured and identified in a two-arm spectrometer  consisting of layers of Si detectors placed on either side of the beam axis in a 2 Tesla magnetic field. Unique features of the PHOBOS spectrometer, such as a close proximity of the sensitive detector layers to the interaction region, little material between the collision vertex and Si layers and high segmentation of Si detectors, permit precise measurements of particles down to very small transverse momenta.

In this paper we present the measured ratios of charged antiparticles to particles for pions, kaons and protons produced near mid-rapidity. These measurements provide information on properties of the system at the chemical freeze-out point and on baryon transport processes. The invariant yields of charged particles at very low transverse momenta are also shown. The results of these measurements are sensitive to large-volume physics phenomena and to the effects of collective radial expansion of the system. An enhancement of the production of very low transverse momentum particles could indicate new long distance scale physics in heavy ion collisions .

## 2 ANTIPARTICLE TO PARTICLE RATIOS

The ratio of negatively to positively charged particles was measured during the RHIC 2001 Run using spectrometer data taken with both polarities of the magnetic field. Event triggering and determination of the collision centrality was provided by two sets of scintillator paddle counters . The rapidity coverage for the ratio measurements extends from about 0.2 to 0.8 for kaons and protons and 0.35 to 1.3 for pions. More details on the tracking, particle identification and details of the ratio measurements can be found elsewhere .

The fully corrected ratios measured within our rapidity acceptance for the 12% most central Au+Au collisions at 200 GeV are : $`<\pi ^{}>/<\pi ^+>=1.025\pm 0.006(stat.)\pm 0.018(syst.)`$, $`<K^{}>/<K^+>=0.95\pm 0.03(stat.)\pm 0.03(syst.)`$ and $`<\overline{p}>/<p>=0.73\pm 0.02(stat.)\pm 0.03(syst.)`$. The $`<\overline{p}>/<p>`$ ratio increases by about 25% over the value at $`\sqrt{s_{_{NN}}}`$ = 130 GeV , indicating a rapidly decreasing net-baryon density near mid-rapidity. The estimated baryochemical potential ($`\mu _B`$) for the system formed in central Au+Au collisions at 200 GeV is $`\mu _B=27\pm 2`$ MeV .

Figure 1 shows preliminary results for the dependence of antiparticle to particle ratios on transverse momentum ($`p_T`$) for the 10% most central collisions at the top RHIC energy. Within our current estimates of the systematic uncertainties, the ratios are independent of $`p_T`$ in the measured rapidity range. The centrality dependence of the measured ratios is shown in Fig. 2. The pion and kaon ratios are consistent with a constant value over the measured centrality range. For the $`<\overline{p}>/<p>`$ ratio, a weak but systematic drop with $`<N_{part}>`$ can be observed. However, within our uncertainties, the ratio is also consistent with a constant value, independent of centrality.

## 3 YIELDS AT VERY LOW TRANSVERSE MOMENTA

To extend our identified particle measurements to lower transverse momenta, we have searched for particles which stop in the fifth Si plane of the spectrometer. The reconstruction procedure is based on the analysis of tracks with large energy depositions in the first five spectrometer planes which are located in the field free region of the spectrometer. To determine the particle mass we check the mass and momentum hypotheses by making cuts on the energy deposited per unit length in every plane and the total deposited energy obtained by summing all the energy depositions. The analysis procedure was tested successfully on samples of simulated low momentum pions, kaons and protons. The measured raw yields were corrected for acceptance and efficiency. The correction factors were obtained by embedding simulated low momentum particles into real data events. Additional corrections, including feed-down from weak decays and contributions from secondary, misidentified and ghost particles, were also applied. Our current estimates of the systematic uncertainties are $`\pm `$20% for pions, $`\pm `$40% for kaons and $`\pm `$50% for protons and include systematic effects in the measured yields and in the estimated correction factors.

Figure 3 shows preliminary results for invariant yields of ($`\pi ^\pm `$), ($`K^\pm `$) and ($`p+\overline{p}`$) measured at mid-rapidity in the transverse momentum ranges from 30 to 50 MeV/c for charged pions, 90 to 130 MeV/c for kaons and 140 to 210 MeV/c for protons and antiprotons for the 15% most central Au+Au collisions at $`\sqrt{s_{_{NN}}}`$ = 200 GeV . Our results are compared in Fig.4 to the predictions of HIJING, RQMD  and single freeze-out model of Broniowski and Florkowski . Also shown in Fig. 4 are PHENIX measurements at 130 GeV . All data and the model predictions shown in Fig.4 have been corrected for feed-down from weak decays. It is evident that our measurements explored the range of very low transverse momenta, below the current measurements from other RHIC experiments (see also ). The yields of ($`p+\overline{p}`$) seem to agree with the extrapolation of the PHENIX measurements at higher $`p_T`$, indicating that the flattening of proton spectra attributed to the collective radial expansion extends down to very low $`p_T`$. Comparison with the model predictions also shows that none of the models presented in Fig.4 can consistently describe the pion, kaon and proton data. Particularly striking differences between the model predictions can be seen for the ($`p+\overline{p}`$) yields where the models differ by a factor of 2 to 6.

## 4 SUMMARY

The PHOBOS Collaboration has measured the yields of identified particles in Au+Au collisions at $`\sqrt{s_{_{NN}}}`$ = 200 GeV. The $`<\overline{p}>/<p>`$ ratio reaches a value of 0.73 indicating that only one quarter of the protons at mid-rapidity have been transported from the beam rapidity regions. The estimated baryochemical potential is a factor of 2 smaller at 200 GeV than at 130 GeV. The antiparticle to particle ratios are consistent with the constant values as a function of transverse momentum and collision centrality.

At very low transverse momenta, we see no enhancement in the particle yields. On the contrary, particle production, particularly for protons and antiprotons, seems to be suppressed, which can result from a rapid transverse expansion of the system. These measurements provide constraints on models and suggest that dynamical processes (rescattering, expansion) play an important role in the description of baryon spectra at low transverse momenta.

Acknowledgments: This work was partially supported by U.S. DOE grants DE-AC02-98CH10886, DE-FG02-93ER40802, DE-FC02-94ER40818, DE-FG02-94ER40865, DE-FG02-99ER41099, and W-31-109-ENG-38 as well as NSF grants 9603486, 9722606 and 0072204. The Polish group was partially supported by KBN grant 2-P03B-10323. The NCU group was partially supported by NSC of Taiwan under contract NSC 89-2112-M-008-024.