text:
# Heavy-light mesons spectrum from the nonperturbative QCD in the einbein field formalism

## 1 Introduction

In spite of a rather long history of discussion and many theoretical attacks, an interest to the properties of heavy-light mesons is still very high. Indeed, heavy-light systems incorporate both properties of the light quarks which are extremely important for the physics of Chiral Symmetry Breaking (CSB), one of the most challenging phenomenon of the nonperturbative QCD, and heavy degrees of freedom which allow application of such profound methods of investigation as Heavy Quark Effective Theory (HQET) or Operator Product Expansion (OPE). The two main sources of information, experiment and lattice simulations, deliver new data on the heavy-light mesonic spectra, including orbital and radial excitations of the $`q\overline{q}`$ pair, which strongly need theoretical identification and description. In the present paper we calculate the masses of several low lying $`D`$ and $`B`$ mesonic states using the method developed in  based on the Hamiltonian approach to the bound states problem in QCD. Starting with the quarks with current masses we intensively use the so-called einbein fields formalism  which not only considerably simplifies calculations, but allows to justify the treatment of the spin-spin and spin-orbit interactions as perturbations, naturally explaining the appearance of the “constituent” quark mass. Results for the spectra and splittings are compared with the experimental and lattice data and a good agreement is found.

We demonstrate that in our model the proper account of the QCD string dynamics gives extra negative contribution into the energy of the orbital excitations and makes the orbitally excited state ($`l=2`$, $`n=0`$ and $`l=3`$, $`n=0`$) share the region of masses with the first radial excitation ($`l=0`$, $`n=1`$). The latter observation may resolve the problem of the resonance $`D(2637)`$ recently claimed by DELPHI Collaboration  and which is under intensive discussion at present (see $`e.g.`$ ). Indeed, in spite of the fact that in our model both above mentioned states lie somewhat higher than the experimentally observed resonance, the latter can be identified with the orbital excitation of charm (with $`l=2`$ or $`l=3`$) rather then with the first radial excitation ($`n=1`$) solving in such a way the problem of an extremely small width of the resonance inconsistent with any estimates of the $`2S`$-state decays .

## 2 Einbein fields formalism and the mesonic Hamiltonian

The einbein fields formalism which has a long history in literature  was introduced as a method allowing to treat the kinematics of relativistic particles, but then it was generalized for the case of spinning particles and strings . One of its applications is the possibility to develop the canonical quantization procedure $`ala`$ Dirac . In its simplest form this formalism can be applied to the point-like scalar relativistic particle described by the action

$$S=_{t_i}^{t_f}L𝑑tL=m\sqrt{\dot{x}^2}\dot{x}_\mu =\frac{x_\mu }{t}$$
(1)
so that the modified form of (1) looks like

$$L=\frac{\mu \dot{x}^2}{2}\frac{m^2}{2\mu },$$
(2)
where the einbein field $`\mu `$ is introduced. Dynamics defined by the equations of motion for Lagrangians (1) and (2) is the same if the einbein field is treated as an independent degree of freedom and the corresponding constraints are introduced and properly treated. For a short review of the einbein fields formalism see $`e.g.`$  and references therein. In the mean time there exists another approach to the einbein fields  which allows to neglect their dependence on the proper time $`t`$ and thus to treat them as variational parameters to be got rid of in the spectrum rather than in the Hamiltonian. Such an approach was used to a success in a set of papers (see $`e.g.`$ ) and proved to be rather accurate reproducing the relativistic spectra with the error within several per cent. In what follows we accept the above accuracy and treat einbeins variationally.

As it was shown in  within the Vacuum Correlators Method (VCM) , writing the gauge invariant Euclidean Green’s function of the $`q\overline{q}`$ meson

$$G_{q\overline{q}}=\mathrm{\Psi }_{q\overline{q}}^+(\overline{x},\overline{y}|A)^+\mathrm{\Psi }_{q\overline{q}}^+(\overline{x},\overline{y}|A)_A\mathrm{\Psi }_{q\overline{q}}(x_1,x_2|A)=\overline{\mathrm{\Psi }}_{\overline{q}}(x_1)\mathrm{\Phi }(x_1,x_2)\mathrm{\Psi }_q(x_2),$$
(3)
with $`\mathrm{\Phi }(x_1,x_2)`$ being the standard path-ordered parallel transporter operator

$$\mathrm{\Phi }(x_1,x_2)=P\mathrm{exp}\left(ig_{x_2}^{x_1}𝑑z_\mu A_\mu \right)$$

and employing Feynman–Schwinger representation for the single particle Green’s functions and the area law for the isolated Wilson loop bounded by the quark and antiquark trajectories, one can arrive at the following result

$$G_{q\overline{q}}=D\mu _1(t_1)D\mu _2(t_2)D\stackrel{}{x}_1D\stackrel{}{x}_2e^{K_1K_2}Tr[\mathrm{\Gamma }^{(f)}(m_1\widehat{D})\mathrm{\Gamma }^{(i)}(m_2\widehat{\overline{D}})\times $$
(4)
$$P_\sigma \mathrm{exp}\left(_0^T\frac{dt_1}{2\mu _1(t_1)}\sigma _{\mu \nu }^{(1)}\frac{\delta }{i\delta s_{\mu \nu }(x_1(t_1))}\right)\mathrm{exp}(_0^T\frac{dt_2}{2\mu _2(t_2)}\sigma _{\mu \nu }^{(2)}\frac{\delta }{i\delta s_{\mu \nu }(x_2(t_2))})\mathrm{exp}(\sigma S)],$$

where $`\mu _1`$ and $`\mu _2`$ are the einbein fields, $`K_i`$ being the kinetic energies of the quarks

$$K_i=_0^T𝑑t_i\left(\frac{m_i^2}{2\mu _i}+\frac{\mu _i}{2}+\frac{\mu _i\dot{\stackrel{}{x}}_i^2}{2}\right),i=1,2,$$
(5)
$`S`$ stands for the minimal area swept by the quark-antiquark trajectories, $`\sigma _{\mu \nu }=\frac{1}{4i}(\gamma _\mu \gamma _\nu \gamma _\nu \gamma _\mu )`$ and $`\delta /\delta s_{\mu \nu }`$ denotes the derivative with respect to the element of the area $`S`$. The form of the kinetic energy (5) immediately follows from the equation (2) if the laboratory gauge $`t=x_0`$ is fixed in the proper time reparametrization group. Equation (4) allows to find the Hamiltonian of the quark-antiquark pair connected with the relativistic string. If one neglects spin-dependent terms in (4), then the corresponding Hamiltonian coincides by the one for the Nambu-Goto string with scalar quarks at the ends which follows from the Lagrangian

$$L(t)=m_1\sqrt{\dot{x}_1^2}m_2\sqrt{\dot{x}_2^2}\sigma _0^1𝑑\beta \sqrt{(\dot{w}w^{})^2\dot{w}^2w^2}x_{10}=x_{20}=t,$$
(6)
where the minimal string is parametrized by the profile function $`w(t,\beta )`$, the dot and the prime denote its derivatives with respect to the proper time $`t`$ and the string coordinate $`\beta `$ respectively. Thus for the centre-of-mass frame one arrives at

$$H=\underset{i=1}{\overset{2}{}}\left(\frac{p_r^2+m_i^2}{2\mu _i}+\frac{\mu _i}{2}\right)+_0^1𝑑\beta \left(\frac{\sigma ^2r^2}{2\nu }+\frac{\nu }{2}\right)+\frac{\stackrel{}{L}^2}{2r^2[\mu _1(1\zeta )^2+\mu _2\zeta ^2+_0^1𝑑\beta \nu (\beta \zeta )^2]},$$

