text:
# 1 Schematic description of the net-works of the chromoelectric flux tubes in case of the tetra-quark (right) with one junction 𝐽 and one anti-junction 𝐽̄ and the penta-quark (left) with two junctions and one anti-junction. Three fluxes belonging in the fundamental 𝟑 color representation are flowing into a junction J and three such fluxes are flowing out of an anti-junction 𝐽̄. Any outgoing/incoming flux in color 𝟑 is equivalent to an outgoing/incoming flux in color 𝟑̄ representation. The fluxes are locally coupled to a singlet via 𝜀_{𝑎𝑏𝑐}Φ_𝑎Φ_𝑏Φ_𝑐 with 𝑎,𝑏,𝑐 the three colors. The drawing shows an idealized case where the length 𝑎 of the flux tubes joining 𝐽 and 𝐽̄ are significantly larger than the radius 𝑏 of tube.

Are there non-strange low-lying penta-quarks and can we understand their width

R.W. Gothea and S. Nussinova,b

aDepartment of Physics & Astronomy, University of South Carolina, Columbia, USA
bSchool of Physics & Astronomy, University of Tel Aviv, Ramat Aviv Tel-Aviv, Israel

Abstract

We argue that the lightest isospin 1/2 partners of the $`Z^+(1530)`$ $`\overline{s}uudd`$ penta-quark predicted by Diakonov, Petrov and Polyakov are not the $`N(1710)`$ mixed anti-decuplet states, but the pure non-strange $`\overline{u}(ud)(ud)`$ and $`\overline{d}(ud)(ud)`$ penta-quark states which may lie as low as $`1200MeV`$. The expected low width of a few $`MeV`$ of such a putative state may explain why it was missed in phase shift analyzes of pion-nucleon scattering.

Exotic tetra-, penta- and hexa-quarks have been discussed in the frame-work of QCD and various approximations thereof for more than three decades. The generally accepted point of view was that such states which certainly should exist are numerous, broad and hence blend in a continues background. The anomalous lightness of the Nambu-Goldstone pions allows most exotics to decay into a stable baryon and pions or just into pions with large phase spaces in all cases. It was also argued that decay rates of non exotics as $`q\overline{q}`$ mesons or $`qqq`$ baryons requiring the creation of an extra and often even specific $`q\overline{q}`$ pair are suppressed  relative to the ”fall-apart” decay of exotics like the $`\mathrm{\Theta }^+\overline{s}uudd\overline{s}u+uddK^+n`$, where all final quarks are already present in the initial state.
If pionic decays are partially blocked by the Zweig rule, suppressing $`s\overline{s}`$ annihilation into the lighter $`u\overline{u}`$ or $`d\overline{d}`$ quarks, the relevant states like the $`I=0`$ $`f(980)`$ and $`I=1`$ $`a(980)`$, which following Jaffe  we take as $`s\overline{s}q\overline{q}`$ 111$`q`$ refers to the lightest quarks $`(u,d)`$, $`Q`$ to heavy quarks $`(c,b)`$ and the intermediate $`s`$ is mentioned explicitly. configurations, can have reasonable $`50100MeV`$ widths despite of the large phase space for pionic decays. It was then suggested that in certain heavier quark systems, such as the $`\overline{c}suud`$ penta-quark  or $`ssuudd`$ $`\mathrm{\Lambda }\mathrm{\Lambda }`$-hexa-quark  states, the favorable strong hyperfine interaction can generate stable multi-quark states which decay only weakly. It has been noted by several authors  that tetra-quarks, particularly the $`cc\overline{q}\overline{q}`$ or $`c\overline{c}q\overline{q}`$, are even more likely to be discovered. This was indeed verified in the Belle experiment , where a remarkable narrow peak was found in the $`J/\mathrm{\Psi }\pi ^+\pi ^{}`$ channel. This state which we believe has the quantum numbers of a $`D^{}\overline{D}`$ in $`s`$-wave, namely $`J^\pi =1^+`$, forms readily in the $`B`$ decay involving two charm quarks in close spatial proximity and with reasonably low relative momenta. Its annihilation into $`J/\mathrm{\Psi }`$ and pions is still suppressed by the need to have the $`c\overline{c}`$ in the $`J/\mathrm{\Psi }`$ configuration causing its remarkable narrowness. This may be the first of a host of other $`Q\overline{Q}q^{}\overline{q}`$ exotica with $`q^{},q\{u,d,s\}`$. If the Babar $`D_s(2317)`$ state  is a $`c\overline{s}(u\overline{u}+d\overline{d})/\sqrt{2}`$ rather than $`p`$-wave $`c\overline{s}`$ state , QCD inequalities  suggest the existence of a third member in the series  starting with $`f(980)s\overline{s}q\overline{q}`$, namely the $`I=0`$ $`c\overline{c}q\overline{q}`$ state with a mass $`m3670MeV`$ and a narrow width decaying into $`\eta _c\eta `$.
A $`\overline{s}uudd`$ penta-quark resonance at $`m(\mathrm{\Theta }^+)=1540MeV`$ has been seen in $`K^+n`$ invariant mass distributions in various real photon experiments off deuterons  and protons  by analyzing the $`K^{}K^+n`$ and $`K_s^0K^+n`$ final states without being swamped by multiple pion complex final states. This resonance has also been reported in $`K^+`$-Xenon scattering  and the possible existence of a very weak evidence in $`K^+d`$ cross section data in the PDG (Particle Data Group compilation)  was noted .
The most remarkable feature of the new state is the narrow width $`\mathrm{\Gamma }_{\mathrm{\Theta }+}<25MeV`$ seen in $`\gamma d`$ and $`\gamma p`$ experiments where this upper bound is given by the experimental resolution. A much more stringent bound of $`\mathrm{\Gamma }<36MeV`$  follows from the lack of a prominent enhancement of the $`K^+d`$ scattering cross section in the relevant $`K^+d`$ momentum interval whose size is fixed by the Fermi motion in the deuteron. An even stronger upper bound of $`\mathrm{\Gamma }<1MeV`$ is derived from a phase shift analysis  and by analyzing along similar lines the $`K^+`$ charge exchange reaction .
Motivated by this new development we would like to suggest schemes where the width can be as small as $`\mathrm{\Gamma }<13MeV`$ and to dispel the belief that, when decays into pseudo scalar particles are allowed and the Zweig rule does not apply, the exotics are always extremely broad. Thus a very different explanation why multi-quark exotics have not been seen before emerges. Their widths are no longer too large but too small, causing the production cross sections for these exotic states scaled by their widths $`\mathrm{\Gamma }`$ to be so tiny that the peaks in the invariant mass distribution escaped detection in earlier lower statistics and/or lower resolution experiments.
A $`Z^+(1530)`$ state with a width of $`515MeV`$, $`I=0`$ and $`J^\pi =\frac{1}{2}^+`$ has been predicted in an extended SU(3)-flavor Skyrme model . It appears as the isospin singlet tip of an anti-decuplet where the corresponding $`I=\frac{1}{2}`$ doublet has been identified with the $`N(1710)`$ $`J^\pi =\frac{1}{2}^+`$ state. Diakonov, Petrov and Polyakov also suggest that the $`\mathrm{\Sigma }(1890)`$ and $`\mathrm{\Xi }(2070)`$ states are the remaining penta-quark states of this anti-decuplet, leaving only the $`Z^+`$ to be discovered.
With the $`Z^+(1530)`$ and the corresponding measured $`\mathrm{\Theta }^+(1540)`$ at hand we should anchor the penta-quark scale at $`1540MeV`$ and look for other even lighter penta-quark states comprising $`u,d`$ and $`s`$ quarks and anti-quarks.
In the anti-decuplet the $`S=0`$ $`N(1710)`$ is obtained from the $`Z^+`$ by the U spin lowering operator, that replaces in $`Z^+\overline{s}uudd`$ either one $`d`$ by a $`s`$ or one $`\overline{s}`$ by a $`\overline{d}`$, yielding

$$|N,I=\frac{1}{2},I_z=\frac{1}{2}>=\sqrt{\frac{1}{3}}|\overline{d}uudd>+\sqrt{\frac{2}{3}}|\overline{s}uuds>.$$
(1)
The matrix element of the SU(3) breaking Lagrangian $`m_s\overline{\mathrm{\Psi }}_s\mathrm{\Psi }_sm_d\overline{\mathrm{\Psi }}_d\mathrm{\Psi }_d`$ in the $`N(1710)`$ state is $`\frac{4}{3}\mathrm{\Delta }_m^{penta}`$ and $`\mathrm{\Delta }_m^{penta}`$ in the case of $`Z^+(1530)`$. $`\mathrm{\Delta }_m^{penta}`$ is an effective $`s`$ and $`d`$ quark mass difference which subsumes also the hyperfine mass splitting effects. The mass differences between these $`I=0`$ and $`I=\frac{1}{2}`$ penta-quark states is then

$$m(N(1710))m(Z^+(1530))=\frac{4}{3}\mathrm{\Delta }_m^{penta}\mathrm{\Delta }_m^{penta}=\frac{1}{3}\mathrm{\Delta }_m^{penta}=180MeV.$$
(2)
This large effective strange versus up or down quark mass difference in the penta-quark system $`\mathrm{\Delta }_m^{penta}=540MeV`$ is more than three times larger than the difference in the standard baryonic decuplet $`\mathrm{\Delta }_m^{baryon}=160MeV`$ that is traditionally identified with the constituent mass difference between the strange $`s`$ and non-strange $`u,d`$ quarks. As we indicate next a large $`\mathrm{\Delta }_m^{penta}=200400MeV`$ could be a better prediction reflecting the large hyperfine splittings in case of the $`q_i\overline{q}_j`$ Nambu-Goldstone pions associated with the spontaneous chiral symmetry breaking which the Skyrme model incorporates.
General arguments based on QCD inequalities motivate

$$\mathrm{\Delta }_m^x=m(X\overline{s})m(X\overline{q})>m_s^cm_q^c=\mathrm{\Delta }_m^c.$$
(3)
This conjecture constituting a stronger variant of Vaffa-Witten’s rigorous result  will be addressed separately in detail. The masses on the left are those of the lowest lying states with given $`J^\pi `$ consisting of the same subsystem $`X`$ with an extra $`\overline{s}`$ or $`\overline{q}`$ quark in the same overall state and $`m_{q,s}^c`$ are the current quark masses. $`\mathrm{\Delta }_m^c`$ is smaller than the standard $`\mathrm{\Delta }_m^{baryon}`$ of $`160MeV`$ in common quark models. It has been noted  that $`\mathrm{\Delta }_m=m(Q\overline{s})m(Q\overline{q})`$ monotonically decreases with $`m_Q`$. This is expected from the $`1/m_Q`$ decrease of the attractive hyperfine interaction.

