text:
# Unconditional Security of the Bennett 1992 quantum key-distribution over lossy and noisy channel

## I Introduction

Quantum key distribution (QKD) is one of the most interesting topics in quantum information processing, which allows the sender (Alice) and the receiver (Bob) to share a secret key with negligibly leakage of information to an eavesdropper (Eve). This task is accomplished by making use of the properties of quantum mechanics, such as no-cloning theorem noclone , indistinguishability of nonorthogonal quantum states, entanglement, and so on peres . The B92 protocol B92 , where only two nonorthogonal states are used, is thought to be secure due to the fact that Eve cannot clone and distinguish the states deterministically.

Although the qualitative reason of the security is simple, it is quite hard to show the quantitative security against so-called coherent attack, which is the most general attack allowed by quantum mechanics. In lc98 ; sp00 , the security of the BB84 protocol BB84  is proven by showing the relationship between QKD and other important protocols in quantum information, such as the entanglement distillation protocol (EDP) EPP  and the Calderbank-Shor-Steane (CSS) quantum error correcting codes CSS . Recently in tki03 , the unconditional security of the B92 protocol has been proved by relating the B92 protocol to an entanglement distillation protocol (EDP) EPP  initiated by local filtering Gisin96 . The assumptions of the proof are that Alice and Bob’s devises are perfect and the quantum channel is loss-free. To study the unconditional security of B92 protocol over a lossy and noisy channel is interesting from the point of view of the quantum information theory, since the channel loss directly affects the security of the B92 protocol. The proof based on the equivalence between the security and EDPs gives us an idea how to distill the EPR state from states that arise from the lossy and noisy channel.

In this paper, we study the unconditional security of the B92 protocol. We assume that Alice has a perfect single photon source, Bob’s polarization measurement is perfect and he has a perfect single photon counter that discriminates between single photon states on one hand and vacuum state or multi-photon states on the other hand. Under these assumptions, we show the unconditional security of the B92 protocol. The proof follows the idea in tki03 . This paper is organized as follows. In Sec. II, we propose the protocol based on EDP with local filtering and CSS codes sp00 , and we reduce this to the B92 protocol. In Sec. III we give the formula for the error estimations for EDP based on CSS code. Finally, in Sec. IV we give the examples of the resulting performance, and the summary and discussion follow.

## II Protocol based on EDP and its reduction to the B92 protocol

In this section, we introduce the protocol that will turn out to be secure and will reduce to the B92 protocol. In this protocol, Alice first prepares a two-qubit nonmaximally entangled states that can be written as

$`|\mathrm{\Psi }_{\mathrm{AB}}={\displaystyle \frac{1}{\sqrt{2}}}\left(|0_z_\mathrm{A}|\phi _0_\mathrm{B}+|1_z_\mathrm{A}|\phi _1_\mathrm{B}\right)`$ (1)
and sends system B to Bob. Here $`|\phi _j\beta |0_x+(1)^j\alpha |1_x`$, ($`j=0,1`$), where $`0<\alpha <1/\sqrt{2}`$, $`\beta \sqrt{1\alpha ^2}`$, and $`\{|0_x,|1_x\}`$ is a basis ($`X`$-basis) of the qubit. $`X`$-basis and $`Z`$-basis is related through the relationship as $`|j_z=[|0_x+(1)^j|1_x]/\sqrt{2}`$.

After receiving the state, Bob first performs the QND measurements that can be described by the POVM (positive operator valued measure) peres

$`Q_\mathrm{s}`$ $`=`$ $`\mathrm{\Pi }_\mathrm{s}`$ (2)
$`Q_\mathrm{v}`$ $`=`$ $`𝟙\mathbb{\Pi }_\mathrm{s}`$ (3)
where $`\mathrm{\Pi }_\mathrm{s}`$ is the projector onto the single photon space and “v” means the detection of vacuum state or multi-photon states. Thanks to this QND measurements, even if Eve sends multi-photon states or vacuum state to Bob, Alice and Bob can work only on the qubit pairs, provided that the outcome is “s”. We also assume that, after this measurement, if the outcome is “s”, Bob performs the “local filtering operation” Gisin96  on qubit B, described by the Kraus operator

$`A_{\mathrm{fil}}{\displaystyle \frac{\gamma }{\beta }}\left(\alpha |0_x_\mathrm{B}0_x|+\beta |1_x_\mathrm{B}1_x|\right),`$ (4)
where $`\gamma (0,1]`$. Note that if the channel is loss-free and noiseless or Eve does nothing, this operation yields the maximally entangled state (EPR state) $`(|0_x_\mathrm{A}|0_x_\mathrm{B}+|1_x_\mathrm{A}|1_x_\mathrm{B})/\sqrt{2}`$ with probability $`2\alpha ^2\gamma ^2`$ due to the fact that the initial state is also written as $`|\mathrm{\Psi }_{\mathrm{AB}}=\beta |0_x_\mathrm{A}|0_x_\mathrm{B}+\alpha |1_x_\mathrm{A}|1_x_\mathrm{B}`$. Thus, the filtering operation with $`\gamma =1`$ is optimum one in this case. When channel is lossy and noisy, or Eve is present, the filtered states may include bit errors, represented by the subspace spanned by $`\{|0_z_\mathrm{A}|1_z_\mathrm{B},|1_z_\mathrm{A}|0_z_\mathrm{B}\}`$, and a phase error, represented by the subspace spanned by $`\{|0_x_\mathrm{A}|1_x_\mathrm{B},|1_x_\mathrm{A}|0_x_\mathrm{B}\}`$. In order to distill the EPR pairs, Alice and Bob run the EDP based on CSS code sp00  and after this EDP Alice and Bob measure each in the $`Z`$-basis of the EPR pairs so that they obtain a secure key. Note that the QND measurement followed by the filtering $`A_{\mathrm{fil}}`$ and $`Z`$-basis measurement is equivalent to a measurement $`_{\mathrm{B92}}^{L,\gamma }`$ described by $`F_0^\gamma =\gamma ^2|\overline{\phi _1}\overline{\phi _1}|/(2\beta ^2)`$, $`F_1^\gamma =\gamma ^2|\overline{\phi _0}\overline{\phi _0}|/(2\beta ^2)`$, $`F_{\mathrm{null}}^\gamma =\mathrm{\Pi }_\mathrm{s}F_0^\gamma F_1^\gamma `$, and $`F_\mathrm{v}=𝟙𝔽_\mathrm{𝟘}^\gamma 𝔽_\mathrm{𝟙}^\gamma 𝔽_{\mathrm{null}}^\gamma `$, where $`|\overline{\phi _j}\alpha |0_x(1)^j\beta |1_x`$. Thus $`_{\mathrm{B92}}^{L,\beta }`$ is the same measurement as in the B92 with perfect photon detector, and $`_{\mathrm{B92}}^{L,\gamma }`$ corresponds to the generalized B92 protocol. For convenience, let $`_{\mathrm{B92}}^\gamma `$ be the measurements described by the set of the POVM elements $`\{F_0^\gamma ,F_1^\gamma ,F_{\mathrm{null}}^\gamma \}`$. The following is the protocol that we want to be shown secure and reduce to the generalized B92 protocol.

Protocol: (1) Alice prepares the state $`|\mathrm{\Psi }_{\mathrm{AB}}`$, and she sends the system B to Bob over a quantum channel. (2) Bob performs the QND measurement described by $`\{Q_\mathrm{s},Q_\mathrm{v}\}`$, and he publicly announces the outcome. They take notes of the outcome. (3) Alice and Bob repeats the step (1) and (2) $`2N`$ times. Let $`2LN`$ be the number of outcomes that resulted in “v”. (4) Alice and Bob discard all the pairs whose outcomes of the QND measurement are “v”. By public discussion, Alice and Bob randomly permute the positions of the remaining $`2(1L)N`$ pairs. (5) For the first $`(1L)N`$ pairs (check pairs), Alice measures system A in the $`Z`$-basis, and Bob performs measurement $`_{\mathrm{B92}}^\gamma `$ on his system. By public discussion, they determine the number $`n_{\mathrm{err}}`$ of errors in which Alice found $`0`$ and Bob’s outcome was $`1`$, or Alice found $`1`$ with Bob’s outcome $`0`$. (6) For the second $`(1L)N`$ pairs (data pairs), Bob performs the filtering $`A_{\mathrm{fil}}`$ on each of his qubits, and announces the total number $`n_{\mathrm{fil}}`$ and the positions of the qubits that have passed the filtering. (7) From $`L`$, $`n_{\mathrm{err}}`$ and $`n_{\mathrm{fil}}`$, they estimate an upper bound for the number of bit errors $`n_{\mathrm{bit}}`$, and an upper bound for the number of phase errors $`n_{\mathrm{ph}}`$, in the $`n_{\mathrm{fil}}`$ pairs. If these bounds are too large, they abort the protocol. (8) They run an EDP that can produce $`n_{\mathrm{key}}`$ nearly perfect EPR pairs if the estimation is correct. (9) Alice and Bob each measures the EPR pairs in $`Z`$-basis to obtain a shared secret key.

In the following, we will reduce our protocol to the B92 protocol. Before this reduction, we should briefly mention what Alice and Bob do experimentally in the B92 protocol. In the B92 protocol, Alice randomly prepares the single photon polarization in the state $`|\phi _0`$ or $`|\phi _1`$ and sends the state to Bob. After receiving the state, Bob performs measurement described by $`_{\mathrm{B92}}^{L,\gamma }`$. When $`\gamma =\beta `$ and Bob detects a single photon state, this measurement effectively performs a polarization measurement on the single photon state in a basis randomly chosen from $`\{|\phi _0,|\overline{\phi }_0\}`$ and $`\{|\phi _1,|\overline{\phi }_1\}`$. In the case where $`\gamma =1`$ and Bob detects a single photon state, this measurement effectively performs a general measurement which describes the optimum unambiguous state discrimination between $`|\phi _0`$ and $`|\phi _1`$ unb .

Thanks to Bob’s perfect detector, if the outcome of the QND measurement is “s”, then Bob’s states are projected onto the single photon space. This projection requires two more steps in our protocol compared to the Protocol 1 in tki03 . The benefit of this projection is that for the states projected onto qubit-pairs Alice and Bob can use the one-way EDP based on CSS code sp00 . In the context of QKD, this EDP allows Alice to perform $`Z`$-basis measurements immediately after she has prepared the state $`|\mathrm{\Psi }_{\mathrm{AB}}`$ and Bob to perform $`Z`$-basis measurements immediately after he has performed the filtering. Now Bob’s measurement is described by $`_{\mathrm{B92}}^\gamma `$. Since Bob performs the QND measurement before the measurement $`_{\mathrm{B92}}^\gamma `$, his measurement is equivalent to the measurement $`_{\mathrm{B92}}^{L,\gamma }`$. As a whole, the protocol has been reduced into the protocol where Alice prepares and sends $`|\phi _0`$ or $`|\phi _1`$ randomly and Bob performs the $`_{\mathrm{B92}}^{L,\gamma }`$ measurement. This completes the reduction of our protocol into the generalized B92 protocol over lossy and noisy quantum channel.

Due to this equivalence, it is sufficient to prove the security of the protocol for the security analysis of the B92. The most important part in the quantitative security analysis of our protocol is to estimate phase errors and bit errors in the step (7). If this estimation is not exponentially reliable, then in the protocol Alice and Bob cannot distill the state sufficiently close to EPR state so that the protocol is not secure lc98 ; sp00 . In the next section, we will give the formula to estimate the bit errors and phase errors.

## III Phase error estimation from the bit error

In this section, we consider the error estimations in the step (7) of the protocol. Thanks to Bob’s perfect detector, the filtered states are qubit states, which allows us to directly apply some inequalities from tki03 .

First, note that bit errors $`n_{\mathrm{bit}}`$ in the data pairs could be determined if Alice measures those pairs in $`Z`$-basis, Bob performs measurement $`_{\mathrm{B92}}^\gamma `$ onto those pairs, and they compare their results. Thus, the bit error rate $`n_{\mathrm{err}}`$ in the check pairs and bit errors $`n_{\mathrm{bit}}`$ in the data pairs are determined by the same measurements. Furthermore thanks to the random permutation in step (4), we can apply the classical probability estimation so that we obtain

$$|n_{\mathrm{bit}}n_{\mathrm{err}}|(1L)Nϵ_1.$$
(5)
For any strategy by Eve, the probability of violating this inequality is asymptotically less than $`\mathrm{exp}\left[(1L)Nϵ_1^2\right]`$.