$$\zeta =\frac{\mu _1+_0^1𝑑\beta \nu \beta }{\mu _1+\mu _2+_0^1𝑑\beta \nu }.$$
(7)
Note that an extra einbein $`\nu `$ following equation (7) has the meaning of the string energy density. Extrema in all three einbeins ($`\mu _1`$, $`\mu _2`$ and $`\nu `$) should be taken in (7). For the case of zero angular momentum $`L=0`$ the last term in (7) vanishes and the extremal value for $`\nu `$ is easily found to be

$$\nu _0=\sigma r,$$
(8)
so that the resulting interaction gives the linearly rising potential $`\sigma r`$. Expansion of the last term in Hamiltonian (7) in powers of $`\sqrt{\sigma }/\mu `$ with the substitution of the extremal value $`\nu _0`$ yields

$$H=H_0+V_{string}$$
(9)
$$H_0=\underset{i=1}{\overset{2}{}}\left(\frac{\stackrel{}{p}^2+m_i^2}{2\mu _i}+\frac{\mu _i}{2}\right)+\sigma r\frac{4}{3}\frac{\alpha _s}{r}C_0,$$
(10)
$$V_{string}\frac{\sigma (\mu _1^2+\mu _2^2\mu _1\mu _2)}{6\mu _1^2\mu _2^2}\frac{\stackrel{}{L}^2}{r}$$
(11)
where we have supplied $`H_0`$ with the constant term $`C_0`$ and the colour Coulomb interaction which comes from the perturbative gluon exchange omitted in (4) for simplicity. Equation (11) gives the first spin-independent correction called the string correction  which is known to reduce the inverse Regge trajectory slope $`dM_L^2/dL=4\sigma `$ specific for the heavy-light limit of the Hamiltonian $`H_0`$, and, when taken at the account at full scale, to lead to the correct string slope $`\pi \sigma `$ .

In what follows we shall consider $`H_0`$ as the zeroth approximation for calculation of the mesonic spectrum. Once the angular momentum $`\stackrel{}{L}`$, the total spin of the quark-antiquark pair $`\stackrel{}{S}`$ and the total momentum $`\stackrel{}{J}=\stackrel{}{L}+\stackrel{}{S}`$ are separately conserved with the Hamiltonian $`H_0`$ then one can specify the mesonic states as terms $`n^{2S+1}L_J`$, where $`n`$ is the radial excitation number.

Let us turn back to the equation (4) and to identify the spin-dependent correction to the zeroth order Hamiltonian $`H_0`$. Following  one finds

$$V_{sd}=\frac{8\pi \kappa }{3\mu _1\mu _2}(\stackrel{}{S}_1\stackrel{}{S}_2)\left|\psi (0)\right|^2\frac{\sigma }{2r}\left(\frac{\stackrel{}{S}_1\stackrel{}{L}}{\mu _1^2}+\frac{\stackrel{}{S}_2\stackrel{}{L}}{\mu _2^2}\right)+\frac{\kappa }{r^3}\left(\frac{1}{2\mu _1}+\frac{1}{\mu _2}\right)\frac{\stackrel{}{S}_1\stackrel{}{L}}{\mu _1}+\frac{\kappa }{r^3}\left(\frac{1}{2\mu _2}+\frac{1}{\mu _1}\right)\frac{\stackrel{}{S}_2\stackrel{}{L}}{\mu _2}$$

$$+\frac{\kappa }{\mu _1\mu _2r^3}\left(3(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})(\stackrel{}{S}_1\stackrel{}{S}_2)\right)+\frac{\kappa ^2}{2\pi \mu ^2r^3}\left(\stackrel{}{S}\stackrel{}{L}\right)(2\mathrm{ln}(\mu r)\gamma _E),\gamma _E=0.57$$
(12)
where $`\kappa =\frac{4}{3}\alpha _s`$ and we have added the term of order $`\alpha _s^2`$ which comes from one-loop calculations and is extensively discussed in literature . Note that up to the last term potential (12) coincides in form with the Eichten–Feinberg–Gromes results , but contains $`\mu `$’s in the denominators instead of masses. In the meantime as clearly seen from (10) $`\mu _i\sqrt{\stackrel{}{p}^2+m_i^2}>m_i`$ or even $`\mu _im_i`$ for light quarks111Let us remind here that these are current quark masses to be denoted $`m_i`$ and to enter equations (1)-(11) (see also Table 1). Einbeins $`\mu `$’s can be viewed as dynamical or “constituent” quark masses that makes it possible to consider (12) as a correction and to justify the expansion made in (11) (see e.g.  for more detailed discussion of this issue).

## 3 Numerical solution

### 3.1 Spectrum of the Hamiltonian $`H_0`$

As mentioned in the introduction einbein fields can considerably simplify the relativistic dynamics. Indeed with the einbeins $`\mu `$’s introduced the kinetic part of the Hamiltonian $`H_0`$ from (10) has a nonrelativistic form, so that the corresponding Schrödinger equation can be written in the reduced dimensionless form (see $`e.g.`$ )

$$\left(\frac{d^2}{d\stackrel{}{x}^2}+|\stackrel{}{x}|\frac{\lambda }{|\stackrel{}{x}|}\right)\chi _\lambda =a(\lambda )\chi _\lambda \lambda =\kappa \left(\frac{2\mu }{\sqrt{\sigma }}\right)^{2/3}$$
(13)
On finding the solutions of equation (13) for $`\chi _\lambda `$ and $`a(\lambda )`$ as functions of the reduced Coulomb interaction strength $`\lambda `$ one has for the extremal values of the einbeins:

$$\mu _1(\lambda )=\sqrt{m_1^2+\mathrm{\Delta }^2(\lambda )}\mu _2(\lambda )=\sqrt{m_2^2+\mathrm{\Delta }^2(\lambda )}\mu (\lambda )=\frac{1}{2}\sqrt{\sigma }\left(\frac{\lambda }{\kappa }\right)^{3/2},$$
(14)
where

$$\mathrm{\Delta }^2(\lambda )=\frac{\sigma \lambda }{3\kappa }\left(a+2\lambda \left|\frac{a}{\lambda }\right|\right).$$

Then the definition of the reduced field $`\mu `$ gives the equation for $`\lambda `$

$$\mu (\lambda )=\frac{\mu _1(\lambda )\mu _2(\lambda )}{\mu _1(\lambda )+\mu _2(\lambda )},$$
(15)
which is the subject to numerical investigation. Unfortunately there is no analytic solution for $`a(\lambda )`$, so one has to generate a selfconsistent solution of both equations, (13) and (15).

Let us give here two more formulae which will be needed for calculating the spin-spin and spin-orbit splittings. For radially excited states one can find

$$\left|\psi (0)\right|^2=\frac{2\mu \sigma }{4\pi }\left(1+\lambda x^2\right),$$
(16)
where

$$r^N=(2\mu \sigma )^{N/3}x^N=(2\mu \sigma )^{N/3}_0^{\mathrm{}}x^{N+2}\left|\chi _\lambda (x)\right|^2,N>32l.$$
(17)
In Table 1 we give the results of numerical calculations with the Hamiltonian $`H_0`$ for $`D`$, $`D_s`$, $`B`$ and $`B_s`$ mesons.

### 3.2 Spin-spin and spin-orbit splittings. Comparison with the experimental and lattice data

In this subsection we calculate the contribution of the $`V_{string}`$ and $`V_{sd}`$ terms. First of all we collect the averaged values of the spin-orbit and spin-tensor interactions between various $`1P`$ and $`1D`$ eigenstates of the Hamiltonian $`H_0`$