$$<\mathrm{\Psi }_H|_{hfi}|\mathrm{\Psi }_H>=\underset{i,j}{}C_{i,j}(H)\frac{(\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j)(\stackrel{}{\lambda }_i\stackrel{}{\lambda }_j)}{m_{q_i}m_{q_j}}$$
(4)
The sum extends over all the quarks and anti-quarks in the hadron $`H`$ and $`C_{i,j}(H)`$ is the relative wave-functions at zero separation of the various $`q_iq_j`$ or $`q_i\overline{q}_j`$ pairs in the hadron $`H`$. The components of $`\stackrel{}{\sigma }`$ are the $`2\times 2`$ Pauli matrices $`\sigma _\alpha `$ with $`\alpha \{1\mathrm{}3\}`$ and those of $`\stackrel{}{\lambda }`$ are the $`3\times 3`$ Gelman matrices $`\lambda _\beta `$ with $`\beta \{1\mathrm{}8\}`$ representing the eight SU(3) color generators in the $`\mathrm{𝟑}`$ and $`\overline{\mathrm{𝟑}}`$ basis of quarks and anti-quarks. For pseudo-scalar hadrons this differences $`\mathrm{\Delta }_m=m(Q\overline{s})m(Q\overline{q})`$ starting for light quarks with $`m(K)m(\pi )360MeV`$ decrease all the way down to $`m(B_s)m(B_u)90MeV`$ for $`b`$ quarks. Nussinov and Shrock have conjectured  that indeed the latter mass difference more correctly represents the current quark mass difference $`\mathrm{\Delta }_m^c`$ in agreement with lattice results . One has to correct for running quark masses  since, due to the increased reduced mass, the $`Q\overline{q}`$ state is smaller than the corresponding $`q\overline{q}`$ state. The empirical strange non-strange $`\mathrm{\Delta }_m^{baryon}`$ values of $`150250MeV`$ can now be explained by this smaller $`\mathrm{\Delta }_m^c`$ and hyperfine interactions using constituent quark masses similar to those in the mesons and the fact that the hyperfine splitting in baryons is reduced relative to that in mesons due to $`<baryon|\stackrel{}{\lambda }_i\stackrel{}{\lambda }_j|baryon>=1/2<meson|\stackrel{}{\lambda }_i\stackrel{}{\overline{\lambda }}_j|meson>`$.
Next we would like to identify the likely non-strange analogue of the $`\mathrm{\Theta }^+(1540)`$ penta-quark and estimate in the same manner $`\mathrm{\Delta }_m^{penta}`$.
In many multiplets like the vector and tensor meson nonets the physical mass eigenstates are adequately represented by the $`I=0`$ non-strange $`|(u\overline{u}+d\overline{d})/\sqrt{2}>`$ and strange $`|s\overline{s}>`$ quark states and are far from the SU(3) flavor singlet $`|(u\overline{u}+d\overline{d}+s\overline{s})/\sqrt{3}>`$ and $`|(u\overline{u}+d\overline{d}2s\overline{s})/\sqrt{6}>`$ octet state. This is due to the fact, manifested in the Zweig rule, that the $`s\overline{s}gluonsu\overline{u}`$ off diagonal matrix elements are smaller than the mass difference between the two strange and two non-strange ($`u`$ or $`d`$) quarks.
An important exception to this are the mesons $`\pi `$, $`K`$, $`\eta `$ and $`\eta ^{}`$ of the pseudo-scalar octet with a strong coupling to the pure glue channel, often attributed to the U(1) axial anomaly. The Skyrme model contains the U(1) axial anomaly and the Skyrmions can be viewed as a coherent $`\pi ,K,\eta `$ fields. Hence the emergence in first order of the SU(3) symmetric anti-decuplet of Diakonov, Petrov and Polyakov can be understood . Note however that a detailed analysis of the radiative widths  suggests that also the $`\eta `$ and $`\eta ^{}`$ are intermediate states between the mixed and pure SU(3) flavor states.
We explore the possibility that the lowest lying partners of the $`\mathrm{\Theta }^+(1540)`$ are essentially the pure non-strange $`\overline{u}(ud)(ud)`$ and $`\overline{d}(ud)(ud)`$ penta-quark states $`P^0`$ and $`P^+`$.
In the penta-quark state we encounter for the first time two di-quarks taken in the simplest model  to be the standard strongly bound $`I=0`$ $`S=0`$ $`(u_1d_1)`$ and $`(u_2d_2)`$ with an relative angular momentum $`L_{12}=1`$. The task of estimating the difference between the analog states $`m(\mathrm{\Theta }^+)m(P)`$ where $`P`$ is obtained by exchanging $`\overline{s}`$ in $`\mathrm{\Theta }^+\overline{s}(u_1d_1)(u_2d_2)`$ by $`\overline{u}`$ or $`\overline{d}`$ seems easier,

$$m(\mathrm{\Theta }^+)m(P)=\mathrm{\Delta }_m^c\underset{q\{u,d\},i\{1,2\}}{}C_{\overline{s},q_i}\frac{(\stackrel{}{\sigma }\stackrel{}{\sigma }_i)(\stackrel{}{\overline{\lambda }}\stackrel{}{\lambda }_i)}{m_sm_q}+\underset{q\{u,d\},i\{1,2\}}{}C_{\overline{q},q_i}\frac{(\stackrel{}{\sigma }\stackrel{}{\sigma }_i)(\stackrel{}{\overline{\lambda }}\stackrel{}{\lambda }_i)}{m_q^2}$$
(5)
where we sum only over the hyperfine interactions between the anti-quark and the four quarks in the penta-quark state because the remaining mutual $`q_iq_j`$ interactions cancels out. Since each of the di-quarks $`u_id_i`$ has $`S=0`$ the expectation value of each of the four hyperfine interactions vanishes. In this lowest order approximation we thus find $`m(\mathrm{\Theta }^+)m(P)\mathrm{\Delta }_m^c100MeV`$. But this result has to be questioned because the penta-quark state with two rigid $`I=0`$ $`S=0`$ $`ud`$ di-quarks is only justified if the remaining anti-quark is a heavy $`c`$ or $`b`$ quark with small mutual $`\overline{Q}q_i`$ hyperfine interactions. This is no longer the case when the anti-quark is a $`\overline{u}`$ or $`\overline{d}`$ because lower energies are obtainable if a $`S=1`$ admixture in the di-quark is allowed while benefiting from the very large $`\overline{q}q`$ hyperfine interaction with $`S_{\overline{q}q}=0`$. The scale for the strength of this hyperfine interaction is set by $`m(K_{\overline{s}q})m(\pi _{\overline{q}q})360MeV`$. This as well as the large mass splitting in the Diakonov, Petrov and Polyakov anti-decuplet widens the likely $`m(\mathrm{\Theta }^+)m(P)`$ range from $`100MeV`$ up to $`350MeV`$.
Such a non-strange penta-quark $`P`$, with a mass of $`12001450MeV`$ and $`I=1/2`$ could be seen in $`\pi ^{}p`$ but not in $`\pi ^+p`$ system. A new resonance cannot be ruled out by $`\pi ^{}p`$ phase shift analyzes, if the resonance is as narrow as $`\mathrm{\Gamma }<5MeV`$ and its cross section small 222The total formation cross section at the maximum of the $`P`$ resonance in the $`\pi ^{}p`$ channel of $`525mb`$ is smaller than the one in the $`K^+n`$ channel at $`W(m_{\mathrm{\Theta }^+})`$ of 37 mb .. The full circle in the Argand diagram is completed within an interval $`\mathrm{\Delta }_W2\mathrm{\Gamma }`$ requiring a high energy resolution and a small step width of less than $`0.5MeV`$ to measure the Argand diagram, not taking into account that even in $`0.5MeV`$ $`W`$ bins the total $`P`$ to $`\pi ^{}p`$ cross section ratio is only $`0.10.3`$. In one of the few high resolution experiments narrow peaks in $`\pi ^{}p`$ channel have indeed been seen in the reaction $`{}_{}{}^{12}C(e,e^{}p\pi ^{})^{11}C`$ at MAMI . The high missing mass resolution of $`\sigma _m=0.27MeV`$ allows to identify two $`4MeV`$ narrow states at an invariant mass of $`1222MeV`$ and $`1236MeV`$ that have been interpreted as bound $`\mathrm{\Delta }^0`$ states. Interestingly a narrow resonance with a mass of $`1225MeV`$ and a width of $`50keV`$ was found to diminish the $`\chi ^2`$ of partial wave analyzes based on the $`\pi N`$ SAID data base .
The most puzzling aspect of the $`\mathrm{\Theta }^+(1540)`$ is the narrow width $`\mathrm{\Gamma }<𝒪(13MeV)`$ established by $`K^+`$ nuclear scattering data. If the $`\mathrm{\Theta }^+(1540)`$ will survive further experimental scrutiny and even more so if lighter narrow analogues exist, then finding a credible scenario for such narrow widths is a challenge to QCD, which should at least qualitatively explain all hadronic data. The penta- and tetra-quark are more complex than ordinary baryons and mesons. In a chromoelectric flux tube model (CFTM) mesons are tubes with the radius b of chromoelectric flux extending between the quark and anti-quark at both ends. It generates a confining linear potential $`V=\sigma r_{q\overline{q}}`$ where $`\sigma `$ is the string tension . This picture is particularly appropriate for states with large angular momenta and length $`r_{q\overline{q}}=ab`$. In this limit linear Regge trajectories with a slope $`\alpha ^{}=2\pi \sigma 1GeV^2`$ arise and in the limit of $`b0`$ a string like model emerges.
In the CFTM baryons consist of three flux tubes which join at a junction $`J`$. The tetra- and penta-quarks have more elaborate color networks with two and three junctions respectively (see Fig.1). The complexity of the penta-quark state suppresses its formation in hadronic collisions due to the specific rearrangements that are required. By detailed balance the system will then also be trapped for a long time while finding the specific path back to the original hadrons. The CFTM offers a concrete realization of this general concept. To allow for the penta-quark state to decay into baryon and meson or for the tetra-quark state to decay into two mesons a junction $`J`$ and the anti-junction $`\overline{J}`$ have to annihilate before the flux tubes can reconnect to generate the final two hadrons.