Since the POVM elements that give the bit error rate on the filtered states $`\mathrm{\Pi }_{\mathrm{bit}}`$ and phase error rate on the filtered states $`\mathrm{\Pi }_{\mathrm{ph}}`$ do not commute POVMelement , the classical probability argument cannot be applied to the estimation of the phase errors from the bit errors. For the estimation, let us perform the gedanken measurements in the $`|i_x_\mathrm{A}|j_x^{}_\mathrm{B}`$ basis $`(i=0,1,j^{}=0,1,\mathrm{v})`$, and let the numbers $`n_{ij^{}}`$ of the pairs found in the state $`|i_x_\mathrm{A}|j_x^{}_\mathrm{B}`$. Here we define $`n_{i\mathrm{v}}`$ as the outcome where Alice obtains the bit value $`i`$ and Bob obtains vacuum or multi-photon states. As far as we consider the measurement outcomes of the POVM elements that are diagonal with respect to $`|i_x_\mathrm{A}|j_x^{}_\mathrm{B}`$ basis, we are allowed to apply the classical argument for the estimation. Thus, we have the following inequalities,

$`\left|\alpha ^2(n_{00}+n_{10})+\beta ^2(n_{01}+n_{11}){\displaystyle \frac{\beta ^2}{\gamma ^2}}n_{\mathrm{fil}}\right|(1L)Nϵ_2`$ (6)
$`\left|\alpha ^2n_{10}+\beta ^2n_{01}{\displaystyle \frac{\beta ^2}{\gamma ^2}}n_{\mathrm{ph}}\right|(1L)Nϵ_3,`$ (7)
which are violated with probability asymptotically less than $`\mathrm{exp}\left[2(1L)Nϵ_2^2\right]`$ and $`\mathrm{exp}\left[2(1L)Nϵ_3^2\right]`$, respectively. Similarly, since neither the noisy channel nor Eve can touch the qubits held by Alice, we have

$$|\alpha ^2N(n_{10}+n_{11}+n_{1\mathrm{v}})|Nϵ_4,$$
(8)
with probability of violation asymptotically less than $`\mathrm{exp}(2Nϵ_4^2)`$.

In order to represent the bit error in the check pairs in terms of gedanken measurements, we introduce the global gedanken measurements in the basis $`|\mathrm{\Gamma }_{ij}`$ $`(i=0,1,j=0,1)`$ and the number of the corresponding measurement outcomes $`m_{ij}`$, where $`|\mathrm{\Gamma }_{ij}(1)^{ij}\beta |i_x_\mathrm{A}|j_x_\mathrm{B}+(1)^{i(j+1)}\alpha |(i+1)_x_\mathrm{A}|(j+1)_x_\mathrm{B}`$. Here the summation is taken modulo 2. Since the POVM element for the bit error is diagonal with respect to the $`|\mathrm{\Gamma }_{ij}`$ basis POVMelement , we have

$`\left|(m_{11}+m_{01})/2{\displaystyle \frac{\beta ^2}{\gamma ^2}}n_{\mathrm{err}}\right|(1L)Nϵ_5,`$ (9)
which is violated with probability asymptotically less than $`\mathrm{exp}(2(1L)Nϵ_5^2)`$.

Since $`\{|\mathrm{\Gamma }_{01},|\mathrm{\Gamma }_{10}\}`$ and $`\{|0_x_\mathrm{A}|1_x_\mathrm{B},|1_x_\mathrm{A}|0_x_\mathrm{B}\}`$ span the same subspace, we can relate $`m_{10}+m_{01}`$ and $`n_{10}+n_{01}`$ by the classical probability estimate as

$$|(m_{10}+m_{01})(n_{10}+n_{01})|(1L)Nϵ_6,$$
(10)
which is violated with probability asymptotically less than $`\mathrm{exp}\left[L(1L)Nϵ_6^2\right]`$. For the estimate of the phase error rate in the data qubit from the bit error rate in the check qubit, we need to relate $`m_{01}/(m_{01}+m_{10})`$ to $`n_{01}/(n_{01}+n_{10})`$. This relation cannot be obtained by classical arguments, since $`|\mathrm{\Gamma }_{01}`$ and $`|0_x_\mathrm{A}|1_x_\mathrm{B}`$ are nonorthogonal. In tki03 , the authors consider the qubits that are invariant under the random permutation in the limit of large number. Then they showed that when we consider the measurement outcomes onto such qubits, we are allowed to regard the qubits as if they arose from the independently and identically distributed quantum source. Their argument can be directly applied to our case, since in step (4) of the protocol we have also applied the random permutation. Thus, we can make use of the exponentially reliable inequality derived in tki03  as

$$\mathrm{sin}^2(\theta _l\theta )ϵ_7\mathrm{sin}^2\varphi _l\mathrm{sin}^2(\theta _l+\theta )+ϵ_8$$
(11)
for $`l=0,1`$, where all the angles are defined in $`[0,\pi /2]`$ by the relations $`n_{11}/(n_{11}+n_{00})=\mathrm{sin}^2\theta _0`$, $`n_{01}/(n_{01}+n_{10})=\mathrm{sin}^2\theta _1`$, $`m_{11}/(m_{11}+m_{00})=\mathrm{sin}^2\varphi _0`$, $`m_{01}/(m_{01}+m_{10})=\mathrm{sin}^2\varphi _1`$, and $`\alpha ^2=\mathrm{sin}^2\theta `$.

The Eqs. (5)–(11), together with the relation $`n_{ij^{}}=N`$ and $`m_{ij}=(1L)N`$ can be used to derive the formula for the estimation of the phase error from bit error. In the following we set $`ϵ_i0`$, $`(i=1,2,\mathrm{}8)`$ by taking the limit $`N\mathrm{}`$ so that Eqs. (5)–(10) are now linear equations.

We first consider the special case for the B92 where bit error rate in the filtered qubits is zero while the phase error rate in the qubits is not zero. Eve can accomplish this by using channel loss. From $`n_{\mathrm{bit}}=0`$, i.e., $`m_{11}=m_{01}=0`$, one can easily find that for $`n_{\mathrm{fil}}/N=2\alpha ^2\gamma ^2(1L)`$