$$\begin{array}{ccc}^1P_1|\stackrel{}{S}_1\stackrel{}{L}|^1P_1=0\hfill & ^1P_1|\stackrel{}{S}_2\stackrel{}{L}|^1P_1=0\hfill & ^1P_1|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^1P_1=\frac{1}{4}\hfill \\ & & \\ ^3P_0|\stackrel{}{S}_1\stackrel{}{L}|^3P_0=1\hfill & ^3P_0|\stackrel{}{S}_2\stackrel{}{L}|^3P_0=1\hfill & ^3P_0|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3P_0=\frac{1}{4}\hfill \\ & & \\ ^3P_1|\stackrel{}{S}_1\stackrel{}{L}|^3P_1=\frac{1}{2}\hfill & ^3P_1|\stackrel{}{S}_2\stackrel{}{L}|^3P_1=\frac{1}{2}\hfill & ^3P_1|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3P_1=\frac{1}{4}\hfill \\ & & \\ ^3P_2|\stackrel{}{S}_1\stackrel{}{L}|^3P_2=\frac{1}{2}\hfill & ^3P_2|\stackrel{}{S}_2\stackrel{}{L}|^3P_2=\frac{1}{2}\hfill & ^3P_2|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3P_2=\frac{1}{20}\hfill \end{array}$$
(18)
$$\begin{array}{ccc}^1D_2|\stackrel{}{S}_1\stackrel{}{L}|^1D_2=0\hfill & ^1D_2|\stackrel{}{S}_2\stackrel{}{L}|^1D_2=0\hfill & ^1D_2|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^1D_2=\frac{1}{4}\hfill \\ & & \\ ^3D_1|\stackrel{}{S}_1\stackrel{}{L}|^3D_1=\frac{3}{2}\hfill & ^3D_1|\stackrel{}{S}_2\stackrel{}{L}|^3D_1=\frac{3}{2}\hfill & ^3D_1|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3D_1=\frac{1}{12}\hfill \\ & & \\ ^3D_2|\stackrel{}{S}_1\stackrel{}{L}|^3D_2=\frac{1}{2}\hfill & ^3D_2|\stackrel{}{S}_2\stackrel{}{L}|^3D_2=\frac{1}{2}\hfill & ^3D_2|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3D_2=\frac{1}{4}\hfill \\ & & \\ ^3D_3|\stackrel{}{S}_1\stackrel{}{L}|^3D_3=1\hfill & ^3D_3|\stackrel{}{S}_2\stackrel{}{L}|^3D_3=1\hfill & ^3D_3|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3D_3=\frac{1}{28}\hfill \end{array}$$
(19)
as well as the transition matrix elements

$$\begin{array}{cc}^1P_1|\stackrel{}{S}_1\stackrel{}{L}|^3P_1=\frac{1}{\sqrt{2}}\hfill & ^1P_1|\stackrel{}{S}_2\stackrel{}{L}|^3P_1=\frac{1}{\sqrt{2}}\hfill \\ & \\ ^1D_2|\stackrel{}{S}_1\stackrel{}{L}|^3D_2=\sqrt{\frac{3}{2}}\hfill & ^1D_2|\stackrel{}{S}_2\stackrel{}{L}|^3D_2=\sqrt{\frac{3}{2}},\hfill \end{array}$$
(20)
which lead to mixing of $`|^1P_1`$, $`|^3P_1`$, and $`|^1D_2`$, $`|^3D_2`$ pairs so that the physical state is subject to the matrix equations of the following type:

$$\left|\begin{array}{cc}E_1E& V_{12}\\ V_{12}^{}& E_2E\end{array}\right|=0$$
(21)
The results of our numerical calculations as well as the experimental and lattice data are given in Tables 2,3. The only fitting parameter we use here is the constant $`C_0`$ which takes the following values:

$$C_0(D)=212MeVC_0(D_s)=124MeVC_0(B)=203MeVC_0(B_s)=124MeV.$$
(22)
As clearly seen from (22) the constant $`C_0`$ practically does not depend on the heavy quark and is completely defined by the properties of the light one.

## 4 Discussion and conclusions

Several comments concerning our choice of parameters are in order here. From Table 1 it is seen that we use practically the same values of the strong coupling constant $`\alpha _s`$ for all states. To justify this action let us note that in the case of a light particle moving in the field of a heavy one the OGE interaction between them depends not on the total mass of the system, but rather on its size. Then the slight variation of the strong coupling is due to small differences in the mesons sizes. The exact value of $`\alpha _s`$ is chosen to be near its frozen value  that emphasizes the nonperturbative nature of the processes responsible for the formation of the heavy-light mesons spectrum. The fits are very weakly sensible to variations of the heavy quarks mass which lead mainly to rescaling of the overall shift constant $`C_0`$, so that the latter remains practically the only fitting parameter as stated in the preceding section.

In conclusion we would like to comment on the possible identification of the resonance $`D(2637)`$ recently claimed by DELPHI Collaboration . It was reported to be extremely narrow, about $`15MeV`$, that leads to a confusion as such a small width is inconsistent with the identification of this state as the first radial excitation $`D_{}^{}{}_{}{}^{}(J^P=1^{})`$. It was found in a number of papers  that in spite of the fact that the mass of this state perfectly coincides with the predictions of quark models, all estimates of its width fail to give such a small value. In the meantime it was found that widths of orbitally excited $`D`$ mesons with quantum numbers $`2^{}`$, $`3^{}`$ could be consistent with the reported value , but then the following two points were put forward as main objections: i) a neighboring slightly more massive state should be observed as well, ii) quark models predict orbitally excited mesons to be at least $`50MeV`$ heavier than needed. Here we would like to comment on the second argument. As one can find from the fit given in Table 2, in our model the first radial and the second orbital excitations share the same region of masses. It is not surprise as the negative mass shift of about $`50MeV`$ for the orbitally excited state, missing in standard quark models, is readily delivered by the string correction given by equation (11). This contribution does not affect radially excited states but is significant for orbital excitations. It was demonstrated in  that the account of the proper dynamics of the QCD string in mesons brings the slope of Regge trajectories to their correct values. Appearance of the term (11) in the meson Hamiltonian and its important contribution into the masses of orbitally excited states is yet another reflection of the general situation that the proper QCD string dynamics is extremely important in description of hadronic properties and thus it should be taken into account. So we conclude that the mysterious $`D(2637)`$ state indeed can be identified with the orbital excitation $`2^{}`$ or $`3^{}`$ rather than with the radial one $`0^{}`$, that resolves the second objective above. This statement may hold true for the corresponding states in the $`B`$ meson spectrum, where DELPHI Collaboration also claims a similar state .

As far as the first one is concerned, our model also predicts another orbitally excited state with the mass $`2728MeV`$, $`i.e.`$ $`65MeV`$ higher than $`D_2^{}`$ given in the table, and its experimental grounds are really not clear.

Our predictions for $`D(2654)`$, $`D(2663)`$ and $`D(2664)`$ lie somewhat higher than the experimentally observed value. We find it to be a reflection of a general lack of the “$`\mu `$-technique” used in this paper, which gives larger errors for higher excited states. The difference of about $`1520MeV`$ between the theoretical predictions and the experimental datum $`2637\pm 6MeV`$ can be explained in this way. Development of a systematical approach to the einbein fields as variational parameters could shed light on the sources of systematical errors and possibly to improve the results presented in this paper.