A necessary but not sufficient condition for the annihilation is that $`J`$ and $`\overline{J}`$ are within a relative distance smaller than the flux tube radius $`b`$. In a simplistic terta-quark model we identify the centers of mass of the di-quarks $`q_1q_2`$ and $`\overline{q^{}}_1\overline{q^{}}_2`$ with the junctions $`J`$ and $`\overline{J}`$. These junctions are treated as scalar particles with the mass of the di-quarks $`m_{dq}650MeV`$. These particles are in the color representation $`\overline{\mathrm{𝟑}}`$ and $`\mathrm{𝟑}`$ with the same linear potential as between the quarks in a meson. The two junctions are in a relative $`p`$-wave. Let $`<r_{J\overline{J}}>`$ in the lowest state be $`a`$. Since the $`l=1`$ wave-function decreases linearly in the small $`r`$ region, the probability $`P(r<b)`$ is approximately $`(b/a)^5`$. For $`\sigma 0.15GeV^2`$ we estimate, using the minimum of the radial potential and the curvature $`d^2V/dr^2`$ at this minimum, that $`r_{min}<r>a`$ is larger than $`0.7fm`$. For $`b`$ in the range of $`0.20.5fm`$ the probability $`P(r<b)`$ is smaller than $`0.0020.2`$.
To estimate the decay rate note that the relative coordinate $`r`$ changes by $`2a`$ during a period $`T`$ exceeding $`a/c0.7fm/c(300MeV)^1`$. During each period $`J\overline{J}`$ annihilation occures with a probability which is smaller than $`P(r<b)`$. Hence the decay rate is $`\mathrm{\Gamma }^{Tetra}<P(r<b)300MeV=0.660MeV`$. In the penta-quark either of the two junctions $`J`$ (see Fig.1) can annihilate with $`\overline{J}`$, therefore we expect roughly twice this width.
The crucial flux tube radius $`b`$ can be determined by pure gluo-dynamics up to small $`1/N_c^2`$ corrections. The large mass of the lightest glue-ball is expected to be above $`m_{gb}1.6GeV`$. This sets the scale for the radius of the flux tube to $`b=1/m_{gb}0.15fm`$. A full QCD lattice simulation, recently performed for $`QQQ`$ baryons with three quarks pinned down at relative distances of $`r_{qq}0.7fm`$, clearly reproduces the $`Y`$ shape of the flux structure shown in an action density contour plot . The corresponding flux junction radius is $`b0.2fm`$. Even in the ground state $`qqq`$ nucleons, where the light quarks move so fast that the short flux tubes get tangled up into a more uniform spherical distribution than in the case of the $`\overline{q}(ud)(ud)`$ penta-quark states, the flux junction radius will remain the same.
On the other hand the distance $`a`$ between $`J`$ and $`\overline{J}`$ in tetra-quarks is related to the distance between the quark and anti-quark in mesons, because both have the same color charges and the same flux tubes. In mesons $`a`$ is essentially the charge radius of the pion $`r_\pi 0.65fm`$. The di-quark masses are $`m_{dq}=650MeV`$ instead of $`350MeV`$ for the light quarks. For a linear potential this scales the size by $`(m_q/m_{dq})^{\frac{1}{3}}`$. However here, the two di-quarks are in a $`p`$-wave. The resulting additional centrifugal barrier effectively doubles the kinetic energy reducing the mass by a factor of about $`2`$, which restores $`<r_{J\overline{J}}>=a`$ back to out initial estimate of about $`0.7fm`$.
Acknowledgments
We like to thank A. Casher for many illuminating discussions of the flux tube model.


format: text

text:
# Heuristic Segmentation of a Nonstationary Time Series

## I Introduction

A stationary time series has statistical properties that do not change under time translation Stratonovich81 . Interestingly, the time series that arise in a large number of phenomena in a broad range of areas—including physiologic systems, economic systems, vehicle traffic systems, and the Internet Ivanov ; Bunde00 ; Goldberger02 ; Stanley02 ; Musha76 ; Leland95 ; Paxson96 ; Crovella97 ; Takayasu ; Abbot99 —are nonstationary. Thus nonstationarity is a property common to both natural and human-influenced phenomena. For this reason, the statistical characterization of the nonstationarities in real-world time series is an important topic in many fields of research and numerous methods of characterizing nonstationary time series have been proposed Method .

One useful approach to quantifying a nonstationary time series is to view it as consisting of a number of time segments that are themselves stationary. The statistical properties of the segments (i) can help us better understand the overall nonstationarity of the time series and (ii) yield practical applications. For example, developing control algorithms for Internet traffic will be easier if we understand the statistical properties of these stationary segments—a behavior that directly corresponds to the “coarser flow rate” of the current network traffic.

In general, it is impossible to obtain the exact segmentation of a nonstationary time series because of the complexity of the calculation. An exact segmentation algorithm requires a computation time that scales as $`O(N^N)`$, where $`N`$ is the number of points in the time series. Hence, such an algorithm is not practical. For this reason, the segmentation of a real-world time series must accomplish a trade-off between the complexity of the calculation and the desired precision of the result.

Bernaola-Galván and co-workers recently proposed a heuristic segmentation algorithm Segment  designed to characterize the stationary durations of heart beat time series. In this algorithm Segment , the calculation cost is reduced by iteratively attempting to segment the time series into only two segments. A decision to cut the times series is made by evaluating a modified Student’s $`t`$-test for the data in the two segments.

The application of this segmentation algorithm reveals that the distribution of the stationary durations in heart beat time series decays as a power law Segment . Intriguingly, a recent analysis of Internet traffic uncovered that the distribution of stationary durations in the fluctuation of the traffic flow density also follows a power law dependence Fukuda01 . Because these signals have their origin in such diverse contexts, these findings suggest that the power law decay of the distribution of the stationary period may be a common occurrence for complex time series. Thus, the correct implementation and interpretation of the results obtained by the segmentation algorithm is essential in understanding the dynamics of the system. In fact, there are many implementation issues concerning the segmentation algorithm of Ref. Segment  that have not yet been addressed explicitly in the literature, especially those concerning the proper estimation of the value of the power-law tail’s exponent in the cumulative distribution of stationary durations.

In this paper we systematically analyze different types of surrogate time series to determine the scope of validity of the segmentation algorithm of Bernaola-Galván and co-workers Segment . In Section II, we briefly explain the segmentation algorithm. In Section III, we present results for the dependence of the exponent of the power law tail on the minimum size of the segments in the distribution of the stationary durations. In Section IV, we consider the effect of the amplitude of the noise and the presence of spike-type noise. In Section V we consider long-ranged temporally-correlated noise. Finally, in Section VI, we summarize our findings.

## II Implementing the Segmentation Algorithm

### II.1 The algorithm

To divide a nonstationary time series into stationary segments Segment , we move a sliding pointer from left to right along the time series and, at each position of the pointer, compute the mean of the subset of the signal to the left of the pointer $`\mu _{\mathrm{left}}`$ and to the right $`\mu _{\mathrm{right}}`$. For two samples of Gaussian distributed random numbers, the statistical significance of the difference between the means of the two samples, $`\mu _1`$ and $`\mu _1`$, is given by the Student’s $`t`$-test statistic prob\_intro

$$t\left|\frac{\mu _1\mu _2}{S_D}\right|,$$
(1)
where

$$S_D=\left(\frac{(N_11)s_1^2+(N_21)s_2^2}{N_1+N_22}\right)^{1/2}\left(\frac{1}{N_1}+\frac{1}{N_2}\right)^{1/2}$$
(2)
is the pooled variance NRC , $`s_1`$, $`s_2`$ are the standard deviations of the two samples, and $`N_1`$ and $`N_2`$ are the number of points in the two samples.

Moving the pointer along our time series from left to right, we calculate $`t`$ as a function of the position in the time series. We use the statistic $`t`$ to quantify the difference between the means of the left-side and right-side time series. Larger $`t`$ means that the values of the mean of both time series are more likely to be significantly different, making point $`t_{\mathrm{max}}`$, with the largest value of $`t`$, a good candidate as a cut point.

We then calculate the statistical significance $`P(t_{\mathrm{max}})`$. Note that $`P(t_{\mathrm{max}})`$ is not the standard Student’s $`t`$-test because we are not comparing independent samples. $`P(t_{\mathrm{max}})`$ is numerically approximated as

$$P(t_{\mathrm{max}})\{1I_{[\nu /(\nu +t_{\mathrm{max}}^2)]}\times (\delta \nu ,\delta )\}^\gamma ,$$
(3)
where $`\gamma =4.19\mathrm{ln}N11.54`$ and $`\delta =0.40`$ are obtained from Monte Carlo simulations Segment , $`N`$ is the size of the time series to be split, $`\nu =N2`$, and $`I_x(a,b)`$ is the incomplete beta function.

If the difference in mean is not statistically significant—i.e., if is smaller than a threshold (typically set to 0.95)—then the time series is not cut. If the difference in means between the left and right part of the time series is statistically significant, then the time series is cut into two segments as long as the means of the two new segments are significantly different from the means of the adjacent segments. If the time series is cut, we continue iterating the above procedure recursively on each segment until the obtained significance value is smaller than the threshold, or the size of the obtained segments is smaller than a minimum size $`\mathrm{}_0`$. We will see that the value of $`\mathrm{}_0`$ is one of the parameters controlling the accuracy of the algorithm.

### II.2 Surrogate Time Series

To investigate the validity of the algorithm, we analyze surrogate time series $`x(t)`$ generated by linking segments with different means. As described in Section I, the cumulative distribution of the stationary durations for some real-world time series is characterized by a power law decay in the tail, so the probability of finding a segment of size larger than $`m`$, i.e., the cumulative distribution of segment sizes in our time series, is

$$P(>m)=\frac{\gamma +1}{m_0^{\gamma +1}}m^\gamma ,\text{for}mm_0,$$
(4)
where $`m_0`$ is the minimum size of a segment.

We generate time series with a power law distribution of segment sizes by the following procedure:

1. draw from the interval \[$`m_0`$,+$`\mathrm{}`$\] a sequence of segment lengths {$`m_i`$} with distribution given by Eq. (4);
2. draw from the interval  a sequence of mean time series values $`\overline{x}_i`$ with uniform probability;
3. draw from the interval \[-$`\sqrt{3}\sigma _ϵ`$,$`\sqrt{3}\sigma _ϵ`$\] a sequence of fluctuation values $`ϵ_i(k_i)`$, for $`k_i=1,\mathrm{},m_i`$, with uniform probability.

The resulting time series is given by

$$x(t)=\overline{x}_i+ϵ_i(k_i),$$
(5)
where $`i`$ is such that

$$\underset{j=0}{\overset{i1}{}}m_j<t\underset{j=0}{\overset{i}{}}m_j,$$
(6)
and $`k_i`$ is such that

$$k_i=t\underset{j=0}{\overset{i1}{}}m_j.$$
(7)
To quantify the level of the noise, we define the ratio

$$R\frac{\sigma _ϵ}{\sigma _{\overline{x}}},$$
(8)
where $`\sigma _{\overline{x}}`$ is the standard deviation of the mean of the segments and $`\sigma _ϵ`$ is the standard deviation of the fluctuations within a segment. For $`\overline{x}`$ uniformly distributed in the interval , $`\sigma _{\overline{x}}=0.3`$.

For each set of parameters ($`R`$, $`\gamma `$, $`m_0`$) we generate ten time series, each with 50,000 data points. Note that knowing a priori the value of $`m_0`$ in a real-world time series is unlikely, but in order to test the algorithm in a consistent way, we consider in the following section $`m_020`$ because that is the resolution limit for the algorithm (see also Appendix A).

Figure 1 displays a surrogate time series, the corresponding stationary durations, and the result of the segmentation algorithm with $`\mathrm{}_0=50`$ and $`m_0=20`$. The segments obtained by the segmentation algorithm do not exactly match the stationary segments in the surrogate time series but the figure strongly suggests that the algorithm provides the correct coarse-grained description of the time series. As expected, the segmentation algorithm cannot extract segments with size $`m<\mathrm{}_0`$.

## III Accuracy of the segmentation algorithm

### III.1 Dependence on $`\mathrm{}_0`$ and $`m_0`$

Figure 2(a) displays the cumulative distribution of segment sizes, which is the probability of finding a segment with size larger than $`\mathrm{}`$ for a surrogate time series split for different values $`\mathrm{}_0`$. The cumulative distributions of the size of the stationary segments cut by the segmentation algorithm for surrogate time series are well fit by a power law decay. For all cases, we find

$$P_{\mathrm{}_0,m_0}(>\mathrm{})\mathrm{}^{\widehat{\gamma }(\mathrm{}_0,m_0)}$$
(9)
with the same exponent value $`\widehat{\gamma }1.0`$ GAMMA , indicating that the segmentation algorithm splits the nonstationary time series into segments with the correct asymptotic statistical properties. However, the range of scales for which we observe a power law decay with $`\widehat{\gamma }\gamma `$ depends strongly on the selection of $`\mathrm{}_0`$.