$$\mathrm{Max}\left(\frac{n_{\mathrm{ph}}}{n_{\mathrm{fil}}}\right)\{\begin{array}{cc}=\frac{\alpha ^2}{\beta ^2\alpha ^2}\frac{L}{1L}\hfill & \text{for }L\beta ^2\alpha ^2\hfill \\ & \\ 1/2\hfill & \text{for }L>\beta ^2\alpha ^2\hfill \end{array}.$$
(12)
The limiting loss rate $`\beta ^2\alpha ^2`$ coincides with the case where Eve employs optimum unambiguous state discrimination measurement unb . In this attack, if she obtains the conclusive outcome, then she sends the corresponding state to Bob, otherwise she does not.

Next we consider the general case. The relationship can be obtained by solving the equations that we have mentioned just above. After straight forward calculations, we find the implicit constraint

$`{\displaystyle \frac{\beta ^2}{\gamma ^2}}\left|n_{\mathrm{fil}}2n_{\mathrm{err}}\right|N\alpha \beta f(x),`$
where

$`f(x)`$ $``$ $`\sqrt{\left(xL_1\right)^2\left(\mathrm{\Delta }+{\displaystyle \frac{L_1}{\alpha ^2\beta ^2}}\right)^2}`$
$`+`$ $`\sqrt{\left(1xL_0\right)^2\left(\beta ^2\alpha ^2\mathrm{\Delta }+{\displaystyle \frac{L_0}{\alpha ^2\beta ^2}}\right)^2}`$
$`\mathrm{\Delta }`$ $``$ $`\left[\left(\beta ^2/\gamma ^2\right)n_{\mathrm{fil}}/N2\alpha ^2\beta ^2\right]/(\beta ^2\alpha ^2)`$
$`x`$ $``$ $`\left(\beta ^2/\gamma ^2\right)2n_{\mathrm{ph}}/N(\beta ^2\alpha ^2)\mathrm{\Delta }`$
$`L_0`$ $``$ $`\alpha ^2n_{1\mathrm{v}}/N+\beta ^2n_{0\mathrm{v}}/N`$
$`L_1`$ $``$ $`\alpha ^2n_{0\mathrm{v}}/N+\beta ^2n_{1\mathrm{v}}/N`$ (13)
Here, $`xL_1`$ and $`1xL_0`$ are equal to $`(n_{10}+n_{01})/N`$ and $`(n_{00}+n_{11})/N`$ respectively so that they must be positive. Combining these positivity constraints with the positivities inside the square roots, we have $`L_1+|\mathrm{\Delta }+\frac{L_1}{\alpha ^2\beta ^2}|x1L_0|\mathrm{\Delta }\alpha ^2+\beta ^2+\frac{L_0}{\alpha ^2\beta ^2}|`$. We can easily confirm that these equations are the same as equations in tki03  where $`\gamma =\beta `$ and the channel has no loss i.e., $`L_0=L_1=0`$. Eq. (III) gives us the possible parameter set $`\{n_{\mathrm{err}},n_{\mathrm{ph}}\}`$ for fixed $`\alpha `$, $`\gamma `$, $`L_0`$, $`L_1`$, $`L`$, and $`n_{\mathrm{fil}}`$. By taking the optimum ratio between $`L_0`$ and $`L_1`$ with respect to $`n_{\mathrm{ph}}`$ for given $`n_{\mathrm{err}}`$, $`n_{\mathrm{fil}}`$ and $`L`$, we obtain an upper bound $`\overline{n}_{\mathrm{ph}}`$ on the number of phase errors $`n_{\mathrm{ph}}`$, as a function of the observed values $`n_{\mathrm{err}}`$, $`n_{\mathrm{fil}}`$ and $`L`$, and the chosen parameters $`\alpha `$ and $`\gamma `$.

## IV Examples and Conclusion

In this section, at first we give the examples of the phase error estimation and resulting performance of the B92 protocol.

To illustrate the accuracy of the phase error estimation and security performance of the B92 protocol, we assign values to the observable data as they would arise from a depolarizing quantum channel with loss, i.e., the state $`\rho `$ evolves into $`L|VV|+(1L)\left[(1p)\rho +_{i=x,y,z}\sigma _i\rho \sigma _i\right]`$, where $`|V`$ is the vacuum state, $`L`$ is the loss rate of quantum channel, $`p`$ is the depolarizing rate, and $`\sigma _i`$ is the Pauli matrix. In Fig. 1, we show the upper bound of the estimated phase error rate, the actual phase error rate, and bit error rate for the case (a) $`p=0.01,L=0`$, and (b) $`p=0.01,L=0.5`$. These figures are independent of $`\gamma `$, since we plot the phase error rates normalized by $`n_{\mathrm{fil}}`$ so that the pre-factor $`\beta ^2/\gamma ^2`$ is cancelled. In each figure, it is seen that when $`|\phi _0|\phi _1|^2`$ becomes smaller, the estimation of the phase errors becomes poorer. On the other hand, larger values of $`|\phi _0|\phi _1|^2`$ make the signal more vulnerable to noise and loss. Furthermore, due to the loss, the estimation in (b) is poorer than that in (a), especially in the region for small $`|\phi _0|\phi _1|^2`$.

Since Alice and Bob apply random permutation in the protocol, the achievable length of the final key is given by $`n_{\mathrm{key}}=(\gamma ^2/\beta ^2)n_{\mathrm{fil}}[1h(n_{\mathrm{bit}}/n_{\mathrm{fil}})h(\overline{n}_{\mathrm{ph}}/n_{\mathrm{fil}})]`$ tki03 ; mh03  , when $`\overline{n}_{\mathrm{ph}}/n_{\mathrm{fil}}1/2`$. Here $`h(p)p\mathrm{log}_2p(1p)\mathrm{log}_2(1p)`$. The positivity of $`n_{\mathrm{key}}`$ means that Alice and Bob generate a secure key. In Fig. 2 (a), we have numerically calculated the key generation rate $`Gn_{\mathrm{key}}/N`$ in the case of $`\gamma =\beta `$ for the fixed value of $`|\phi _0|\phi _1|^2`$, $`L`$ and $`p`$, and plot $`G`$ optimized over $`|\phi _0|\phi _1|^2`$ as a function of $`p`$. The optimum value of $`|\phi _0|\phi _1|^2`$ is shown in Fig. 2 (b). Since the optimization of the ratio between $`L_0`$ and $`L_1`$ is complicated, we have optimized the ratio numerically. From Fig. 2 (a), it is seen that the B92 protocol as described here is secure up to $`p0.034`$ (in the case of $`L=0`$), $`p0.023`$ (in the case of $`L=0.2`$), and $`p0.012`$ (in the case of $`L=0.5`$). The use of the general B92 measurement, i.e., $`\gamma =1`$ yields the higher secret key gain, however this does not change the cutoff of the gain nor the optimum angle in Fig. 2 (b). This is again because of the fact that the cutoff is determined by the error rate normalized by $`n_{\mathrm{fil}}`$.

In the BB84 protocol or 6-state protocol b98 , thanks to the randomly chosen conjugate bases, phase error estimation is more precise than that in the B92 protocol. The poor estimation in the B92 makes this protocol more vulnerable to eavesdropping. Actually, it is known that the BB84 protocol and 6-state protocol are secure up to $`p0.165`$ sp00  and $`p0.1905`$ l01  corresponding to bit error rates of $`11\%`$ and $`12.7\%`$, respectively.

In summary, we have shown that the security proof of the B92 over loss-free channel tki03  can be adapted to lossy and noisy channels assuming Bob’s detector that discriminates between single photon states on one hand and vacuum state or multi-photon states on the other hand. We have derived the general formula for the estimation of the phase errors from bit errors. We have also shown that the boundary loss rate at which Eve can induce phase errors freely without inducing bit errors coincides with the case where Eve employs the optimum unambiguous state discrimination measurement. Finally, we have shown the examples of the security assuming the quantum channel is the depolarizing quantum channel with loss.

It is an interesting problem to compare our results to the case where Eve performs only individual attack. For the B92 with loss, the security is estimated in tki03-2 . However in that paper, it is assumed that Alice and Bob employ error discarding protocol so that we cannot compare our result to their result directly. We leave this problem for future studies. Another interesting problem is to consider the case where Bob’s detector is not perfect. In this case, some filtered states may remain vacuum state or multi-photon states so that our proof cannot be directly applied. The security analysis of the B92 with coherent state is another interesting problem. We hope our study is helpful for solving these problems and we leave them as future studies.

## V Acknowledgement

K.T appreciates the warmth and hospitality of the Quantum Information Theory Group, Universität Erlangen-Nürnberg in Germany during his stay. We thank to Marcos Curty, Peter van Loock, Philippe Raynal, and Aska Dolinska for helpful discussions. This work has been supported in part by the Academic Frontiers Student Exchange Promotion Program under Grant No. 60 by the Ministry of Education, Culture, Sports, Science and Technology, Japan, and by the DFG under the Emmy-Noether programme.


format: text

text:
# Occurrence and Stability of Apsidal Resonance in Multiple Planetary Systems

## 1 INTRODUCTION

The detecting of extrasolar planetary systems reveals fruitful results during the past years. More than 100 extrasolar planets has been inferred by the Doppler radial velocity measurements to the solar-type stars (California and Carnegie Planet Search, 2003), among them 10 multiple-planet systems are confirmed. For a multiple-planet system, the dynamical stability of the system under planetary interaction is an important issue concerning the dynamical evolution as well as the possible existence of habitable zone of the system.

There are many effects which can affect the stability of a multiple-planet system. For the orbits of planets with small or modest eccentricities and inclinations, mean motion resonances between planets can sometimes lead to stable configurations. Another effect is the secular resonance between the planets. An apsidal resonance occurs when the relative apsidal longitudes of the two orbits $`\mathrm{\Delta }\varpi `$ librates about $`0`$ (aligned resonance) and $`\pi `$ (anti-aligned resonance) during the evolution. Due to the aligned apsidal resonance, the two planets on elliptic orbits can greatly reduce the possibility of close encounters, thus it is believed that aligned apsidal resonance can stabilize the interacting planets.

For the 10 multiple planetary systems observed to date (GJ876, 47UMa, HD82943, HD12661, HD168443, HD37124, HD38529, HD74156, Upsilon Andromedae, and 55 Cancri), the best fit orbital parameters inferred from the radial velocity observations show 6 pairs of planets could be in apsidal resonance: HD82943(Goździewski & Maciejewski, 2001), Upsilon Andromedae c and d(Chiang, Tabachnik, & Tremaine 2001), GJ876(Lee & Peale, 2002), 47UMa(Laughlin, Chambers, & Fischer 2002), HD12661(Goździewski & Maciejewski, 2003; Lee & Peale, 2003), and 55 Cancri b and c(Ji et al., 2003). Ubiquitous as it is, the apsidal resonance phenomenon is worthy to be studied in detail. In this paper, we are interested when the apsidal resonance occurs and whether it really leads to a stable configuration between planets, since an anti-aligned resonance could lead to close encounters between orbits with large eccentricities .

For the occurrence of apsidal resonance, Laughlin et al. (2002) gives a criterion based on the Laplace-Lagrange solution of the secular perturbation system. In this paper, the criterion is represented to a more explicit form in section II. The probability that the apsidal resonance happens is also derived according to the criteria. In section III, we study the stability of orbits in apsidal resonances by calculating the largest Lyapunov exponents of the orbits with a general three-body model. The behavior of orbits in migration are studied with a torqued three-body model. The conclusions and the applications of the criteria to the other multiple planetary systems are discussed in the final section.

## 2 LOCATIONS OF APSIDAL RESONANCE

In this section we derive the explicit criteria under which the apsidal secular resonance may occur. To that aim, the linear secular perturbation theory is employed for two interacting planets under the attraction of the host star. For a planetary system with two planets, hereafter we denote all the quantities of the host star, the inner and outer planets with subscripts “0”, “1” and “2”, respectively. So the three bodies have masses $`m_0,m_1,m_2`$, respectively, where $`m_1,m_2m_0`$. In the present study we address the coplanar problem only, so the two planets are on the orbits with osculating orbital elements $`(a_1,e_1,\varpi _1,M_1)`$ and $`(a_2,e_2,\varpi _2,M_2)`$, respectively, where $`a,e,\varpi ,M`$ are the semi-major axis, eccentricity, longitude of pericenter and mean anomaly of the orbit, respectively. We adopt the commonly used unit system, i.e., the mass unit is the solar mass, the length unit is 1AU, and the time unit is 1yr/($`2\pi `$).

### 2.1 Linear secular perturbation theory revisited

We start with the linear secular theory following Murray and Dermott (1999). For the coplanar case, the disturbing function for the motions of planets $`m_1`$ and $`m_2`$ are given as:

$$\begin{array}{cc}R_1=& n_1a_1^2[\frac{1}{2}A_{11}e_1^2+A_{12}e_1e_2\mathrm{cos}(\varpi _1\varpi _2)]\hfill \\ R_2=& n_2a_2^2[\frac{1}{2}A_{22}e_2^2+A_{21}e_1e_2\mathrm{cos}(\varpi _1\varpi _2)]\hfill \end{array}$$
(1)
where $`n_1,n_2`$ are the mean motion of planets $`m_1`$ and $`m_2`$, respectively; $`A_{ij}`$ are elements of matrix given by

$$\left(\begin{array}{cc}A_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)=\left(\begin{array}{cc}c_1& c_0c_1\\ c_0c_2& c_2\end{array}\right)$$
(2)
where $`c_k>0(k=0,1,2)`$ are functions of $`a_1,a_2,m_0,m_1,m_2`$ defined as

$$\begin{array}{cc}\hfill c_0=& b_{\frac{3}{2}}^{(2)}(\alpha )/b_{\frac{3}{2}}^{(1)}(\alpha )\frac{5}{4}\alpha (1\frac{1}{8}\alpha ^2)\hfill \\ \hfill c_1=& \frac{1}{4}n_1\frac{m_2}{m_0+m_1}\alpha ^2b_{\frac{3}{2}}^{(1)}(\alpha )\hfill \\ \hfill c_2=& \frac{1}{4}n_2\frac{m_1}{m_0+m_2}\alpha b_{\frac{3}{2}}^{(1)}(\alpha )\hfill \end{array}$$
(3)
with $`b_{\frac{3}{2}}^{(i)}(\alpha )(i=1,2)`$ being the Laplace coefficients, and $`\alpha =a_1/a_2<1`$. Fig.1 shows the approximation of $`c_0`$ by the above formula. Quantitatively, the error of approximation is less than $`1\%`$ for $`\alpha <0.66`$, and less than $`5\%`$ for $`\alpha <0.91`$. So the approximation is quite well for the study of the planetary system. Moreover, we define

$$\xi =\frac{c_2}{c_1}=\frac{1}{\alpha }\frac{n_2m_1(m_0+m_1)}{n_1m_2(m_0+m_2)}q\alpha ^{1/2},$$
(4)
with $`q=m_1/m_2`$ and the terms with orders of $`O(\frac{m_1}{m_0}),O(\frac{m_2}{m_0})`$ or higher are neglected in the above approximation, since $`\frac{m_1}{m_0},\frac{m_2}{m_0}10^3`$ in the planetary systems. Denote $`g_1,g_2`$ as the two eigenvalues of matrix (2), and the corresponding eigenvectors are $`S_i\left(\begin{array}{c}\mathrm{cos}\theta _i\\ \mathrm{sin}\theta _i\end{array}\right)`$ , where $`\theta _i(\pi /2,\pi /2)`$ and

$$\mathrm{cos}\theta _i=\frac{c_0c_1}{\sqrt{(c_1g_i)^2+c_0^2c_1^2}},\mathrm{sin}\theta _i=\frac{c_1g_i}{\sqrt{(c_1g_i)^2+c_0^2c_1^2}},(i=1,2)$$
(5)
with

$$\begin{array}{c}g_1=\frac{1}{2}[(c_1+c_2)+\sqrt{(c_1c_2)^2+4c_0^2c_1c_2}],\hfill \\ g_2=\frac{1}{2}[(c_1+c_2)\sqrt{(c_1c_2)^2+4c_0^2c_1c_2}].\hfill \end{array}$$
(6)
Define

$$\begin{array}{cc}\rho _1\mathrm{tan}\theta _1\hfill & =\frac{1}{2c_0}(1\xi \sqrt{(1\xi )^2+4c_0^2\xi }),\hfill \\ \rho _2\mathrm{tan}\theta _2\hfill & =\frac{1}{2c_0}(1\xi +\sqrt{(1\xi )^2+4c_0^2\xi }).\hfill \end{array}$$
(7)
So $`\rho _1<0`$ and $`\rho _2>0`$, or $`\frac{\pi }{2}<\theta _1<0<\theta _2<\frac{\pi }{2}`$. The scaling factor $`S_i`$ (i=1,2) can be expressed in terms of initial eccentricities $`e_{10}`$, $`e_{20}`$ and $`\mathrm{\Delta }\varpi _0=\varpi _{20}\varpi _{10}`$:

$$\begin{array}{c}S_1=\frac{[\rho _2^2e_{10}^22\rho _2e_{10}e_{20}\mathrm{cos}\mathrm{\Delta }\varpi _0+e_{20}^2]^{1/2}}{|\rho _1\rho _2|\mathrm{cos}\theta _1}\frac{F}{|\rho _1\rho _2|\mathrm{cos}\theta _1},\hfill \\ S_2=\frac{[\rho _1^2e_{10}^22\rho _1e_{10}e_{20}\mathrm{cos}\mathrm{\Delta }\varpi _0+e_{20}^2]^{1/2}}{|\rho _1\rho _2|\mathrm{cos}\theta _2}\frac{G}{|\rho _1\rho _2|\mathrm{cos}\theta _2}.\hfill \end{array}$$
(8)
The secular system with disturbing functions (1) is integrable and the solutions can be written as:

$$\begin{array}{cc}e_1& =\frac{1}{|\rho _1\rho _2|}[F^2+2FG\mathrm{cos}\mathrm{\Delta }\psi +G^2]^{1/2}\hfill \\ e_2& =\frac{1}{|\rho _1\rho _2|}[\rho _1^2F^2+2\rho _1\rho _2FG\mathrm{cos}\mathrm{\Delta }\psi +\rho _2^2G^2]^{1/2},\hfill \\ e_1e_2\mathrm{sin}\mathrm{\Delta }\varpi & =\frac{1}{\rho _1\rho _2}FG\mathrm{sin}\mathrm{\Delta }\psi \hfill \\ e_1e_2\mathrm{cos}\mathrm{\Delta }\varpi & =\frac{1}{(\rho _1\rho _2)^2}[\rho _1F^2+(\rho _1+\rho _2)FG\mathrm{cos}\mathrm{\Delta }\psi +\rho _2G^2].\hfill \end{array}$$
(9)
where $`\mathrm{\Delta }\psi =\psi _2\psi _1=(g_2t+\beta _2)(g_1t+\beta _1)`$, with $`t`$ the time and $`\beta _1,\beta _2`$ given by:

$$\begin{array}{cc}\mathrm{sin}\beta _1=\frac{1}{F}(h_{10}\rho _2h_{20}),\hfill & \mathrm{sin}\beta _2=\frac{1}{G}(h_{10}\rho _1h_{20}),\hfill \\ \mathrm{cos}\beta _1=\frac{1}{F}(k_{10}\rho _2k_{20}),\hfill & \mathrm{cos}\beta _2=\frac{1}{G}(k_{10}\rho _1k_{20}).\hfill \end{array}$$
(10)
where $`h_{i0}=e_{i0}\mathrm{sin}\varpi _{i0}`$, $`k_{i0}=e_{i0}\mathrm{cos}\varpi _{i0}`$, $`(i=1,2)`$. From (9), it’s easy to verify that the evolution of $`e_1,e_2`$ obeys an integral:

$$\frac{e_1^2}{A_{12}}+\frac{e_2^2}{A_{21}}=D,$$
(11)
where $`D`$ is a constant which depends only on the initial parameters. Moreover, from (8) and (9), the maximum of $`e_1`$ and minimum of $`e_2`$ occur at $`\mathrm{cos}\mathrm{\Delta }\psi =1`$ (as $`\rho _1<0`$), with values:

$$e_{1max}=\frac{F+G}{|\rho _1\rho _2|},e_{2min}=\frac{|\rho _1F+\rho _2G|}{|\rho _1\rho _2|}.$$
(12)
The minimum of $`e_1`$ and maximum of $`e_2`$ occur at $`\mathrm{cos}\mathrm{\Delta }\psi =1`$, with

$$e_{1min}=\frac{|FG|}{|\rho _1\rho _2|},e_{2max}=\frac{\rho _2G\rho _1F}{|\rho _1\rho _2|}.$$
(13)
Thus, we can obtain the maximum excursions of $`e_1`$ and $`e_2`$ for any given $`e_{10},e_{20},\mathrm{\Delta }\varpi _0`$ as follows:

$$\mathrm{\Delta }e_1=\frac{(F+G)|FG|}{|\rho _1\rho _2|},\mathrm{\Delta }e_2=\frac{(\rho _2G\rho _1F)|\rho _1F+\rho _2G|}{|\rho _1\rho _2|}$$
(14)
### 2.2 The explicit criteria

With the help of the last equation of (9), the criterion in Laughlin et al. (2002) for the apsidal resonance can be expressed as

$$S=\left|\frac{(\rho _1+\rho _2)FG}{\rho _1F^2+\rho _2G^2}\right|<1$$
(15)
Since when $`S<1`$, the values of $`\mathrm{\Delta }\varpi `$ can not reach $`\pi /2`$ or $`3\pi /2`$ (thus $`\mathrm{cos}\mathrm{\Delta }\varpi 0`$), so it must librate about $`0`$ or $`\pi `$. On the contrary, when $`S>1`$, it is possible that $`\mathrm{\Delta }\varpi `$ will reach $`\pi /2`$ or $`3\pi /2`$ , thus it will circulate in $`[0,2\pi ]`$.

The equation (15) equivalent to, after some algebra manipulations:

$$\frac{F}{G}>\mathrm{max}(1,\rho _1/\rho _2),\mathrm{or},0<\frac{F}{G}<\mathrm{min}(1,\rho _1/\rho _2),$$
(16)
In view of (8), the above relations are equivalent to :

$$\frac{e_{20}}{e_{10}}<\frac{2\rho _1\rho _2}{\rho _1+\rho _2}\mathrm{cos}\mathrm{\Delta }\varpi _0,$$
(17)
or,

$$\frac{e_{20}}{e_{10}}>\frac{\rho _1+\rho _2}{2}\frac{1}{\mathrm{cos}\mathrm{\Delta }\varpi _0}>0.$$
(18)
By substituting (3) (4) and (7) into the above expressions, we finally obtain:

$$\frac{e_{20}}{e_{10}}<\frac{5}{2}\frac{q\alpha ^{3/2}(1\frac{1}{8}\alpha ^2)}{1q\alpha ^{1/2}}\mathrm{cos}\mathrm{\Delta }\varpi _0,$$
(19)
or,

$$\frac{e_{20}}{e_{10}}>\frac{2}{5}\frac{1q\alpha ^{1/2}}{\alpha (1\frac{1}{8}\alpha ^2)}\frac{1}{\mathrm{cos}\mathrm{\Delta }\varpi _0}>0.$$
(20)
These are the explicit criteria for the occurrence of the apsidal secular resonance. Equations (17) and (18) are obtained with the linear secular perturbation theory, while to get (19) and (20), we use the approximations of $`c_0`$ and $`\xi `$ in (3) and (4).

We call the libration region defined in (19) the down-libration region, and that defined in (20) the up-libration region. Whether the down-libration or up-libration is the aligned or anti-aligned resonance depends on the sign of $`(1q\alpha ^{1/2})`$. For $`q\alpha ^{1/2}<1`$, down-libration occurs only when $`\pi /2<\mathrm{\Delta }\varpi _0<3\pi /2`$ (thus it is the anti-aligned resonance) and up-libration occurs when $`\mathrm{\Delta }\varpi _0<\pi /2`$, $`\mathrm{\Delta }\varpi _0>3\pi /2`$ (the aligned resonance). This is the case of HD12661 system. If $`q\alpha ^{1/2}>1`$, the conclusions are reversed, which is the case of 47UMa system. For the critical case $`q\alpha ^{1/2}=1`$, according to (19) and (20), all the orbits are in libration except those with $`\mathrm{\Delta }\varpi _0=\pi /2,3\pi /2`$. Fig.2 shows a typical phase space and the evolution of an orbit in the linear secular perturbation system (1). The jumps of $`\mathrm{\Delta }\varpi `$ are due to the cross of origin in the $`(e_i\mathrm{cos}\varpi _i,e_i\mathrm{sin}\varpi _i)`$ plane. Fig.3 shows the resonance region in the $`e_{20}\mathrm{\Delta }\varpi _0`$ plane defined by (19)(20) with different $`e_{10}`$. The parameters $`\alpha `$ and $`q`$ are taken from the planetary systems HD12661 and 47 UMa(listed in Table 1-2). The boxes around the present configuration dots show the uncertainties of the elements (also listed in Table 1-2).

There are two limit cases for the criteria (19) and (20):

(1) $`\alpha 0`$. The minimum $`e_{20}`$ for up-libration and maximum $`e_{20}`$ for down-libration can be obtained by setting $`\mathrm{cos}\mathrm{\Delta }\varpi =1`$ in (19) and (20). When $`\alpha 0`$, the minimum $`e_{20}`$ for up-libration tends to very large and maximum $`e_{20}`$ for down-libration tends to zero. So when the two planets are far away, the libration regions in $`e_{20}\mathrm{\Delta }\varpi _0`$ plane can be negligible small.

(2) $`e_{i0}0`$. This happens when one of the planets is in a near-circular orbit, and it is just the case discussed in Malhotra (2002), namely for two planets initially in nearly circular orbits, an impulse perturbation may impart a finite eccentricity to one planet’s orbit. When $`e_{10}0`$, criterion (20) is always fulfilled if $`\pi /2<\mathrm{\Delta }\varpi _0<\pi /2`$ for $`q\alpha ^{1/2}<1`$ or $`\pi /2<\mathrm{\Delta }\varpi _0<3\pi /2`$ for $`q\alpha ^{1/2}>1`$. The boundary curves are the limits of those with $`e_{10}`$ tends to zero in Fig.3a or Fig.3b, with both the minimum $`e_{20}`$ for up-libration and the maximum $`e_{20}`$ for down-libration tend to zero. Thus half of the $`e_{20}\mathrm{\Delta }\varpi _0`$ plane is the possible resonance region with up-libration. So the probability of these two planets captured into apsidal resonance tends to $`50\%`$. This conclusion agrees with that of Malhotra (2002). It’s possible that these two planets are captured in anti-aligned apsidal resonance, which depends on the sign of $`1q\alpha ^{1/2}`$. Similarly, when $`e_{20}0`$, criterion (19) is always fulfilled, thus half of $`e_{20}\mathrm{\Delta }\varpi _0`$ plane are possible resonance region with down-libration, and the probability of the two planets in apsidal resonance (either aligned or anti-aligned resonance) tends to $`50\%`$.

To compare the above libration regions obtained by the criteria with those calculated from the secular perturbation system, we integrate the orbits of the secular perturbation system for the HD12661 system. Fig.4 shows the diagrams of orbits in $`e_2\mathrm{\Delta }\varpi `$ plane, the initial values of the studied orbits have $`e_{10}=0.1`$, with $`\mathrm{\Delta }\varpi _0=0`$ for Fig.4a and $`\mathrm{\Delta }\varpi _0=\pi `$ for Fig.4b. According to the criteria(19) and (20), $`e_{20}>0.022`$ for the aligned resonance at $`\mathrm{\Delta }\varpi _0=0`$, and $`e_{20}<0.375`$ for the anti-aligned resonance at $`\mathrm{\Delta }\varpi _0=\pi `$, which coincide with those given from the secular perturbation system.

Fig.5 shows the contours of the excursions of $`e_1`$,$`e_2`$ for different $`e_{20}`$ and $`\mathrm{\Delta }\varpi _0`$. As we can see, generally the orbits in apsidal resonance have relatively smaller $`\mathrm{\Delta }e_1`$, $`\mathrm{\Delta }e_2`$, especially the turning points of the contour curves lie in one of libration boundary curves. Thus from the linear secular perturbation theory, the orbits in apsidal resonance, either in aligned resonance or anti-aligned resonance, are more stable than those in non-resonance region.

### 2.3 Area of the libration region

From the above criteria, we can calculate the probability that the two planets fall in apsidal resonance in the space of initial orbital elements . We define the probability as the area of the libration region in the $`e_{20}\mathrm{\Delta }\varpi _0`$ plane for a given $`e_{10}`$. For the down-libration case, according to (19), it’s possible that the peak value of the $`e_{20}\mathrm{\Delta }\varpi _0`$ curve can be above unit for larger $`e_{10}`$ (as the $`e_{10}=0.35`$ case in Fig.3a, and the $`e_{10}=0.50`$ case in Fig.3b). We set $`\mathrm{\Delta }\varpi _d`$ as the half width of the down-libration region where the boundary curve reaches $`e_{20}=1`$. Fig.3a and Fig.3b show $`\mathrm{\Delta }\varpi _d`$ for the $`e_{10}=0.35`$ and $`e_{10}=0.50`$ curves, respectively. Define

$$Q_d=\frac{5}{2}\frac{e_{10}q\alpha ^{3/2}(1\frac{1}{8}\alpha ^2)}{|1q\alpha ^{1/2}|},$$
(21)
then,

$$\mathrm{\Delta }\varpi _d=\{\begin{array}{cc}\mathrm{arccos}(1/Q_d),\hfill & \mathrm{if}Q_d>1\hfill \\ 0,\hfill & \mathrm{if}Q_d1\hfill \end{array}$$
(22)
and the area ratio of the down-libration area to the total area of $`e_{20}\mathrm{\Delta }\varpi _0`$ plane is,

$$\begin{array}{cc}P_d\hfill & =\frac{1}{\pi }[\mathrm{\Delta }\varpi _d+_{\mathrm{\Delta }\varpi _a}^{\mathrm{\Delta }\varpi _a+\frac{\pi }{2}\mathrm{\Delta }\varpi _d}Q_d\mathrm{cos}\mathrm{\Delta }\varpi d\mathrm{\Delta }\varpi ]\hfill \\ & =\frac{1}{\pi }[\mathrm{\Delta }\varpi _d+Q_d(1\mathrm{sin}\mathrm{\Delta }\varpi _d)]\hfill \end{array}$$
(23)
where the lower integration limit is the beginning point of the down-libration region in $`\mathrm{\Delta }\varpi _0`$-axis ($`\mathrm{\Delta }\varpi _a=\frac{\pi }{2}`$ for $`q\alpha ^{1/2}<1`$ and $`\mathrm{\Delta }\varpi _a=\frac{3\pi }{2}`$ for $`q\alpha ^{1/2}>1`$). Fig.3a and Fig.3b show the $`\mathrm{\Delta }\varpi _a`$ for the $`e_{10}=0.35`$ and $`e_{10}=0.50`$ curves, respectively. As one can see, the area ratio for the down-libration increases linearly with $`e_{10}`$ when $`e_{10}`$ is small, since in the interval one has $`\mathrm{\Delta }\varpi _d=0`$ in equation(23). However, for larger $`e_{10}`$, the increase of area ratio is no longer linear since $`\mathrm{\Delta }\varpi _d0`$.

Similarly, for the up-libration resonance, if we set $`\mathrm{\Delta }\varpi _u`$ as the half wide of the up-libration region when the boundary curve meets $`e_{20}=1`$(see Fig3) and define

$$Q_u=\frac{2}{5}\frac{e_{10}|1q\alpha ^{1/2}|}{\alpha (1\frac{1}{8}\alpha ^2)},$$
(24)
so,

$$\mathrm{\Delta }\varpi _u=\mathrm{arccos}(Q_u),$$
(25)
and the area ratio of the up-libration area to the total area of $`e_{20}\mathrm{\Delta }\varpi _0`$ plane is,

$$\begin{array}{cc}P_u\hfill & =\frac{1}{\pi }[\mathrm{\Delta }\varpi _u_{\mathrm{\Delta }\varpi _b}^{\mathrm{\Delta }\varpi _b+\mathrm{\Delta }\varpi _u}Q_u\frac{1}{\mathrm{cos}\mathrm{\Delta }\varpi }𝑑\mathrm{\Delta }\varpi ]\hfill \\ & =\frac{1}{\pi }[\mathrm{\Delta }\varpi _uQ_u\mathrm{ln}(\frac{1+\mathrm{sin}\mathrm{\Delta }\varpi _u}{\mathrm{cos}\mathrm{\Delta }\varpi _u})]\hfill \end{array}$$
(26)
where the lower integration limit $`\mathrm{\Delta }\varpi _b`$ is the center of the up-libration region ($`\mathrm{\Delta }\varpi _b=0`$ for $`q\alpha ^{1/2}<1`$ and $`\mathrm{\Delta }\varpi _b=\pi `$ for $`q\alpha ^{1/2}>1`$, see Fig.3).

In the early evolution of planetary systems, both $`q`$ and $`\alpha `$ may vary due to planetary formation and migration. We fix $`e_{10}=0.35`$ for the HD12661 systems, and see the variation of libration area ratio with $`\alpha `$ or $`q`$. Fig.6 shows the variation of libration area ratios with $`\alpha `$ and $`q`$. In Fig.6a $`q1.46`$ is fixed as the observed value. One can see the ratios increase with $`\alpha `$ before they reach the maximum (unit) at $`\alpha _{max}=q^20.47`$, which is the critical case, and then decrease. In Fig.6b $`\alpha 0.32`$ is fixed as the observed values. Again the curves reach the maximum at $`q=\alpha ^{1/2}1.76`$, the critical case.

## 3 STABILITY OR ORBITS IN RESONANCE

Since the linear secular perturbation theory is an approximation to the real three-body system, the above criteria obtained from the linear perturbation theory has its limitation. To apply the linear criteria to the predicting of the apsidal secular resonance, we integrate the orbits in a general three-body (co-planar) system, where the longitudes of the ascending nodes and the inclinations of the two planet orbits are assumed to be zero ($`\mathrm{\Omega }_1=\mathrm{\Omega }_2=0`$, $`i_1=i_2=0`$) in the paper. We adopt the RKF7(8) (Runge-Kutta-Fehlberg) integrator with adaptive step-sizes to integrate the orbits. Generally the step is alternated between $`0.00625`$yr and $`0.0125`$yr, so there are 80-160 steps in a period of planet orbit with a semi-major axis of 1AU, and the final error of the Hamiltonian of the three-body system after $`10M`$ years’ evolution is less than $`10^9`$.

We define an index to indicate whether an orbit is in libration region or not. Choose a serial of discrete time during the evolution of orbits (for example, every 12.5 years), give an index $`I_n`$ for each time so that $`I_n=0`$ if at that time $`\pi /2<\mathrm{\Delta }\varpi <\pi /2`$ and $`I_n=1`$ if $`\pi /2<\mathrm{\Delta }\varpi <3\pi /2`$. Then the average values of $`I_n`$ over very large $`n`$, denoted by $`<I_n>`$, shows roughly the character of the orbit during the studied period of time according to

$$\mathrm{Index}=<I_n>\{\begin{array}{cc}0& \mathrm{aligned}\mathrm{libration}\\ 0.5& \mathrm{circulation}\\ 1& \mathrm{anti}\mathrm{aligned}\mathrm{libration}\\ \mathrm{others}& \mathrm{mixed}\end{array}$$
(27)
Fig.7 shows the index for the HD12661 planet system for $`10\times 51`$ orbits in the interval $`[0,0.5]\times [0,\pi ]`$ of the $`e_{20}\mathrm{\Delta }\varpi _0`$ plane(the same for the following calculations), the initial eccentricity is $`e_{10}=0.10`$ in Fig.7a and $`e_{10}=0.35`$ in Fig.7b. As one can see, most of the orbits in the libration region predicted by linear criteria (19) and (20) are in real libration region when $`e_{20}`$ is small. The discrepancies between the linear system and the three-body one mainly occur on larger $`e_{20}`$ and the boundary between the libration and circulation region.

Next we want to see whether the orbits, either in apsidal resonance or not, may have different stability in the general co-planar three-body model. The stability of an orbit in a Hamiltonian system is related with the topology (regular or chaotic) of the phase space, so we calculate the largest Lyapunov Characteristic Exponent (LCE) to indicate whether the corresponding orbit is in a regular or chaotic region. The LCE at finite time $`\chi (t)`$ is calculated for few orbits up to $`t=10`$Myr (denoted as $`\chi _7`$), and for most orbits up to $`t=1`$Myr (denoted by $`\chi _6`$). Fig.8 shows the LCEs of four orbits in the four different kinds of region in the phase space. For curve (a), $`\chi (t)`$ decrease linearly with $`t`$, thus the orbit corresponding to curve (a) has zero LCE, and is in a regular region. Curve (b) shows a very small but non-zero LCE, so the orbit corresponding to (b) is in a very weak chaotic region. Both curves (a) and (b) have $`\chi _610^5`$yr<sup>-1</sup>, $`\chi _710^6`$yr<sup>-1</sup>. Curve (c) tends to a constant value, with $`\chi _6\chi _710^4`$yr<sup>-1</sup>, so the orbit corresponding to (c) should be in a strong chaotic region. The orbit corresponding to curve (d) is unstable with the outer planet escape before $`10^7`$Myr, and $`\chi (t)`$ for such an orbit is generally great than $`10^2`$yr<sup>-1</sup> before escape. Thus by calculating $`\chi _6`$ (unit: yr<sup>-1</sup>), we can tell at least three different kinds of orbits:

$$\chi _6\{\begin{array}{cc}(10^5,10^4]& \mathrm{regular}\mathrm{or}\mathrm{weak}\mathrm{chaos}\\ (10^4,10^2]& \mathrm{strong}\mathrm{chaos}\\ (10^2,10^1]& \mathrm{unstable}\end{array}$$
(28)
We calculate $`\chi _6`$ for the HD12661 system with initial eccentricity $`e_{10}=0.10`$ in one run and $`e_{10}=0.35`$ in another run, and the other initial parameters are taken as the observed values. Fig.9 show the results. The boundary curves for the corresponding $`e_{10}`$ are also plotted in the diagram. We find for small $`e_{10}`$, the orbits, whether in aligned resonance, anti-aligned resonance or non-resonance regions, do not show much difference about the LCEs. So in this case, whether an orbit is in apsidal resonance do not affect much on its stability. However, for larger initial $`e_{10}`$, orbits in aligned resonance seem to be more stable since they have much lower LCE as compared with those in the anti-aligned resonance or circulation regions with same $`e_{20}`$. This example shows for larger $`e_{10},e_{20}`$, planets in aligned resonance regions would be relatively more stable. From Fig.9b, we can also see the present configuration of HD12661b and HD12661c is in the boundary of a chaotic region, if $`\mathrm{\Delta }\varpi _0`$ is set to $`130.7^0`$, which is symmetric with the observed value $`130.7^0`$. This conclusion has been obtained by Kiseleva-Eggleton et al. (2002).

Finally, we study the role of apsidal resonance on the stability of the orbits when the planets are migrating. In the early stage of planet evolution, the protoplanets and the stellar-disk might be coexisting and interacting, thus planet migration might happen due to nebular tides(see, e.g., Ward 1997). We adopt the torqued three-body model as in Laughlin et al. (2002), and for the sake of simplicity, we consider the case that only the outer planet experiences an azimuthal torque due to the planet-disk interaction. We take the azimuthal acceleration as $`f_2=2\times 10^6`$AU<sup>2</sup>yr<sup>-1</sup> (which is smaller than that used in Laughlin et al. (2002), because here we concern the qualitative evolutions only), and study both the forward ($`t>0`$) and backward ($`t<0`$) evolutions of orbits under this acceleration. We calculate the orbits for the HD12661 system with initial $`e_{10}=0.35`$, and $`M_{10},M_{20}`$ are random chosen. All the other initial parameters of orbits are taken from table 1. The evolution time span is $`50000`$ years. We find in all the studies cases(both forward and backward), the semi-major axis of the inner planet do not have secular changes, just as the results in Laughlin et al. (2002). Fig.10a shows the index for the orbits during the evolution and Fig.10b shows the $`a_2`$ at the final time for the forward case. Some of the orbits which were initially in apsidal resonance region become mixed, though, according to Fig.6a, the libration regions are enlarged due to the increase of $`\alpha `$ (from 0.32 to approximately 0.41) by planet migration. For orbits in the aligned resonance, they are more stable and with smaller final $`a_2`$, while for orbits in anti-aligned resonance with larger $`e_{20}`$, or in circulation region, they tend to be more unstable due to close encounters and thus have larger final $`a_2`$, and most of them will escape soon in the following evolutions. Fig.11 shows the backward case. We find the conclusions are more or less similar. Though the libration regions shrink due to the decrease of $`\alpha `$ by the planet migration (from 0.32 to approximately 0.23) in this case, the configurations of apsidal resonance are generally preserved during the migration. Thus most of planet systems which are observed in apsidal resonance now may be in resonance before the migration begins. Again, the orbits in aligned resonance seem to be more stable in the sense of having modest final $`a_2`$, which may be the results of less close encounters between the two planets.

## 4 CONCLUSIONS AND DISCUSSIONS

In this paper,we have studied the occurrence as well as the stability of the apsidal resonance. The apsidal resonance occurs when the equations (19)(20) are fulfilled. We find the occurrence of apsidal resonance depends only on the mass ratio $`q=m_1/m_2`$, semi-major axis ratio $`\alpha =a_1/a_2`$(in secular systems, $`a_1,a_2`$ are constants), initial eccentricity ratio $`e_{10}/e_{20}`$ and the relative apsidal longitude $`\varpi _{20}\varpi _{10}`$ of the two planets. The criteria are based on the Laplace-Lagrange secular solution of linear perturbation theory. Based on these criteria, the ratio of librating to non-librating orbits in the $`e_{20}\mathrm{\Delta }\varpi _0`$ plane can be obtained analytically, which is given in (23)(26). We also find for two planets on the orbits with large eccentricities, they can be in a stable configuration only when they are in aligned apsidal resonance. When the planets are migrated under the planet-disk interactions, more than half of the studied orbits preserves the configurations of apsidal resonance.

The linear secular perturbation theory is applicable only when the two planets are not in a lower order mean motion resonance. Since in lower order resonances, the variations of $`\varpi _1`$ and $`\varpi _2`$ are not guided by the secular dynamics, but the resonance angles. For example, for the $`j:(kj)`$ resonance, if the two resonance angles $`\theta _1=j\lambda _1+(kj)\lambda k\varpi _1`$, $`\theta _2=j\lambda _1+(kj)\lambda k\varpi _2`$ librate around $`0`$ or $`\pi `$, then the relative longitude of pericenter $`\mathrm{\Delta }\varpi =(\theta _1\theta _2)/k`$ must librate around either $`0`$ or $`\pi `$. Thus we think in this case, the lower order mean motion resonance and the apsidal resonance are not independent, and the former one guide the dynamics. In the case that only $`\theta _1`$ librates, we believe that the apsidal resonance is very difficult to occur, since in this case $`\varpi _1`$ is guided by $`\theta _1`$, thus it can not have similar variations with $`\varpi _2`$.

Beauge et al. (2002) find that, there may exist some asymmetric stationary solutions in mean motion resonance region, where both the resonant angles and $`\mathrm{\Delta }\varpi `$ are constants with values different from 0 or $`\pi `$. We think such kind of solutions are due to the mean motion resonance and can only exist in the resonance regions, since such kinds of apsidal resonance solutions with $`\mathrm{\Delta }\varpi `$ librates about constants with values different from 0 or $`\pi `$ can not be found in the linear secular perturbation system.

For the ten known multiple-planet systems (see Table 8 of Fischer et al. (2003) for a list of the elements), the situation of whether apsidal resonance happens between their planets can be classified roughly into three groups(see Table 3 for the extensions of $`\mathrm{\Delta }\varpi _0`$ when apsidal resonance would occur for the observed $`q,\alpha ,\frac{e_{20}}{e_{10}}`$):

(1) planets both in apsidal resonance and mean motion resonance. 55 Cancri b and c (Fig.12a), GJ876 b and c , HD82943 b and c are in apsidal resonance since they are in mean resonances 3:1, 2:1 and 2:1 respectively. In fact, HD82943 b and c can in apsidal resonances without in the mean motion resonance, while GJ876 b and c are near the boundary of libration according to the linear secular dynamics(Fig.13a).

(2) planets in apsidal resonances far away from lower order mean motion resonances. HD12661 b and c, 47 UMa b and c, Ups And c and d(Fig.12b) are in this type. They are in apsidal resonance without the existence of any strong mean resonances. Moreover, HD12661 b and c seems to be in anti-aligned apsidal resonance, which is in the boundary of a chaotic region.

(3) planets not in apsidal resonance either due to the negligible small libration region in the $`e_{20}\mathrm{\Delta }\varpi _0`$ plane or without suitable $`\mathrm{\Delta }\varpi _0`$. The two planets in HD38529, HD168443, HD74156 are not in apsidal resonance, since the libration regions in the $`e_{20}\mathrm{\Delta }\varpi _0`$ are negligible small for the observed $`q`$,$`\alpha `$ and $`e_2/e_1`$(Table 3). HD37124 b and c have large aligned libration area for the present parameters, and their eccentricities are not small, but they are not in resonance due to the present values of $`\mathrm{\Delta }\varpi _0`$ if $`\mathrm{\Omega }_b=\mathrm{\Omega }_c=0`$ are assumed (Fig.13b).

However, due to the unknown of the inclinations and the longitude of ascending nodes in the orbital fit from the observation data, it is still too early to make conclusions for some planet systems whether the planets are in apsidal resonances. For example, for the HD12661 b and c, they are believed to be in the anti-aligned libration now. This is achieved by assuming $`\mathrm{\Omega }_1=\mathrm{\Omega }_2=0`$, and they are in the boundary of a chaotic region. Alternative choices of $`\mathrm{\Omega }_i`$ (i=1,2) may change the conclusion. Especially if we chose $`\mathrm{\Omega }_1\mathrm{\Omega }_2180^0`$, the two planets will in the aligned apsidal resonance , which is stable according to Fig.9b. Similarly, HD37124 b and c could in apsidal resonance if suitable parameters $`\mathrm{\Omega }`$ are observed.

Similar secular resonance may also happen due to nearly the same averaging precessing rate of the ascending nodes between the two planets in mutually inclined orbits. Since the inclination perturbations are isolated from the eccentricity ones in the linear secular perturbation theory, we will address this problem in a separate paper (Zhou and Sun, in preparation).

We would like to thank the anonymous referee for his valuable suggestions. This work is supported by the Natural Science Foundation of China (No. 10233020), the Special Funds for Major State Basic Research Project (G200077303), and a grant from the Department of Education of China for Ph.D candidate training(20020284011).


text:
# Quasi-one-dimensional Bose gases with large scattering length

## Abstract

Bose gases confined in highly-elongated harmonic traps are investigated over a wide range of interaction strengths using quantum Monte Carlo techniques. We find that the properties of a Bose gas under tight transverse confinement are well reproduced by a 1d model Hamiltonian with contact interactions. We point out the existence of a unitary regime, where the properties of the quasi-1d Bose gas become independent of the actual value of the 3d scattering length $`a_{3d}`$. In this unitary regime, the energy of the system is well described by a hard rod equation of state. We investigate the stability of quasi-1d Bose gases with positive and negative $`a_{3d}`$.

In recent years the study of quasi-1d quantum Bose gases has attracted a great deal of interest. Intriguing properties of quasi-1d gases, such as the exact mapping between interacting bosons and non-interacting fermions, have been predicted Girardeau ; Olshanii ; Petrov . A bosonic gas that behaves as if it consisted of spinless fermions, a so-called Tonks-Girardeau (TG) gas, cannot be described within mean-field theory since it exhibits strong correlations; instead, a many-body framework is called for. While experimental evidence of quasi-1d behavior has been reported for bosonic atomic gases under highly-elongated harmonic confinement EXP1D , TG gases have not been observed yet. It has been suggested, however, that TG gases can be realized experimentally for either low atomic densities or strong atom-atom interaction strengths. The 3d $`s`$-wave scattering length $`a_{3d}`$, and hence the strength of atom-atom interactions, can be tuned to essentially any value, including zero and $`\pm \mathrm{}`$, by utilizing a magnetic atom-atom Feshbach resonance EXPFR1 ; EXPFR2 .

Utilizing a two-body Feshbach resonance, 3d degenerate gases with large scattering length $`a_{3d}`$ have been studied experimentally and theoretically. For $`a_{3d}\pm \mathrm{}`$, it is predicted that the behavior of the strongly-correlated gas is independent of $`a_{3d}`$ Pethick . For homogeneous 3d Bose gases, this unitary regime can most likely not be reached experimentally since three-body recombination is expected to set in when $`a_{3d}`$ becomes comparable to the average interparticle distance. Three-body recombination leads to cluster formation, and hence “destroys” the gas-like state. The situation is different for Fermi gases, for which the unitary regime has already been achieved experimentally EXPFR2 . In this case, the Fermi pressure stabilizes the system even for large $`|a_{3d}|`$. It has been predicted that three-body recombination processes are suppressed for strongly interacting 1d Bose gases gangardt . These studies raise the question whether a highly-elongated inhomogeneous Bose gas, that is, an inhomogeneous quasi-1d Bose gas, is stable as $`a_{3d}\pm \mathrm{}`$.

This Letter investigates the properties of a quasi-1d Bose gas at zero temperature over a wide range of interaction strengths within a microscopic, highly accurate many-body framework. We find that the system i) is well described by a 1d model Hamiltonian with contact interactions and renormalized coupling constant Olshanii  for any value of the 3d scattering length $`a_{3d}`$; ii) behaves like a TG gas for a critical positive value of $`a_{3d}`$; iii) reaches a unitary regime for large values of $`|a_{3d}|`$, where the properties of the quasi-1d Bose gas become independent of the actual value of $`a_{3d}`$ and are similar to those of a hard-rod gas; and iv) becomes unstable against cluster formation for a critical negative value of $`a_{3d}`$.

