text:
# Robustness of Quadratic Solitons with Periodic Gain

## Abstract

We address the robustness of quadratic solitons with periodic non-conservative perturbations. We find the evolution equations for guiding-center solitons under conditions for second-harmonic generation in the presence of periodic multi-band loss and gain. Under proper conditions, a robust guiding-center soliton formation is revealed.

Multicolor optical soliton formation mediated by cascading of quadratic nonlinearities has been demonstrated experimentally during the last few years in a variety of geometries and frequency-mixing processes, in settings for spatial, temporal and spatio-temporal trapping of light.1 Soliton signals exist in particular in the process of second-harmonic generation (SHG) that is addressed here, where solitons form in waveguides and in bulk crystals by the mutual trapping between the fundamental frequency and second-harmonic waves. Multidimensional soliton families exist above a threshold light intensity for all values of the phase-mismatch between the waves, and most of such solitons have been shown to be dynamically stable under propagation with the equations that model the ideal light evolution under conditions of focused and pulsed SHG. Adiabatic soliton decay and amplification in the presence of weak loss or gain, have been also studied.2-4

In this Letter we address the robustness of quadratic solitons against strong, but periodic non-conservative perturbations. To start the program we consider soliton formation in the presence of multi-frequency losses and large, but rapidly-varying periodic gain. Our goal is to derive the corresponding guiding-center evolution equations and to expose the robustness of the existing solitons under proper conditions. We believe that the results reported bear a generic fundamental interest to the robustness of quadratic soliton formation in structures with periodic non-conservative perturbations. Moreover, they might find direct applications in reduced models of multi-color laser systems with intracavity frequency generation, including self-frequency doubling schemes, operating in the solitonic regime.5,6

Here the focus is on solitons formed in one-dimensional structures under conditions for non-critical type I SHG, but the analysis can be extended to different physical settings. The evolution of the slowly-varying envelopes of the light waves in the presence of multi-frequency band loss and periodic gain can be described by the reduced equations7

$`i{\displaystyle \frac{a_1}{\xi }}{\displaystyle \frac{\alpha _1}{2}}{\displaystyle \frac{^2a_1}{s^2}}+a_1^{}a_2\mathrm{exp}(i\beta \xi )=i\mathrm{\Gamma }_1(\xi )a_1,`$ (1)
$`i{\displaystyle \frac{a_2}{\xi }}{\displaystyle \frac{\alpha _2}{2}}{\displaystyle \frac{^2a_2}{s^2}}+a_1^2\mathrm{exp}(i\beta \xi )=i\mathrm{\Gamma }_2(\xi )a_2,`$ (2)
where $`a_1`$ and $`a_2`$ are the normalized amplitudes of the fundamental frequency (FF) and second-harmonic (SH) waves. In the case of spatial solitons, $`\alpha _1=1`$, and $`\alpha _2=k_1/k_20.5`$, where $`k_\nu `$, with $`\nu =1,2`$ are the linear wave numbers at both frequencies. In the case of temporal solitons, $`\alpha _\nu `$ stand for the ratio between the group-velocity-dispersions existing at both frequencies. The transverse and longitudinal coordinates are normalized to the beam or pulse width and to the diffraction or dispersion length at the fundamental frequency, respectively. The parameter $`\beta `$ is the scaled phase mismatch, and $`\mathrm{\Gamma }_\nu (\xi )`$ stand for periodic gain and loss. Let $`\mathrm{\Gamma }_\nu (\xi )=\mathrm{\Gamma }_\nu ^{(0)}+R_\nu (\xi )`$, where $`\mathrm{\Gamma }_\nu ^0`$ are the average gain or loss at the FF and SH frequencies, and $`R_\nu (\xi )`$ are periodic functions, with period $`L1`$, and zero mean.

To derive the evolution equations of the guiding-center solitons in the presence of the periodic gain, we first use the approach originally introduced by Mollenauer, Evangelides and Haus for the case of Kerr solitons propagating in optical fibers.8 Let $`A_1=a_1`$ and $`A_2=a_2\mathrm{exp}(i\beta \xi )`$. The explicit periodic gain can be removed from equations (1)-(2) by making the transformations $`A_\nu (\xi ,s)=U_\nu (\xi ,s)\mathrm{exp}[_0^\xi R_\nu (\xi ^{})\text{d}\xi ^{}]`$. Substitution in (1)-(2) leads to the evolution equations

$`i{\displaystyle \frac{U_1}{\xi }}{\displaystyle \frac{\alpha _1}{2}}{\displaystyle \frac{^2U_1}{s^2}}i\mathrm{\Gamma }_1^{(0)}U_1+ϵ_2(\xi )U_1^{}U_2=0,`$ (3)
$`i{\displaystyle \frac{U_2}{\xi }}{\displaystyle \frac{\alpha _2}{2}}{\displaystyle \frac{^2U_2}{s^2}}i\mathrm{\Gamma }_2^{(0)}U_2\beta U_2+ϵ_3(\xi )U_1^2=0,`$ (4)
where the resulting longitudinally-varying nonlinear coefficients write

$`ϵ_2(\xi )=\mathrm{exp}\left[{\displaystyle _0^\xi }R_2(\xi ^{})\text{d}\xi ^{}\right],`$ (5)
$`ϵ_3(\xi )=\mathrm{exp}\left[{\displaystyle _0^\xi }\left(2R_1(\xi ^{})R_2(\xi ^{})\right)\text{d}\xi ^{}\right].`$ (6)
To proceed further one now assumes that the wave evolution over a period of the map is mostly dictated by the gain and loss, which are responsible for fast amplitude oscillations of the fields $`A_\nu `$, in addition to the residual effects induced by the nonlinearity. Therefore, $`U_\nu `$ vary slowly over a period of the map, so that one can average equations (3)-(4) to approximately get

$`i{\displaystyle \frac{U_1}{\xi }}{\displaystyle \frac{\alpha _1}{2}}{\displaystyle \frac{^2U_1}{s^2}}i\mathrm{\Gamma }_1^{(0)}U_1+\eta _2U_1^{}U_20,`$ (7)
$`i{\displaystyle \frac{U_2}{\xi }}{\displaystyle \frac{\alpha _2}{2}}{\displaystyle \frac{^2U_2}{s^2}}i\mathrm{\Gamma }_2^{(0)}U_2\beta U_2+\eta _3U_1^20,`$ (8)
where

$$\eta _{2,3}=<ϵ_{2,3}(\xi )>=\frac{1}{L}_0^Lϵ_{2,3}(\xi )\text{d}\xi ,$$
(9)
are the averaged nonlinear coefficients over a period of the map. Therefore, the evolution equations for the slowly-varying averaged fields $`U_{\nu ,0}(\xi ,s)`$, defined as

$$U_{\nu ,0}(\xi ,s)=<\mathrm{exp}[_0^\xi R_\nu (\xi ^{})\text{d}\xi ^{}]>U_\nu (\xi ,s),$$
(10)
write

$`i{\displaystyle \frac{U_{1,0}}{\xi }}{\displaystyle \frac{\alpha _1}{2}}{\displaystyle \frac{^2U_{1,0}}{s^2}}i\mathrm{\Gamma }_1^{(0)}U_{1,0}+U_{1,0}^{}U_{2,0}0,`$ (11)
$`i{\displaystyle \frac{U_{2,0}}{\xi }}{\displaystyle \frac{\alpha _2}{2}}{\displaystyle \frac{^2U_{2,0}}{s^2}}i\mathrm{\Gamma }_2^{(0)}U_{2,0}\beta U_{2,0}+\eta U_{1,0}^20,`$ (12)
where $`\eta =\eta _2\eta _3/\eta _1^2`$, with $`\eta _1=<ϵ_1(\xi )>`$, and

$$ϵ_1(\xi )=\mathrm{exp}\left[_0^\xi R_1(\xi ^{})\text{d}\xi ^{}\right].$$
(13)
One thus concludes that under the conditions where the approximations used to derive (11)-(12) hold, the guiding-center solitons are given by those existing without gain and loss when $`\mathrm{\Gamma }_\nu ^{(0)}=0`$,9 or otherwise with compensated gain and loss,10 but with properly renormalized amplitudes. To expose the value and dependencies of the renormalization coefficient $`\eta `$, consider the illustrative map $`R_\nu (\xi )=G_\nu H(\xi )`$, where $`H(\xi )`$ is the Heaviside function