For $`\mathrm{}`$ greater than about $`5\mathrm{}_0`$, the tails of the distributions are close to power law decays with $`\widehat{\gamma }\gamma =1.0`$. Moreover, all $`P(\mathrm{})`$ follow the original curve for $`\mathrm{}>\text{1,000}`$, i.e., the algorithm correctly identifies the large segments independently of the selection of $`\mathrm{}_0`$. For $`\mathrm{}<5\mathrm{}_0`$, the distributions do not follow power law decays. The origin of this behavior lies in the fact that (i) for $`\mathrm{}_0=20=O(m_0)`$, there are not enough data points to reliably perform the Student’s $`t`$-test, so one cannot reasonably expect any statistical procedure to be able to extract those short segments, and (ii) for $`\mathrm{}_0m_0`$, one fails to extract many of the stationary segments. The reason for the latter is that the value of $`\mathrm{}_0`$ is in this case considerably larger than the size of the shortest segments in the time series, so the algorithm is forced to merge a number of short segments into longer ones with size greater than $`\mathrm{}_0`$. This process gives rise to an excess of segments with sizes between $`\mathrm{}_0`$ and $`2\mathrm{}_0`$.

Table 1 shows the mean and standard deviation of the estimated exponent value $`\widehat{\gamma }`$ calculated from surrogate time series for several values of $`m_0`$ and $`\mathrm{}_0`$ (see Appendix A for details on how to to estimate $`\widehat{\gamma }`$). Our results indicate that $`\widehat{\gamma }`$ depends on both $`m_0`$ and $`\mathrm{}_0`$: If $`\mathrm{}_0m_0`$, $`\widehat{\gamma }`$ overestimates $`\gamma `$, while if $`\mathrm{}_00(m_0)`$, the algorithm correctly estimates the value of the exponent $`\gamma `$; cf. Fig. 2(b).

### III.2 Dependence on $`\gamma `$

Next, we focus on the dependency of the accuracy of the segmentation algorithm on the value of $`\gamma `$. Figure 3(a) displays the cumulative distribution of segment sizes for surrogate time series generated with $`\gamma =2.0`$. A challenge for the segmentation algorithm is that Eq. (4) indicates that the probability of finding segments with size shorter than $`4m_0`$ is 90% for $`\gamma =2.0`$, which can lead to the “aggregation” of several consecutive segments of small size into a single longer segment.

As we found for $`\gamma =1`$, the tails of the distributions follow power law decays for large $`\mathrm{}`$, showing that the algorithm yields segments with the proper statistical properties. In Fig. 3(b), we show the dependence of $`\widehat{\gamma }`$ on $`\mathrm{}_0`$ and $`m_0`$ for $`\gamma =2.0`$. We find that for $`\mathrm{}_0/m_0<4`$ the algorithm extracts segments with distributions of lengths that decay in the tail as power laws with exponents that are quite close to $`2.0`$ .

In Tables 2 and 3, we report the values of $`\widehat{\gamma }`$ for $`\gamma =2.0`$ and $`\gamma =3.0`$, respectively. We find that for small $`m_0`$ and large $`\gamma `$, one over-estimates $`\gamma `$. Thus, we surmise that for $`\gamma >3.0`$, it becomes impractical to estimate $`\gamma `$ accurately, except for extremely long time series. This fact is not as serious a limitation as one may think because for large $`\gamma `$ it is always difficult to judge whether a distribution decays in the tail as an exponential or as a power law with a large exponent.

## IV Robustness of the algorithm with regard to noise

### IV.1 Amplitude of fluctuations around a segment’s mean

Another factor that may affect the performance of the segmentation algorithm of Bernaola-Galván and co-workers is the amplitude of the fluctuations within a segment. It is plausible that greater noise amplitudes will increase the difficulty in identifying the boundaries of neighboring segments. Thus, we next analyze the effect of the amplitude of the noise for surrogate time series.

Figure 4 demonstrates that for large R, the segmentation algorithm yields few short segments. This result arises from the concatenation of neighboring segments with means that become statistically indistinguishable due to the large value of $`\sigma _ϵ`$.

We show in Fig. 5 the cumulative distributions of segment sizes for $`\gamma =1.0`$, $`m_0=20`$, and for different values of $`R`$. For large $`\mathrm{}`$ and $`R3`$, we find $`\widehat{\gamma }\gamma `$, while for $`R=4`$ the algorithm becomes ineffective at extracting the stationary segments in the time series. It is visually apparent that for $`R=4.0`$ the fluctuations within a segment are so much larger than the jumps between the means of the stationary segments that the segmentation becomes totally unable to parse the different segments; cf. Fig 4d.

In Fig. 6, we show $`\widehat{\gamma }`$ for different values of $`R`$. For $`\gamma =1.0`$, we estimate $`\widehat{\gamma }\gamma `$ for $`0<R<4`$. In contrast, for $`\gamma =2.0`$, $`\widehat{\gamma }\gamma `$ only for $`0<R<1.5`$. When $`R`$ increases, the segmentation algorithm is unable to cut the segments because the greater amplitude of the fluctuations inside a segment decreases the significance of the differences between regions of the time series. This effect yields very large segments, which results in very small estimates of $`\widehat{\gamma }`$. This effect is even stronger for $`\gamma =3.0`$, for which we find $`\widehat{\gamma }\gamma `$ only for $`0<R<0.6`$.

### IV.2 Spike noise

Next we analyze the effect of spike noise on the performance of the segmentation algorithm. We generate surrogate time series as before and then for each $`t`$ make, with probability $`p`$, $`x(t)=2`$. The effect of this procedure is illustrated in Fig.7 for four distinct values of $`p`$. The figure also suggests that the segmentation algorithm yields a good coarse-grained description of the surrogate time series for $`p`$ as large as 0.1. This result suggests that the algorithm is robust to the existence of uncorrelated spike noise in the data (Fig. 8).

## V Correlated noise

In this Section, we investigate the effect of long-range correlations in the fluctuations around the segment’s mean on the performance of the segmentation algorithm. This study is particularly important because real-world time series often display long-range power law decaying correlations.

### V.1 Segmentation of correlated noise with no segments

We generate a temporally-correlated noise whose power spectrum decays as $`S(f)f^\beta `$ Makse96 . The surrogate time series consists of 60,000 points, with mean 0. Figure 9 displays the cumulative distribution of segment lengths. The curves show the different noise correlations: $`\beta =0.3`$, 0.5, and 1.0 ($`1/f`$). Note that we have confirmed only one long segment for $`\beta =0.0`$, as one expects. On the other hand, for small $`\beta `$, we still observe a longer segment, although the curves decay rapidly for large $`\beta `$. Thus, the segmentation algorithm divides a correlated noise for large $`\beta `$ into many small pieces of stationary segments, because the given time series is already nonstationary for this case. Also, it is hard to determine the functional form of the plots for all cases. Even if the plots are assumed to follow power laws, the relationship between $`\beta `$ and $`\widehat{\gamma }`$ is quantitatively obscure.

### V.2 Segmentation of correlated noise with segments

Here, in order to show the ability of the identification of segments, we analyze surrogate time series concatenating segments with long-range correlated random variables whose power spectrum is followed by a power law with an exponent $`\beta `$. We show in Figs. 10(a)-(c) typical surrogate time series generated with $`\gamma =1.0`$, $`m_0=20`$, $`R=1.0`$, $`\mathrm{}_0=20`$, and different temporal correlations: (a) $`\beta =0.1`$, (b) $`\beta =0.3`$, and (c) $`\beta =0.5`$. The figure suggests that the segmentation algorithm can correctly parse the short segments but that long segments get cut multiple times, especially for $`\beta =0.5`$. This result is to be expected because the strong correlations in the noise lead to marked changes in the mean.

In Figs. 10(d)-(f), we show the cumulative distributions of segment sizes for $`\beta =0.1,0.3,\text{and }0.5`$. The data confirm quantitatively the visual impression gained from Figs. 10(a)-(c), i.e., that longer segments get cut multiple times. In particular, for $`\beta =0.5`$, the distributions clearly deviate from the power law, independent of the selection of $`\mathrm{}_0`$. However, this fact should not be seen as a shortcoming of the algorithm; for large $`\beta `$, correlated noise for a long segment is already nonstationary, so that the algorithm is indeed cutting a nonstationary signal into stationary durations.

## VI Discussion

In this paper we analyzed nonstationary surrogate time series with different statistical properties in order to investigate the validity of the segmentation algorithm of Bernaola-Galván and co-workers Segment . Our results demonstrate that this heuristic segmentation algorithm can be extremely effective in determining the stationary regions in a time series provided that a few conditions are fulfilled: First, one must have enough data points in the time series to yield a large number of segment lengths, otherwise one will not be able to reach the aymptotic regime of the tail of the distribution of segment sizes.

Second, the ratio of the amplitude of the fluctuations within a segment to the typical jump between the means of the stationary segments must be relatively small (less than about $`0.6`$) in order for one to trust the output of the segmentation algorithm. This concern contrast with the case of spike noise in the data which affects the performance of the segmentation algorithm only weakly.

Finally, if there are long-range temporal correlations of the fluctuations around the mean of the segment, then the segmentation algorithm correctly cuts the time series into the stationary segments for small $`\beta `$. However, for $`\beta >0.3`$, the fluctuations inside long segments become nonstationary, which results in the algorithm detecting many “stationary” durations inside these long segments.

Our analysis provides a number of clear guidelines for using effectively the segmentation algorithm of Bernaola-Galván et al. Segment :

1. One must perform the segmentation for a number of different values of $`\mathrm{}_0`$ in order to identify the region for which the tails of the distributions of segment sizes reach the asymptotic scaling behavior. (Note: If $`\widehat{\gamma }`$ is large, than the estimation error can be quite considerable, especially if $`m_0`$ is small.)
2. One must calculate the ratio between the standard deviation of the mean value of the segment and the standard deviation of the fluctuations within a segment after performing the segmentation. If the $`R`$ than 0.6, then there is the possibility that $`\widehat{\gamma }`$ is considerably under-estimating the true value of $`\gamma `$.

#### Acknowledgments

We thank S. Havlin, P. Ch. Ivanov, and especially P. Bernaola-Galván for discussions. We thank NIH/NCRR (P41 RR13622) and NSF for support. L.A.N.A. acknowledges a Searle Leadership Fund Award for support.

## Appendix A Performance of the algorithm for fixed segment sizes

As a special case, we analyze the time series whose stationary duration is fixed in order to discuss the minimum resolution of the segment algorithm. We concatenate segments with constant lenght $`m`$ with alternating means of 0.0 and 1.0. We then add fluctuations to those segments with a standard deviation of 0.29. We define the fraction of successfully split segments

$$Q\frac{\text{Number of segments correctly identified by the algorithm}}{\text{Number of segments in surrogate data}},$$
(10)
where $`Q=1.0`$ corresponds to perfect segmentmentation.