Our study is based on the 3d Hamiltonian $`H_{3d}`$,

$`H_{3d}={\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle \frac{\mathrm{}^2}{2m}}_i^2+{\displaystyle \frac{m}{2}}\left(\omega _\rho ^2\rho _i^2+\omega _z^2z_i^2\right)\right]+{\displaystyle \underset{i<j}{}}V(r_{ij}),`$ (1)
which describes $`N`$ spin-polarized mass $`m`$ bosons under highly-elongated confinement with $`\omega _z=\lambda \omega _\rho `$, where $`\lambda 1`$. The coordinates $`\rho _i=\sqrt{x_i^2+y_i^2}`$ and $`z_i`$ denote, respectively, the radial and axial position of the $`i`$-th particle, $`r_{ij}=|\stackrel{}{r}_i\stackrel{}{r}_j|`$ denotes the interparticle distance between atom $`i`$ and $`j`$, and $`V(r)`$ denotes the two-body interatomic potential. We consider two different potentials: i) a purely repulsive hard-sphere (HS) potential, $`V^{HS}(r)=\mathrm{}`$ for $`r<a_{3d}`$, and zero otherwise; and ii) a more realistic short-range (SR) potential, which can support two-body bound states, $`V^{SR}(r)=V_0/\text{cosh}^2(r/r_0)`$, with well depth $`V_0`$. For $`V^{HS}`$, the $`s`$-wave scattering length $`a_{3d}`$ coincides with the range of the potential. For $`V^{SR}`$, in contrast, $`r_0`$ determines the range of the potential, while the scattering length $`a_{3d}`$ is a function of $`r_0`$ and $`V_0`$. In our calculations, $`r_0`$ is fixed at a value much smaller than the transverse oscillator length, $`r_0=0.1a_\rho `$, where $`a_\rho =\sqrt{\mathrm{}/m\omega _\rho }`$. To simulate the behavior of $`a_{3d}`$ near a field-dependent Feshbach resonance, we vary the well depth $`V_0`$ and consequently the scattering length $`a_{3d}`$. Importantly, $`a_{3d}`$ diverges for particular values of $`V_0`$ (the inset of Fig. 1 shows one such divergence). At each divergence, a new two-body $`s`$-wave bound state is pulled in.