The authors would like to thank Yu.A.Simonov for constant interest to the work and useful discussions, R.N.Faustov for reading the manuscript and for valuable comments, and B.L.G.Bakker for providing the numerical solutions for the reduced Schrödinger equation. Financial support of RFFI grants 00-02-17836 and 00-15-96786 and INTAS-RFFI grant IR-97-232 is gratefully acknowledged.


format: text

text:
# Dispersion-managed soliton in optical fibers with zero average dispersion

## Figure captions:

Fig.1. Time-dependence of first, second, third and fourth iteration (curves 1,2,3,4 respectively).

Fig.2. First iteration $`A^{(1)}(t)`$ (curve 1) versus eq. $`(\text{Dispersion-managed soliton in optical fibers with zero average dispersion})`$ (curve 2).


text:
# Sensitivity of the Radial Velocity Technique in Detecting Outer Planets

## 1 Introduction

Radial velocity surveys of nearby stars have been employed in the search for extra-solar planets for nearly two decades. An important element of these surveys is the sensitivity; i.e. what is the minimum-mass planet that can be detected at a given orbital period? This issue boils down to a question about how well one can detect periodic signals in a finitely sampled data set.

There are two main regimes of analysis: the “short-period” regime, characterized by $`\tau T_0`$ and the “long-period” regime with $`\tau T_0`$. Here $`\tau `$ is the orbital period and $`T_0`$ is the duration of the survey. The planetary objects found to date through the RV technique are short period objects, with periods less than a year or so (e.g. Mayor & Queloz (1995); Marcy & Butler (1996); Fischer et al. (1998)). The sensitivity of RV searches is very well understood in this short orbital period regime and compact analytical formulae exist (Lomb (1976); Scargle (1982); Horne & Baliunas (1986); Nelson & Angel (1998); Cumming, Marcy, & Butler 1999).

According to current theoretical prejudice, however, one expects giant planets to be primarily formed in the colder regions of the proto-planetary nebula, and thus we expect such objects to possess periods in the range of many years to centuries (Boss (1995)). Unfortunately, it is well known that the sensitivity of RV technique is considerably worse for planets with large orbital periods. Earlier analyses in the long-period regime have primarily concentrated on estimating the sensitivity through simulations (e.g. Nelson & Angel (1998), Cumming et al. (1999)).

This paper has two goals: to provide analytical insight into the sensitivity of the RV technique in the long-period regime, and to address the issues of detection and detection efficiency. We note that almost all of the previous work has concentrated on setting upper limits and not addressed the issue of detection. It is also worth noting that much of the previous work has been based on the periodogram method (e.g. Scargle (1982); Horne & Baliunas (1986); Walker et al. (1995); Cumming et al. (1999)). Following, Nelson & Angel (1998) , we argue at various points in the paper why a Least Squares approach is preferable to a periodogram approach. Essentially the Least Squares approach, in contrast to the traditional periodogram approach, offers the most general method and needs no modifications in the long-period regime or in the sparse-data regime.

The plan of the paper is as follows. We summarize the Least Squares approach and derive the basic equations in §2. In §3 we provide analytical estimates for obtaining a detection in the absence of any signal in the short and long-period regimes. In §4 we carry out simulations and obtain estimates for minimum detectable signals in the presence of noise. We conclude in §5.

## 2 Radial Velocity Technique: Basic Equations

For simplicity, we will assume circular orbits throughout this discussion. A planet in a circular orbit undergoes acceleration, and because the linear momentum of the system must be conserved, the star undergoes a reflex acceleration. It is this acceleration that directly informs us of the presence of the companion. However, the observable is the velocity:

$$v(t)=A\mathrm{sin}(2\pi t/\tau +\varphi )+\gamma $$
(1)
where the amplitude, $`A`$ is given by

$$A=\left(\frac{2\pi GM}{\tau }\right)^{\frac{1}{3}}\frac{M_p\mathrm{sin}i}{M_{}}.$$
(2)
Here, $`\gamma `$ is the radial velocity of the planetary system, $`\varphi `$ is the orbital phase, $`\tau `$ is the orbital period, $`M_p`$ is the planet mass, $`M_{}`$ is the stellar mass, $`M=M_p+M_{}M_{}`$ (as $`M_pM_{}`$), and $`i`$ is the inclination of the orbit with respect to the plane of the sky.

The sensitivity of an RV experiment is essentially defined as the minimum mass planetary companion that can be detected at a given period. For a given measurement of $`A`$, we invert Equation 2 to obtain a relation between $`M_p`$ and $`\tau `$:

$$M_p\mathrm{sin}i=A\left(\frac{M_{}^2}{2\pi G}\right)^{\frac{1}{3}}\tau ^{\frac{1}{3}}.$$
(3)
We see that $`M_p\tau ^{\frac{1}{3}}`$, which means that our sensitivity decreases as we look for longer-period orbits. Thus, it is not surprising that the first detection of an extra-solar planet was a planet with very short orbital period (Mayor & Queloz (1995)). There is an additional effect which makes it difficult to identify planets with long periods. With only a baseline of a few years, the RV surveys can only observe some fraction of $`A`$ for orbital periods longer than the observing baseline. In this (long-period) regime, the sensitivity depends critically on the orbital phase. The most sensitive mass estimates are obtained when the acceleration in the radial direction attains the extreme value. In contrast, when the radial acceleration is close to zero, the lack of curvature in the orbit makes the signal difficult to distinguish from the unknown systemic velocity of the system. These two effects combine to explain why all the planetary detections to date have orbital periods of less than 3 years (Mayor & Queloz (1995); Marcy & Butler (1996); Butler & Marcy (1996); Butler et al. (1997); Noyes et al. (1997); Cochran et al. (1997); Marcy et al. (1998); Fischer et al. (1998); Mayor et al. (1998); Marcy et al. (1999); Queloz et al. (1999)).

### 2.1 Least Squares Fitting of Sinusoids

The basic RV analysis consists of fitting the observations to the model specified in Equation 1 and then inferring the mass of the planet\[s\] from the fitted values through Equation 3. As noted by several authors (e.g. Scargle (1982); Horne & Baliunas (1986); Nelson & Angel (1998)) the most optimal fitting method is obtained by using the technique of Least Squares. To enable the use of linear Least Squares fitting techniques we proceed by deriving a linear model equation from Equation 1:

$$v(t)=v_c\mathrm{cos}(\omega t_i)+v_s\mathrm{sin}(\omega t_i)+\gamma ;$$
(4)
here $`v_c=A\mathrm{sin}\varphi `$, $`v_s=A\mathrm{cos}\varphi `$ and $`\omega =2\pi /\tau `$.

There appears to be considerable discussion and debate about the $`\gamma `$ term in the literature (e.g. Walker et al. (1995); Cumming et al. (1999)). The origin of this debate can be traced to two issues. First, it is an assumption of periodogram analyses that the signal is a pure sinusoidal wave (i.e. $`\gamma =0`$). Indeed, the classic periodogram problem is the sine wave (with no offset) buried in zero-mean noise. The second issue is that astronomers are rarely interested in determining the precise radial velocity of a star and its planetary system–i.e. the systemic motion. For this reason, $`\gamma `$ is usually seen as a “nuisance” parameter. In some periodogram analyses, $`\gamma `$ is first determined from the mean of the data and in others $`\gamma `$ is determined in conjunction with $`v_c`$ and $`v_s`$. Cumming et al. (1999) refer to these as “fixed” and “floating” mean methods, respectively while Walker et al. (1995) call them “correlated” and “uncorrelated”. However, our view is that such methods for including $`\gamma `$ (and the additional extension of including a constant acceleration term; see §4.3) only make the periodogram closer to the Least Squares method. Thus, on philosophical grounds of generality we prefer the Least Squares approach. We justify this choice in greater detail below.