$$H(\xi )=\{\begin{array}{ccc}1,& & nL<\xi <(n+1/2)L,\hfill \\ 1,& & (n+1/2)L<\xi <(n+1)L.\hfill \end{array}$$
(14)
Substitution in the above expressions and performing all the averages yields

$$\eta =\frac{F[G_2L/2]F[(2G_1G_2)L/2]}{F^2[G_1L/2]},$$
(15)
where

$$F[\mu ]=\frac{1\mathrm{exp}(\mu )}{\mu }.$$
(16)
Of particular interest are the weak maps with $`\mu 1`$, corresponding to $`G_\nu L1`$. In such cases, one might expand (15) in powers of the corresponding parameters to get, at order $`L^2`$

$$\eta 1+\frac{(G_1G_2)^2L^2}{48}.$$
(17)
Note that when $`G_1=G_2`$, equation (15) gives $`\eta =1`$, regardless the value of the map period $`L`$. However, this does not mean that the guiding-center evolution is given by (11)-(12) with $`\eta =1`$ at all orders of $`L`$, because the derivation of (11)-(12), hence the governing equations themselves, are only intended to hold when $`L1`$. To elucidate the applicability limits of (11)-(12), it is worth deriving the guiding-center evolution equations using more mathematically systematic approaches.11,12 Next we outline the outcome of the asymptotic expansion method developed by Kivshar and co-workers for similar problems but for Kerr solitons.12 Such approach has been recently employed to obtain guiding-center evolution equations of light signals propagating in quasi-phase-matched quadratic structures.13

Assuming perfect periodicity of the functions $`R_\nu (\xi )`$ as above, one can express the periodic gain and all fields as Fourier series of the form

$`R_\nu (\xi )={\displaystyle \underset{n0}{}}d_{\nu ,n}\mathrm{exp}(in\lambda \xi ),`$ (18)
$`A_\nu (\xi ,s)={\displaystyle \underset{n}{}}U_{\nu ,n}(\xi ,s)\mathrm{exp}(in\lambda \xi ),`$ (19)
where $`\lambda =2\pi /L`$ and the Foruier coefficients $`U_{\nu ,n}(\xi ,s)`$ vary much slower than the corresponding carrier $`\mathrm{exp}(in\lambda \xi )`$. For the map (14) considered here one has

$$d_{\nu ,2n}=0,d_{\nu ,(2n+1)}=\frac{2G_\nu }{i\pi (2n+1)}.$$
(20)
Assuming the amplification maps to vary rapidly over a diffraction/dispersion lenght, i.e. $`L1`$, so that the spatial frequency $`\lambda `$ is large enough, one can expand the harmonic amplitudes $`U_{\nu ,n0}`$, as power series of $`1/\lambda `$. Namely,

$$U_{\nu ,n0}=\underset{m>0}{}\frac{F_{\nu ,n}^{(m)}}{\lambda ^m}.$$
(21)
Substitution in the governing equations and matching leading-order contributions leads to

$$F_{\nu ,n}^{(1)}=i\frac{d_{\nu ,n}}{n}U_{\nu ,0}.$$
(22)
The guiding-center evolution is obtained by solving recursively for the higher-order contributions $`F_{\nu ,n}^{(m>1)}`$ and substituting the corresponding expressions in the evolution equations for the average fields $`U_{\nu ,0}`$. At order $`1/\lambda ^2`$, one exactly finds (11)-(12), with (17).

Equations (11), (12) and (17) are the central result of this paper. First, they reveal that when $`G_1=G_2`$, at order $`L^2`$ the guiding-center evolution equations are identical to those without the periodic gain, regardless the value of the amplitudes $`G_\nu `$. This is a remarkable result, that emphasizes the robustness of quadratic solitons with the type of perturbations considered. Second, the typical values of the gain coefficient that is obtained in Erbium-doped lithium niobate or potassium titanyl phosphate around the third telecommunication window centered at 1.55 $`\mu `$m fall in the range 0-2 dB/cm.5,14 In the case of spatial solitons, this yields values of the normalized gain coefficient of the order of $`G_\nu 10^1`$.3 With such values and letting $`L10^1`$, one always obtains negligible corrections of the order of $`[\eta 1]10^5`$.

To confirm that under proper conditions (11)-(12) hold, we solved (1)-(2) numerically for different input conditions and maps. Figure 1 shows typical examples of the outcome, when $`\mathrm{\Gamma }_\nu ^{(0)}=0`$. The plots correspond to the phase-mismatch $`\beta =3`$, but analogous results were obtained for other values. To emphasize that guiding-center solitons form with gigantic values of the gain-loss amplitude, provided that the period of the map is small enough so that the guiding-center approach is justified, we display results for a map with $`G_1=20`$, $`G_2=0`$ and $`L=0.2`$. Figures 1(a)-(b) show the propagation of a guiding-center soliton excited by the input $`A_\nu (0)=\stackrel{~}{A}_\nu (0)/\rho _\nu `$, where $`\stackrel{~}{A}_\nu (0)`$ are the corresponding stationary solitons existing without gain and loss with energy flow $`\stackrel{~}{I}=\stackrel{~}{I}_1+\stackrel{~}{I}_2=(|\stackrel{~}{A}_1|^2+|\stackrel{~}{A}_2|^2)\text{d}s=30`$,9 while $`\rho _1=(\eta _2\eta _3)^{1/2}`$, and $`\rho _2=\eta _2`$, are the renormalization factors dictated by (7)-(8). For the map considered one has $`\rho _10.5`$, so that in the renormalized input the FF energy is enhanced by a factor of four. Figures 1(c)-(d) show the excitation of a guiding-center soliton in the same map but with only FF input light carrying the same energy flow as above in the Gaussian shape $`A_1(0)=(2/\pi )^{1/4}I_1^{1/2}\mathrm{exp}(s^2)`$. Here $`I_1=\stackrel{~}{I}_1/\rho _1^2`$.

Figure 2 shows the detailed evolution of the field amplitudes in the case shown in Figs. 1(a)-(b). Because of the large excursions of the amplitudes, the actual evolution differs slightly from that predicted by equations (11)-(12). Such departures are responsible, e.g., of the small resonance peaks appearing in the Fourier spectra of the FF and SH evolving signals, as shown in Figure 3.

When $`\mathrm{\Gamma }_\nu ^{(0)}0`$, but the net loss at one frequency is compensated for by the presence of a net gain at the other frequency band, the system (11)-(12) allows stationary soliton solutions,10 which are chirped. However, an interesting result revealed by the guiding-center evolution equations (11)-(12) is that, in the absence of a net gain or loss (i.e., $`\mathrm{\Gamma }_\nu ^{(0)}=0`$), on average the guiding-center quadratic solitons are chirpless. Accordingly, in contrast to what is found with solitons in other periodic systems, e.g., with dispersion-managed solitons propagating in optical fibers,15 one concludes that with chirpless input conditions the guiding-center quadratic solitons with periodic gain are best excited when the first domain has the whole nominal length. The outcome of our numerical simulations confirms that such is indeed the case, as shown in Figure 4. The plot shows the evolution over a period of the map of the quantities

$$C_\nu (\xi )=\frac{i}{4}\frac{s(A_\nu A_{\nu ,s}^{}A_\nu ^{}A_{\nu ,s})\text{d}s}{s^2|A_\nu |^2\text{d}s},$$
(23)
where $`A_{\nu ,s}`$ stands for $`A_\nu /s`$. The evolution displayed corresponds to a guiding-center soliton similar to that shown in Figs. 1(a)-(b), but for $`G_1=G_2=10`$.

Analogous results than those shown above were obtained for a variety of values of the phase-mismatch $`\beta `$, the input energy flow $`I`$, the signal shape, and the map amplitude. Naturally, this is so provided that the input conditions, the gain amplitude, and the period $`L`$ make a guiding-center approach justified.11-12 Otherwise, e.g., when $`L`$, $`G_\nu `$ or the input $`I`$ are too large, the guiding-center evolution fails, yielding a totally new scenario. In particular, under such conditions resonance phenomena that make the wave propagation unstable can occur. Figure 5 shows a typical example. The plot shows the unstable evolution of the renormalized solution with $`\stackrel{~}{I}=30`$ in a map with $`G_1=2`$, with $`L=1`$.

In conclusion, the evolution equations for guiding-center quadratic solitons propagating in structures with multi-frequency losses and rapidly-varying periodic gain have been presented. Under proper conditions, robust multicolor quadratic soliton formation has been revealed. Results might find applications to reduced models of multicolor laser systems with intracavity frequency generation, including self-frequency doubling schemes, operating in the solitonic regime. Extension of the analysis to general maps that include the details of the laser structures, including the pump light-matter interaction,5,14 is worth investigating.

This work was supported by the Generalitat de Catalunya. Ole Bang acknowledges support from the Danish Technical Research Council through Talent Grant No. 9800400.


format: text

text:
# Microscopic substructure effects in potential-model descriptions of the 7Be(𝑝,𝛾)⁸B reaction

## Abstract

The spectroscopic factor arises from short-range effects in the nuclear wave function. On the other hand, cross sections for external capture reactions, such as 7Be$`(p,\gamma )^8`$B at low energies, depend primarily on the asymptotic normalization of the bound-state wave function – a long-range property. We investigate the relationship between potential models and the full many-body problem to illustrate how microscopic substructure effects arise naturally in the relevant transition matrix element and can be (in part) accounted for by a spectroscopic factor.

The 7Be(p,$`\gamma `$)8B reaction at solar energies ($`E_{cm}`$ $``$ 20 keV) plays an important role in the ‘solar neutrino puzzle’ since the neutrino event rate in the existing chlorine and water Cerenkov detectors is either dominated by or almost entirely due to the high-energy neutrinos produced in the subsequent $`\beta `$ decay of 8B. Direct measurements of the 7Be(p,$`\gamma `$)8B rate, however, are extremely difficult, since the absolute cross sections diminish exponentially at low energies. Thus, theoretical studies of this process become very valuable. In addition to microscopic theories, such as the nuclear shell model or cluster models, potential models provide a popular framework for such investigations. In the latter approach, an average potential is used to generate single-particle wave functions and to calculate various observables; microscopic substructure effects are (partly) accounted for through the use of spectroscopic factors. This strategy, however, has recently been called into question by Csótó Csoto00 , who argued that spectroscopic factors should not be included in potential-model calculations of the 7Be$`(p,\gamma )^8`$B cross section.

Csótó’s argument is based on the observation that the 7Be$`(p,\gamma )^8`$B reaction at low energies depends only on the asymptotic normalization of the 8B bound-state wave function Muk99 , whereas the spectroscopic factor arises from the short-ranged properties of the wave function. However, the unstated assumption underlying this line of reasoning, namely that the short-range and long-range parts of the wave function are independent, is incorrect. For example, the bound state is normalized by $`d^{3N}𝒓|\psi (𝒓)|^2=1`$, a global condition which illustrates that a change in the wave function at any radius affects the normalization integral and thus the value of the wave function everywhere, including in the tail region. Moreover, there exists a more subtle connection between the asymptotic normalization and the short-range behavior of the wave function based on the bound-state Lippmann-Schwinger equation Jennings00 . In fact, the asymptotic normalization can be given in terms of an integral over the interior of the nucleus Jennings00 .

The arguments presented in ref. Csoto00 , as well as the connection between short-range and long-range parts of the wave function that was outlined in ref. Jennings00 , have motivated us to reconsider the proper treatment of microscopic nuclear structure effects in potential-model descriptions of external capture reactions. In this paper we study how potential models are related to the full nuclear many-body problem. We set up a perturbative scheme for the many-body system which, in the lowest-order approximation, gives the potential model and discuss higher-order corrections coming from the microscopic nuclear structure. The spectroscopic factor will be shown to arise naturally in this framework.

The eight-body wave functions for the 8B ground state and the proton-7Be scattering state are written as $`|\psi _8^B(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6,𝝆_7)`$ and $`|\psi _8^k(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6,𝝆_7)`$, respectively. Here $`k`$ is the asymptotic value of the momentum of the scattering particle relative to the 7Be system, and we have introduced the Jacobi coordinates

$`𝝆_j{\displaystyle \frac{_{i=1}^jm_i𝒓_i}{_{i=1}^jm_i}}𝒓_{j+1}\text{for }j\text{=1,}\mathrm{}\text{,8},`$ (1)
with $`𝒓_9=0`$. The dependence of the wave functions on the center-of-mass coordinate, $`𝑹=𝝆_8`$, can be factored out in the form of a plane wave and will be suppressed throughout this paper.

The eigenstates of the seven-body 7Be Hamiltonian, $`|\psi _7^n(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6)`$, form a complete set of wave functions when both bound and scattering states are included. Using this complete set, the eight-body bound-state wave function for 8B can be written as

$`|\psi _8^B(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6,𝝆_7)={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}|\psi _7^n(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6)\varphi _B^n(𝝆_7),`$ (2)
where the summation symbol indicates a sum over the bound states and an integral over the continuum states. The expansion coefficients,

$`\varphi _B^n(𝝆_7)=\psi _7^n(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6)|\psi _8^B(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6,𝝆_7),`$ (3)
are obtained by inverting eq. (2). Note that $`\varphi _B^1(𝝆_7)`$ is the spectroscopic amplitude as defined in ref. Csoto00 . Its norm is simply the spectroscopic factor, $`S=𝑑𝝆_7|\varphi _B^1(𝝆_7)|^2`$ and we will show below that $`\varphi _B^1(𝝆_7)`$ can be identified as a one-body wave function.

The 8B Hamiltonian, $`H_8`$, can be written as $`H_8=H_7\frac{_{𝝆7}^2}{2\mu _7}+_{i=1}^7V(|𝒓_i𝒓_8|)`$, where $`H_7`$ is the Hamiltonian for the 7Be nucleus and $`\mu _7=m_8_{i=1}^7m_i/_{i=1}^8m_i`$ denotes the reduced mass of the eight-body system. Since $`|\psi _7^n`$ and $`|\psi _8^B`$ are eigenstates of $`H_7`$ and $`H_8`$, respectively, we have:

$`\psi _7^n|H_8|\psi _8^B`$ $`=`$ $`E_8^B\psi _7^n|\psi _8^B`$ (4)
$`=`$ $`E_7^n\psi _7^n|\psi _8^B+\psi _7^n|{\displaystyle \frac{_{𝝆7}^2}{2\mu _7}}+{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _8^B.`$
Inserting the expansion given in eq. (2) for $`|\psi _8^B`$, we obtain a set of exact, Schrödinger-like, coupled equations for the coefficients:

$`(E_8^BE_7^n)\varphi _B^n(𝝆_7)={\displaystyle \frac{_{𝝆7}^2}{2\mu _7}}\varphi _B^n(𝝆_7)+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\psi _7^n|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^m\varphi _B^m(𝝆_7).`$ (5)
The normalization condition for the $`\varphi _B^n(𝝆_7)`$ follows from eq. (2):

$`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle d^3𝝆_7|\varphi _B^n(𝝆_7)|^2}=1.`$ (6)
The spectroscopic factor is thus $`S=1_{n=2}^{\mathrm{}}d^3𝝆_7|\varphi _B^n(𝝆_7)|^2`$ and differs from one because of the $`n1`$ terms in the sum.

From eq. (5) we can extract information on the long-range behavior of the 8B ground state wave function. The matrix element $`\psi _7^n|_{i=1}^7V(|𝒓_i𝒓_8|)|\psi _7^m`$, for $`n`$ or $`m`$ corresponding to a 7Be bound state, has a range on the order of the potential convoluted with the bound-state wave function. Thus it falls off rapidly with increasing $`𝝆_7`$. When both $`\psi _7^n`$ and $`\psi _7^m`$ describe continuum states, the matrix element has a long-ranged, but infinitesimal, tail. Solving eq. (5) outside the range of the potential, we find that the $`\varphi _B^n(𝝆_7)`$ decouple and fall off exponentially, $`\varphi _B^n(𝝆_7)\mathrm{exp}[\sqrt{2\mu _7(E_7^nE_8^B)}\rho _7]`$, where $`\rho _7=|𝝆_7|`$. Since $`(E_7^nE_8^B)`$ is smallest for $`n=1`$, the long-range behavior of the 8B wave function is dominated by $`\varphi _B^1(𝝆_7)`$ and the associated asymptotic normalization is given by $`A_n=lim_{\rho _7\mathrm{}}\varphi _B^1(𝝆_7)\rho _7\mathrm{exp}[+\sqrt{2\mu _7(E_7^nE_8^B)}\rho _7]`$. Coulomb effects have been ignored in these considerations.

In order to study the relationship between potential models and the full many-body problem, we consider the above set of coupled differential equations for arbitrary distances. However, we introduce a well-defined approximation scheme. The off-diagonal matrix elements $`V_{nm}\psi _7^n|_{i=1}^7V(|𝒓_i𝒓_8|)|\psi _7^m`$, $`nm`$, will be treated as a perturbation. In zeroth order we take $`V_{nm}=0`$ for $`nm`$, that is, only the diagonal terms survive in eq. (5). The equations for the coefficients $`\varphi _B^n(𝝆_7)`$ decouple and take the form of Hartree equations for the eighth particle moving in the mean field of the other nucleons. The Hartree potential, $`V_{nm}`$, is different for each $`\varphi _B^n(𝝆_7)`$ since the two-body potential is convoluted over different seven-body configurations. Assuming that the ground state of 8B is built on the ground state of 7Be, only $`\varphi _B^1(𝝆_7)`$ can occur in the zeroth-order expansion of $`|\psi _8^B`$, that is, $`|\psi _8^B_0=|\psi _7^1\varphi _B^1(𝝆_7)`$. Coefficients $`\varphi _B^n(𝝆_7)`$, with $`n>1`$, are associated with excited states of the 7Be system and therefore contribute (to lowest order) only to higher-energy states of 8B. Since the 8B states are orthonormal, it follows directly that the $`\varphi _B^n(𝝆_7)`$ which correspond to bound states are normalized to one and those associated with continuum states are normalized to a delta function.

From the considerations above, it becomes apparent that each $`\varphi _B^n(𝝆_7)`$ is normalized and satisfies a single-particle, Schrödinger-like, wave equation with a mean field of the Hartree type. Therefore, the $`\varphi _B^n(𝝆_7)`$ can be identified with potential-model wave functions. We recall that $`\varphi _B^1(𝝆_7)`$ is, by definition, the spectroscopic amplitude, and thus find that the potential-model wave function associated with the ground state of 8B yields the lowest-order approximation to the spectroscopic amplitude.

Since, in lowest order, $`\varphi _B^n(𝝆_7)=0`$ for $`n>1`$, it follows from eq. (5) that there is no first-order correction to the coefficient $`\varphi _B^1(𝝆_7)`$; the next correction to $`\varphi _B^1(𝝆_7)`$ is of second order in the off-diagonal coupling. Nevertheless, coefficients $`\varphi _B^n(𝝆_7)`$ with $`n>1`$ contribute to the 8B ground state in first order, thus the norm of $`\varphi _B^1(𝝆_7)`$ changes from $`S=1`$ to $`S<1`$. The deviation of the spectroscopic factor from one is of second order in the off-diagonal matrix elements $`V_{nm}`$.

The 7Be-p scattering state $`|\psi _8^S`$ can be expanded similarly to the 8B bound state:

$`|\psi _8^S(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6,𝝆_7)={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}|\psi _7^n(𝝆_1,𝝆_2,𝝆_3,𝝆_4,𝝆_5,𝝆_6)\varphi _S^n(𝝆_7).`$ (7)
Proceeding as before, we obtain a set of coupled equations for the expansion coefficients $`\varphi _S^n(𝝆_7)`$:

$`(E_8^SE_7^n)\varphi _S^n(𝝆_7)={\displaystyle \frac{_{𝝆7}^2}{2\mu _7}}\varphi _S^n(𝝆_7)+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\psi _7^n|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^m\varphi _S^m(𝝆_7)`$ (8)
A perturbative treatment of the coupling matrix elements $`V_{nm}`$ is possible here as well: In lowest order the equations decouple and only $`\varphi _S^1(𝝆_7)`$ contributes to $`|\psi _8^S`$. The $`\varphi _S^n(𝝆_7)`$ with $`n>1`$ contribute to the scattering state in first order and the next correction to $`\varphi _S^1(𝝆_7)`$ is of second order in the $`V_{nm}`$. There is no spectroscopic factor associated with the scattering wave function, since $`|\psi _8^S`$ is normalized by its value at asymptotically large distances, rather than by an integral.

Considering the long-range behavior of $`|\psi _8^S`$, we infer the following from eq. (8): If the incident proton has a kinetic energy less than that of the first excited 7Be state ($`E_7^2`$=0.429 MeV Tilley99 ), $`\varphi _S^1(𝝆_7)`$ becomes oscillatory for large $`𝝆_7`$, while the $`\varphi _S^n(𝝆_7)`$ with $`n>1`$ are exponentially damped. Hence, for low reaction energies, the $`n=1`$ term dominates the long-range behavior of the scattering wave function: $`|\psi _8^S|\psi _7^1\varphi _S^1(𝝆_7)`$. At higher energies, $`\varphi _S^n(𝝆_7)`$ with $`n>1`$ display oscillatory behavior and thus a coupled-channels approach to the problem may become necessary.

The dominant contribution to the direct capture process 7Be(p,$`\gamma `$)8B arises from an electric dipole transition between the 7Be-p scattering state and the 8B ground state Kim87 . The relevant matrix element is given by:

$`\psi _8^B|{\displaystyle \underset{i=1}{\overset{8}{}}}e_i𝒓_i|\psi _8^S`$ $`=`$ $`\psi _8^B|{\displaystyle \underset{i=1}{\overset{7}{}}}ϵ_i𝝆_i|\psi _8^S`$ (10)
$`=`$ $`ϵ_7{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle d^3𝝆_7\varphi _B^n(𝝆_7)𝝆_7\varphi _S^n(𝝆_7)}`$
$`+{\displaystyle \underset{\genfrac{}{}{0pt}{}{n=1}{m=1}}{\overset{\mathrm{}}{}}}\psi _7^n|{\displaystyle \underset{i=1}{\overset{6}{}}}ϵ_i𝝆_i|\psi _7^m{\displaystyle d^3𝝆_7\varphi _B^n(𝝆_7)\varphi _S^m(𝝆_7)},`$
where the “Jacobi charges”, $`ϵ_i`$ ($`i=1,\mathrm{},8`$), are defined such that $`_{i=1}^8e_i𝒓_i=_{i=1}^8ϵ_i𝝆_i`$ and $`ϵ_8`$ equals $`Q`$, the total charge of the system. Since the final-state wave function, $`|\psi _8^B`$, is anti-symmetric and the dipole operator, $`_{i=1}^8e_i𝒓_i`$, is symmetric, it is not necessary to anti-symmetrize the initial-state wave function, $`|\psi _8^S`$. Without loss of generality, we can thus take the eighth nucleon to be the incident particle. Note also that the diagonal matrix elements in the last line, $`\psi _7^n|_{i=1}^6ϵ_i𝝆_i|\psi _7^n`$, vanish since $`_{i=1}^6ϵ_i𝝆_i`$ has odd parity and we can choose the 7Be wave functions, $`|\psi _7^n`$, to be states of good parity.

To lowest order in the perturbative expansion for the wave functions only terms with $`n=m=1`$ survive. Thus we obtain:

$`\psi _8^B|{\displaystyle \underset{i=1}{\overset{8}{}}}e_i𝒓_𝒊|\psi _8^Sϵ_7{\displaystyle d^3𝝆_7\varphi _B^1(𝝆_7)𝝆_7\varphi _S^1(𝝆_7)},`$ (11)
where the Jacobi charge

$`ϵ_7=\mu _7\left({\displaystyle \frac{_{i=1}^7e_i}{_{i=1}^7m_i}}{\displaystyle \frac{e_8}{m_8}}\right),`$ (12)
accounts for the fact that the photon couples to both the initial nucleus and to the incident proton ChristyDuck61 . The lowest-order term of the perturbative expansion corresponds to the one-body approximation as given, for example, in ref. ChristyDuck61 .

The first-order corrections to the matrix element of the dipole operator originate entirely from the $`n=1`$, $`m1`$ and $`n1`$, $`m=1`$ terms in the second sum of eq. (10). Note that $`\psi _7^1|_{i=1}^6ϵ_i𝝆_i|\psi _7^m`$, for $`m>1`$, contains dipole excitations of the 7Be core, that is, the transition operator acts on the internal coordinates of the seven-body system. These first-order effects are short-ranged and will be suppressed at low energies, for which the reaction takes place predominately in the tail of wave function.

While there are no first-order corrections to the first sum in eq. (10), there are several second-order effects: First, there are $`n1`$ contributions, which will be suppressed at low energies since $`\varphi _B^n(𝝆_7)`$ and $`\varphi _S^n(𝝆_7)`$ are exponentially damped for large $`𝝆_7`$. Second, there are corrections to both $`\varphi _B^1(𝝆_7)`$ and $`\varphi _S^1(𝝆_7)`$ which are due to the off-diagonal coupling. Since $`\varphi _B^1(𝝆_7)`$ and $`\varphi _S^1(𝝆_7)`$ depend on only one coordinate, these modifications can be accommodated by adjusting the one-body potential in a manner similar to the construction of an optical-model potential. Third, $`\varphi _B^1(𝝆_7)`$ is no longer normalized to one but to the spectroscopic factor, as follows from eq. (6). In addition, there are second-order corrections to the second sum in eq. (10), such as $`n1`$, $`m1`$ contributions. These, again, correspond to dipole excitations of the 7Be core and are thus short-range effects, which are negligible at small energies.

From the above considerations it follows that the spectroscopic factor should be explicitly included in potential-model descriptions of external capture reactions. When the off-diagonal matrix elements $`V_{nm}`$ are treated perturbatively, the spectroscopic factor arises naturally in the normalization of the bound-state wave function at second order and is present at all higher orders. For low-energy reactions first-order contributions and second-order effects not contained in the spectroscopic factor are either small or may be included by an appropriate modification of the one-body potential. At higher energies these first- and second-order corrections may become important and it will no longer be sufficient to accommodate many-body correlations through the use of spectroscopic factors alone.

The arguments for including the spectroscopic factor can be taken beyond perturbation theory. At low energies, eq. (11) will be valid more generally since the neglected terms are associated with short-range effects while the reaction occurs at large distances. Thus we have to evaluate an integral which contains two functions of one variable, a scattering wave function, $`\varphi _S^1(𝝆_7)`$, and a bound-state wave function, $`\varphi _B^1(𝝆_7)`$. Formally, we can construct Schrödinger-like equations for these coefficients using Green functions. To illustrate this, we rewrite the coupled differential equations for the expansion coefficients $`\varphi _\alpha ^n(𝝆_7)`$ (eqs. (5) and (8)) as follows:

$`{\displaystyle \frac{_{𝝆7}^2}{2\mu _7}}\varphi _\alpha ^1(𝝆_7)`$ $`+`$ $`\psi _7^1|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^1\varphi _\alpha ^1(𝝆_7)(E_8^\alpha E_7^1)\varphi _\alpha ^1(𝝆_7)`$ (13)
$`=`$ $`{\displaystyle \underset{m=2}{}}\psi _7^1|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^m\varphi _\alpha ^m(𝝆_7)`$
$`{\displaystyle \frac{_{𝝆7}^2}{2\mu _7}}\varphi _\alpha ^n(𝝆_7)`$ $`+`$ $`{\displaystyle \underset{m=2}{\overset{\mathrm{}}{}}}\psi _7^n|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^m\varphi _\alpha ^m(𝝆_7)(E_8^\alpha E_7^n)\varphi _\alpha ^n(𝝆_7)`$ (14)
$`=`$ $`\psi _7^n|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^1\varphi _\alpha ^1(𝝆_7)`$
Here $`\alpha =B`$ or $`S`$, for the bound and scattering states, respectively. The latter equation is valid for $`n>1`$ and can be expressed as:

$`{\displaystyle \underset{m=2}{\overset{\mathrm{}}{}}}D_{nm}^{(\alpha )}(𝝆_7)\varphi _\alpha ^m(𝝆_7)`$ $`=`$ $`U_n^{(\alpha )}(𝝆_7),`$ (15)
where

$`D_{nm}^{(\alpha )}(𝝆_7)`$ $`=`$ $`\left[{\displaystyle \frac{_{𝝆7}^2}{2\mu _7}}(E_8^\alpha E_7^1)\right]\delta _{nm}+\psi _7^n|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^m,`$
$`U_n^{(\alpha )}(𝝆_7)`$ $`=`$ $`\psi _7^n|{\displaystyle \underset{i=1}{\overset{7}{}}}V(|𝒓_i𝒓_8|)|\psi _7^1\varphi _\alpha ^1(𝝆_7).`$ (16)
Let $`G`$ be the Green function matrix which solves $`_{k=2}^{\mathrm{}}D_{nk}^{(\alpha )}(𝝆_7)G_{(\alpha )}^{km}(𝝆_7,𝝆_7^{})`$ $`=\delta (𝝆_7𝝆_7^{})\delta _{nm}`$. It then follow that

$`\varphi _\alpha ^n(𝝆_7)`$ $`=`$ $`{\displaystyle d^3𝝆_7^{}\underset{k=2}{\overset{\mathrm{}}{}}G_{(\alpha )}^{nk}(𝝆_7,𝝆_7^{})\psi _7^k|\underset{i=1}{\overset{6}{}}V(|𝒓_i𝒓_8|)|\psi _7^1\varphi _\alpha ^1(𝝆_7^{})}`$ (17)
is a solution of eq. (14). Inserting this result into the right-hand side of eq. (13) yields an exact one-body, Schrödinger-like, wave equation for $`\varphi _\alpha ^1(𝝆_7)`$. For the bound-state case, that is, for $`\alpha =B`$, this wave function is normalized to the spectroscopic factor.

To summarize, we have studied the relationship between potential-model descriptions of the capture reaction 7Be(p,$`\gamma `$)8B and the full many-body problem. To this end, we set up a perturbative scheme for the many-body system and identified the potential-model wave function as the zeroth-order approximation to the spectroscopic amplitude. We showed that the lowest-order term in the perturbative expansion of the relevant transition matrix element corresponds to the standard one-body approximation; the spectroscopic factor was found to arise naturally in the normalization of the bound-state wave function at second order. At low energies, where first-order effects and additional second-order corrections are small, the spectroscopic factor allows one to incorporate the most important many-body correlations in potential-model descriptions of external capture reactions and should thus be included explicitly in the calculations. At higher energies, where contributions from the interior of the nucleus become important, additional microscopic nuclear structure effects arise and including the spectroscopic factor will no longer be sufficient.

###### Acknowledgements.

This work was supported by funding from the Natural Sciences and Engineering Research Council of Canada.


text:
# Hydrogen Lyman-𝛼 Absorption Predictions by Boltzmann Models of the Heliosphere

## 1 Introduction

The basic heliospheric structure created by the interaction between the fully ionized solar wind and the partially ionized local interstellar medium (LISM) can to first order be modeled by considering plasma interactions alone, and therefore many models both new and old have been presented that ignore the neutral LISM particles entirely (e.g., Holzer, 1989; Steinolfson, Pizzo, & Holzer, 1994; Wang & Belcher, 1998). However, the charge exchange mechanism, whereby an interstellar neutral loses its electron to a proton, allows the neutrals to take part in the solar-wind/LISM interaction in many important ways (see review by Zank, 1999).

For example, many interstellar neutrals that penetrate the heliopause are ionized by charge exchange with solar wind protons and are then picked up by the outflowing solar wind. These “pick-up ions” dominate the internal energy of the solar wind beyond about 10 AU (Burlaga et al., 1994; Gloeckler et al., 1994), and the acceleration of some of these ions at the termination shock is believed to be the origin of anomalous cosmic rays (Fisk, Kozlovsky, & Ramaty, 1974; Klecker, 1995).

Heliospheric models that treat the neutral gas and plasma in a self-consistent manner predict somewhat different properties than plasma-only models, such as shorter distances to the termination shock, heliopause, and bow shock (Baranov & Malama, 1993, 1995; Pauls, Zank, & Williams, 1995; Zank et al., 1996). These models also suggest that neutral hydrogen in the heliosphere should be very hot, with temperatures of order 20,000–40,000 K. The importance of this prediction is that this high temperature gas should produce H I Lyman-$`\alpha `$ absorption broad enough to be separable from the interstellar absorption observed toward nearby stars, meaning hydrogen in the outer heliosphere can potentially be directly detectable. This absorption has in fact been observed using observations by the Hubble Space Telescope (HST). Linsky & Wood (1996) detected excess Lyman-$`\alpha `$ absorption in HST spectra of the nearby stars $`\alpha `$ Cen A and B. They and Gayley et al. (1997) both demonstrated that the properties of this excess absorption are consistent with a heliospheric origin.

In upwind directions, such as that toward $`\alpha `$ Cen, the heliospheric H I column density is dominated by compressed, heated, and decelerated material just outside the heliopause, which constitutes the so-called “hydrogen wall.” This material accounts for most of the non-LISM absorption observed toward $`\alpha `$ Cen. An additional detection of heliospheric H I absorption only $`12^{}`$ from the upwind direction was provided by HST observations of 36 Oph (Wood, Linsky, & Zank, 2000). For downwind lines of sight the H I density is much lower than in the hydrogen wall, but the sightline through the heated heliospheric H I is longer, potentially allowing heliospheric Lyman-$`\alpha `$ absorption to be observed downwind as well as upwind (Williams et al., 1997). Izmodenov, Lallement, & Malama (1999) report a detection of heliospheric absorption along a downwind line of sight toward the star Sirius.

In addition to detections of hot H I surrounding the Sun, there have also been many reported detections of analogous “astrospheric” material surrounding other Sun-like stars (Wood, Alexander, & Linsky, 1996; Dring et al., 1997; Wood & Linsky, 1998). Wood & Linsky (1998) have illustrated how these observations might potentially be used to infer properties of the winds of these stars and their interstellar environments.

Models that self-consistently treat the neutrals and plasma are essential to help interpret these observations, but such models are a complex theoretical and computational problem. The fundamental difficulty is that neutrals in the heliosphere are far from equilibrium, and their velocity distribution at a given location can be highly non-Maxwellian. One approach is to treat the plasma as one fluid and the neutral hydrogen as three separate fluids, one for each distinct environment in which charge exchange occurs (Zank et al., 1996). Gayley et al. (1997) used such four-fluid models to reproduce the Lyman-$`\alpha `$ absorption observed toward $`\alpha `$ Cen. However, the four-fluid treatment of the problem remains an approximation and may not be entirely successful in reproducing the actual H I velocity distribution (Baranov, Izmodenov, & Malama, 1998), or the expected Lyman-$`\alpha `$ absorption for a given line of sight.

Lipatov, Zank, & Pauls (1998) adopted a particle-mesh method for solving the neutral Boltzmann equation, which should yield accurate particle distribution functions, and Müller, Zank, & Lipatov (2000) further developed this code (see also Müller & Zank, 1999). In this paper, we use heliospheric models created by this Boltzmann code to predict heliospheric Lyman-$`\alpha `$ absorption profiles that can be compared with the HST observations. Models are created assuming different Mach numbers for the inflowing LISM to see which best matches the observed amount of absorption. For the $`\alpha `$ Cen line of sight, a similar analysis has been carried out by Gayley et al. (1997) but using four-fluid rather than Boltzmann models, and assuming different LISM parameters. We compare the absorption predicted by the models with lines of sight observed by HST, including $`\alpha `$ Cen and the aforementioned 36 Oph and Sirius lines of sight, which all show evidence for heliospheric absorption (Wood et al., 2000; Izmodenov et al., 1999).

## 2 The Observations

By now, at least twenty or so useful lines of sight through the LISM have been observed by HST (for a partial list, see Linsky, 1998). Evidence for heliospheric absorption has only been found for three of them, but some of the other lines of sight with low LISM column densities might still be useful for providing upper limits for heliospheric absorption. Thus, in addition to the three lines of sight that suggest heliospheric absorption (36 Oph A, $`\alpha `$ Cen, and Sirius) we also work with three other lines of sight (31 Com, $`\beta `$ Cas, and $`ϵ`$ Eri) that sample different directions through the heliosphere.

In Figure 1, we display the Lyman-$`\alpha `$ lines of these six stars. The $`\theta `$ values shown in the figure are the angles from the upwind direction, which range from the nearly upwind line of sight toward 36 Oph A ($`\theta =12^{}`$) to the nearly downwind line of sight toward $`ϵ`$ Eri ($`\theta =148^{}`$). All of the spectra except that of 36 Oph A were taken using the Goddard High Resolution Spectrograph (GHRS) instrument. Of these spectra, all but that of Sirius utilized the high resolution Echelle-A grating. The Sirius spectrum was taken using the lower resolution G140M grating, explaining the larger binsize of those data. The 36 Oph A data were obtained by the Space Telescope Imaging Spectrograph (STIS), which replaced GHRS in 1997. The resolution of that spectrum, which was observed with STIS’s E140H grating, is comparable to that of the Echelle-A GHRS spectra.

The spectra in Figure 1 are plotted on a heliocentric velocity scale. All show broad, saturated H I absorption near line center, and narrower deuterium (D I) absorption about 80 km s-1 blueward of the H I absorption. Also plotted are the assumed intrinsic stellar Lyman-$`\alpha `$ lines (solid lines) and the best estimates for the interstellar absorption (dotted lines), which fit the data only for 31 Com and $`\beta `$ Cas. These LISM absorption estimates are based on previously published work, which we now discuss briefly.

As mentioned in §1, the $`\alpha `$ Cen data provided the first evidence for heliospheric absorption. Both members of the $`\alpha `$ Cen binary system were observed. The existence of two independent data sets for the same line of sight allowed for a very nice consistency check, and Linsky & Wood (1996) demonstrated that the results of their analysis were the same for both the $`\alpha `$ Cen A and $`\alpha `$ Cen B data. It is the $`\alpha `$ Cen B spectrum that we display in Figure 1, since it has somewhat higher signal-to-noise (S/N) than the $`\alpha `$ Cen A data.

Using the D I line to define the central velocity and temperature of the interstellar material, Linsky & Wood (1996) found they could not fit the $`\alpha `$ Cen H I profile with LISM absorption alone. Excess H I absorption existed on both sides of the Lyman-$`\alpha `$ line (see Fig. 1), especially on the red side. Heliospheric absorption is believed to be responsible for the excess absorption on the red side, and stellar astrospheric absorption may be responsible for the less prominent excess absorption on the blue side (Gayley et al., 1997).

The LISM H I column density toward $`\alpha `$ Cen could not be constrained well at all, so two models of the absorption were presented: one with a deuterium-to-hydrogen (D/H) ratio similar to previously measured values for the local interstellar cloud, $`\mathrm{D}/\mathrm{H}1.5\times 10^5`$ (Linsky, 1998), and one with a much smaller D/H value ($`\mathrm{D}/\mathrm{H}=6\times 10^6`$). In Figure 1, we show the model with the generally accepted D/H value.

The analysis of the 36 Oph data proved to be very similar to that of $`\alpha `$ Cen, with the apparent existence of excess Lyman-$`\alpha `$ absorption on both the red and blue sides of the line (see Fig. 1). Wood et al. (2000) argue that only contributions of both the heliospheric and astrospheric absorption can account for all of the excess absorption.

The Sirius data were first presented by Bertin et al. (1995a, b). The Sirius Lyman-$`\alpha `$ spectrum has lower resolution and lower S/N than the other spectra in Figure 1. Furthermore, unlike the other lines of sight, there are two interstellar components instead of just one. Nevertheless, Bertin et al. (1995a, b) found that if they constrained the H I absorption by assuming the LISM temperature and local D/H ratio reported by Linsky et al. (1993) toward Capella (T=7000 K and $`\mathrm{D}/\mathrm{H}=1.65\times 10^5`$, respectively), they could not account for excess H I absorption on both the blue and red sides of the line (see Fig. 1). They interpreted the larger excess on the blue side to be due to an optically thick stellar wind, and the red side excess to be due to absorption from an evaporative interface between the local cloud and surrounding hot ISM.

Izmodenov et al. (1999) presented a different interpretation of these data. They proposed that the red excess is due to heliospheric absorption and the blue excess is due to analogous astrospheric absorption. They supported this interpretation using a self-consistent kinetic/gasdynamic model, which suggests that a combination of heliospheric and astrospheric material can produce the required amount of absorption. The estimated LISM absorption shown in Figure 1 is from Izmodenov et al. (1999).

The 31 Com, $`\beta `$ Cas, and $`ϵ`$ Eri data shown in Figure 1 were first analyzed by Dring et al. (1997). The intrinsic stellar Lyman-$`\alpha `$ lines we assume in the figure are somewhat different than those derived by Dring et al. (1997), but the LISM parameters used to compute the interstellar absorption profiles are nearly identical. The 31 Com and $`\beta `$ Cas spectra can be fitted beautifully by single interstellar absorption components, but the $`ϵ`$ Eri spectrum is more complex. The H I absorption is significantly blueshifted relative to the D I absorption, indicating a substantial amount of excess H I absorption on the blue side of the line (see Fig. 1). Dring et al. (1997) interpreted this excess absorption to be astrospheric in origin.

## 3 The Models

We construct many models of the heliosphere using the Boltzmann technique described by Lipatov et al. (1998), Müller & Zank (1999), and Müller et al. (2000), and we also present a four-fluid model for comparison (see Table 1). Lyman-$`\alpha `$ absorption is computed for each model for comparison with the HST data. The models assume standard values for the fully ionized solar wind at 1 AU: $`\mathrm{n}(\mathrm{H}^+)=5.0`$ cm-3, $`v=400`$ km s-1, and $`T=10^5`$ K. For the inflowing partially ionized LISM, we assume a velocity of 26 km s-1 and a neutral hydrogen density of $`\mathrm{n}(\mathrm{H})=0.14`$ cm-3, which are consistent with observations (Linsky et al., 2000; Lallement et al., 1995; Quémerais et al., 1994; Gloeckler et al., 1993).

Other input parameters are varied, and the parameters assumed for the various models are listed in Table 1. Observations suggest a likely value of $`\mathrm{T}8000`$ K for the interstellar temperature (Wood & Linsky, 1998; Witte, Banaszkiewicz, & Rosenbauer, 1996) and a value of $`\mathrm{n}(\mathrm{H}^+)0.1`$ cm-3 for the interstellar proton density (Wood & Linsky, 1997), so these are the numbers assumed for most of our models. In Models 8 and 9, we also experiment with smaller values that are still plausible assumptions for the real LISM (see Table 1).

However, the parameter we vary the most is a parameter, $`\alpha `$, relating the proton pressure to the total LISM plasma pressure, $`P=\alpha n(H^+)kT`$. For the known velocity of the inflowing LISM, the Mach number M is related to $`\alpha `$ by $`\mathrm{M}=220(\alpha T)^{0.5}`$, with $`T`$ measured in K. For a pure hydrogen plasma, $`\alpha =2`$, based on the assumption that the electrons and the protons of the plasma have equal number densities and temperatures. The existence of significant cosmic ray and magnetic field pressure in the LISM could potentially increase $`\alpha `$ well above 2. Neither of these additional sources of pressure is well constrained by observation, so we experiment with a wide range of $`\alpha `$’s (see Table 1).

In Figure 2, we display the temperature and density distributions for protons and neutral hydrogen in the heliosphere for Model 1. Panel 2a shows the locations of the heliospheric boundaries, and streamlines in panel 2b indicate the paths of protons in the model. The LISM plasma is heated at the bow shock (BS) and diverted around the heliosphere (see Fig. 2a-b). The tangential discontinuity of the heliopause (HP) separates the LISM plasma from the solar wind. In Model 1, the distances from the Sun to the HP and to the BS in the upwind direction are 100 AU and 260 AU, respectively. Inside the heliosphere, the cool, fast solar wind goes through a termination shock (TS), which is at 62 AU in Model 1. The solar wind is heated to temperatures in excess of $`10^6`$ K and diverted tailward. The neutral temperatures and densities are computed by taking moments of the particle distributions computed by the Boltzmann code (see Müller et al., 2000). In the region between HP and BS, the overdensity in neutral hydrogen referred to as the hydrogen wall is clearly visible in Figure 2d. Downwind of this hydrogen wall, neutral hydrogen is depleted but hot ($`10^5`$ K, see Fig. 2c). Note that if $`\alpha >6`$ then $`\mathrm{M}<1`$ (assuming $`\mathrm{T}=8000`$ K), in which case there will be no bow shock (Models 4–7). A more in-depth description of Models 1–10 and their properties will be provided in a future paper (Müller, Zank, & Wood 2000, in preparation).

The LISM protons are shock heated at the bow shock, and if there is no bow shock they are still heated by adiabatic compression at the heliopause. The interstellar neutrals are then heated by charge exchange with these protons, creating the high temperatures seen in Figure 2c. Charge exchange with outflowing solar wind protons also occurs, resulting in a flux of hot neutrals back through the heliopause and bow shock. Further charge exchanges can deposit additional energy both in the hydrogen wall and in front of the bow shock, and such processes can have measurable effects on the heliospheric structure (Zank et al., 1996). In general, higher heliospheric hydrogen temperatures and densities are observed for models with lower values of $`\alpha `$, which is not surprising since a lower $`\alpha `$ corresponds with a higher Mach number for the inflowing LISM and therefore more heating. Because higher temperatures and column densities will produce broader Lyman-$`\alpha `$ absorption profiles, models with higher Mach numbers will also generally produce more absorption.

The Boltzmann models provide us with velocity distributions throughout the heliosphere that we can use to compute the absorption profiles. In Figure 3, radial velocity distributions are displayed for neutral hydrogen particles in the upwind and downwind directions for Model 1. In order to create these distributions, particles are summed within $`100\times 100`$ AU boxes centered on the locations indicated in the figure. Poissonian error bars ($`N^{0.5}`$) are displayed for each velocity bin. The bin size of the histogram is 6 km s-1. Gaussians have been fitted to the distributions (dashed lines), and the poor fits illustrate the non-Maxwellian character of the distributions.

The complexity of the distributions is due in part to the fact that the neutral particles are created by charge exchange with protons throughout the heliosphere with very different temperatures and flow velocities. The main peak consists mostly of neutral particles created by charge exchange with heated LISM protons in the hydrogen wall. Following the nomenclature used in Zank et al. (1996), these are called “component 1” neutrals. The distributions also contain some particles created by charge exchange with the very hot, decelerated solar wind protons in between the termination shock and heliopause, which for the upwind case in Figure 3 are seen after flowing back across the heliopause. These particles, the “component 2” neutrals, are particularly prevalent in the far wings of the distributions. The component 2 neutrals are far more abundant in downwind directions, explaining why the downwind distribution is broader and has more extended wings than the upwind distribution. Finally, there are a few particles in the distributions created by charge exchange with solar wind protons near the Sun, thereby creating a small population of neutrals with speeds of about $`+400`$ km s-1 (see Fig. 3), which are referred to as “component 3” neutrals. At the upwind location shown in Figure 3, the peak of the distribution lies at about $`16`$ km s-1, which represents a significant deceleration from the $`26`$ km s-1 speed of the undisturbed LISM. At the downwind location shown in the figure, the average speed is about $`+28`$ km s-1, slightly faster than the $`+26`$ km s-1 LISM flow speed.

In Figure 4, we show the H I Lyman-$`\alpha `$ absorption predicted by Model 1 for an upwind and downwind line of sight. In computing such absorption profiles, we first divide the radial path through the heliosphere into many segments. We assign an H I density to each segment based on the number of particles found within a $`4\times 5`$ AU box containing the segment. To get viable statistics for the distribution of radial velocities for the segment, we must use a larger area of influence to determine the distribution histogram. We use a circle with a radius of 30 AU containing the segment for this purpose.

Velocities can be linearly mapped onto wavelengths using the familiar relation

$$\lambda =\lambda _0\left(1\frac{v}{c}\right)^1,$$
(1)
where $`\lambda _0=1215.6701`$ Å is the rest wavelength of Lyman-$`\alpha `$. Thus, the velocity distribution defines a line profile function, $`\varphi _\lambda `$, which is normalized so that

$$_0^{\mathrm{}}\varphi _\lambda 𝑑\lambda =1.$$
(2)
The opacity profile in frequency space is (in cgs units)

$$\tau _\nu =0.02654fN\varphi _\nu ,$$
(3)
where $`f`$ is the oscillator absorption strength and $`N`$ is the column density (i.e., the density times the pathlength of the line segment in question). For Lyman-$`\alpha `$, $`f=0.4164`$ (Morton, 1991). Since $`\varphi _\nu =(\lambda ^2/c)\varphi _\lambda `$, the opacity profile in wavelength space is

$$\tau _\lambda =\frac{0.02654fN\lambda ^2}{c}\varphi _\lambda .$$
(4)
We compute these opacity profiles for each segment of the line of sight and then add them all up. The absorption profile is then simply $`I_\lambda =I_0\mathrm{exp}(\tau _\lambda )`$, where $`I_0`$ is the assumed background Lyman-$`\alpha `$ flux. In Figure 4, we simply assume $`I_0=1`$ for all wavelengths.

In the procedure described above, we have assumed that the line profile is regulated solely by Doppler motions. This is a very good approximation here since the velocity distributions are much broader than the core of the natural line broadening profile, and column densities are too low for the Lorentzian wings of the natural profile to become evident. We have also ignored the fine structure of Lyman-$`\alpha `$, but this has only a very minor effect since the two fine structure components are separated by only 1.3 km s-1.

The H I properties displayed in Figures 2c-d were computed by taking moments of the velocity distribution functions. In Figure 4, we show absorption profiles computed directly from these moments (dashed lines). These profiles, however, are significantly different from those computed directly from the distributions, indicating that flow velocities and temperatures computed from moments do not define the characteristics of the distribution well enough to be used to accurately determine absorption profiles. Accurate absorption profiles can be computed only from the distributions themselves.

## 4 Comparing the Models with the Data

In Figures 5 and 6 we show the heliospheric absorption predicted for the lines of sight observed by HST for Models 1–7. The heliospheric absorption is combined with the LISM absorption (dotted lines) and then convolved with the instrumental line spread function before being plotted. For the sake of clarity, we display four of the models in Figure 5 and four in Figure 6, with Model 4 being shown in both. We have zoomed in on the red side of the absorption profiles since that is where most of the heliospheric absorption is located, astrospheric absorption being the more likely explanation for any excess absorption on the blue side (see §2).

The models certainly have no trouble producing observable absorption. In fact, the models tend to predict too much absorption, especially downwind. This downwind discrepancy may be even more dramatic than the figures suggest, because our models only extend 1000 AU in that direction, which may be far enough for lines of sight with $`\theta 125^{}`$ but is probably not far enough to incorporate all the heliospheric H I further downwind. Thus, we will underestimate the amount of heliospheric absorption toward Sirius and $`ϵ`$ Eri.

The absorption exhibits some interesting nonlinear behavior, which is particularly pronounced in the downwind direction. As expected, the amount of absorption decreases with increasing $`\alpha `$ between $`\alpha =2`$ and $`\alpha =5`$, but the absorption increases between $`\alpha =5`$ and $`\alpha =9.6`$ before decreasing once again between $`\alpha =9.6`$ and $`\alpha =18`$ (see, e.g., the $`ϵ`$ Eri panels of Figs. 5-6). Similar behavior is observed upwind, although it is not as dramatic. The transition between decreasing absorption and increasing absorption at $`\alpha =5`$ is presumably associated with the boundary between a supersonic and a subsonic LISM wind, which is at $`\alpha =6`$.

Unfortunately, none of the models does a very good job in reproducing the observations. Some of the apparent disagreement is illusory, however, because it can be fixed by tweaking either the assumed stellar Lyman-$`\alpha `$ profile or the assumed LISM absorption, or both.

In Figure 7, we show how this is done for two examples: the $`\alpha =2.0`$ model for the 36 Oph line of sight and the $`\alpha =12.5`$ model for the 31 Com line of sight. In both cases, we allow the assumed LISM absorption parameters to be altered within the uncertainties derived in the initial analyses of these data, which are determined partially from the analyses of other LISM lines (Dring et al., 1997; Wood et al., 2000). We also alter the shapes of the assumed stellar profiles to try to make the models fit the data (solid lines: original, dashed lines: altered profiles).

However, there are limits to the shapes a reasonable profile can have. This is especially true for 31 Com. The 31 Com Lyman-$`\alpha `$ emission line is very broad, meaning we cannot allow the slope of the line to become very large within the absorption region. This limits how effective profile changes are in helping us to force the heliospheric absorption model to fit the data. The narrower 36 Oph emission line allows us more leeway in altering the assumed profile. Thus, the changes from the original profiles to the new profiles in Figure 7 are much larger for 36 Oph than for 31 Com.

Another problem for 31 Com is that there is no evidence for excess absorption on the blue side of the line that would suggest the presence of astrospheric absorption, meaning that we cannot completely ignore the quality of fit on that side of the line, in contrast with the 36 Oph example. The quality of the 31 Com fit in Figure 7 is noticeably worse than that of the simpler LISM-only fit in Figure 1. Nevertheless, we conclude that the 31 Com fit in Figure 7 is just good enough to claim that the $`\alpha =12.5`$ model (Model 6) is barely consistent with the 31 Com observations at $`\theta =73^{}`$. Based on Figure 7, we also claim that the $`\alpha =2.0`$ model (Model 1) is consistent with the 36 Oph observations at $`\theta =12^{}`$.

We have performed an analysis like that in Figure 7 for all combinations of models and data, and in Figure 8 we show the best fits that can be obtained by such alterations of stellar profile and LISM absorption for Models 1–7. Many of the models that did not appear to agree well with the data in Figures 5–6 do fit the data well in Figure 8, but many still do not. In general, discrepancies near the base of the absorption (e.g., at $`15<v<25`$ km s-1 for $`\alpha `$ Cen) are much harder to remedy than discrepancies farther in the wing of the absorption (e.g., at $`v>25`$ km s-1 for $`\alpha `$ Cen). Large discrepancies near the base can generally only be fixed by introducing unreasonable fine structure into the assumed stellar Lyman-$`\alpha `$ profile.

Based on the results displayed in Figure 8, we provide in Table 1 our evaluation of which observed lines of sight are inconsistent with which models. None of the models discussed to this point (Models 1–7) are consistent with every line of sight. The low $`\alpha `$ models do better upwind and the high $`\alpha `$ models do better downwind. In general, the problem with the models that do not fit the data is that they predict too much absorption, the exception being the 36 Oph line of sight for which most of the models predict too little absorption. Model 7 is the only model that does not predict way too much absorption along the downwind line of sight to $`ϵ`$ Eri.

For Models 8 and 9 in Table 1, we repeated the $`\alpha =2.0`$ model but with a lower LISM proton density and temperature, respectively, in order to see what effect this has on the predicted absorption (see §3). The n(H+) and T values assumed for these models are not the most likely values based on our current understanding of the LISM, but are still within the realm of possibility (see Table 1). Figure 9 compares the absorption predicted by all the $`\alpha =2.0`$ models. The figure shows that lowering the proton density decreases the amount of absorption downwind but has little effect elsewhere. Decreasing the temperature decreases the amount of absorption in all directions.

The changes in absorption indicated in Figure 9 are fairly modest, however, and certainly not enough to change the results shown in Figure 8, which suggest that $`\alpha =2.0`$ models are consistent only with the 36 Oph line of sight. If lower proton densities and/or temperatures were assumed for the models with other assumed $`\alpha `$’s, it might allow some of them to become consistent with some of the observed lines of sight. The $`\alpha =12.5`$ model could be made consistent with the $`ϵ`$ Eri data, for example (see Fig. 8), but it would worsen the disagreement with the 36 Oph data.

Thus, we conclude that small changes in n(H+) and T will not solve the overall disagreement between the current Boltzmann models and the data. The models predict too much absorption downwind and sidewind relative to the amount of absorption observed upwind. Given the non-linear response of heliospheric models to different boundary conditions, we cannot completely exclude better agreement between data and models in another area of parameter space, perhaps one where the solar wind parameters are different than those that we assume. Moreover, barring unaccounted systematic errors in the assumed stellar profiles and the ISM absorption component, the heliospheric models employed in this study do not yet include all possibilities of interaction in the heliosphere. All of the models neglect the magnetic fields of both the solar wind and interstellar medium, the latter being very poorly constrained observationally. The Boltzmann kinetic models neglect neutral-neutral collisions, as well as proton-neutral collisions that are not accompanied by charge exchange, and neutral depletion due to photoionization close to the Sun.

Ironically, the one model listed in Table 1 that is consistent with all lines of sight is the four-fluid model (Model 10), which uses a less sophisticated treatment of the heliospheric H I than the Boltzmann models. Figure 10 shows the comparison of predicted and observed absorption for Model 10, both before and after the assumed stellar profile and ISM parameters are tweaked to maximize the quality of fit (as was done in Fig. 8 for Models 1–7). The agreement is good enough to claim that the model is consistent with all six lines of sight. Note that the absorption predicted by the four-fluid Model 10 is vastly different from that predicted by the Boltzmann Model 1 (see Fig. 10 and $`\alpha =2.0`$ model in Fig. 5), despite the fact that their input parameters are identical (see Table 1). This illustrates just how sensitive the predicted absorption is to how the H I velocity distributions are treated in the models.

There are several ways to interpret the results illustrated in Figure 10. It should be understood that the four-fluid model and the Boltzmann/kinetic models represent somewhat different physics for the neutral H distributions. The underlying assumption inherent in the four-fluid model is that some scattering of neutral H is present to ensure that the distributions of the individual components, if not the full distribution, is approximately Maxwellian. This process is absent in the kinetic models (see Zank et al. 1996 and Williams et al. 1997 for discussions concerning the four-fluid model assumptions). Figure 10 indicates that a good fit to the Lyman-$`\alpha `$ data along the six sightlines requires that the heliospheric neutral H distribution be very like that produced by the four-fluid model. That Boltzmann Model 1 and the four-fluid Model 10 produce quite different absorption profiles suggests the possibility that additional physical processes with respect to H may need to be incorporated in the Boltzmann codes. Unfortunately, we have not yet explored parameter space in sufficient detail to rule out the possibility that a Boltzmann model does in fact exist that results in acceptable agreement with the Lyman-$`\alpha `$ data, assuming reasonable input parameters. In addition, other factors not included, explicitly or implicitly, in either the Boltzmann or four-fluid models (such as magnetic fields) may still be required to fully address these issues.

## 5 Summary

We have compared H I Lyman-$`\alpha `$ absorption profiles computed using fully self-consistent kinetic/hydrodynamic models of the heliosphere with profiles observed in UV spectra obtained by HST. We make this comparison for six different lines of sight through the heliosphere toward six nearby stars. Our results are summarized as follows:

Most of our models use a Boltzmann particle code for the neutrals, allowing us to estimate neutral velocity distributions throughout the heliosphere. These velocity distributions are far from Maxwellian, having extended wings that lead to broad Lyman-$`\alpha `$ absorption profiles.

Average temperatures and flow velocities can be computed by taking moments of the velocity distributions, and Lyman-$`\alpha `$ absorption profiles can be computed from these moments. However, we find that such profiles do not generally agree well with more accurate profiles computed directly from the distributions. Thus, the moments of the distributions apparently do not provide enough information about their shapes to yield accurate absorption profiles.

The amount of absorption predicted by the models generally increases as the assumed Mach number of the interstellar wind is increased (i.e., the $`\alpha `$ parameter is decreased). However, there are some interesting non-linear effects near the transition between subsonic and supersonic models (i.e., near $`M=1`$).

The Boltzmann models tend to predict too much absorption in sidewind and downwind directions, and too little absorption upwind. This problem is exacerbated at higher Mach numbers. Varying the assumed interstellar temperature and proton density does not seem to help.

In contrast to the Boltzmann models, a four-fluid model that uses a less sophisticated, multi-fluid treatment of the neutrals is consistent with the data. The absorption predicted by this model is very different from that predicted by a Boltzmann model with the same input parameters, emphasizing how crucial the treatment of the neutrals is for accurately predicting Lyman-$`\alpha `$ absorption profiles.

Possible reasons for the lack of success for the Boltzmann models in reproducing the data include the neglect of collisions not involving charge exchange in those models, or perhaps an insufficient exploration of parameter space. There may also be other physical processes that have not been considered in either the Boltzmann or four-fluid models which could be important and affect our results, such as processes involving magnetic fields.

Support for this work was provided by NASA grant NAG5-9041 to the University of Colorado. GPZ and HRM acknowledge the partial support of NASA grant NAG5-6469.


text:
# Beat-Frequency Models of Kilohertz QPOs

## 1. Introduction

Soon after the launch of the Rossi X-ray Timing Explorer (RXTE) in late 1995, observations with it of neutron star low-mass X-ray binaries revealed kilohertz quasi-periodic brightness oscillations (QPOs) in the accretion-powered emission from many of these sources (see, e.g., van der Klis 2000 for a review). These oscillations have high frequencies (up to $``$1300 Hz \[van Straaten et al. 2000\], the highest-frequency astrophysical oscillations ever observed), high amplitudes (up to 15% rms in the 2-60 keV band of the Proportional Counter Array on RXTE), and high coherences (with quality factors $`Q\nu /\mathrm{FWHM}>100`$ in many cases), and often appear as two (but no more) high-frequency oscillations in a single power density spectrum. Beat-frequency models (BFMs) were quickly proposed for this phenomenon (Strohmayer et al. 1996; Miller, Lamb, & Psaltis 1998). In these models, the higher-frequency of the two oscillations is attributed to the orbital frequency of gas at some special radius near the star, and the lower-frequency oscillation is a beat between this orbital frequency and the stellar spin frequency. These models are consistent with many of the trends evident in the early data, including the approximate constancy of the frequency difference between the two simultaneous kilohertz oscillations and the close match in four sources of this frequency difference with the stellar spin frequency inferred from brightness oscillations during thermonuclear X-ray bursts.

It has, however, been established recently that in several sources the frequency difference is not constant, and indeed can vary by more than 50 Hz. Moreover, this variation is systematic: the higher the lower peak frequency, the lower the frequency difference. The explicability of this behavior in the beat-frequency picture has direct bearing on some of the most important inferences drawn from the kilohertz QPOs. For example, only in BFMs is the leveling-off of the frequency of both kilohertz QPOs possibly observed from 4U 1820–30 (Zhang et al. 1998) a signature of the presence of the innermost stable circular orbit (ISCO), a crucial prediction of strong-gravity general relativity. Hence only in BFMs can one infer a high gravitational mass $`M>2.1M_{}`$ for this source, which constrains strongly the equation of state of the high-density matter in the core of neutron stars.

The changing separation frequency observed in several sources provided part of the motivation for the development of other models of the kilohertz QPOs, in particular the relativistic precession models (e.g., Stella & Vietri 1998). In these models the close match between the separation frequency and the spin frequency inferred from burst brightness oscillations is a coincidence, but they do predict the qualitative effect of a separation frequency that drops with increasing kilohertz QPO frequency.

Here we discuss the beat-frequency model in light of these new developments and contrast it with alternate pictures. In § 2 we describe the observational trends that motivated the development of beat-frequency models. We then elaborate on these models, in particular the sonic-point beat-frequency model. In § 3 we discuss the evidence for a changing difference frequency in the four sources Sco X-1, 4U 1608–52, 4U 1728–34, and 4U 1735–44. We show that an aspect of the sonic-point beat-frequency model, included in the dynamics but originally omitted from the frequency estimates, naturally accommodates the changing difference frequency and can quantitatively fit the data. Finally, in § 4 we contrast some of the predictions of the beat-frequency model with the predictions of the relativistic precession model, and discuss analysis that might be done with current data to help discriminate between the two interpretations. We also explore the impact of future observations, both with the upcoming generation of high spectral resolution satellites (such as Chandra and XMM) and with longer-term projects such as Constellation-X and a hypothetical high-area follow-on to RXTE.

## 2. Motivation for Beat-Frequency Models

As discussed in the introduction, soon after the discovery of kilohertz QPOs it was established that these oscillations have (1) high frequency, (2) high amplitude, and (3) high coherence, and that there are always two or fewer kilohertz QPOs in a given power density spectrum. In addition, the separation frequency appeared consistent with constant in many sources and close to the spin frequency inferred from burst brightness oscillations in the four sources where this could be tested. It is now known (see § 3) that in several, and perhaps all, sources, the separation frequency is not constant (although it is still close to the inferred spin frequency), and in fact decreases systematically with increasing lower peak frequency. In § 3 we discuss how this new result may be interpreted within the beat-frequency model.

The high frequency indicates that the source of the brightness oscillations is close to the neutron star. A natural candidate for these oscillations is the orbital frequency at some special radius. Given that the burst oscillation frequency is most convincingly interpreted as the stellar spin frequency or its first overtone (see, e.g., Strohmayer & Markwardt 1999), the close match of the separation frequency with the inferred stellar spin frequency suggests a sideband relation between the two simultaneous kilohertz QPOs. The lack of a third kilohertz QPO suggests a beat-frequency relation, because most other mechanisms will produce both an upper sideband and a lower sideband. For example, amplitude modulation of one frequency by another will produce two sidebands of equal strength. For all of these reasons, beat-frequency mechanisms are natural to propose for the kilohertz QPO phenomenon.

The qualitative match of beat-frequency expectations with the data is not enough to accept this model: it is also important to establish a reason why a particular radius in the accretion disk would be selected. The high coherence of the oscillations requires, in this model, that the radial range from which the oscillations are generated has a fractional width less than $`1/Q`$, or less than $``$1% in the most stringent cases. Moreover, the high amplitudes observed in some sources indicates that the luminosity cannot be generated in just the small range of radii where the frequency is determined.

For these reasons, the sonic-point beat-frequency model was proposed (Miller, Lamb, & Psaltis 1998). In this model, the special radius at which the frequencies are generated is where the inward radial velocity of the gas in the disk increases rapidly with decreasing radius. As described in Miller, Lamb, & Psaltis (1998), this rapid increase in inward radial velocity is usually caused by radiation drag, which removes angular momentum from the accreting gas and hence causes it to spiral inwards. If for some reason radiation drag is not effective in a particular instance, as might happen if the accretion rate is so high that the optical depth from the stellar surface is large, then a similar rapid increase in velocity will nonetheless occur near the ISCO, where the inspiral of gas will open up simply because of gravitational effects. The radial velocity typically goes from subsonic to supersonic in this transition, and it is therefore convenient to label it $`r_{\mathrm{sonic}}`$, the sonic point. However, in this model it is only the rapid change in velocity, and not the specific fact that the gas crosses a sonic point, that is important.

In this model we assume that there are dense clumps of gas orbiting at many radii near the star, and as we show in Miller, Lamb, & Psaltis (1998) only the clumps near $`r_{\mathrm{sonic}}`$ will produce sharp brightness oscillations. As gas streams from the clump onto the star, it produces a bright pattern of impact that rotates with the clump (see Miller, Lamb, & Psaltis 1998), and hence a distant observer sees a modulation of the flux from the system as the bright impact spot passes into and out of view. This occurs at approximately the orbital frequency of the clumps at $`r_{\mathrm{sonic}}`$. If the intensity of radiation from the stellar surface is also modulated at the stellar spin frequency, as it will be if a weak stellar magnetic field funnels extra matter towards one or two magnetic poles, then this variation in intensity modifies the mass accretion rate from the clumps, and hence changes the total luminosity from the impact spots. As shown in Miller, Lamb, & Psaltis (1998), this creates an observed luminosity modulation from the system at approximately the beat frequency between the orbital frequency of the clumps and the stellar spin frequency. In this model the luminosity in the oscillations is generated at the stellar surface, where most of the gravitational energy is released, and hence the oscillations can have high amplitudes.

We now demonstrate that in this model the slow inward drift of the clumps, included in the dynamical calculations but not the frequency estimates of Miller, Lamb, & Psaltis (1998), generally produces a separation frequency between the pair of kilohertz QPO peaks that is less than the stellar spin frequency.

## 3. Changing Difference Frequency

### 3.1. Observational Evidence

By late 1996, observations of Sco X-1 by van der Klis and colleagues (van der Klis et al. 1997; see also Méndez & van der Klis 2000) demonstrated that the separation frequency between the kilohertz QPO peaks in this source is definitely not constant, and in fact that the separation clearly decreases with increasing lower peak frequency. This is contrary to the simplest expectations of the beat frequency model, in which the clumps generating the frequency are orbiting at a constant radius and hence the separation is extremely close to the stellar spin frequency.

Many explanations for this behavior were discussed in the community. Most of them centered on the apparently unique nature of Sco X-1 in this respect and hinged on other of its properties such as its very high (near-Eddington) luminosity at the point when the change in the difference frequency was most pronounced. These explanations included jet models (van der Klis et al. 1997) and other unpublished ideas that involved effects such as the finite thickness of the accretion disk. These models all had different physical bases and even different mathematical models for the changing difference frequency, but were all able to fit the relatively smooth $`\mathrm{\Delta }\nu `$ vs. $`\nu _{\mathrm{lower}}`$ curve for Sco X-1.

In the last year or so, however, careful analyses by Méndez and colleagues have shown that Sco X-1 is not unique in having a changing difference frequency. Specifically, the low-luminosity sources 4U 1608–52 (Méndez et al. 1998), 4U 1728–34 (Méndez & van der Klis 1999), and 4U 1735–44 (Ford et al. 1998) also have difference frequencies that drop with increasing lower peak frequency, and the separation frequency in 4U 1636–536 is slightly but significantly less than the spin frequency inferred from burst brightness oscillations (Méndez, van der Klis, & van Paradijs 1998). More generally, Psaltis et al. (1998) showed that, within the observational uncertainties, all of the kilohertz QPO sources have separation frequencies that are consistent with this trend. This is therefore clearly not an effect that requires a high luminosity, and more general explanations must be sought.

### 3.2. Explanation Within the Beat-Frequency Model

As described in detail in Lamb & Miller (2000), and as we now discuss, in the sonic-point model the inward drift of the clumps tends to increase the frequency of the lower kilohertz peak and decrease the frequency of the upper kilohertz peak. This produces a difference frequency that is generally less than the stellar spin frequency, in agreement with observations.

If a clump were to orbit in a perfect circle, then the inspiral time and phase traversed during inspiral of gas from the clump to the stellar surface would be the same regardless of when the gas separated from the clump. In this case, consider two successive “beats”, i.e., two successive maxima of the mass flow rate from the clump. These occur a time $`\mathrm{\Delta }T=1/(\nu _{\mathrm{orb}}\nu _{\mathrm{spin}})`$ apart. Because the inspiral time to the surface is the same for each of the two beats, the arrival time at the surface is separated by exactly that same $`\mathrm{\Delta }T`$, and the observed frequency of the beat is therefore just $`1/\mathrm{\Delta }T=\nu _{\mathrm{orb}}\nu _{\mathrm{spin}}`$.

Now consider clumps that are drifting in slowly. Then the inspiral time decreases as time progresses. If we again think of two successive beats, suppose that the inspiral time for the second beat is a small time $`\delta t`$ less than the inspiral time for the first beat. Then the interval between releases of gas from the clump in the first and second beats is still $`\mathrm{\Delta }T=1/(\nu _{\mathrm{orb}}\nu _{\mathrm{spin}})`$, but the arrival times at the surface are now separated by $`\mathrm{\Delta }T\delta t`$. Therefore the observed frequency, which depends on surface arrival times because the surface is where the luminosity is generated, is $`1/(\mathrm{\Delta }T\delta t)`$, which is higher than it would have been without the inward drift. This is somewhat analogous to Doppler blueshifting. Further analysis and comparison with numerical simulations shows (Lamb & Miller 2000) that the upper peak frequency is also affected. In fact, the inward drift decreases the upper peak frequency, typically by a fractional amount that is $``$50% of the fractional amount by which the lower peak frequency is increased compared to the circular orbit approximation. Both the increase in the lower peak frequency and the decrease in the upper peak frequency act to decrease the frequency separation.

As described in Lamb & Miller (2000), the magnitude of the change in the frequencies depends on the details of the inspiral, and in particular on how the total phase and time of the inspiral depend on the radial location of the clump. However, much of the physics is captured in a simple model in which the frequency change depends only on the ratio of two radial velocities: the velocity $`v_{\mathrm{clump}}`$ of the clump near the sonic point and a characteristic velocity $`v_{\mathrm{gas}}`$ of the gas when it leaves the clump (see Lamb & Miller 2000 for more details). The observed frequencies are then approximately $`\nu _{\mathrm{low}}=\nu _{\mathrm{beat}}/(1v_{\mathrm{clump}}/v_{\mathrm{gas}})`$ and $`\nu _{\mathrm{high}}=\nu _{\mathrm{orb}}(1\frac{1}{2}v_{\mathrm{clump}}/v_{\mathrm{gas}})`$, where $`\nu _{\mathrm{beat}}\nu _{\mathrm{orb}}\nu _{\mathrm{spin}}`$ and $`\nu _{\mathrm{orb}}`$ is the orbital frequency at the sonic point. Numerical simulations show that the characteristic velocity $`v_{\mathrm{gas}}`$ of the gas is fairly insensitive to radius, and hence for simplicity we assume that $`v_{\mathrm{gas}}`$ is independent of radius.

What remains is to produce a physical model for $`v_{\mathrm{clump}}`$ as a function of radius. The radial velocity of a clump with angular momentum $`J`$ under the influence of a torque $`N`$ is $`v_r=N/(J/r)`$. In a Schwarzschild spacetime, for clumps whose mass does not depend on the distance from the star,

$$J/r\frac{1}{2}\frac{r6M}{\left(r3M\right)^{3/2}},$$
(1)
where we use geometrized units in which $`Gc=1`$. The singularity at $`r=6M`$ in the equation for the radial velocity is only apparent. In reality, the velocity deviates from this simple formula near the ISCO. Numerical simulations show (Lamb & Miller 2000) that for plausible torques the radial velocity of the clump can be small enough that even near the ISCO the inspiral time of the clumps is large, and hence the coherence of the brightness oscillations can be as high as observed. The simulations also show that for small enough torques the separation frequency can remain almost constant even near the ISCO. Therefore, in this model the observations of 4U 1820–30, which indicate an approximately constant separation frequency even when the kilohertz QPO frequencies may be leveling off, are still consistent with the asymptotic frequency being approximately the orbital frequency at the ISCO (Lamb & Miller 2000).

As a particular model of the torque, assume that the clumps interact only weakly with each other and with any background gas in the accretion disk, so that the primary torque is exerted by the stellar magnetic field. The general form of the torque at radius $`r`$ is then $`NB^2Ar`$, where $`A`$ is the cross-sectional area of the clump and $`B`$ is the strength of the stellar magnetic field at the location of the clump. Assume also that the clump mass is approximately independent of the distance from the star, and that the clumps are in pressure equilibrium with the external magnetic field. Suppose that the clumps are kept at the local Compton temperature, which is roughly independent of the radius. Then $`nkTB^2r^6`$ for a dipolar field, meaning that the diameter of the clumps scales as $`r^2`$ and the cross-sectional area scales as $`r^4`$. Then the torque scales as $`Nr^1`$.

This model can then be fit to the data, using the relations between radial velocity and upper and lower QPO frequency given above. For a source such as Sco X-1 which has an unknown spin frequency, there are three parameters: the stellar gravitational mass $`M`$, the spin frequency $`\nu _{\mathrm{spin}}`$, and the coefficient of the torque. For a source such as 4U 1728–34 with a known $`\nu _{\mathrm{spin}}`$=364 Hz, there are two free parameters. Figure 1 shows the results of the fitting. It is encouraging that one physical model can fit both of these sources as well as it does, especially given the high quality of the data for Sco X-1 (this model also fits 4U 1608–52 and 4U 1735–44).

## 4. Tests of Models With Current and Future Data

With the immense archive of data accumulated with RXTE, there are a number of further analyses that could in principle be performed to discriminate between models or test predictions of models. In this section we discuss such tests, particularly as they might discriminate between beat-frequency models and relativistic precession models, and also look to the future to consider the qualitatively new types of observations that may be available with the upcoming generation of X-ray satellites.

### 4.1. Tests With Current Data

An important difference between the predictions of the beat-frequency model and the relativistic precession model concerns expectations for the value of the stellar spin frequency. In the beat-frequency model it is expected that the stellar spin frequency is close to (although slightly larger than) the separation frequency, and is also very close to the burst oscillation frequency or half the burst frequency, depending on the symmetry. In contrast, as recently emphasized by Psaltis & Norman (2000), in the simplest version of the relativistic precession model the data on the lower frequency horizontal branch type brightness oscillations combined with the kilohertz oscillations requires that the stellar spin frequency exceed $``$500 Hz in all sources, and be especially large (up to $``$ 900 Hz) in the Z sources. Simple estimates of corrections due to fluid viscosity (Psaltis & Norman 2000) suggest that these corrections would increase further the required spin frequency. The predictions for $`\nu _{\mathrm{spin}}`$ therefore differ greatly in the two pictures, and if a clear spin frequency is detected in the persistent, accretion-powered emission then this will have strong bearing on the modeling. In addition, if it happens that the spin frequency is not equal to the burst oscillation frequency or half of it, this will overturn the apparently compelling arguments that burst brightness oscillations can only be explained via a rotational modulation mechanism.

The two models also differ in their predictions of the frequencies of other, weak brightness oscillations. In the relativistic precession model, there is an upper sideband of the orbital frequency expected (Psaltis & Norman 2000). The amplitude of this oscillation need not be as great as the amplitude at the lower sideband; for example, if the underlying perturbations of the accretion disk cut off sharply at frequencies exceeding the observed orbital frequency then the upper sideband could be weaker. Nonetheless, the oscillation is expected to be there at some level. In the beat frequency model there is an oscillation expected at twice the observed lower peak frequency. This is expected for several reasons. For example, the waveform generated by the beat is not expected to be perfectly sinusoidal, and hence overtones will be generated. In addition, if the radiation pattern rotating with the star, which helps produce the beat, is itself not perfectly sinusoidal, then it will generate weak overtones of the lower peak frequency.

The predicted frequencies discussed above depend only on the observed frequencies. If the observed frequencies are $`\nu _{\mathrm{low}}`$ and $`\nu _{\mathrm{high}}`$, then the predicted frequency in the relativistic precession model is $`2\nu _{\mathrm{high}}\nu _{\mathrm{low}}`$ and the predicted frequency in the beat frequency model is $`2\nu _{\mathrm{low}}`$. The frequencies will change with time for a given source, but these predictions are ideally suited for testing by use of the shift-and-add method first used by Méndez et al. (1998) to discover a weak brightness oscillation in 4U 1608–52.

### 4.2. Tests With Future Data and Satellites

The launch of a new generation of high spectral resolution satellites such as Chandra and XMM opens up new ways to probe the strong gravity and dense matter of neutron stars. For example, suppose that simultaneous observations of a kilohertz QPO source are performed with RXTE and XMM. If an Fe K$`\alpha `$ profile from the inner edge of the nearly-circular flow is detected and characterized with XMM, this profile gives the radius of that inner edge in units of the gravitational mass of the neutron star. If at the same time a pair of kilohertz QPOs is detected with RXTE, then the frequency of the upper peak, combined with the radius in units of the gravitational mass, yields the mass of the star. This is true independent of where the inner edge is; that is, it need not be at the innermost stable circular orbit. Therefore, such simultaneous measurements could provide a clean way to measure the gravitational masses of neutron stars in low-mass X-ray binaries.

Future, high-area missions such as Constellation-X or a hypothetical high-area follow-up timing mission to RXTE could provide even more information, by allowing characterization of the waveforms of the brightness oscillations seen in accretion-powered and burst-powered emission. An example is shown in the left panel of Figure 2, which shows theoretical waveforms for burst brightness oscillations. The curve for the more compact star is broader, as is expected due to the extra gravitational light deflection. In addition, the curve for the larger (less compact) star is more asymmetric, due to the greater Doppler shifts from the surface rotation velocity. The waveform therefore encodes information about both the mass and radius of the star, meaning that repeated observations of bursts from a single source can constrain the mass and radius tightly, with consequent constraints on the equation of state of the high-density matter in the core of neutron stars.

The right panel of Figure 2 shows an example of the constraints possible with Constellation-X (outer contour, at 1$`\sigma `$) and a hypothetical 10 m<sup>2</sup> future timing instrument (inner contour, at 1$`\sigma `$). Clearly, waveform fitting can in principle yield very precise information about the mass and radius of individual neutron stars, and therefore about the equation of state of matter at supranuclear densities.

High-area timing missions also will potentially allow us to do qualitatively new types of tests of the strong-gravity predictions of general relativity. For example, we have just discussed two independent ways of estimating the gravitational mass of a neutron star: by combining Fe K$`\alpha `$ profiles with kilohertz QPOs, and by waveform fitting of burst brightness oscillations. If either of these is successful for an individual source, then we can predict precisely the orbital frequency at the innermost stable circular orbit. With high-area timing instruments we expect to see more cases like 4U 1820–30, in which there is leveling-off of the QPO frequency and thus evidence for the ISCO. The match of this asymptotic frequency with the frequency predicted from the gravitational mass and general relativity will provide us with unprecedented quantitative tests of general relativity in strong gravity. In conclusion, therefore, the continued qualitative and quantitative agreement of the beat-frequency model with observations of kilohertz QPOs has not only yielded important new constraints on the equation of state of the dense matter in the core of neutron stars and, possibly, the first direct evidence for unstable orbits around neutron stars, a key prediction of general relativity. It also indicates strongly that future observations of these sources, especially with high-area timing missions, will allow us to continue to make qualitative leaps in our observational understanding of strong gravity and dense matter.

## ACKNOWLEDGEMENTS

We thank Michiel van der Klis and Tod Strohmayer for discussions about possible information from high-area timing missions. This research was supported in part by NASA ATP grant number NRA-98-03-ATP-028.

## REFERENCES

Ford, E. C., van der Klis, M., van Paradijs, J., Méndez, M., Wijnands, R. A. D., & Kaaret, P. 1998, ApJ, 508, L155

Lamb, F. K., & Miller, M. C. 2000, in preparation

Méndez, M., & van der Klis, M. 1999, ApJ, 517, L71

Méndez, M., & van der Klis, M. 2000, MNRAS, in press (astro-ph/0006243)

Méndez, M., van der Klis, M., & van Paradijs, J. 1998, ApJ, 506, L117

Méndez, M., van der Klis, M., Wijnands, R. A. D., Ford, E. C., van Paradijs, J., & Vaughan, B. A. 1998, ApJ, 505, L23

Miller, M. C., Lamb, F. K., & Psaltis, D. 1998, ApJ, 508, 791

Psaltis, D. et al. 1998, ApJ, 501, L95

Psaltis, D., & Norman, C. 2000, ApJ, submitted (astro-ph/0001391)

Stella, L. & Vietri, M. 1998, ApJ, 492, L59

Strohmayer, T. E., & Markwardt, C. B. 1999, ApJ, 516, L81

Strohmayer, T. E., Zhang, W., Swank, J. H., Smale, A., Titarchuk, L., & Day, C. 1996, ApJ, 469, L9

van der Klis, M. 2000, ARA&A, in press (astro-ph/0001167)

van der Klis, M., Wijnands, R. A. D., Horne, K., & Chen, W. 1997, ApJ, 481, L97

van Straaten, S., Ford, E. C., van der Klis, M., Méndez, M., & Kaaret, P. 2000, ApJ, submitted (astro-ph/0001480)

Zhang, W., Smale, A. P., Strohmayer, T. E., & Swank, J. H. 1998, ApJ, 500, L171