Our numerical calculations (see below) are performed for the well depths $`V_0`$ shown in the inset of Fig. 1.

We consider situations where the bosonic gas described by Eq. (1) is in the 1d regime for any value of the 3d scattering length $`a_{3d}`$, which implies $`N\lambda 1`$. If the range of $`V(r)`$ is much smaller than $`a_\rho `$, it is predicted Olshanii  that the properties of the 3d gas are well described by the 1d contact-interaction Hamiltonian $`H_{1d}`$,

$`H_{1d}={\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2}{z_i^2}}+{\displaystyle \frac{m}{2}}\omega _z^2z_i^2\right)+g_{1d}{\displaystyle \underset{i<j}{}}\delta (z_{ij})+`$
$`N\mathrm{}\omega _\rho ,`$ (2)
where $`g_{1d}`$ is an effective coupling constant. The limit $`g_{1d}0`$ corresponds to the weakly-interacting mean-field regime, while $`g_{1d}\mathrm{}`$ corresponds to the strongly-interacting TG regime. For positive $`g_{1d}`$ and vanishing confinement ($`\omega _\rho =\omega _z=0`$), Eq. (2) reduces to the Lieb-Liniger (LL) Hamiltonian Lieb , whose gas-like properties have been studied in detail. For negative $`g_{1d}`$ and $`\omega _\rho =\omega _z=0`$, Eq. (2) supports cluster-like bound states McGuire ; little is known about gas-like states in this case.