Essentially the basic fact is that $`\gamma `$ is needed to represent the physical model correctly. $`\gamma `$ may be an uninteresting parameter, but it is as unknown as $`v_c`$, $`v_s`$ and $`\tau `$. $`\gamma `$ can be dropped from Equation 1 only if it can be demonstrated that it is not covariant with the remaining three parameters, $`v_c`$, $`v_s`$ and $`\tau `$. As shown below, this is not the case and thus one must solve for all the parameters simultaneously. In astrometry, the problem is considerably worse with the position and proper motion being covariant with the parameters of a potential planetary orbit. As with the RV case, one must solve for all unknown parameters simultaneously rather than sequentially, and the periodogram analyses would then include “floating means”and “first derivatives”.

Only in the short-period regime can $`\gamma `$ potentially get decoupled from the other parameters. For this to happen, we need dense sampling over a number of cycles. In the long-period regime, as noted below, the cross-talk never disappears and one must solve for $`\gamma `$, regardless of the density of sampling. The covariance of $`\gamma `$ with orbital parameters is easily seen here since in this regime, we only measure a portion of the orbit and in this limit, the orbit can be approximated by a linear term (constant velocity) and curvature (acceleration). The first term is covariant with the systemic motion. In the astrometry context, Black & Scargle (1982) were the first to recognize the consequences of this covariance.

The parameters of our data set are as follows: the duration of the survey is $`T_0`$, $`v^{}(t_j)`$ is the measured RV at epoch $`t_j`$, and $`n_0`$ is the number of measured epochs. With no loss of generality, we let our time go from $`t=T_0/2`$ to $`t=T_0/2`$. This device allows us to simplify the normal equations (see below). To find the three unknowns, $`v_c`$, $`v_s`$ and $`\gamma `$, we minimize

$$X^2=\underset{i=0}{\overset{n_01}{}}[v^{}(t_i)v_c\mathrm{cos}(\omega t_i)v_s\mathrm{sin}(\omega t_i)\gamma ]^2$$
(5)
with respect to $`v_c`$, $`v_s`$, and $`\gamma `$. This yields a matrix equation for the three unknowns:

$$\left(\begin{array}{ccc}_{i=0}^{n_01}\mathrm{cos}^2(\omega t_i)& _{i=0}^{n_01}\mathrm{cos}(\omega t_i)\mathrm{sin}(\omega t_i)& _{i=0}^{n_01}\mathrm{cos}(\omega t_i)\\ _{i=0}^{n_01}\mathrm{cos}(\omega t_i)\mathrm{sin}(\omega t_i)& _{i=0}^{n_01}\mathrm{sin}^2(\omega t_i)& _{i=0}^{n_01}\mathrm{sin}(\omega t_i)\\ _{i=0}^{n_01}\mathrm{cos}(\omega t_i)& _{i=0}^{n_01}\mathrm{sin}(\omega t_i)& n_0\end{array}\right)\times \left(\begin{array}{c}v_c\\ v_s\\ \gamma \end{array}\right)$$

$$=\left(\begin{array}{c}_{i=0}^{n_01}v^{}(t_i)\mathrm{cos}(\omega t_i)\\ _{i=0}^{n_01}v^{}(t_i)\mathrm{sin}(\omega t_i)\\ _{i=0}^{n_01}v^{}(t_i)\end{array}\right).$$
(6)
In the short-period regime where $`\tau T_0`$, we observe many cycles and thus, under most reasonable sampling schemes, the sinusoidal summations, $`\mathrm{sin}(\omega t)`$, $`\mathrm{cos}(\omega t)`$, and $`\mathrm{sin}(\omega t)\times \mathrm{cos}(\omega t)`$ will average to zero, while $`\mathrm{sin}^2(\omega t)`$ and $`\mathrm{cos}^2(\omega t)`$ will average to 1/2. Thus, in the short-period regime the matrix in Equation 6 becomes diagonal, and the fit parameters are given by

$$v_c=\frac{2}{n_0}\underset{i=0}{\overset{n_01}{}}v^{}(t_i)\mathrm{cos}(\omega t_i),$$
(7)
$$v_s=\frac{2}{n_0}\underset{i=0}{\overset{n_01}{}}v^{}(t_i)\mathrm{sin}(\omega t_i),$$
(8)
$$\gamma =\frac{1}{n_0}\underset{i=0}{\overset{n_01}{}}v^{}(t_i).$$
(9)
In the long-period regime, $`\tau >T_0`$, the matrix is not diagonal since most of the terms do not average to zero. However, by design (and assuming a reasonable sampling scheme), the sinusoidal summations involving only one power of $`\mathrm{sin}(\omega t`$) will average to zero yielding the following normal equations:

$$v_c\underset{i=0}{\overset{n_01}{}}\mathrm{cos}^2(\omega t_i)+\gamma \underset{i=0}{\overset{n_01}{}}\mathrm{cos}(\omega t_i)=\underset{i=0}{\overset{n_01}{}}v^{}(t_i)\mathrm{cos}(\omega t_i)$$
(10)
$$v_s\underset{i=0}{\overset{n_01}{}}\mathrm{sin}^2(\omega t_i)=\underset{i=0}{\overset{n_01}{}}v^{}(t_i)\mathrm{sin}(\omega t_i)$$
(11)
$$v_c\underset{i=0}{\overset{n_01}{}}\mathrm{cos}(\omega t_i)+\gamma n_0=\underset{i=0}{\overset{n_01}{}}v^{}(t_i).$$
(12)
## 3 Type I Errors

Type I errors describe the probability that a high amplitude will be obtained even when no signal is present in the data. We assess the probability of type I errors, assuming that our data set consists of Gaussian noise with a mean of zero and a standard deviation of $`\sigma _0`$.

Short-Period Regime. From Equations 79 we note that for most reasonable sampling schemes, $`v_c`$, $`v_s`$ and $`\gamma `$ are simply sums of Gaussian variables and thus from the Gaussian addition theorem, all three derived parameters are also Gaussian variables. Specifically, $`v_c`$ and $`v_s`$ obey a Gaussian distribution with a mean of zero and a standard deviation of

$$\sigma =\sqrt{\frac{2}{n_0}}\sigma _0,$$
(13)
where $`n_0`$ is the number of measurements. Denoting by $`V_{1s}`$ the value of $`|v_c|`$ (or $`|v_s|`$) that is exceeded in 1% of cases, we note that

$$V_{1s}=2.61\sigma =\frac{3.69}{\sqrt{n_0}}\sigma _0.$$
(14)
From Equation 9 we see that $`\gamma `$ follows a Gaussian distribution with zero mean and a variance of $`\sigma ^2/n_0`$. Denoting $`\mathrm{\Gamma }=|\gamma |`$ we note that the 99th percentile value of $`\mathrm{\Gamma }`$ is $`\mathrm{\Gamma }_{1s}=2.61\sigma _0/\sqrt{n_0}`$. Thus $`\mathrm{\Gamma }_{1s}=0.71V_{1s}`$.

Because we are interested in the fitted amplitude, we will now examine the combined statistics of $`v_c`$ and $`v_s`$. The probability density function for $`v_c`$ and $`v_s`$ is

$$\rho (v_c,v_s)dv_cdv_s=\frac{1}{2\pi \sigma ^2}e^{v_c^2/2\sigma ^2}e^{v_s^2/2\sigma ^2}dv_cdv_s.$$
(15)
Denoting by $`K`$ the square of the inferred amplitude, $`K=v_c^2+v_s^2`$, we note that the probability density function of $`K`$ is an exponential and find the probability to be