We plot in Fig. 11 $`Q`$ as a function of $`m`$ for different values of $`\mathrm{}_0`$ . For $`m<20`$, the segmentation algorithm does not yield the correct segments in the surrogate time series even though the segment’s means are quite different. This result suggests that the resolution of the segmentation algorithm is $`20`$. We also find that for $`\mathrm{}_0=5`$, the algorithm splits the time series into too many segments. This result suggests that for optimal performance $`\mathrm{}_010`$.

## Appendix B Estimation of $`\widehat{\gamma }`$

To estimate the exponent $`\widehat{\gamma }`$ of the power laws, we calculate the power law exponent at a segment size $`\mathrm{}_n`$ from a small region around $`\mathrm{}_n`$

$$\widehat{\gamma }(\mathrm{}_n)=\frac{\mathrm{log}(P(\mathrm{}_{n+\tau }))\mathrm{log}(P(\mathrm{}_n))}{\mathrm{log}(\mathrm{}_{n+\tau })\mathrm{log}(\mathrm{}_n)},$$
(11)
where we set the region to $`\tau =5`$. Figure 12 shows the behavior of the power law exponent over $`\mathrm{}`$ corresponding to Fig. 2(a); the black lines are the behavior of the exponent at segment size $`\mathrm{}`$, and the broken lines are the behavior of the exponent of the power law in the distribution of the surrogate time series for comparison.

The range of $`\mathrm{}_n`$ for the calculation is between a starting point (where the exponential decay disappears) and 10,000. For $`\mathrm{}_0=20`$, the value of the exponent decreases for $`\mathrm{}<100`$, which corresponds to the exponential decay as shown in Fig. 2(a). However, the value of the exponent fluctuates around 1.0 for $`\mathrm{}>100`$, indicating that the algorithm can reproduce the correct statistical behavior. For $`\mathrm{}_0=50`$, the value of the exponent is close to 1.0 for $`50<\mathrm{}<\text{10,000}`$. Moreover, for $`\mathrm{}_0=100`$ and $`400`$, the curves decay quickly for the smaller $`\mathrm{}`$, and the value of the exponent tends to be overestimated.


text:
# Production of Long-Lived Ultracold Li2 Molecules from a Fermi gas

## Abstract

We create weakly-bound Li2 molecules from a degenerate two component Fermi gas by sweeping a magnetic field across a Feshbach resonance. The atom-molecule transfer efficiency can reach 85 $`\%`$ and is studied as a function of magnetic field and initial temperature. The bosonic molecules remain trapped for $`0.5`$s and their temperature is within a factor of 2 from the Bose-Einstein condensation temperature. A thermodynamical model reproduces qualitatively the experimental findings.

Feshbach resonances constitute a unique tool to tune the microscopic interactions in ultracold bosonic and fermionic gases Feshbach ; Naturespecial ; Inouye98 . These resonances arise when the total energy of a pair of colliding atoms matches the energy of a bound state of another hyperfine manifold, leading to the resonant occupation of this state during the collision. Thus, by means of an external magnetic field, one is able to change the magnitude and sign of the scattering length. In bosonic samples, the collapse of Bose-Einstein condensates (BEC) for negative scattering length Jilacollapse , soliton formation Khaykovich02 ; Truscott02 , and coherent oscillations between an atomic condensate and molecules Donley02  have been observed. For fermions with attractive interaction, the superfluid transition temperature is predicted to be maximum near a Feshbach resonance Nozieres85 ; Randeria93 ; Holland01 ; Ohashi02 ; Milstein02 .

In this letter we give another striking example of the control of interactions in a Fermi gas. We perform time-dependent experiments near a Feshbach resonance to produce in a reversible manner ultracold and trapped molecules from a quantum degenerate fermionic 6Li gas. The production efficiency exceeds 80% and the observed molecule lifetime reaches half a second. The phase-space density of these bosonic molecules is on the order of one, the highest value reported so far. Using a similar method with fermionic 40K atoms, the JILA group recently reported molecule production with a lifetime of 1 ms, atom to molecule conversion efficiency of 50% and direct measurement of the molecular binding energy Regal03 . Molecules have also recently been produced from 87Rb and 133Cs condensates Grimm03 ; Durr03  and also from a cold 133Cs cloud Chin03 . Molecule formation has also been achieved through one-photon or two-photon photoassociation but with a considerably lower phase-space density than reported here Wynar00 . Our work paves the way to Bose-Einstein condensation of molecules and to the study of the crossover between the regime of molecular BEC and the regime of superfluid BCS pairing in Fermi systems Nozieres85 ; Randeria93 ; Ohashi02 ; Milstein02 .

In 6Li a broad ($`100`$ Gauss) Feshbach resonance exists between the two Zeeman sublevels of the hyperfine ground state : $`|1/2,1/2`$ and $`|1/2,1/2`$ at a magnetic field of 810 Gauss, see Fig. 1 Gupta03 ; Bourdel03 . In a recent study of this Feshbach resonance, we have reported an anomalous negative value for the gas interaction energy between 700 Gauss and 810 gauss, i.e. below resonance, where the scattering length $`a`$ is positive Bourdel03 . We suggested that this could be due to the presence of weakly bound molecules confined simultaneously with the cloud of ultracold fermions, and a recent theoretical paper explains our results Kokkelmans03 .

The method used here to produce molecules is illustrated in Fig.1 and was suggested for Bose gases in Timmermans98 ; Abeelen99 ; Julienne . It consists in scanning over a Feshbach resonance from the region of attractive interaction ($`a<0`$) to region 2 in Fig. 1 where $`a`$ is large and positive, and where a weakly bound molecular state exists with energy $`E_b=\mathrm{}^2/ma^2`$ (where $`m`$ is the atomic mass). It is thus energetically favorable to populate this bound state. Having twice the polarisability of individual atoms, the molecules are expected to be confined in an optical dipole trap with the same oscillations frequencies as atoms and twice the trap depth.

Experimentally, using the setup described in Bourdel03 , we prepare $`N_1=N_{}+N_{}=\mathrm{1.5\hspace{0.17em}10}^5`$ atoms in a $`50`$($`10`$)$`\%`$ mixture of $`|1/2,1/2`$ and $`|1/2,1/2`$ in a crossed Nd:YAG dipole trap at $`B=1060G`$ where the scattering length between $`|1/2,1/2`$ and $`|1/2,1/2`$ is large and negative $`a=150`$nm (position 1 in Fig. 1). For a power of $`1.53`$W (resp. $`4.3`$W) in the horizontal (resp. vertical) trapping beam propagating along Ox (resp. Oy), oscillation frequencies are $`\omega _x/2\pi =2.9(3)`$kHz, $`\omega _y/2\pi =5.9(5)`$kHz and $`\omega _z/2\pi =6.5(6)`$kHz. By evaporation in the optical trap, the temperature $`T`$ of the gas mixture can be tuned between $`0.2`$ and $`0.5T_F`$ where $`T_F`$ is the Fermi temperature defined by $`T_F=\mathrm{}\overline{\omega }(3N_1)^{1/3}/k_B`$ and $`\overline{\omega }=(\omega _x\omega _y\omega _z)^{1/3}`$.

The Fermi gas quantum degeneracy is measured through analysis of absorption images after abrupt ($`20\mu `$s) switch-off of the magnetic field and time of flight expansion of the cloud. Free 6Li atoms are detected at zero magnetic field, for which $`a0`$, using laser light tuned to the D2 line, position 4 in Fig. 1. Having prepared the degenerate Fermi gas at position 1, we sweep the magnetic field in 50 ms to position 2 where $`a`$ is large and positive, the region where a weakly-bound molecular state exists. Absorption images indeed reveal that the number of atoms $`N_2`$ as counted after $`B`$ switch-off to position 4 as explained above, has dramatically decreased to $``$20$`\%`$ of the initial number $`N_1`$. After letting the system at position 2 for a variable wait time $`t_\mathrm{w}`$ between $`0`$ and 5 s, we sweep the magnetic field back to position 1 in 50 ms and count the number of atoms, $`N_3`$, after this round-trip through resonance.

Surprisingly, when position 2 is at $`689`$ Gauss and $`t_\mathrm{w}=0`$, we find $`N_3N_1`$, indicating that no loss has occurred in 100 ms. This proves that the decrease in detected atom number at position 2 is not due to losses. We are thus led to conclude that: (i) at position 2, 80$`\%`$ of the atoms were in a molecular bound state which is not detected using light resonant with the atoms at position 4, (ii) the atom-to-molecule formation process is reversible. $`(N_3N_2)/2`$ represents thus the number of molecules at position 2.

In comparison to other recent experiments on bosons and fermions Donley02 ; Durr03 ; Regal03 ; Grimm03 , our 50 ms scan over the Feshbach resonance is very slow and, in particular, much longer than the collision time in the cloud. This allows the system to relax towards thermal equilibrium between atoms and molecules during the magnetic field ramp.

The relative fraction of atoms bound in a molecule, $`(N_3N_2)/N_3`$, and the temperature of the atoms at position 2 and 3 are plotted in Fig. 2 as a function of the magnetic field at position 2. The initial parameters at $`1060`$G for this experiment are $`T=4.7\mu `$K, $`T_F=11\mu `$K, and $`N_1=\mathrm{8\hspace{0.17em}10}^4`$. The frequencies of the trap are $`\omega _x/2\pi =2.2(2)`$kHz, $`\omega _y/2\pi =4.6(4)`$kHz, $`\omega _z/2\pi =5.1(5)`$kHz. One observes that very few molecules are detected above $`0.77`$kG. When decreasing the magnetic field, the fraction of molecules increases up to $`60\%`$. Below $`650`$G, losses become important and the molecular fraction decreases, an effect that we will study below.

A key parameter in our detection scheme is the switch off time ($`20\mu `$s) of the Feshbach magnetic field from position 2. Indeed, the molecules may or may not dissociate in this process. Molecules will dissociate if the relative rate of change of the binding energy $`\text{d}E_b/E_b\text{d}t`$ is greater than the oscillation pulsation $`E_b/\mathrm{}`$. So, the dissociation of molecules is governed by the parameter $`\alpha =\mathrm{}\text{d}E_b/E_{b}^{}{}_{}{}^{2}\text{d}t`$. If $`\alpha 1`$, they dissociate and appear as free atoms in the absorption image. Our detection scheme is unable to prove molecule formation. On the contrary, if the binding energy is sufficiently large at position 2 ($`\alpha 1`$), the molecular state follows adiabatically the change of $`B`$ to a deeper bound state which is off-resonant with the detection laser. In this case, $`N_3>N_2`$ is the signature of molecule formation.

For $`B=700`$G, one has $`\alpha 1`$. Reducing the rate $`\text{d}B/\text{d}t`$ by a factor of $`500`$, we observe a shift of the threshold for molecule detection to $`800`$G, as expected from the expression of $`\alpha `$ and the position of the Feshbach resonance at $`810`$G Bourdel03 . As a second consequence of the role of the switch-off time, we attribute the decrease of the molecular fraction above $`720`$G in Fig. 2a to the crossover between the adiabatic and non adiabatic regimes.