Below, we solve the many-body Schrödinger equation for the 1d Hamiltonian with confinement, Eq. (2), for positive and negative $`g_{1d}`$, and relate our results to those for the 3d Hamiltonian, Eq. (1). Considering two bosons in a highly-elongated geometry that interact through a regularized zero-range pseudo-potential, Olshanii Olshanii  shows that the effective 1d coupling constant $`g_{1d}`$ can be expressed in terms of the known 3d scattering length $`a_{3d}`$,

$`g_{1d}={\displaystyle \frac{2\mathrm{}^2a_{3d}}{ma_\rho ^2}}{\displaystyle \frac{1}{1Aa_{3d}/a_\rho }},`$ (3)
where $`A=|\zeta (1/2)|/\sqrt{2}=1.0326`$. Alternatively, $`g_{1d}`$ can be expressed through the effective 1d scattering length $`a_{1d}`$, $`g_{1d}=2\mathrm{}^2/(ma_{1d})`$ Olshanii , where

$`a_{1d}=a_\rho \left({\displaystyle \frac{a_\rho }{a_{3d}}}A\right).`$ (4)
For $`a_{3d}>0`$ (but $`a_{3d}a_\rho `$), $`g_{1d}`$ approaches the unrenormalized coupling constant, $`g_{1d}^0=2\mathrm{}^2a_{3d}/(ma_\rho ^2)`$, which is obtained by averaging the 3d coupling constant $`g_{3d}=4\pi \mathrm{}^2a_{3d}/m`$ over the transverse oscillator ground state (see, e.g., Ref. Petrov ). For these values of $`a_{3d}`$, the 1d Hamiltonian given in Eq. (2) with $`g_{1d}`$ replaced by $`g_{1d}^0`$, describes a quasi-1d system accurately DMC .

For $`a_{3d}a_\rho `$ and for negative $`a_{3d}`$, in contrast, the confinement induced renormalization becomes important, and the effective 1d coupling constant $`g_{1d}`$ and scattering length $`a_{1d}`$, Eqs. (3) and (4), have to be used. Figure 1 shows $`g_{1d}`$ (dashed line) and $`a_{1d}`$ (solid line) as a function of $`a_{3d}`$. At the critical value $`a_{3d}^c=0.9684a_\rho `$ (indicated by a vertical arrow in Fig. 1), $`g_{1d}`$ diverges while $`a_{1d}`$ goes through zero. At the 3d resonance, that is for $`a_{3d}\pm \mathrm{}`$, $`g_{1d}`$ and $`a_{1d}`$ each reach an asymptotic value ($`g_{1d}=1.9368a_\rho \mathrm{}\omega _\rho `$ and $`a_{1d}=1.0326a_\rho `$, respectively, indicated by horizontal arrows in Fig. 1). Tuning $`a_{3d}`$ to large values hence allows a unitary quasi-1d regime, where $`g_{1d}`$ and $`a_{1d}`$ are independent of $`a_{3d}`$, to be entered.