$$P(K<K_0)=1e^{K_0/2\sigma ^2}.$$
(16)
Thus, the squared velocity amplitude that is exceeded by pure fluctuations in only 1% of the cases, is

$$K_{1s}=2\sigma ^2\mathrm{ln}\left(\frac{1}{10.99}\right)=\frac{18.42\sigma _0^2}{n_0}.$$
(17)
Long-Period Regime. In this regime, the expressions for $`v_c`$, $`v_s`$ and $`\gamma `$ are not as simple as those in the short-period regime, and so we resort to simulations. To this end, we create a synthetic data set consisting of Gaussian noise with $`\sigma _0=3`$ m s-1, which is the best accuracy thus far obtained for radial velocity measurements of this type (Butler et al. (1996)). We sample the synthetic data at one month intervals for $`T_0=12`$ years. We explore fitted periods from 5 years to 100 years, choosing the interval between sampled periods so as to result in a 1 radian decrease in the total number of orbital cycles over the the length of the observations,

$$\mathrm{\Delta }\tau =\frac{\tau ^2}{2\pi T_0}.$$
(18)
Thus the sequence of the periods which we consider is 60, 64, 69, 74, 80, 87, 95, 105, 117, 132, 152, 177, 212, 261, 337, 463, 699, and 1239 months.

For each period, $`\tau `$, we simulate $`N=1000`$ data sets and carry out the Least Squares fit. We set $`K_1`$ equal to the 10th highest $`K`$ that arises. Clearly, 99% of the fitted $`K`$’s will lie below this value. A plot of $`K_1`$ versus $`\tau `$ is shown in Figure 1. As Figure 1 indicates, our simulations are in excellent agreement with the expected value of $`K_{1s}`$ in the short period regime (Equation 17).

Thus far, our analysis has followed that of Nelson & Angel (1998), but at this point the two analyses diverge. Nelson & Angel (1998) assume that $`K_1`$ in the long-period regime is a power law and from their simulations find $`K_1\tau ^{3.64}`$. The value of this exponent has no natural explanation and indeed it was this unusual value that motivated the present analysis.

The principal goal of our analysis was to understand the behavior of $`K_1`$ in the long-period regime and to derive an analytical expression for $`K_1`$ that is valid for all fitted periods. We begin by examining the covariance matrix of $`v_c`$, $`v_s`$, and $`\gamma `$. At short periods, we find these three fitted parameters are uncorrelated; this is expected from Equations 79.

However, the situation is quite different in the long-period regime as can be seen from the plots in Figure 2. From this figure we conclude the following:
1. $`v_c`$ and $`v_s`$ are correlated.
2. $`v_c`$ and $`\gamma `$ are anti-correlated.
3. $`v_s`$ and $`\gamma `$ are uncorrelated.
The corollary to this result is that phase becomes important in the long-period regime. Indeed, as can be seen from Figure 2, $`\varphi \pm 90^{}`$ is preferred. Nelson & Angel (1998) and Cumming et al. (1999) recognized that the phase would become non-random in the long period regime. However, neither they nor others have explored the full implications of this fact. Below we look into this issue in more detail.

We now provide a simple (physical) explanation of the results displayed in Figure 2. When $`\tau T_0`$, sinusoids with random phase can be fit to the Gaussian data. Clearly, the amplitude of these sinusoids cannot be significantly bigger than the vertical scale of the data (which is approximately $`\sigma _0=3\mathrm{m}\mathrm{s}^1`$). However, when $`\tau >T_0`$, this is no longer necessarily the case. The maximum value of $`K`$ is obtained when we fit a cosine to random data since a cosine is flat around $`t=0`$; we remind the reader of the choice of our time baseline, $`[T_0/2,T_0/2]`$. Recall that for small $`t`$, $`\mathrm{cos}t1t^2/2`$. Thus, the size of the fitted cosine is limited by how much it deviates from a constant in the range from $`t=[0,T_0/2]`$. This deviation is given by

$$1\mathrm{cos}\left(\frac{\pi T_0}{\tau }\right).$$
(19)
The amplitude of this cosine is then chosen so that this deviation from a constant will be roughly equal to the vertical scale of the noise ($`\sigma _0`$). So, Equation 19 tells us the fraction of the total fitted amplitude, and thus the actual fitted amplitude is

$$V_{c_1}=\frac{2V_{1s}}{(1\mathrm{cos}\left[\frac{\pi T_0}{\tau }\right])};\mathrm{for}\tau T_0$$
(20)
here, $`V_{c_1}`$ is the value of $`|v_c|`$ which will be exceeded in 1% of Least Squares fits to Gaussian noise, $`V_{1s}`$ is the corresponding $`V_{c_1}`$ in the short-period regime (Equation 14), and 2 is a normalization factor (since the above deviation was peak-to-peak).

Now, we consider the behavior of $`v_s`$ in the long-period regime. Here, there is a slight difference: because the full amplitude of a sine centered around $`t=0`$ can be observed in the interval $`[\tau /4,\tau /4]`$, the “long-period regime for $`v_s`$” does not actually begin until $`\tau >2T_0`$. However, once we are in the long-period regime, our analysis is much the same. In particular, we find the amplitude of a sine wave such that in the interval from $`t=[T_0/2,T_0/2]`$, the sine wave does not exceed the vertical scale of the noise. Following the same lines as the analysis for $`v_c`$, it is clear that this amplitude must be

$$V_{s_1}=\frac{V_{1s}}{\mathrm{sin}\left(\frac{\pi T_0}{\tau }\right)};\mathrm{for}\tau >2T_0$$
(21)
This function is plotted with the simulated data in Figure 4. Not surprisingly, $`V_{s_1}`$ is considerably smaller than $`V_{c_1}`$.

Understanding the behavior of $`\gamma `$ is not difficult in light of our understanding of $`v_c`$. Since the mean value of the simulated data is zero, our fitted function must also, in the mean, be zero. The anti-correlation between $`\gamma `$ and $`v_c`$ in Figure 2 confirms this expectation. From the above discussion, we know that the fitted signal is primarily a cosine with an amplitude given by Equation 20. We must choose a $`\gamma `$ that translates the fitted cosine such that the result is “centered” around $`y=0`$. This translation is accomplished by subtracting the fitted amplitude $`A`$ (which produces a cosine wave whose maximum value is zero) and then adding back $`\sigma _0/2`$ in order to properly center the fitted signal. Thus, we expect

$$\mathrm{\Gamma }_1=\mathrm{\Gamma }_{1s}\left\{\frac{3}{1\mathrm{cos}\left[\frac{\pi T_0}{\tau }\right]}\frac{1}{2}\right\}.$$
(22)
Note that we expect no correlation between $`v_s`$ and $`\gamma `$, because $`v_s`$ is already centered around zero, and requires no translation (or physically, a sinusoidal orbit with $`\varphi =0^{}`$ is non-degenerate with constant velocity motion).

Equipped with the behaviors of $`V_{c_1}`$ and $`V_{s_1}`$, we are now in a position to write down an analytical expression for $`K_1`$ ($`V_{c_1}^2+V_{s_1}^2`$) at all periods:

$$K_1=\{\begin{array}{cc}K_{1s}=\frac{18.42\sigma _0^2}{n_0}\hfill & \text{for }\tau T_0\hfill \\ \frac{4K_{1s}}{(1\mathrm{cos}\left[\frac{\pi T_0}{\tau }\right])^2}\hfill & \text{for }T_0<\tau <2T_0\hfill \\ \frac{5K_{1s}}{4}\left\{\frac{4}{(1\mathrm{cos}\left[\frac{\pi T_0}{\tau }\right])^2}+\frac{1}{\mathrm{sin}^2\left(\frac{\pi T_0}{\tau }\right)}\right\}\hfill & \text{for }\tau >2T_0\hfill \end{array}.$$
(23)
Here, the factor of 5/4 is chosen to enforce continuity of $`K_1`$ for all periods. Since $`V_cV_s`$ for $`\tau >T_0`$, to a very good approximation, we can express Equation 23 as

$$K_1=\frac{4K_{1s}}{(1\mathrm{cos}\left[\frac{\pi T_0}{\tau }\right])^2};\mathrm{for}\tau >T_0.$$
(24)
We prefer this compact formula over the more exact formula (Equation 23). As demonstrated by Figures 1 and 35, the analytical expressions derived above provide excellent fits to the simulated data. Thus, we have fulfilled our original objective of determining an analytical understanding of $`K_1`$.

We now turn to the importance of the phase term. As can be seen from Figure 2, in the short period regime, the distribution of $`v_c`$ and $`v_s`$ is cylindrically symmetric. However, in the long-period regime the distribution becomes highly elliptical, and the Least Squares fit overwhelmingly prefers orbits with $`\varphi =\pm 90^{}`$. An orbit sampled at a phase of $`\pm 90^{}`$ has the smallest radial acceleration and thus for a given set of measurements yields the largest value of $`K`$.

We have so far considered only the statistics of $`K`$, which is perfectly reasonable when the $`v_c`$$`v_s`$ distribution is cylindrically symmetric. However, in the long-period regime this distribution is elliptical, in which case the exclusive use of the radial parameter $`K`$ is bound not to be optimal. Thus in the long-period regime, we need to look at both the $`v_c`$ and $`v_s`$, or equivalently the phase and amplitude of the fitted parameters.

To this end, we define ellipses in the $`v_c`$$`v_s`$ plane, $`ϵ_1`$, such that 1% of simulated fitted pairs lie outside this ellipse. The parameters of these ellipses are easily obtained since we have analytical expressions for $`V_{c_1}`$ and $`V_{s_1}`$ (Equations 20 and 21). Further discussion and exploitation of this “amplitude-phase” ($`K`$-$`\varphi `$) analysis is postponed to Section 4.

Uneven and Sparse Sampling. We re-iterate that the above results have been derived under the assumption of dense, even sampling. Here, we examine the validity of these results in the case of random sampling, paying particular attention to the regime where $`n_0`$ is small.

As discussed in Section 2.1, Equations 1417 are only valid if the chosen sampling scheme preserves the independence of $`v_c`$, $`v_s`$, and $`\gamma `$. If we have a large number of samples over a number of cycles, then even if our sampling is completely random, we expect that the dense sampling of the sinusoidal summations will preserve the diagonality of Equation 6 (and hence Equations 79 will hold). However, for sparse sampling (i.e. small $`n_0`$), our choice of sampling scheme becomes important. In particular, if we sample the data evenly in the interval $`[T_0/2,T_0/2]`$, then Equation 6 remains diagonal, and the independence of $`v_c`$, $`v_s`$, and $`\gamma `$ is preserved. In contrast, if the data are sparsely and unevenly sampled, then the summations of $`\mathrm{sin}(\omega t)`$, $`\mathrm{cos}(\omega t)`$, and $`\mathrm{sin}(\omega t)\times \mathrm{cos}(\omega t)`$ may no longer average to zero, and $`v_c`$, $`v_s`$, and $`\gamma `$ will be covariant. Physically, this corresponds to the fact that for a small number of randomly spaced samples, the sampled sinusoid might look like a straight line (e.g. if all the samples happen to lie at the zero-crossings of the sine-wave). Thus, for sparse and uneven sampling, Equations 1417 do not accurately describe the sensitivity of the Least Squares fit.

We verify these assertions through simulations: for $`n_0=\{5,\mathrm{},25\}`$, we simulate $`N=1000`$ data sets and do Least Squares fits to sinusoids with $`\tau =0.2T_0`$. The sampling times are drawn from a distribution given by $`t_jϵ[\mathrm{\Delta }t\times (jR),\mathrm{\Delta }t\times (j+R)]`$. Here, $`\mathrm{\Delta }t`$ denotes the sampling interval, given by $`T_0/n_0`$, and $`R`$ is a parameter describing the unevenness of the sampling. $`R=0`$ gives an even sampling scheme, and $`R=1/2`$ gives random sampling. The results of these simulations are shown in Figure 6. We see that for $`n_010`$, one suffers substantial losses in sensitivity if the data are unevenly sampled. This phenomenon was noted by Cumming et al. (1999) . However, the basic reason why we see this phenomenon is because Equation 17 is no longer valid. This discussion re-emphasizes of the generality and robustness of the Least Squares approach over the periodogram approach.

We now repeat the above simulations using fitted sinusoids with $`\tau =10T_0`$, to explore the effects of sparse and uneven sampling in the long-period regime. The results of this simulation are shown in Figure 7. From this figure we see that in the long-period regime, the evenness of the sampling scheme is not nearly as important as it is in the short-period regime. This is to be expected, because our analysis of the behavior of $`K_1`$ in the long-period regime takes into account the covariance between $`v_c`$, $`v_s`$, and $`\gamma `$, and thus we expect it to be applicable to sparsely and unevenly sampled data.

We end this section by acknowledging that our treatment of the statistics is not accurate in the low $`n_0`$ regime. At the very least one needs to be aware of the loss of the degrees of freedom because three parameters are obtained from the data. The correct value of the degrees of freedom will depend on what statistic is being estimated. This issue is not central to the main goal of this paper and we intend to investigate this issue in a later paper.

## 4 Type II Errors

So far we have computed the probability of detecting an apparent signal generated purely by noise. In the language of inference, we have discussed Type I probabilities. Our analysis was centered on obtaining the statistics of $`K`$, the squared amplitude from the fitted parameters $`v_c`$ and $`v_s`$. In particular, we estimated $`K_1`$, the 99th percentile of $`K`$ (see Equation 24). However, following this analysis, we noted that in the long-period regime, $`v_c`$ and $`v_s`$ form an elliptical distribution, and thus merely looking at the statistics of $`K`$, a radial parameter, was not optimal.

We now consider Type II probabilities – the probability of failing to detect a genuine signal due to contamination by noise. The goal of this section is to understand the statistics of $`K`$ in the presence of both signal and noise. To this end, we simulate a data set that consists of signal and noise:

$$v^{}(t_i)=\sqrt{K}\mathrm{sin}(\frac{2\pi t_i}{\tau }+\varphi )+N(t_i),$$
(25)
where $`\sqrt{K}`$ is the amplitude of the signal, and $`N(t_i)`$ is the Gaussian noise. We let $`\varphi `$ be drawn from a uniform distribution in the interval from $`[0,2\pi ]`$, an appropriate assumption for circular orbits. We choose an initial signal amplitude of $`\sigma _0/2`$, and then do $`N=1000`$ Least Squares fits (with the same parameters as in Section 3).

### 4.1 Amplitude-Only Analysis

Following the path used for Type I errors, we will first base our analysis on $`K`$, i.e. we will ignore the issue of the elliptical nature of the $`v_c`$$`v_s`$ distribution. At each period, we determine how many of the fitted amplitudes lie below $`K_1`$. If it is less than $`0.01N`$, then we have found the value of $`K`$ such that 99% of the fitted amplitudes lie above $`K_1`$. We call this $`K_{99}`$. If, however, more than 1% of the fitted amplitudes lie below $`K_1`$, we increment the signal amplitude by $`K^{0.6}/20`$ until we find $`K_{99}`$. A plot of $`K_{99}/K_1`$ versus period is shown in Figure 8.

We note that our choice of 99% confidence for $`K`$ is arbitrary and also very conservative. Most observers would be keen to make a discovery rather than set stringent upper limits. Thus one may wish to consider $`K_{50}`$ (or $`K_{90}`$), which gives the signal amplitude necessary to be detected 50% (or 90%) of the time. Plots of $`K_{90}`$ and $`K_{50}`$ versus period are shown in Figure 9 and Figure 10, respectively.

### 4.2 Amplitude-Phase Analysis

In Section 3, we showed that in the long-period regime, the $`v_c`$$`v_s`$ distribution is elliptical. In view of this, the phase of the signal in Equation 25 is critical. This point is best understood in the results of the simulations displayed in Figure 11 for five cases: (a) no signal, (b)–(d), a signal with amplitude $`\sqrt{K_1}`$ (at that period) and $`\varphi =0^{},45^{},90^{}`$, respectively, and (e) a signal with amplitude of $`\sqrt{K_1}`$ and the phase being randomly (for each simulation) chosen from the range $`[0,2\pi ]`$ (uniform probability density function).

In Figure 11, the amplitude-only analysis would consider all points that lie inside the circle of radius $`K_1`$ to be indistinguishable from those produced by noise. However, in the framework of the $`K`$-$`\varphi `$ analysis, we would consider all points within the ellipse to be indistinguishable from those produced from noise. The superiority of the $`K`$-$`\varphi `$ analysis is evident given that the size of the ellipse is smaller than that of the circle. This superiority is reflected numerically as follows: the fraction of simulations in which the signal is reliably recovered by the $`K`$-$`\varphi `$ analysis is (100%, 100%, 33%, 95%) for cases (b–e) whereas for the amplitude-only analysis the corresponding fractions are (25%, 35%, 33%, 34%).

A comparison of Figures 810 shows the increasing gain of the $`K`$-$`\varphi `$ analysis as the tolerance for committing type II error is increased. When the Type II probability is set to 99%, the $`K`$-$`\varphi `$ analysis lowers the $`K`$ threshold by 20% whereas if we are willing to accept a 50% type II error probability then the $`K`$ threshold is lowered by nearly a factor of 30!

### 4.3 Linear Analysis of Long Term Trends

Above we illustrated the superiority of an amplitude-phase analysis over an amplitude-only analysis. However, some previous authors (Walker et al. (1995); Cumming et al. (1999)) have analyzed long-term trends in the data (i.e. apparent signals with $`\tau T_0`$) by examining the slope of the best-fit straight line, denoted by $`a`$. This best-fit line is subtracted from the data, and an amplitude-only analysis is subsequently performed on the residuals. Here, we discuss this “slope analysis” in the context of the results of Sections 4.1 and 4.2.

To determine the sensitivity of slope analysis, we simulate 1000 data sets with Gaussian noise of zero mean and $`\sigma _0=3\mathrm{m}\mathrm{s}^1`$, sampled at one month intervals for $`T_0=12`$ years. For each simulated data set, we perform a Least Squares fit to a straight line, given by $`at+b`$. The 99th percentile slope, denoted by $`A_1`$, is determined by finding the value of $`|a|`$ such that 99% of simulated data sets yield $`|a|<A_1`$. Next, we determine the type II errors, by injecting a sinusoidal signal into the simulated data (Equation 25) and finding the necessary squared signal amplitude ($`K`$) such that the signal can be reliably recovered. This analysis is the same as that carried out in Section 4.1 except that instead of using the amplitude of fitted sinusoids, we use the slope of fitted straight lines to set confidence limits. Thus, at each sampled period, we can distinguish a real detection from one produced by noise if the fitted $`|a|`$ is greater than $`A_1`$.

The results of these simulations are shown in Figures 810. At sampled periods where the amplitude-only analysis is more sensitive than the slope analysis, only the former is shown. We note that although the slope analysis may be more sensitive than the amplitude-only analysis (in the long-period regime), the amplitude-phase analysis still yields a significant improvement. Moreover, the discrepancy between the amplitude-phase and slope analyses is largest for $`\tau 3T_0`$. We examine this discrepancy further by performing an analysis analogous to that shown in Figure 11e. In particular, for fitted periods of $`\tau =2T_0`$ and $`\tau =10T_0`$, we inject simulated data sets consisting of sinusoids of random phase and $`K=K_1`$ plus Gaussian noise. The percentage of signals recovered by the amplitude-only, amplitude-phase, and slope analyses (respectively) is (57%, 87%, 59%) for $`\tau =2T_0`$, and (34%, 95% 91%) for $`\tau =10T_0`$.

We can understand the relative sensitivities of the $`K`$-$`\varphi `$ and slope analyses by noting that the slope analysis is an explicitly linear technique, and thus it throws away any information that is contained in the curvature of the sine or cosine components of the fitted sinusoid. The $`K`$-$`\varphi `$ analysis, in contrast, utilizes all of the information contained in the linear and curvature components of the fitted sinusoid. These curvature components will tend to zero for $`\tau T_0`$, and thus in this regime the sensitivity of the slope analysis will approach that of amplitude-phase analysis. However, for $`\tau 23T_0`$, these curvature terms are significant, and the $`K`$-$`\varphi `$ analysis yields a substantial improvement over a slope analysis.

The above discussion (in our view) re-emphasizes the generality of the Least Squares approach. While Walker et el. (1995) and Cumming et al. (1999) introduce a modified periodogram to account for a slope in the data, the Least Squares approach can be applied without modification, as it has been in §4.2.

## 5 Conclusions

We have developed an analytical understanding of the sensitivity of the radial velocity technique of planet detection in the regime where the orbital period is longer than the total baseline of the observations. We also examined the sensitivity in the short-period regime, paying particular attention to the case of sparsely sampled data. Moreover, we have illustrated the benefits of Least Squares fitting over the equivalent, but more complicated technique of periodogram analysis; while the periodogram must be modified to deal with long-period signals, or with sparsely sampled data, the Least Squares approach can be applied in its basic form.

We have also discovered the potentially exciting new result that in the long-period regime one obtains additional information from the phase. Analyses of the RV data to date (e.g. Walker et al. (1995); Nelson & Angel (1998); Cumming et al. (1999)) have been based either on amplitude-only analyses, or on amplitude analyses supplemented by slope analyses. However, as dramatically demonstrated by Figures 9 and 10, $`K`$-$`\varphi `$ analysis has significant advantages over these previous analyses, especially in the interesting regime where $`\tau 23T_0`$.

Thus, we propose using a confidence test based on amplitude and phase, characterized by ellipses in the $`v_c`$$`v_s`$ plane, $`ϵ_1`$. The analytical expressions we have developed for $`V_{c_1}`$ and $`V_{s_1}`$ provide the analytical behavior of $`ϵ_1`$ as a function of period. For real RV data, like those mentioned above, we suggest that the fitted $`v_c`$ and $`v_s`$ (at each period) be compared to the $`ϵ_1`$ ellipses. Points lying outside these ellipses would signify the detection of periodic signals.

We have applied this $`K`$-$`\varphi `$ analysis technique to the RV data of Walker et al. (1995). While Nelson & Angel (1998) examined this data with an amplitude-only analysis and reported several marginal detections, our technique yielded several clear detections of periodic signals. Although we cannot say whether these periodicities represent planetary signals (as opposed to stellar cycles or periodic systematic errors), it is clear that the use of amplitude and phase allows previously unknown periodicities to be detected. In the future, we hope to apply this $`K`$-$`\varphi `$ analysis to more comprehensive RV surveys, and to detect long-period companions.