The temperature of the atoms at position 2 and 3 is plotted in Fig. 2b. A decrease of the magnetic field increases the binding energy of the molecules. Therefore, the molecular fraction increases and so does the temperature. However, sweeping back the magnetic field to the initial position, the temperature returns close to its initial value, demonstrating the reversibility of the molecular formation process. According to Petrov03 , in our experiment the rate of 3-body recombination to the weakly bound state is greater than $`10^4`$Hz (for $`B>665`$G). Moreover, the atom-atom, atom-molecule and molecule-molecule Petrovprivate  collision rates are also large compared to the rate at which the resonance is scanned. Therefore the system remains close to thermal equilibrium. During the molecule formation process, the released binding energy must be converted into kinetic energy, therefore the increase of temperature in Fig. 2b, as $`a`$ decreases, is expected.

Interestingly, the efficiency of molecule formation strongly depends on the trap depth and reaches values up to $`80\%`$ in Fig. 3, the highest efficiency reported so far. The trap depth is slowly (800 ms) reduced at $`B=1060`$G, before the magnetic field sweep through resonance. The gas temperature and the Fermi temperature decrease by evaporation and adiabatic cooling due to the reduction of the trap oscillation frequencies. At $`689`$G, $`a=78`$nm, $`E_b/k_\mathrm{B}=12\mu `$K. When reducing the trap depth, $`T`$ and $`T_F`$ become small compared to $`E_b`$ and the molecular fraction increases towards 1. Other experiments show that both evaporation (reduction of the ratio $`T/T_\mathrm{F}`$), and adiabatic cooling (reduction of $`T`$ and $`T_\mathrm{F}`$ by the same factor) effectively lead to an increase of the molecular fraction. For instance, at constant $`T_F=11\mu `$K, as the temperature is decreased from $`0.5T_\mathrm{F}`$ to $`0.3T_\mathrm{F}`$, we find that the molecular fraction at $`689`$G increases from $`40\%`$ to $`65\%`$. Similarly, at constant $`T/T_\mathrm{F}=0.3`$, as $`T`$ varies from $`4.5\mu `$K to $`1.8\mu `$K, it increases from $`50\%`$ to $`85\%`$.

Our results can be understood within a thermodynamical model assuming that atoms and molecules are in equilibrium in a trap during the magnetic field sweep Kokkelmans03 . We assume that the sweep is adiabatic, so the entropy is conserved and it is calculated for the initial situation where no molecules are present at $`B=1060`$G. The atom-atom, atom-molecule, molecule-molecule interactions were either neglected (solid line in Fig. 2), or included via a mean field contribution proportional to $`a`$ Petrov03 ; Petrovprivate , in the regime where it is allowed (dotted line). The molecular fraction is then determined by a set of two parameters $`T/T_\mathrm{F}`$ and $`E_b/T_\mathrm{F}`$. As the binding energy gets more negative, the fraction of molecules increases since the molecular state becomes energetically more favorable. Moreover, the temperature increases as each molecule formation event transfers the binding energy into kinetic energy. Both models agree fairly well with our data in the region where both the detection of molecules is efficient and losses are unimportant ($`660B720`$G). In the interacting case, a perfect agreement is achieved if we allow for a 10$`\%`$ mismatch in the quality of the spin mixture. In the strongly interacting region where $`na^31`$, the non-interacting model is shown with a dashed line because we do not expect our model of separated gases of fermionic atoms and bosonic molecules to be valid. Those models also account qualitatively for the increase in the molecular fraction with the evaporation in Fig. 3. Indeed, as the trap depth is lowered, $`T`$ decreases, the ratio $`E_b/T`$ increases and the bound state will then have a higher occupancy.

In order to investigate prospects for the formation of a molecular BEC, two questions immediately arise: what are the lifetime and degeneracy parameter of the produced molecular cloud ? The decay of the molecules in the optical trap is found by measuring the quantity $`(N_3N_2)/2`$ as a function of the wait time $`t_\mathrm{w}`$ at position 2 before returning to 3 for the detection process. Fig. 4 presents the number of remaining molecules as a function of the time elapsed at point 2 in a trap with frequencies $`\omega _x/2\pi =0.95(10)`$kHz, $`\omega _y/2\pi =2.0(2)`$kHz, $`\omega _z/2\pi =2.2(2)`$kHz. In Fig. 4 a and b, the initial temperatures at position 1 are the same, $`1.1\mu `$K $`=0.3T_\mathrm{F}`$. Fig 4a corresponds to $`B=689G`$, $`a=78`$nm, $`E_b/k_B=12\mu `$K, whereas Fig. 4b corresponds to $`B=636G`$, $`a=35`$nm, and $`E_b/k_B=60\mu `$K. The initial fraction of unpaired atoms (not engaged in a molecule) is $`17\%`$ in (a) and $`25\%`$ in (b). For a change in $`a`$ by the modest factor of $`2.2`$, the two samples exhibit strikingly different lifetimes. Decays are well fitted by 2-body + 1-body loss processes with initial lifetimes of $`500`$ms in (a) ($`689`$G) and $`20`$ms in (b) whereas the 1-body decay exceeds 4 seconds.

A strong decrease in the lifetime of molecules with decreasing $`a`$ is expected assuming that they undergo collisional relaxation to deep bound states Petrovprivate . However, the observed factor of 25 in lifetimes is surprisingly large. At $`689`$G, we can evaluate the 2-body loss rate to be $`G=2.4_{1.6}^{+3.2}\mathrm{\hspace{0.17em}10}^{13}`$cm3.s-1. At $`636`$G, we have initially $`N_30.3N_1`$. These fast losses indicate that very likely some evaporation is involved in the limited lifetime and that the system has not reached thermal equilibrium. Therefore, the 2-body loss rate cannot be safely estimated.

Let us now evaluate the phase space density of the trapped molecules. In Fig. 2, at $`689`$G, there are $`\mathrm{1.8\hspace{0.17em}10}^4`$ molecules, confined with $`\mathrm{3.3\hspace{0.17em}10}^4`$ atoms. The measured atom temperature is $`T_{at}=6.7\mu `$K. Assuming thermal equilibrium between atoms and molecules, $`T_{at}=T_m`$ for magnetic fields above $`650`$G, we obtain the peak density of molecules $`n_m\mathrm{4\hspace{0.17em}10}^{13}`$cm-3. Then the critical temperature for molecule condensation is reduced due to the interactions and is $`3.5\mu `$K giving $`T/T_C2`$. Since $`n_ma_m^3\mathrm{5\hspace{0.17em}10}^31`$ with $`a_m=0.6a=47`$nm (Petrovprivate ), the molecular gas is in the dilute regime and the mean distance between molecules is larger than the typical size of a molecule $`a`$. In fact, for all data between $`675`$G and $`750`$G, the phase space density of the molecules is not far from the condensation point. Since, the lifetime of the molecules, $`500`$ms, is long in this region when compared to the molecule-molecule elastic collision time $`1/(n_m\mathrm{\hspace{0.17em}8}\pi a_m^2v_m)3\mu `$s, (for $`a_m=47`$nm), it should be possible to evaporate the molecules further to reach the Bose-Einstein condensation threshold.

In summary, we have produced long-lived and trapped Li2 dimers. The atom-molecule conversion efficiency can approach $`1`$ when the Fermi quantum degeneracy is strong. The lifetime of the trapped molecules strongly depends on the scattering length. Current research concentrates on evaporation towards molecular BEC. Prospects for producing superfluid Fermi mixtures and for investigating the transition between molecular condensates and superfluid Fermi gases are promising Petrovprivate ; Ohashi02 ; Milstein02 ; Carr03 .

We are grateful to D. S. Petrov, S. Stringari, Y .Castin, R. Combescot and J. Dalibard for useful discussions. This work was supported by CNRS, Collège de France, and Région Ile de France. S. Kokkelmans acknowledges a Marie Curie grant from the E.U. under contract number MCFI-2002-00968. Laboratoire Kastler Brossel is Unité de recherche de l’Ecole Normale Supérieure et de l’Université Pierre et Marie Curie, associée au CNRS.


text:
# Sonic-Point and Spin-Resonance Model of the Kilohertz QPO Pairs

## 1 Introduction

Recent discoveries made with RXTE have clarified the relationship between the pairs of kilohertz quasi-periodic oscillations (kilohertz QPOs) observed in the accretion-powered X-ray emission of neutron stars in low-mass binary systems (see van der Klis 2000; Lamb 2002), the high-frequency oscillations seen during thermonuclear X-ray bursts (see Strohmayer 2001; Strohmayer & Bildsten 2003), and the spin frequency of the star. Key developments include the discovery of burst oscillations and kilohertz QPOs in the 401 Hz accretion-powered X-ray pulsar SAX J1808.4$``$3658 (Chakrabarty et al. 2003; Wijnands et al. 2003) and burst oscillations in the 314 Hz accretion-powered X-ray pulsar XTE J1814$``$338 (Markwardt & Swank 2003; Markwardt, Strohmayer, & Swank 2003b).

SAX J1808.4$``$3658 was discovered and identified as an X-ray burst source during its September 1996 outburst (in ’t Zand et al. 1998). During its April 1998 outburst, a 401 Hz periodic oscillation (“pulsation”) was detected in its persistent emission (Wijnands & van der Klis 1998; Chakrabarty & Morgan 1998), but neither kilohertz QPOs nor burst oscillations were detected then. Subsequent analysis of an X-ray burst recorded by BeppoSAX during the September 1996 outburst yielded a marginal detection of an oscillation at $`400\pm 2`$ Hz (in ’t Zand et al. 2001). SAX J1808.4$``$3658 was extensively observed throughout its October 2002 outburst and periodic oscillations, burst oscillations, and kilohertz QPOs were all detected with high confidence (Chakrabarty et al. 2003; Wijnands et al. 2003), making SAX J1808.4$``$3658 a Rosetta stone for understanding these phenomena.

The results of Chakrabarty et al. (2003) establish conclusively that the spin frequency of the neutron star in the SAX J1808.4$``$3658 system is 401 Hz: no oscillation at half this frequency was detected, and the upper limit on the rms amplitude of any such oscillation was determined to be 0.014%, 300 times smaller than the amplitude of the 401 Hz oscillation. This excludes the possibility that the spin frequency is 200.5 Hz.

Wijnands et al. (2003) discovered a pair of kilohertz QPOs in SAX J1808.4$``$3658. Their frequency separation $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ is within 3% of $`\nu _{\mathrm{spin}}/2`$, demonstrating that $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ is related to the spin of the star and that the system somehow generates a frequency difference approximately equal to half the spin frequency, a result unexplained by existing models.

Chakrabarty et al. (2003) observed four X-ray bursts during the October 2002 outburst ofSAX J1808.4$``$3658. All showed a $``$401 Hz oscillation, demonstrating that in this star, $`\nu _{\mathrm{burst}}\nu _{\mathrm{spin}}`$. Markwardt et al. (2003b) have analyzed three X-ray bursts from XTE J1814$``$338. All showed nearly coherent $``$314 Hz oscillations, demonstrating that $`\nu _{\mathrm{burst}}\nu _{\mathrm{spin}}`$ in this star, also. These results make compelling the previously strong evidence (Strohmayer et al. 1996, 1997, 1998a, 1998b; Strohmayer & Markwardt 1999; Muno et al. 2000; Muno et al. 2002; Strohmayer & Markwardt 2002) that burst oscillations are generated by the spin of the star and suggest that $`\nu _{\mathrm{burst}}`$ is close to $`\nu _{\mathrm{spin}}`$ in all the neutron stars in which burst oscillations have been observed. If so, then $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ is $`\nu _{\mathrm{spin}}`$ in the stars with $`\nu _{\mathrm{spin}}360`$ Hz but is $`\nu _{\mathrm{spin}}/2`$ in the stars with $`\nu _{\mathrm{spin}}400`$ Hz. A successful model must therefore explain why some stars have $`\mathrm{\Delta }\nu _{\mathrm{QPO}}\nu _{\mathrm{spin}}`$, whereas others have $`\mathrm{\Delta }\nu _{\mathrm{QPO}}\nu _{\mathrm{spin}}/2`$.