We solve the Schrödinger equation for the many-body 3d Hamiltonian, Eq. (1), numerically using the diffusion quantum Monte Carlo (DMC) technique. Our interest is in the lowest gas-like many-body state. For $`V^{HS}`$, the lowest gas-like state of $`H_{3d}`$ coincides with the many-body ground state; to describe this state we can hence use the “standard” DMC technique (see, e.g., Boronat ). For $`V^{SR}`$, in contrast, the many-body ground state is a cluster-like bound state. To describe the lowest-lying gas-like state, i.e., an excited state, we hence use the fixed-node DMC (FN-DMC) method Reynolds . For a given many-body nodal surface, the FN-DMC method allows the Schrödinger equation to be solved for approximate excited states. For $`N=2`$, we obtain the exact nodal surface by direct diagonalization. Assuming that the scattering properties between each atom pair are unaltered by the presence of other atoms, we then use the $`N=2`$ nodal surface to construct the many-body nodal surface. For dilute quasi-1d gases we expect that the FN-DMC approach as implemented here results in highly accurate many-body energies.

Figure 2 shows the resulting 3d energy per particle, $`E/N\mathrm{}\omega _\rho `$, as a function of $`a_{3d}`$ for $`N=5`$ under quasi-1d confinement, $`\lambda =0.01`$, for the hard-sphere and the short-range two-body potential, $`V^{HS}`$ (diamonds) and $`V^{SR}`$ (asterisks), respectively. For small $`a_{3d}/a_\rho `$, the energies for these two two-body potentials agree within the statistical uncertainty. For $`a_{3d}a_\rho `$, however, clear discrepancies are visible. The DMC energies for $`V^{SR}`$ cross the TG energy per particle (indicated by a dashed horizontal line), $`E/N\mathrm{}\omega _\rho =\mathrm{}\omega _\rho \lambda N/2`$, very close to the value $`a_{3d}^c=0.9684a_\rho `$ \[indicated by a vertical arrow in Fig. 2(b)\],

while the energies for $`V^{HS}`$ cross the TG energy per particle at a notably smaller value of $`a_{3d}`$.

To compare our results obtained for the 3d Hamiltonian, $`H_{3d}`$, with those for the 1d Hamiltonian, $`H_{1d}`$, we also solve the Schrödinger equation for $`H_{1d}`$, Eq. (2). For positive values of the coupling constant $`g_{1d}`$ we calculate the many-body ground state energy by the exact DMC method. For $`g_{1d}<0`$, however, the 1d Hamiltonian supports many-body bound states, and, as in the 3d case, we use the FN-DMC method to describe the lowest-lying gas-like state. For $`N=2`$ and $`g_{1d}<0`$, the first excited state of the Schrödinger equation for $`H_{1d}`$, Eq. (2), has a node at $`z_{12}=a_{1d}`$. To solve the many-body Schrödinger equation for $`H_{1d}`$ with negative $`g_{1d}`$ for the lowest gas-like state by the FN-DMC technique, we parameterize our many-body nodal surface in terms of $`a_{1d}`$. This many-body nodal surface is expected to be good if the density of the gas is low.

In addition to the 3d energy per particle, Fig. 2 shows the resulting 1d energy per particle obtained by solving the Schrödinger equation for $`H_{1d}`$, Eq. (2), for the renormalized coupling constant $`g_{1d}`$ \[squares, Eq. (3)\], and the unrenormalized coupling constant $`g_{1d}^0`$ (pluses), respectively. The 1d energies calculated using the two different coupling constants agree well for small $`a_{3d}`$, while clear discrepancies become apparent for larger $`a_{3d}`$. Importantly, the 1d energies calculated using the renormalized 1d coupling constant $`g_{1d}`$ agree well with the 3d energies calculated using the short-range potential $`V^{SR}`$ (asterisks) up to very large values of the 3d scattering length $`a_{3d}`$, and also for negative $`a_{3d}`$. In contrast, at large $`a_{3d}`$ the 1d energies deviate clearly from the 3d energies calculated using the hard-sphere potential $`V^{HS}`$ (diamonds). We conclude that the renormalization of the effective 1d coupling constant $`g_{1d}`$ and the 1d scattering length $`a_{1d}`$ are crucial to reproduce the results of the 3d Hamiltonian $`H_{3d}`$ when $`a_{3d}a_\rho `$ and when $`a_{3d}`$ negative. Small deviations between the 1d energies calculated using the renormalized 1d coupling constant $`g_{1d}`$ and the 3d energies calculated using the short-range potential $`V^{SR}`$ remain; we attribute these to the finite range of $`V^{SR}`$.

If the size of the cloud is much larger than the harmonic oscillator length $`a_z`$, where $`a_z=\sqrt{\mathrm{}/m\omega _z}`$, it has been shown that the properties of the 1d LL Hamiltonian $`H_{1d}`$, $`g_{1d}>0`$, are well described by a simple equation of state using the local density approximation (LDA) Dunjko . For $`g_{1d}<0`$, we instead apply the equation of state for 1d hard-rods (HRs). Recall that the many-body nodal surface of the lowest-lying gas-like state of $`H_{1d}`$ with $`g_{1d}<0`$ is well parameterized by $`a_{1d}`$. For $`z_{ij}>a_{1d}`$, the corresponding wave function coincides with that of $`N`$ 1d hard rods of size $`a_{1d}`$. For small values of the 1d gas parameter, $`n_{1d}a_{1d}1`$, where $`n_{1d}`$ denotes the linear density, we hence expect that the lowest-lying gas-like state of the 1d many-body Hamiltonian with $`g_{1d}<0`$ is well described by a system of hard rods of size $`a_{1d}`$. The exact energy per particle of the uniform hard rod system is given by $`E/N=(\pi ^2\mathrm{}^2n_{1d}^2/6m)/(1n_{1d}a_{1d})^2`$ Girardeau . For trapped systems with $`N\lambda 1`$ we obtain the expansion

$`E/N\mathrm{}\omega _\rho =\mathrm{}\omega _\rho {\displaystyle \frac{N\lambda }{2}}\left(1+{\displaystyle \frac{128\sqrt{2}}{45\pi ^2}}\sqrt{N\lambda }{\displaystyle \frac{a_{1d}}{a_\rho }}+\mathrm{}\right).`$ (5)
The first term corresponds to the energy per particle in the TG regime. Similarly, the linear density in the center of the cloud is to lowest order given by the TG result, $`n_{1d}=\sqrt{2N\lambda }/(\pi a_\rho )`$. In the unitary limit, that is for $`a_{1d}/a_\rho =1.0326`$, expression (5) becomes independent of $`a_{3d}`$, and depends only on $`N\lambda `$.

Lines in Fig. 2 show the resulting 1d energies per particle for the LL equation of state ($`g_{1d}>0`$) as well as the HR equation of state ($`g_{1d}<0`$) calculated within the LDA. Remarkably, the LDA energies nearly coincide with the 1d many-body DMC energies (pluses and squares, respectively); finite-size effects play a role only for $`a_{3d}a_\rho `$. Our calculations establish for the first time that a simple treatment, i.e., a HR equation of state treated within the LDA, describes trapped quasi-1d gases with negative coupling constant $`g_{1d}`$ well over a wide range of the 3d scattering length $`a_{3d}`$. For $`a_{3d}^{}0`$, that is, for large $`a_{1d}`$, the HR equation of state using the LDA, cannot properly describe trapped quasi-1d gases, which are expected to become unstable against formation of cluster-like many-body bound states for $`a_{1d}1/n_{1d}`$. We hence investigate the regime with negative $`a_{3d}`$ in more detail within a many-body framework.

By comparing with results for the 3d Hamiltonian, Eq. (1), we have shown above that the 1d Hamiltonian $`H_{1d}`$, Eq. (2), provides an excellent description of quasi-1d gases. We hence base our stability analysis of quasi-1d gases with large effective 1d scattering length $`a_{1d}`$ on $`H_{1d}`$. We solve the many-body Schrödinger equation for $`H_{1d}`$ using the variational quantum Monte Carlo (VMC) method. Our variational $`N`$-particle Bijl-Jastrow-type wave function consists of one- and two-body terms. The one-body terms are written as a function of a single variational parameter $`\alpha `$, which determines the size of the atomic gas. The two-body term is parameterized by $`a_{1d}`$, and explicitly accounts for correlations.

Figure 3 shows the resulting VMC energy per particle for $`N=5`$ and $`\lambda =0.01`$ as a function of the variational parameter $`\alpha `$ (Gaussian width) for four different $`a_{1d}`$. For $`a_{1d}/a_\rho =1`$ and $`2`$, Fig. 3 shows a local minimum at $`\alpha _{min}a_z`$. The minimum VMC energy nearly coincides with the essentially exact DMC energy, which suggests that our variational wave function provides a highly accurate description of the quasi-1d many-body system. The energy barrier at $`\alpha 0.2a_z`$ decreases with increasing $`a_{1d}`$, and disappears for $`a_{1d}/a_\rho 3`$. We interpret this vanishing of the energy barrier as an indication of instability stability . For small $`a_{1d}`$, the energy barrier separates the gas-like state from cluster-like bound states. For larger $`a_{1d}`$, this energy barrier disappears and the gas-like state becomes unstable against cluster formation.

We additionally performed variational calculations for larger $`N`$ and different $`\lambda `$. We find that the onset of instability of quasi-1d Bose gases can be described by the product of the 1d scattering length $`a_{1d}`$ and the linear density at the trap center $`n_{1d}`$. To be specific, our many-body calculations suggest that a quasi-1d gas is stable for $`a_{1d}n_{1d}0.35`$, and becomes unstable for $`a_{1d}n_{1d}0.35`$. Our analysis suggests that the quasi-1d unitary regime can be reached experimentally. By tuning the 3d scattering length, it is further possible to investigate the onset of instability. By reducing $`\lambda `$ one should be able to stabilize relatively large quasi-1d systems.

In conclusion, we investigated the energetics

of a Bose gas under highly-elongated harmonic confinement over a wide range of the 3d scattering length. We find that the quasi-1d gas can be described by a many-body 1d model Hamiltonian with zero-range interactions and renormalized coupling constant. For $`a_{3d}\pm \mathrm{}`$, the quasi-1d gas enters a unitary regime, where all properties of the system are independent of $`a_{3d}`$. In the vicinity of the unitary regime, the quasi-1d system behaves like a gas of HRs. For negative $`a_{3d}`$, quasi-1d gases become unstable against cluster formation for a critical value of $`a_{1d}n_{1d}`$.

GEA and SG acknowledge support by the Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR). DB acknowledges support by the NSF (grant 0331529) and by the BEC Center at the University of Trento, and BEG by the NSF through a grant to ITAMP.


text:
# X-ray Bursts from the Accreting Millisecond Pulsar XTE J1814-338

## 1 Introduction

Recently, Markwardt & Swank (2003) reported the discovery of the fifth accreting millisecond pulsar, XTE J1814-338 (hereafter J1814). This pulsar has a 314.36 Hz spin frequency, resides in a binary with an orbital operiod of 4.275 hours, and has a minimum companion mass of $`0.15M_{}`$ (Markwardt et al. 2003). This is the widest, and most massive binary of the 5 known accreting millisecond pulsars. The system parameters of J1814 closely match those of the “classical”, non-pulsing low mass X-ray binaries (LMXBs), for example the 3.8 hour binary 4U 1636-53. These systems are not generally observed as persistent pulsars, but the neutron star spins are known from the detection of oscillations during thermonuclear bursts (see Strohmayer & Bildsten 2003 for a review of burst oscillations). Evidence is emerging that the accreting millisecond pulsars must have stronger large scale magnetic fields than the non-pulsing LMXBs (Galloway et al. 2002; Cumming, Zweibel & Bildsten 2001). The fundamental quantity controlling the field strength may be the long term accretion rate, if the accretion can “bury” the field. These pulsars are all transients and/or subluminous compared to other LMXBs, suggesting a lower long term accretion rate. Magnetic fields also play a key role in the frequency stability of burst oscillations (Cumming, Zweibel & Bildsten 2001). Detailed studies of J1814 and its bursts may thus provide insight into the magnetic fields of LMXBs, and why some do not show persistent pulsations.

The recent discoveries of coherent oscillations during a superburst from 4U 1636-53 (Strohmayer & Markwardt 2002), and burst oscillations from the millisecond pulsar SAX J1808.4-3658 (hereafter J1808, Chakrabarty et al. 2003), have conclusively established that the frequency of burst oscillations are directly related to the spin of the neutron star. Since the discovery of J1814 in June, 2003 extensive observing with RXTE has been undertaken. At the time of this writing a total of 27 thermonuclear bursts have been observed. Pulse trains at the 314.36 Hz spin frequency have been detected in all of these bursts, making this source the second for which burst oscillations have been observed at the known spin frequency of the neutron star. In this Letter we report on a study of these first bursts observed from J1814. As we show below, the burst oscillations are phase locked to the persistent pulsations. We find good evidence for frequency drifting near the onset of only two of the bursts. Most interestingly, we detect the first significant harmonic structure in burst oscillation waveforms. We use spectroscopy of the bursts to obtain a distance constraint. We conclude with some discussion of the implications of our findings.

## 2 Burst Time Profiles and Spectroscopy

Some basic properties of a representative sample of the bursts are given in Table 1. Almost all the bursts last typically 2 minutes and are characterized by a rise time (as defined by the time it takes the photon count rate to grow from 10% to 100% of the peak value) of 5 to 10 seconds, and $`1/e`$ decay times, as measured from the peak onward, of 20 to 35 seconds. Two different temporal components are present: a bright and fast component with peak rates above 600 c s-1PCU-1 and durations of at most 30 seconds, and a slow component that is responsible for the overall duration. The weak bursts appear to lack the fast component. The timescale of the longer component suggests that the neutron star is accreting hydrogen-rich material at a rate where stable hydrogen burning is either not occurring (implying $`\dot{m}<9\times 10^2`$ g cm-2s-1) or the burning rate is too slow to catch up with the supply of fresh fuel (implying $`\dot{m}>2\times 10^3`$ g cm-2s-1; Fujimoto, Hanawa & Miyaji 1981, Fushiki & Lamb 1987; Cumming & Bildsten 2000). The peak count rate, corrected for pre-burst persistent levels and evaluated with 0.5 second time resolution, ranges between $``$0.4 to 1.4 kcts s-1PCU-1. Burst 27, the last burst observed, is qualitatively different from the others in that its rise time is shorter, and it reached a higher peak flux. These differences suggest a higher fraction of helium fuel for this burst. If stable hydrogen burning were taking place prior to this burst, then the lower $`\dot{m}`$ at this epoch should have allowed for more of the accumulating hydrogen to burn, perhaps resulting in a higher helium fraction.

The background-subtracted spectral data could be well fitted with single temperature black body radiation and a negligible absorbing column, $`n_\mathrm{H}`$. The evolution of the temperature and black body radius are similar for most bursts, except that the peak temperatures of the fainter bursts, and burst 27, are $`1.7`$ and 3 keV, while all others are in the range from 2.3 to 2.5 keV. None of the bursts convincingly show spectroscopic evidence of photospheric radius expansion. Therefore all bursts, with the exception of burst 27, appear to be sub-Eddington. This is also consistent with the presence of burst oscillations throughout the peaks of these bursts (Muno et al. 2002). Again, burst 27 is different in that it’s burst oscillations show a gap across the peak (see Figure 1, top right). Such gaps are seen in many radius expansion bursts with oscillations (Muno et al. 2002), suggesting that burst 27 reached the Eddington limit. As noted above, this burst also had the highest peak color temperature, $`>3`$ keV, of any of the bursts. If this burst was indeed Eddington limited then we get a distance contraint based on the peak bolometric flux of 2.64$`\times 10^8`$ erg s-1cm-2. Assuming an Eddington limit of 2$`\times 10^{38}`$ erg s-1 (for a canonical neutron star and a hydrogen dominated photosphere) we get a distance of $`8`$ kpc. The uncertainty on this distance is estimated to be roughly 20% (Kuulkers et al. 2003).

### 2.1 Burst Oscillations

For all the bursts we have 125 $`\mu `$-sec time resolution PCA event mode data. We first barycentered the data using the JPL DE405 solar system ephemeris and the source position determined from PCA scans (Markwardt & Swank 2003). Our timing analysis is performed using the $`Z_n^2`$ statistic (Buccheri et al. 1983; Strohmayer & Markwardt 2002). We first computed dynamic power spectra for each burst. We used 4 s intervals to compute $`Z_1^2`$, and we started a new interval every 0.25 seconds. Significant oscillations were detected in each burst in the vicinity of the 314 Hz pulsar frequency. The pulse trains during most of the bursts lasted from 30 - 100 seconds. These are longer than the pulse trains observed from most burst oscillation sources, and allow the time evolution of the pulsed amplitude to be explored in some detail. Most bursts show complex, episodic modulations in the amplitude of $`35\%`$ (rms) on timescales ranging from 7 - 15 s. We will explore this interesting behavior in more detail in a sequel. The dynamic spectra indicate that the frequency during most bursts is stable, showing little evidence for the drifting seen in other sources (Strohmayer & Markwardt 1999; Muno et al. 2002).

We performed a phase coherent timing analysis for each burst. We used a linear frequency evolution model, $`\nu (t)=\nu _0+\delta \nu t`$ (ie. the phase is quadratic in time), to model the phase of the oscillations. We tracked the phases in successive time bins throughout each burst and found the values of $`\nu _0`$ and $`\delta \nu `$ which gave the best fit to the data. The best fit is found by minimizing $`\chi ^2=_{i=1}^N(\varphi _i(\nu _0,\delta \nu )\mu )^2/\sigma _{\varphi _i}^2`$, where $`\mu =(1/N)_{i=1}^N\varphi _i(\nu _0,\delta \nu )`$, $`\varphi _i`$ is the phase deduced for the ith bin using the model parameters and all X-ray events in the bin, and $`N`$ is the number of phase bins.

For each burst the linear frequency model gave an acceptable fit to the phase timing data, and in no case was there a need for a linear drift term, $`\delta \nu `$, larger than the local change in frequency due to orbital motion. Figure 1 shows dynamic power spectra (top) and phase timing residuals (bottom) for bursts 11 and 27. These were the two bursts which showed significant episodes of frequency drift. In both bursts the frequency evolution shows a spin-down of a few tenths of Hz near the onset of the burst, followed by a statistically weaker recovery (spin up). This behavior appears contrary to that found for J1808, which showed a rapid, few Hz spin up during the rise of a burst (Chakrabarty et al. 2003), although the observed spin up in J1808 might conceivably be related to the recovery phase of the evolution seen in J1814. Nevertheless, the magnitude of the frequency change was much larger in J1808 ($`4`$ Hz). The observation of spin downs in J1814 seems atypical of bursters in general, since most observed drifts are spin ups (see Muno et al. 2002).

We compared the best fitting values for $`\nu _0`$ for bursts 1 - 12 with the predicted frequency of the persistent pulsations based on the orbital solution of Markwardt et al. (2003). The results for these 12 bursts are shown in Figure 2. The frequencies of the oscillations during each burst are consistent with the predicted, orbitally modulated, pulsar frequency. Not only are the frequencies consistent, but the phase of the burst oscillations is strongly locked to that of the persistent pulsations. In Figure 3 we compare the phase folded pulse profile from the first 56 ksecs of non-burst data (thick solid curve) with the co-added pulse profiles from bursts 2 and 3. The burst data were folded with the same ephemeris as the persistent pulse data. Both the fundamental and 1st harmonic components of the profile are aligned. Typical 1$`\sigma `$ upper limits to any average phase offsets during bursts for the fundamental and 1st harmonic components are $`0.7`$ and $`1.2\%`$, respectively. These limits are smaller than the corresponding phase offset for bursts from J1808 of $`11\%`$ (Chakrabarty et al. 2003).

We computed the average amplitudes of the burst oscillations from the phase folded pulse profiles (see Table 1). We define the amplitude as $`a=(I_{max}I_{min})/(I_{max}+I_{min})`$, where $`I`$ is the countrate of the folded profile. For these estimates we have not subtracted off the preburst, persistent emission. Doing so would increase the amplitudes modestly. The amplitudes are in the range from 14 - 18 %, except for burst 27, which showed systematicalluy weaker oscillations. This is larger than the amplitude of the persistent pulsations, which ranges from about 9 - 12 % when measured in the same way after first subtracting off the detector background. This confirms that the black body burst flux must be modulated at the spin frequency. If it were not the amplitude during the bursts would be less than 1%.

The burst oscillations also show a significant harmonic component in their pulse profiles. We demonstrate this in Figure 4, which shows the co-added pulse profile from bursts 2 and 3, as well as a model fit using two Fourier components. The solid curve is the best model, while the dashed and dotted lines show the fundamental and first harmonic components. The harmonic component is strongly required, with an amplitude of $`36\pm 2.6`$ counts s-1. Co-adding bursts 2 and 3 in phase gives $`Z_1^2=95.2`$ at the first harmonic. Since $`Z_1^2`$ is distributed as $`\chi ^2`$ with two degrees of freedom, the value of 95.2 is a highly significant detection. The harmonic content varies significantly from burst to burst. When expressed as the ratio of pulsed amplitude at the first harmonic to that at the fundamental, it ranges from essentially 0 – 0.27. We note that this is always less than the harmonic content of the persistent pulsations, for which the same ratio is $`0.33`$. In terms of countrate the first harmonic component for the persistent pulsations around the epoch of bursts 2 and 3 is about 4.5 counts s-1, which is much to small to account for the harmonic content during these bursts.

## 3 Discussion and Summary

J1814 is now the second pulsar to show burst oscillations at the spin frequency of the neutron star. The large amplitude of the modulations confirms that the burst flux is spin modulated, as first suggested by Strohmayer et al. (1996). As in J1808 (Chakrabarty et al. 2003), the phase of the burst oscillations in J1814 are linked to the phase of the persistent pulsations. It seems likely that this phasing results from the magnetic field of the neutron star. In persistent pulsars the magnetic field must be strong enough to channel the accretion flow and produce an asymmetric accretion luminosity. Such a field would then likely also influence the surface distribution of nuclear fuel, resulting in a thermonuclear flux distribution during bursts which tracks the magnetic field geometry. A strong field will also enforce co-rotation of the surface layers, minimizing shearing and thus frequency drifting. The fact that the burst oscillation sources which do not show persistent pulsations also have the frequency drifts with the longest relaxation time seems consistent with this idea (Strohmayer & Markwardt 1999).

Chakrabarty et al. (2003) found evidence for frequency drifting at the onset of a burst from J1808. Although most of the bursts from J1814 that we have examined show no strong evidence of drifting, bursts 11 and 27 are the exceptions. These bursts show evidence for phase drifting near burst onset (Figure 2, lower panels). For these bursts both the phase residuals and the dynamic power spectra suggest that the oscillation frequency was initially higher than the spin frequency and then decreased below the spin frequency before recovering again. For both bursts the evolution takes place in an $`5`$ s interval during the burst rise.

Perhaps the most interesting result is the detection of significant harmonic content in the burst oscillations. Two important conclusions are that the harmonic content must be intrinsic to the modulation of the thermonuclear flux, and secondly, the harmonic content of the burst oscillations, expressed as the ratio of power at the first harmonic to that at the fundamental, is always less than the harmonic content of the persistent pulsations. As noted by several authors (Strohmayer et al. 1998; Miller & Lamb 1998; Weinberg, Miller & Lamb 2001; Muno et al. 2002), the harmonic content of burst oscillations can in principle be used to constrain the neutron star radius. This is because of the relativistic effects caused by the high surface velocities of these pulsars. The bursts from J1814 with the highest first harmonic content have a harmonic to fundamental ratio greater than 0.25. In contrast, upper limits on this ratio in other burst oscillation sources can be as small as 0.05 (Muno et al. 2002). For a 10 km neutron star, $`v_{surf}0.07c`$ for the spin frequency of J1814. At this velocity such a large harmonic ratio is extremely difficult to produce if the emission is isotropic from a single bright region on the neutron star surface (Weinberg, Miller & Lamb 2001; Muno, Ozel & Chakrabarty 2002). Our results suggest that substantial beaming of the radiation likely occurs if the modulation is from a single bright region on the surface. Alternatively, a geometry with two antipodal bright regions is possible. Detailed modeling of the energy dependent pulse profiles could provide constraints on the stellar compactness, radius and perhaps the magnetic field strength and geometry.

## 4 Figure Captions