Here we propose a modification of the original sonic-point beat-frequency model (Miller, Lamb, & Psaltis, hereafter MLP98; Lamb & Miller 2001) that can explain within a single framework why $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ is close to $`\nu _{\mathrm{spin}}`$ for some stars but close to $`\nu _{\mathrm{spin}}/2`$ for others. In § 2 we discuss the general implications of the new observations for kilohertz QPO mechanisms and in § 3 we propose specific mechanisms for generating the upper and lower kilohertz QPOs. We discuss our conclusions in § 4.

## 2 General Inferences from the Observations

Taken together, the new and previous observations provide strong hints about the mechanisms that generate the kilohertz QPO pairs:

1. It appears highly likely that the frequency of one of the two kilohertz QPOs reflects the orbital frequency of gas in the inner disk. The kilohertz QPOs have frequencies similar to those of orbital motion near neutron stars. The frequencies of the kilohertz QPOs vary by hundreds of Hertz on time scales as short as minutes (see, e.g., Méndez et al. 1999; van der Klis 2000). This is possible if their frequencies are related to orbital motion at a radius that varies (Lamb 2002). These frequency variations exclude mechanisms like the one proposed by Abramowicz & Kluzniak (2001) to explain the high-frequency QPOs observed in black hole candidates and by Kluzniak et al. (2003) to explain the kilohertz QPOs observed in SAX J1808.4$``$3658, because this type of mechanism requires a low-order resonance between the geodesic frequencies of test particles orbiting at a fixed radius and the resonance disappears when the two frequencies change.

2. The star’s spin is somehow involved in producing the frequency separation of the two kilohertz QPOs in a pair. This involvement is clear in SAX J1808.4$``$3658, where $`\mathrm{\Delta }\nu _{\mathrm{QPO}}\nu _{\mathrm{spin}}/2`$. It is strongly indicated in the other kilohertz QPO sources because in all cases where both $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ and $`\nu _{\mathrm{burst}}`$ have been measured, $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ is approximately commensurable with the stellar spin frequency inferred from $`\nu _{\mathrm{burst}}`$ (see van der Klis 2000; Lamb 2002). In four cases, $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ is consistent with $`\nu _{\mathrm{spin}}`$ or with half $`\nu _{\mathrm{spin}}`$ to better than 5% (4U 1608$``$52: $`\nu _{\mathrm{spin}}=319`$ Hz, $`\mathrm{\Delta }\nu _{\mathrm{QPO}}=301.3\pm 7.9`$ Hz \[Méndez et al. 1997; Méndez 2000, private communication\]; 4U 1702$``$429: $`\nu _{\mathrm{spin}}=329`$ Hz, $`\mathrm{\Delta }\nu _{\mathrm{QPO}}=333\pm 5`$ Hz \[Markwardt, Strohmayer, & Swank 1999\]; KS 1731$``$260: $`\nu _{\mathrm{spin}}=524`$ Hz, $`\mathrm{\Delta }\nu _{\mathrm{QPO}}=260\pm 10`$ Hz \[Smith, Morgan, & Bradt 1997; Wijnands & van der Klis 1997\]; and SAX J1808.4$``$3658: $`\nu _{\mathrm{spin}}=401`$ Hz, $`\mathrm{\Delta }\nu _{\mathrm{QPO}}=195\pm 6`$ Hz \[Chakrabarty et al. 2003; Wijnands et al. 2003\]). In all other known cases, the largest value of $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ is approximately consistent with either $`\nu _{\mathrm{spin}}`$ or $`\nu _{\mathrm{spin}}/2`$.

3. A mechanism that produces a single sideband is indicated. Most mechanisms that modulate the X-ray brightness at two frequencies (such as amplitude modulation) would generate at least two strong sidebands. Although weak single and double sidebands have been detected close to the frequency of the lower kilohertz QPO (Jonker, Méndez, & van der Klis 2000), at most two strong kilohertz QPOs are observed in a given system (van der Klis 2000; Méndez & van der Klis 2000). This suggests that the frequency of one QPO is the primary frequency while the other is generated by a single-sideband mechanism. Beat-frequency mechanisms naturally produce a single sideband. Because one QPO frequency is almost certainly an orbital frequency, the most natural mechanism would be one in which the second frequency is generated by a beat with the star’s spin frequency or with another orbital frequency.

4. A successful model should explain why stars with high spin frequencies (as inferred from their burst oscillations) tend to have $`\mathrm{\Delta }\nu _{\mathrm{QPO}}\nu _{\mathrm{spin}}/2`$, whereas stars with low spin frequencies tend to have $`\mathrm{\Delta }\nu _{\mathrm{QPO}}\nu _{\mathrm{spin}}`$.

All existing beat-frequency models predict that if the fundamental frequencies are (1) an orbital frequency $`\nu _{\mathrm{orb}}`$ and (2) the frequency produced by a pattern with $`n`$-fold symmetry rotating with the star, then the principal beat frequency will be $`n(\nu _{\mathrm{orb}}\nu _{\mathrm{spin}})`$ (Lamb et al. 1985). The sonic-point beat-frequency model does predict the appearance of weaker QPOs at a variety of other frequencies (MLP98). One of these could by chance differ from the fundamental orbital frequency by about $`\nu _{\mathrm{spin}}/2`$, but this would be an ad hoc explanation and would not explain a strong QPO at this frequency. Consequently, when evidence arose that in many cases $`\nu _{\mathrm{QPO1}}\nu _{\mathrm{QPO2}}\nu _{\mathrm{burst}}/2`$, it was suggested that perhaps in these stars $`\nu _{\mathrm{burst}}2\nu _{\mathrm{spin}}`$ (MLP98; see also Strohmayer et al. 1998b, Muno et al. 2001). The compelling evidence that $`\nu _{\mathrm{QPO1}}\nu _{\mathrm{QPO2}}\nu _{\mathrm{spin}}/2`$ in SAX J1808.4$``$3658 requires a model in which this frequency arises naturally. At the same time, beat-frequency models have many advantages, as discussed above. We propose a beat-frequency model in which $`\nu _{\mathrm{QPO1}}`$ is either approximately $`\nu _{\mathrm{QPO2}}\nu _{\mathrm{spin}}`$ or approximately $`\nu _{\mathrm{QPO2}}\nu _{\mathrm{spin}}/2`$.

## 3 The Sonic-Point and Spin-Resonance Model

We first show that the star’s magnetic and radiation fields will excite vertical motion in the disk at the “spin-resonance” radius where $`\nu _{\mathrm{orb}}\nu _{\mathrm{spin}}`$ equals the vertical epicyclic frequency $`\nu _\psi `$, and that the orbital frequency at this radius is approximately $`\nu _{\mathrm{spin}}/2`$. We propose that the upper kilohertz QPO is generated by clumps orbiting at the sonic point via the same mechanism as in the original sonic-point beat-frequency model (MLP98). Indeed, there is new evidence supporting this identification. We then show that interaction of the radiation pattern generated by clumps orbiting at the sonic-point radius $`r_{sp}`$ with the vertical motion of gas at the spin-resonance radius can produce a second QPO at either $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}/2`$ or $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}`$. We propose that this second QPO is the lower kilohertz QPO.

### 3.1 Generation of motions with frequency $`\nu _{\mathrm{spin}}/2`$

The magnetic and radiation fields of a neutron star are not perfectly symmetric about its rotation axis. Consequently, as the star spins and the gas in the disk orbits, each element of gas experiences periodically varying radial and vertical forces. If the frequency of this forcing is close to a natural frequency, the response of the element of gas may be large. The two most basic natural frequencies of an element of gas—aside from its orbital frequency—are its radial and vertical epicyclic frequencies. Vertical epicyclic motion appears more likely to generate a QPO, for two reasons. First, oscillation of an element of gas in the vertical direction can generate an X-ray signal by, for example, scattering radiation produced elsewhere in the system into or out of the observer’s line of sight. In contrast, oscillation of an element of gas in the radial direction does not offer such possibilities. Second, vertical motion of gas in an accretion disk may be weakly damped (Marković & Lamb 1998) whereas radial motion is strongly damped (Marković & Lamb 2000).

Perturbation of the disk flow by a magnetic or radiation field that rotates with the star will preferentially excite vertical motion in the disk at the spin-resonance radius $`r_{sr}`$ defined implicitly by the resonance condition

$$\nu _{\mathrm{spin}}\nu _{\mathrm{orb}}(r_{sr})=\nu _\psi (r_{sr}),$$
(1)
where $`\nu _\psi (r)`$ is the vertical epicyclic frequency at radius $`r`$. Here $`\nu _{\mathrm{spin}}\nu _{\mathrm{orb}}(r_{sr})`$ is the forcing frequency experienced by an element of gas orbiting at $`r_{sr}`$. Numerical simulations (see Fig. 1) show that the vertical displacement of gas is much greater at the resonant radius than at any other radius. In a Newtonian $`1/r`$ gravitational potential, $`\nu _\psi (r)=\nu _{\mathrm{orb}}(r)`$. In general relativity, $`\nu _\psi (r)`$ is not exactly equal to $`\nu _{\mathrm{orb}}(r)`$, but the difference is $`<2`$ Hz at the radii of interest (where $`\nu _{\mathrm{orb}}<300`$ Hz). Consequently, at the resonance radius where vertical motion is preferentially excited, $`\nu _{\mathrm{orb}}\nu _\psi \nu _{\mathrm{spin}}/2`$, i.e., at this radius the orbital and vertical frequencies are both approximately half the star’s spin frequency.

Each individual element of gas at the spin-resonance radius orbits the star with frequency $`\nu _{\mathrm{orb}}(r_{\mathrm{sr}})\nu _{\mathrm{spin}}/2`$. However, the pattern created by the elements of gas that are above the plane moves around the star with frequency $`\nu _{\mathrm{pattern}}=\nu _{\mathrm{spin}}`$. Stated differently, the azimuth where the gas is highest above the disk advances at the rate $`\nu _{\mathrm{spin}}`$, even though no gas moves at this rate. The reason for this is that the response of the gas is almost in phase with the periodic driving force. For example, if the force on a particular element of gas is maximal at the time-dependent azimuthal position $`\varphi _{\mathrm{max}}(t)=\varphi _0+\nu _{\mathrm{spin}}t`$ and when the gas is furthest above the disk plane, then as the forcing field rotates, the gas in the disk at $`r_{sr}`$ will be at the peak of its vertical excursion at the azimuth $`\varphi =\varphi _{\mathrm{max}}(t)`$.

So far we have considered the effect of forces acting only on one (the ‘upper”) surface of the disk. The second harmonic of $`\nu _{\mathrm{spin}}`$ is weak in the persistent X-ray emission of all known millisecond accretion-powered pulsars and, where both have been determined, $`\nu _{\mathrm{burst}}`$ is close to $`\nu _{\mathrm{spin}}`$ rather than $`2\nu _{\mathrm{spin}}`$, indicating that we observe only one (the “upper”) magnetic pole of the star and its X-ray beam as the star rotates. The star’s opposite pole and its X-ray beam interact with clumps when they are below the lower surface of the disk with a phase that will further amplify their vertical motions.

### 3.2 Generation of the upper kilohertz QPO

We propose that the upper kilohertz QPO is related to the orbital motion of clumps at a special radius near the neutron star. As shown in MLP98, the radius $`r_{sp}`$ of the sonic point in the disk is a strong candidate for this special radius. The sonic transition at $`r_{sp}`$ is a strong-field general relativistic effect and is usually produced by the force on the gas in the disk exerted by radiation from the star. At $`r_{sp}`$, the inward radial velocity in the accretion disk increases sharply. The sharpness of this change in the radial velocity matches the observed narrow width of the upper kilohertz QPO in power density spectra. In addition, radiation and magnetic forces are likely to enhance clumping at the sonic radius over clumping further out in the disk. As in the original sonic-point model, we assume the frequency $`\nu _{\mathrm{QPO2}}`$ of the upper kilohertz QPO is generated by mapping of accretion streams from clumps orbiting at $`r_{sp}`$ onto the surface of the star. This mechanism for generating the upper kilohertz QPO is supported by the observed short-term anticorrelation of the frequency of the upper kilohertz QPO and the X-ray flux in 4U 1608$``$52 (Yu, van der Klis, & Jonker 2001; Yu & van der Klis 2002).

As discussed in Lamb & Miller (2001), the inward drift of the clumps at $`r_{sp}`$ causes the footpoints of the streams on the stellar surface to revolve around the star with an angular velocity slightly different from $`\nu _{\mathrm{orb}}(r_{sp})`$. The radiation pattern they produce rotates with frequency $`\nu _{\mathrm{QPO2}}\nu _{\mathrm{orb}}(r_{sp})\delta \nu _{\mathrm{rad}}`$, where $`\delta \nu _{\mathrm{rad}}`$ is typically positive and increases in magnitude from nearly zero to a few tens of Hertz as $`r_{sp}`$ decreases.

### 3.3 Generation of the lower kilohertz QPO

Scattering of radiation from the neutron star by the clumps orbiting in the disk close to the sonic radius will create a radiation pattern on the surface of the disk that rotates with frequency $`\nu _{\mathrm{orb}}(r_{\mathrm{sp}})`$. This pattern will preferentially illuminate or shadow the gas orbiting at the spin-resonance radius, when it is above the disk plane. The rotating radiation pattern will be scattered much more strongly by clumps orbiting at the spin-resonance radius than by clumps elsewhere because clumps at this radius are making large vertical excursions, unlike clumps elsewhere in the disk.

Generation of a lower kilohertz QPO with frequency $`\nu _{\mathrm{QPO1}}\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}/2`$.—Suppose first that the gas in the disk near the spin-resonance radius is highly clumped. When illuminated, each clump orbiting at $`r_{sr}`$ scatters radiation in all directions. In effect, each clump redirects the radiation propagating outward from the sonic radius in the modest solid angle that it subtends (as seen from the sonic radius) into all directions. From the point of view of a distant observer, each individual clump looks like a light bulb that is blinking on and off with a frequency equal to $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{orb}}(r_{sr})\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}/2`$. If there are only a modest number of clumps at $`r_{sr}`$, the effect is somewhat like what one would see if a movie marquee had only a few light bulbs, so the blinking of the individual bulbs dominates the time variation. Because the radiation is scattered in all directions, an observer does not have to be close to the disk plane to see the X-ray flux modulation. The left-hand panel of Figure 2 shows the power spectrum of the flux variation generated in a simulation in which five randomly-positioned clumps scatter the radiation pattern coming from the sonic radius. The peak at $`\nu _{\mathrm{orb}}(r_{\mathrm{sp}})\nu _{\mathrm{spin}}/2`$ is clearly dominant.

Generation of a lower kilohertz QPO with frequency $`\nu _{\mathrm{QPO1}}\nu _{\mathrm{orb}}(r_{\mathrm{sp}})\nu _{\mathrm{spin}}`$.—Suppose now that the gas in the disk near the spin-resonance radius is less highly clumped. There may be a larger number of smaller clumps or the flow may even be relatively smooth. As before, each element of gas is oscillating vertically with frequency $`\nu _{\mathrm{spin}}/2`$. Together they form a pattern of raised fluid elements that rotates around the star with frequency $`\nu _{\mathrm{spin}}`$. Because a large number of fluid elements are scattering radiation to the observer at any given moment, their individual contributions blend together, so the dominant time variation has frequency $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}`$. The effect is somewhat like what one sees when viewing a movie marquee with hundreds of light bulbs blinking in phase so that the pattern formed by the lighted bulbs moves around the marquee. In this case the brightness variation produced by the pattern of lighted bulbs dominates the brightness variation produced by the individual bulbs. The right-hand panel of Figure 2 shows the power spectrum of the flux variation generated in a simulation in which 500 randomly-positioned clumps scatter the radiation pattern coming from the sonic radius. The peak at $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}`$ is clearly dominant.

Variation of the kilohertz QPO frequency difference.—In the model proposed here, the frequency of the upper kilohertz QPO is close to but less than $`\nu _{\mathrm{orb}}(r_{sp})`$ and the difference varies by as much as several tens of Hertz. However, the frequency of the lower kilohertz QPO is very close to either $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}`$ or $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}/2`$, depending on whether the flow at $`r_{sp}`$ is smooth or clumpy. Consequently, the frequency separation between the two QPOs is not exactly $`\nu _{\mathrm{spin}}`$ (or $`\nu _{\mathrm{spin}}/2`$). Table 1 shows the results of fits of this model to the observed changing frequency separation in several sources (see van der Klis et al. 1997; Méndez et al. 1998a; Méndez et al. 1998b; Ford et al. 1998; Psaltis et al. 1998; Méndez & van der Klis 1999; Jonker, Méndez, & van der Klis 2002), adopting the model of Lamb & Miller (2001; note that in that work the frequency of the lower kilohertz QPO also varies by a few tens of Hertz, whereas here only the frequency of the upper kilohertz QPO varies). The fit is reasonable for all these sources except Sco X-1, which has an unusually complicated $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$$`\nu _{\mathrm{QPO2}}`$ relation.

## 4 Discussion and Conclusions

The sonic-point and spin-resonance model proposed here appears able to explain the observed frequencies of the two kilohertz QPOs, the variation of their frequencies by hundreds of Hertz, the closeness of their frequency separation $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ to $`\nu _{\mathrm{spin}}`$ in some sources and to $`\nu _{\mathrm{spin}}/2`$ in others, and the variation of $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ with $`\nu _{\mathrm{QPO2}}`$ (accretion rate) in a given source.

The difference between “fast” and “slow” rotators.—In the model proposed here, $`\nu _{\mathrm{QPO1}}\nu _{\mathrm{orb}}(r_{\mathrm{sp}})\nu _{\mathrm{spin}}/2`$ if the resonance with the stellar spin occurs where the gas in the disk is highly clumped, whereas $`\nu _{\mathrm{QPO1}}\nu _{\mathrm{orb}}(r_{\mathrm{sp}})\nu _{\mathrm{spin}}`$ if the resonance occurs where the gas in the disk is relatively smooth. Magnetic forces may cause the gas in the accretion disk to become more clumped as it approaches the neutron star (MLP98; Lamb & Miller 2001). Consequently, the parameters that may be most important in determining whether the flow at the spin resonance radius $`r_{sr}`$ is clumpy or smooth are the star’s spin frequency and magnetic field. For a given stellar magnetic field, the flow is likely to be more clumpy if the star is spinning rapidly and $`r_{sr}`$ is therefore close to the star. For a given spin rate, the flow is likely to be more clumpy if the star’s magnetic field is stronger.

The four sources with $`\nu _{\mathrm{spin}}>400`$ Hz and measurable frequency separations have $`\mathrm{\Delta }\nu _{\mathrm{QPO}}\nu _{\mathrm{spin}}/2`$ whereas the three sources with $`\nu _{\mathrm{spin}}<400`$ Hz have $`\mathrm{\Delta }\nu _{\mathrm{QPO}}\nu _{\mathrm{spin}}`$. With such a small sample, one cannot make any definite statements, but the apparent trend is consistent with the sonic-point and spin-resonance model. The model suggests that if kilohertz QPOs are detected in the recently-discovered 185 Hz, 191 Hz, and 314 Hz accretion-powered X-ray pulsars XTE J0929$``$314 (Galloway et al. 2002), XTE J1807$``$294 (Markwardt, Smith, & Swank 2003a), and XTE J1814$``$338 (Markwardt & Swank 2003), their frequency separations should be approximately equal to the respective spin frequencies. The 435 Hz spin frequency of XTE J1751$``$305 (Markwardt et al. 2002) is high enough that $`\mathrm{\Delta }\nu _{\mathrm{QPO}}`$ could be either approximately 435 Hz or approximately 217 Hz; QPOs at both frequencies might even be detectable.

Other oscillation frequencies.—We expect the response at the spin-resonance radius defined by $`\nu _{\mathrm{orb}}\nu _{\mathrm{spin}}/2`$ to be largest, but there may be other, weaker motions at other low-order resonances. In principle, weak beat-frequency QPOs could be generated by responses at the radii where $`\nu _{\mathrm{orb}}/\nu _{\mathrm{spin}}`$ is close to ratios of small integers, such as 1/3, 2/3, and so on. The strength of such secondary QPOs would depend on the smoothness of the flow and the rates at which vertical motions are damped at these radii. If the flow is clumpy, so that the primary beat is with the local orbital frequency, then several beat frequencies could be observed. If the flow is instead smooth, so that the primary beat is with the pattern frequency, then because the pattern frequency is always the spin frequency, no other significant QPOs would be expected. In either case, the stronger the damping is, the weaker the other QPOs would be relative to the one at the fundamental spin-resonance radius.

We note that there is no known reason why the mechanism for producing a lower kilohertz QPO proposed in the original sonic-point beat-frequency model would not operate. Apparently this mechanism does not produce a strong QPO in the fast rotators, but it might produce a weak QPO in these sources. If it operates in the slow rotators, it would produce a QPO near $`\nu _{\mathrm{orb}}(r_{sp})\nu _{\mathrm{spin}}`$ that might appear as a sideband to the lower kilohertz QPO.

Although the model outlined here is qualitatively consistent with the basic properties of the kilohertz QPOs, many aspects of the model require further exploration and development. For example, what if any sidebands or additional QPOs are to be expected? How large are the vertical excursions that would be expected in a model that includes more of the physics of accretion disks? Is the model quantitatively consistent with the observational upper limits on X-ray oscillations at $`\nu _{\mathrm{spin}}`$? It will be important to pursue these questions, but in the meantime resonant excitation of gas in the disk appears to be worth exploring as a mechanism for generating the lower kilohertz QPO.

We thank Deepto Chakrabarty and Michiel van der Klis for helpful discussions. This work was supported in part by NSF grant AST 0098399 and NASA grant NAG5-12030 at Illinois, and by NSF grant AST 0098436 at Maryland.