text:
# Nonradial Oscillations of Neutron Stars with a Solid Crust

## 1 Introduction

After the discovery of the $`r`$ mode instability in neutron stars driven by gravitational radiation reaction (Andersson 1998; Friedman & Morsink 1998), nonradial oscillations of relativistic neutron stars has attracted much wider interest than before in astrophysics (see a recent review, e.g., by Andersson & Kokkotas 2001). The relativistic formulation of nonradial oscillations of $`fluid`$ neutron stars was first given by Thorne and his collaborators (e.g., Thorne & Campolattaro 1967; Price & Thorne 1969; Thorne 1969a,b), and later on extended to the case of neutron stars with a solid crust in their interior (Schumaker & Thorne 1983, see also Finn 1990). Since one of the main concerns of these studies was gravitational waves generated by the stellar pulsations (e.g., Thorne 1969a), it was essential to include the perturbations in the metric $`g_{\alpha \beta }`$ in order to obtain a consistent description of the gravitational waves. However, it was the metric perturbations that made it extremely difficult to treat both analytically and numerically the oscillation modes in relativistic stars.

McDermott, van Horn, & Scholl (1983) introduced a relativistic version of the Cowling approximation, in which all the Eulerian metric perturbations $`\delta g_{\alpha \beta }`$ are neglected in the relativistic oscillation equations derived by Thorne & Campolattaro (1967). Under the relativistic Cowling approximation, they calculated the $`p`$-, $`f`$-, and $`g`$\- modes of relativistic $`fluid`$ neutron star models to examine how the oscillation modes depend on the model properties. The relativistic Cowling approximation employed by McDermott, van Horn, & Scholl (1983) was shown to be good enough to calculate the $`p`$-modes for non-rotating polytropic stars by Lindblom & Splinter (1990), who compared the fully relativistic oscillation frequencies to those obtained in the relativistic Cowling approximation. Assuming slow rotation, Yoshida & Kojima (1997) have also carried out similar computations for the $`f`$\- and $`p`$-modes of polytropic models and they confirmed the good applicability of the relativistic Cowling approximation. Quite recently, Yoshida & Lee (2002) have shown in the relativistic Cowling approximation the existence of the relativistic $`r`$ modes111We should note that lots of studies about the general relativistic $`r`$-mode have been done so far (e.g., Kojima 1998; Kojima & Hosonuma 2000; Lockitch, Andersson & Friedman 2001; Yoshida 2001; Ruoff & Kokkotas 2001). with $`l=m`$, which are regarded as a counter part of the Newtonian $`r`$ modes. We believe that the relativistic Cowling approximation is quite useful to investigate the oscillation modes of relativistic stars, although we understand that under the approximation we cannot discuss inherently relativistic oscillation modes like $`w`$-modes (see, e.g., Andersson, Kojima, & Kokkotas 1996).

In their Newtonian calculations, McDermott, van Horn, & Hansen (1988) have shown that cold neutron stars with a solid crust can support a rich variety of oscillation modes (see also Strohmayer 1991; Lee & Strohmayer 1996; Yoshida & Lee 2001). The purposes of this paper are to calculate the various oscillation modes of relativistic neutron stars with a solid crust in the relativistic Cowling approximation, and to discuss the effects of general relativity on the modal property of the oscillation modes. We regard this paper as an extension of the studies by McDermott, van Horn, & Scholl (1983) and McDermott, van Horn, & Hansen (1988). In §2, we derive relativistic oscillation equations for the solid crust in the relativistic Cowling approximation following the formulation developed by Schumaker & Thorne (1983) and Finn (1990). Numerical results are discussed in §3 for nonradial modes of three relativistic neutron star models with a solid crust. §4 is devoted to discussion and conclusion. In this paper, we use units in which $`c=G=1`$, where $`c`$ and $`G`$ denote the velocity of light and the gravitational constant, respectively.

## 2 Formulation

For modal analysis of neutron stars with a solid crust, we solve general relativistic pulsation equations derived under the relativistic Cowling approximation, in which all the metric perturbations in the matter equations are neglected. We assume that the solid crust is in a strain-free state in the equilibrium unperturbed state and the strain in the crust are generated by small amplitude perturbations superposed on the unperturbed state (see, e.g., Aki & Richards 1980). For the background unperturbed state, it is therefore possible to assume a static and spherical symmetric state, for which the geometry is given by the following line element:

$`(ds^2)_{(0)}=g_{(0)\alpha \beta }dx^\alpha dx^\beta =e^{2\nu (r)}dt^2+e^{2\lambda (r)}dr^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\phi ^2,`$ (1)
where the subscript $`(0)`$ is used to refer to the unperturbed state. The four-velocity of the stellar material in the equilibrium is given by

$$u_{(0)}^\alpha =e^{\nu (r)}t^\alpha ,$$
(2)
where $`t^\alpha `$ denotes the timelike Killing vector of the unperturbed spacetime.

In the relativistic Cowling approximation, the basic equations for pulsations are obtained from the energy and momentum conservation laws:

$`u^\alpha _\beta T_\alpha ^\beta =0,(\mathrm{energy}\mathrm{conservation}\mathrm{law})`$ (3)
$`q_\gamma ^\alpha _\beta T_\alpha ^\beta =0,(\mathrm{momentum}\mathrm{conservation}\mathrm{law})`$ (4)
where $`_\alpha `$ is the covariant derivative associated with the metric, $`T_\beta ^\alpha `$ is the energy-momentum tensor, and $`q_\beta ^\alpha `$ is the projection tensor with respect to the fluid four-velocity $`u^\alpha `$, which is defined by

$`q_\beta ^\alpha =\delta _\beta ^\alpha +u^\alpha u_\beta ,`$ (5)
where $`\delta _\alpha ^\beta `$ denotes the Kronecker delta.

For relativistic analysis of vibration of the solid crust in neutron stars, we follow the formulation given by Schumaker & Thorne (1983) and Finn (1990). If we define the rate of shear $`\sigma _{\alpha \beta }`$ and the rate of expansion $`\theta `$ as

$$\sigma _{\alpha \beta }=\frac{1}{2}q_\alpha ^\gamma q_\beta ^\delta (_\gamma u_\delta +_\delta u_\gamma )\frac{1}{3}q_{\alpha \beta }\theta ,$$
(6)
$$\theta =q_\alpha ^\beta _\beta u^\alpha =_\alpha u^\alpha ,$$
(7)
the shear strain tensor $`\mathrm{\Sigma }_{\alpha \beta }`$ is determined, in terms of $`\sigma _{\alpha \beta }`$ and $`\theta `$, as a solution of the differential equation given by

$$L_u\mathrm{\Sigma }_{\alpha \beta }=\frac{2}{3}\theta \mathrm{\Sigma }_{\alpha \beta }+\sigma _{\alpha \beta },$$
(8)
where $`L_u`$ is the Lie derivative along the four-velocity field of the matter (Carter & Quintana 1972). Once equation (8) is solved to give the shear strain tensor $`\mathrm{\Sigma }_{\alpha \beta }`$ in terms of $`\sigma _{\alpha \beta }`$ and $`\theta `$ (see Schumaker & Thorne 1983, Finn 1990), the total stress-energy tensor is given by

$$T_{\alpha \beta }=\rho u_\alpha u_\beta +pq_{\alpha \beta }2\mu \mathrm{\Sigma }_{\alpha \beta },$$
(9)
where $`\mu `$ stands for the isotropic shear modulus, and we have assumed a Hookean relationship between the shear strain and stress tensors. Here, $`\rho `$ and $`p`$ mean the mass-energy density and the isotropic pressure, respectively. Note that $`\mathrm{\Sigma }_{(0)\alpha \beta }=0`$ since we have assumed the strain-free state in the equilibrium of the star.

To discuss nonradial oscillations of a star, we employ a Lagrangian perturbation formalism (see, e.g., Friedman & Schutz 1975; Friedman 1978), in which a Lagrangian displacement vector is introduced to connect fluid elements in the equilibrium state to the corresponding elements in the perturbed state. In this formalism, the Lagrangian change $`\mathrm{\Delta }Q`$ in a quantity is related to its Eulerian change $`\delta Q`$ by

$`\mathrm{\Delta }Q=\delta Q+L_\zeta Q,`$ (10)
where $`L_\zeta `$ denotes the Lie derivative along the displacement vector $`\zeta ^\alpha `$. If we apply this formalism to oscillations of relativistic stars, we have, for example,

$`\mathrm{\Delta }g_{\alpha \beta }=\delta g_{\alpha \beta }+_\alpha \zeta _\beta +_\beta \zeta _\alpha ,`$ (11)
and

$`\mathrm{\Delta }u^\alpha ={\displaystyle \frac{1}{2}}u^\alpha u^\beta u^\gamma \mathrm{\Delta }g_{\beta \gamma }.`$ (12)
The Eulerian perturbation in the velocity field $`\delta \widehat{u}^\beta q_{(0)}^{}{}_{\beta }{}^{\alpha }\delta u^\beta `$ is then given by

$`\delta \widehat{u}^\alpha =q_{(0)}^{}{}_{\beta }{}^{\alpha }(L_{u_{(0)}}\zeta )^\beta .`$ (13)
Notice that $`\delta g_{\alpha \beta }=0`$ and hence $`\delta \widehat{u}^\alpha =\delta u^\alpha `$ in the relativistic Cowling approximation. Because the background state is in a hydrostatic equilibrium, the time dependence of all the perturbed quantities can be given by $`e^{i\sigma t}`$, where $`\sigma `$ is a constant frequency measured by an inertial observer at the spatial infinity. The relation between the Lagrangian displacement $`\zeta ^\alpha `$ and the velocity perturbation $`\delta \widehat{u}^\alpha `$ is then given as an algebraic equation:

$`\delta \widehat{u}^\alpha =i\sigma e^{\nu (r)}\zeta ^\alpha .`$ (14)
Note that the gauge freedom in $`\zeta ^\alpha `$ has been used to demand the relation $`u_\alpha \zeta ^\alpha =0`$. The velocity field of a neutron star in a perturbed state may be written as

$`u^\alpha =u_{(0)}^\alpha +\eta \delta \widehat{u}^\alpha +O(\eta ^2),`$ (15)
where $`\eta `$ is a small expansion parameter introduced for our convenience. The mass-energy density $`\rho `$ and the pressure $`p`$ in the perturbed state are also given by

$$\rho =\rho _{(0)}+\eta (\rho _{(0)}+p_{(0)})\left(\frac{\mathrm{\Gamma }\delta p}{p_{(0)}}\zeta ^\alpha A_\alpha \right)+O(\eta ^2),$$
(16)
$$p=p_{(0)}+\eta \delta p+O(\eta ^2),$$
(17)
where $`A_\alpha `$ is the relativistic Schwarzschild discriminant defined as

$$A_\alpha =\frac{1}{\rho _{(0)}+p_{(0)}}_\alpha \rho _{(0)}\frac{1}{\mathrm{\Gamma }p_{(0)}}_\alpha p_{(0)},$$
(18)
and $`\mathrm{\Gamma }`$ is the adiabatic index defined as

$`\mathrm{\Gamma }={\displaystyle \frac{\rho _{(0)}+p_{(0)}}{p_{(0)}}}\left({\displaystyle \frac{p_{(0)}}{\rho _{(0)}}}\right)_{ad}.`$ (19)
To derive equation (16), we have used the adiabatic condition for the perturbations:

$`\mathrm{\Delta }p={\displaystyle \frac{\mathrm{\Gamma }p_{(0)}}{\rho _{(0)}+p_{(0)}}}\mathrm{\Delta }\rho .`$ (20)
Assuming spherical symmetry of the background spacetime, we may expand the perturbations in terms of appropriate tensor spherical harmonic functions. The Lagrangian displacement, $`\zeta ^k`$ and the pressure perturbation, $`\delta p/(\rho _{(0)}+p_{(0)})`$ can be written as

$$\zeta ^r=rS_l(r)Y_l^m(\theta ,\phi )e^{i\sigma t},$$
(21)
$$\zeta ^\theta =\left(H_l(r)\frac{Y_l^m(\theta ,\phi )}{\theta }T_l^{}(r)\frac{1}{\mathrm{sin}\theta }\frac{Y_l^{}^m(\theta ,\phi )}{\phi }\right)e^{i\sigma t},$$
(22)
$$\zeta ^\phi =\frac{1}{\mathrm{sin}^2\theta }\left(H_l(r)\frac{Y_l^m(\theta ,\phi )}{\phi }+T_l^{}(r)\mathrm{sin}\theta \frac{Y_l^{}^m(\theta ,\phi )}{\theta }\right)e^{i\sigma t},$$
(23)
$$\frac{\delta p}{\rho _{(0)}+p_{(0)}}=\delta U_l(r)Y_l^m(\theta ,\phi )e^{i\sigma t},$$
(24)
where $`Y_l^m`$ is a spherical harmonic function (Regge & Wheeler 1957; Thorne 1980).

By substituting equations (9), (15), (16), and (17), together with the perturbations defined by equations (21) to (24), into equations (3) and (4), and collecting the terms proportional to the small expansion parameter $`\eta `$, we obtain a system of perturbation equations for the solid crust:

$`r{\displaystyle \frac{dz_1}{dr}}=\left(1+{\displaystyle \frac{2\alpha _2}{\alpha _3}}+U_2\right)z_1+{\displaystyle \frac{1}{\alpha _3}}z_2+{\displaystyle \frac{\alpha _2}{\alpha _3}}l(l+1)z_3,`$ (25)
$`r{\displaystyle \frac{dz_2}{dr}}`$ $`=`$ $`\left\{(3U_2+U_1e^{2\lambda }c_1\overline{\sigma }^2)V_1+{\displaystyle \frac{4\alpha _1}{\alpha _3}}(3\alpha _2+2\alpha _1)\right\}z_1+\left(V_24{\displaystyle \frac{\alpha _1}{\alpha _3}}\right)z_2`$ (26)
$`+\left\{V_12\alpha _1\left(1+{\displaystyle \frac{2\alpha _2}{\alpha _3}}\right)\right\}l(l+1)z_3+e^{2\lambda }l(l+1)z_4,`$
$`r{\displaystyle \frac{dz_3}{dr}}=e^{2\lambda }z_1+{\displaystyle \frac{e^{2\lambda }}{\alpha _1}}z_4,`$ (27)
$`r{\displaystyle \frac{dz_4}{dr}}`$ $`=`$ $`\left(V_1+6\mathrm{\Gamma }{\displaystyle \frac{\alpha _1}{\alpha _3}}\right)z_1{\displaystyle \frac{\alpha _2}{\alpha _3}}z_2`$ (28)
$`\left\{c_1\overline{\sigma }^2V_1+2\alpha _1{\displaystyle \frac{2\alpha _1}{\alpha _3}}(\alpha _2+\alpha _3)l(l+1)\right\}z_3(3+U_2V_2)z_4,`$
$`r{\displaystyle \frac{dz_5}{dr}}={\displaystyle \frac{e^{2\lambda }}{\alpha _1}}z_6,`$ (29)
$`r{\displaystyle \frac{dz_6}{dr}}=(3+U_2V_2)z_6\left\{c_1\overline{\sigma }^2V_1\alpha _1(l^{}+1)(l^{}2)\right\}z_5,`$ (30)
where the dependent variables $`z_1`$ to $`z_6`$ are defined as

$`z_1=S_l(r),`$ (31)
$`z_2=2\alpha _1e^\lambda {\displaystyle \frac{d}{dr}}\left(re^\lambda S_l(r)\right)+\left(\mathrm{\Gamma }{\displaystyle \frac{2}{3}}\alpha _1\right)\left\{{\displaystyle \frac{e^\lambda }{r^2}}{\displaystyle \frac{d}{dr}}\left(r^3e^\lambda S_l(r)\right)l(l+1)H_l(r)\right\},`$ (32)
$`z_3=H_l(r),`$ (33)
$`z_4=\alpha _1\left(e^{2\lambda }r{\displaystyle \frac{dH_l(r)}{dr}}+S_l(r)\right),`$ (34)
$`z_5=T_l^{}(r),`$ (35)
$`z_6=\alpha _1e^{2\lambda }r{\displaystyle \frac{dT_l^{}(r)}{dr}}`$ (36)
and the various quantities which appear in the coefficients are

$`\alpha _1={\displaystyle \frac{\mu }{p_{(0)}}},\alpha _2=\mathrm{\Gamma }{\displaystyle \frac{2}{3}}\alpha _1,\alpha _3=\mathrm{\Gamma }+{\displaystyle \frac{4}{3}}\alpha _1,`$ (37)
$`V_1=\left(1+{\displaystyle \frac{\rho _{(0)}}{\epsilon p_{(0)}}}\right)r{\displaystyle \frac{d\nu }{dr}},V_2={\displaystyle \frac{\rho _{(0)}}{\epsilon p_{(0)}}}r{\displaystyle \frac{d\nu }{dr}},U_1=\left({\displaystyle \frac{d\nu }{dr}}\right)^1{\displaystyle \frac{d}{dr}}\left(r{\displaystyle \frac{d\nu }{dr}}\right),U_2=r{\displaystyle \frac{d\lambda }{dr}},`$ (38)
$`c_1=\epsilon {\displaystyle \frac{M}{R^3}}re^{2\epsilon \nu }\left({\displaystyle \frac{d\nu }{dr}}\right)^1,`$ (39)
$`e^{2\lambda }=\left(1\epsilon {\displaystyle \frac{2M(r)}{r}}\right)^1,M(r)={\displaystyle _0^r}4\pi r^2\rho _{(0)}𝑑r,`$ (40)
$`r{\displaystyle \frac{d\nu }{dr}}=\epsilon e^{2\lambda }\left(4\pi r^2\epsilon p_{(0)}+{\displaystyle \frac{M(r)}{r}}\right),r{\displaystyle \frac{d\lambda }{dr}}=\epsilon e^{2\lambda }\left(4\pi r^2\rho _{(0)}{\displaystyle \frac{M(r)}{r}}\right).`$ (41)
Here, $`M=M(R)`$ and $`R`$ are the mass and the radius of the star, and $`\overline{\sigma }=\sigma \sqrt{R^3/M}`$ is the dimensionless frequency. In equations (38)-(41), we have introduced a factor $`\epsilon `$ to indicate the hidden factor $`1/c^2`$ that represents the strength of the general relativistic degree. Fully relativistic oscillation equations for the solid crust are regained for $`\epsilon =1`$, and the Newtonian oscillation equations are obtained in the limit of $`\epsilon 0`$, which are the same as those derived by McDermott, van Horn, & Hansen (1988). We can see that equations (25)-(28) are decoupled from equations (29)-(30) because the background spacetime is spherical symmetric. The former describe spheroidal (or polar parity) oscillations, and the latter toroidal (or axial parity) ones.

In the fluid regions, the oscillation equations to be solved are given by

$`r{\displaystyle \frac{dy_1}{dr}}=\left(3{\displaystyle \frac{V_1}{\mathrm{\Gamma }}}+U_2\right)y_1\left({\displaystyle \frac{V_1}{\mathrm{\Gamma }}}{\displaystyle \frac{l(l+1)}{c_1\overline{\sigma }^2}}\right)y_2,`$ (42)
$`r{\displaystyle \frac{dy_2}{dr}}=(e^{2\lambda }c_1\overline{\sigma }^2+rA_r)y_1(U+rA_r)y_2,`$ (43)
where

$`y_1=S_l(r),y_2=\left(r{\displaystyle \frac{d\nu }{dr}}\right)^1\delta U_l(r)=c_1\overline{\sigma }^2H_l.`$ (44)
Note that equations (42) and (43), which come from equations (25)-(28) in the limit of $`\mu 0`$, are essentially the same as those derived by McDermott, van Horn, & Scholl (1983) for spheroidal oscillation modes. Equations (29)-(30), which describe toroidal torsional oscillations, become trivial in the limit of $`\mu 0`$ because there is no restoring force that is responsible for the toroidal torsional oscillations. To calculate the spheroidal modes, equations (25)-(28) are integrated in the solid crust, and equations (42)- (43) in the fluid core and the surface ocean. For the toroidal torsional modes, equations (29)-(30) should be solved only in the solid crust.

For the spheroidal modes the outer boundary condition is given at the stellar surface ($`r=R`$) by $`\mathrm{\Delta }p=0`$, which reduces to

$`y_1y_2=0,`$ (45)
and the inner boundary condition is the regularity condition at the stellar center given by

$`c_1\overline{\sigma }^2y_1ly_2=0.`$ (46)
The jump conditions at the interfaces between the fluid and the solid regions are given by the continuity conditions, across the perturbed interface, of the stress $`P_\alpha `$ on the interface, where

$`P_\alpha =(pq_{\alpha \beta }2\mu \mathrm{\Sigma }_{\alpha \beta })N^\beta ,`$ (47)
and $`N^\beta `$ is the unit vector normal to the perturbed interface (see, e.g., Finn 1990). The normal form of the core/crust or crust/ocean interface is given by $`d(r\eta \zeta ^r)=n_\alpha dx^\alpha `$ for the radial coordinate of the perturbed interface given by $`r=R_i+\eta \zeta ^r`$, where $`R_i`$ is the radius of the core/crust or crust/ocean interface in the unperturbed state. Thus, the corresponding unit normal one-form $`N_\alpha `$ have the components as follows:

$`N_t`$ $`=`$ $`e^\lambda \eta i\sigma rS_l(r)Y_l^m(\theta ,\phi )e^{i\sigma t}+O(\eta ^2),`$
$`N_r`$ $`=`$ $`e^\lambda +O(\eta ^2),`$ (48)
$`N_\theta `$ $`=`$ $`e^\lambda \eta rS_l(r){\displaystyle \frac{Y_l^m(\theta ,\phi )}{\theta }}e^{i\sigma t}+O(\eta ^2),`$
$`N_\phi `$ $`=`$ $`e^\lambda \eta rS_l(r){\displaystyle \frac{Y_l^m(\theta ,\phi )}{\phi }}e^{i\sigma t}+O(\eta ^2),`$
The jump conditions for the spheroidal modes at the interface are then given by

$`y_1=z_1,V_1(y_1y_2)=z_2,z_4=0.`$ (49)
For the toroidal torsional modes, the boundary conditions are applied at the top and the bottom of the solid crust, and they are

$`z_6=0.`$ (50)
By imposing the boundary conditions given above, we can solve our basic equations as an eigenvalue problem with respect to the eigenvalue $`\overline{\sigma }`$. Here, we employ a Henyey type relaxation method to obtain numerical solutions to our basic equations (see, e.g., Unno et al. 1989).

## 3 Numerical Results

Neutron star models that we use in this paper are the same as those used in the modal analysis by McDermott, Van Horn, & Hansen (1988). These models are taken from the evolutionary sequences for cooling neutron stars calculated by Richardson et al. (1982), where the envelope structure is constructed by following Gudmundsson, Pethick & Epstein (1983). They are composed of a fluid core, a solid crust and a surface fluid ocean. The interior temperature is finite and is not constant as a function of the radial distance $`r`$. The models are not barotropic and the Schwarzschild discriminant $`|A|`$ has finite values in the interior of the star. The models we use are called NS05T7, NS05T8, and NS13T8, and their physical properties such as the total mass $`M`$, the radius $`R`$, the central density $`\rho _c`$, the central temperature $`T_c`$ and the relativistic factor $`GM/c^2R`$ are summarized in Table 1 (for other quantities, see McDermott et al. 1988).

To classify the various oscillation modes of the three component models, we use almost the same nomenclature as that employed by McDermott et al. (1988). We let $`p_k`$ refer to the acoustic modes of the $`k`$ th overtone. The internal gravity modes confined in the surface ocean are denoted as $`g_k^s`$ and those in the fluid core as $`g_k^c`$, where $`k`$ indicates the overtone number. The eigenfrequencies of the $`f`$ modes are usually found between the $`p_1`$ and $`g_1`$ modes. Associated with the fluid-solid interfaces in the models, there are two interfacial modes, which we denote as $`i_1`$ and $`i_2`$ such that $`\sigma (i_1)\sigma (i_2)`$. The $`s_k`$ modes are spheroidal shear dominated modes of the $`k`$ th overtone, the amplitudes of which are strongly confined in the solid crust. The $`t_k`$ modes are toroidal shear dominated modes of the $`k`$ th overtone propagating only in the solid crust. Note that the $`p_k`$-, $`g_k`$-, $`f`$-, $`i_{1(2)}`$, and $`s_k`$ modes are classified as spheroidal modes while the $`t_k`$ modes as toroidal modes.

In the Newtonian limit of $`\epsilon 0`$, we calculate various oscillation modes with $`l=2`$ for the three neutron star models NS05T7, NS05T8, and NS13T8, and tabulate the Newtonian eigenfrequencies, which we denote as $`\overline{\sigma }_N`$, in Tables 2 through 4. These tables confirm that the Newtonian frequencies obtained in this paper are in good agreement with those computed by McDermott et al. (1988). Assuming $`\epsilon =1`$, we also compute the corresponding relativistic oscillation modes with $`l=2`$, and tabulate the relativistic eigenfrequencies, which we denote as $`\overline{\sigma }_{GR}`$, in Tables 2 through 4, in which the relative frequency differences defined as $`\delta \sigma /\sigma (\sigma _{GR}\sigma _N)/\sigma _{GR}`$ are also given. The relative difference $`\delta \sigma /\sigma `$ takes various values for different oscillation modes, and the value of $`\delta \sigma /\sigma `$ of a given mode is also dependent on the physical properties of the neutron star models. We find $`|\delta \sigma /\sigma |<0.1`$ for most of the oscillation modes we calculate. The elastic $`s_k`$ and $`t_k`$ modes with $`k1`$ have particularly small values of $`|\delta \sigma /\sigma |`$, which are at most a few percent. However, there are some exceptions, and the values of $`\delta \sigma /\sigma `$ become as large as $`0.45`$ for the $`t_0`$ modes and $`1`$ for the $`g^s`$\- and $`g^c`$-modes for the model NS13T8, which is the most compact one among the three. If we consider the gravitational redshift as one of the general relativistic effects that are responsible for the difference from the Newtonian oscillation frequency, we may have $`\stackrel{~}{\sigma }_{GR}=e^{\nu (r)}\sigma _N`$ and hence $`(\stackrel{~}{\sigma }_{GR}\sigma _N)/\stackrel{~}{\sigma }_{GR}=1e^{\nu (r)}`$, which is negative since $`\nu (r)<0`$. Although this estimation gives a consistent result for the surface $`g^s`$ modes confined in the surface ocean and for the core $`g^c`$ modes trapped in the fluid core, it does not always give good estimations of the difference as suggested by the tables. In fact, we have positive $`\delta \sigma /\sigma `$ for some oscillation modes. This is understandable, however, since the eigenfunctions usually possess large radial extent in the interior of stars, and the eigenfrequencies are dependent on the equilibrium quantities which contain terms and factors which are absent in the Newtonian oscillation equations.

For the various oscillation modes of NS13T8, we display the displacement vector $`\zeta ^\alpha `$ versus the radial coordinates $`r/R`$ or $`\mathrm{log}(1r/R)`$ in Figures 1 through 8, where the amplitude normalization of the eigenfunctions is given by $`y_1=1`$ at $`r=R`$ for the spheroidal modes and $`z_5=1`$ at $`r=R_{out}`$ for the toroidal modes with $`R_{out}`$ being the radius of the ocean-crust interface. Notice that for the $`g^c`$ modes, the normalization condition $`H_l(R_{in})=1`$ is adapted because the modes are very insensitive to the behavior of the eigenfunctions in the surface ocean, where $`R_{in}`$ is the radius at the bottom of the crust. The eigenfunctions obtained in the Newtonian calculations are also plotted for comparison in each of the figures. From the figures 1 through 8, we can see that the basic properties of the eigenfunctions of the relativistic modes are not very much different from those of the corresponding Newtonian modes, except for the interfacial modes (see the next paragraph). In this sense, we can say that general relativity does not bring about any essential changes in the modal properties represented by the eigenfunctions.

As shown by Figures 3 and 4, the relativistic eigenfunctions of the interfacial modes $`i_1`$ and $`i_2`$ are at first sight quite different from the Newtonian eigenfunctions. Since the eigenfunctions of the relativistic $`i_{1(2)}`$ mode are rather similar to those of the Newtonian $`i_{2(1)}`$ mode, it is tempting to make a guess that the modal properties have been exchanged between the two modes when we move from Newtonian dynamics to relativistic dynamics. In their Newtonian calculation, McDermott et al. (1988) showed that the properties of the $`i`$ modes are very sensitive to the changes in the physical quantities near the interfaces, and mentioned a numerical experiment in which they computed the two interfacial modes by artificially reducing the bulk modulus $`\mu `$ and found an avoided crossing between the two modes through which the modal characters are interchanged with each other. Inspired by this report, we have carried out a numerical experiment, calculating the two interfacial modes as a function of $`\epsilon `$ for the model NS13T8. Plotting the frequencies of the two modes versus $`\epsilon `$ in Figure 9, we confirm the occurrence of an avoided crossing between the two modes at $`\epsilon 0.67`$, through which the modal characters of the two have been exchanged with each other. As $`\epsilon `$ increases from $`\epsilon =0`$ (Newtonian limit), the effective local shear modulus $`\mu `$ in the solid crust decreases because of the gravitational redshift effects, as a result of which the two interfacial modes experience the avoided crossing when we move from Newtonian dynamics to relativistic dynamics.

## 4 Conclusion

In this paper, we have calculated a variety of oscillation modes of relativistic neutron stars with a solid crust in the relativistic Cowling approximation, in which all metric perturbations are ignored. We find $`|\delta \sigma /\sigma |`$ is less than $`0.1`$ for most of the oscillation modes with $`l=2`$, although there are some exceptions such as the surface gravity modes, the core gravity modes, and the nodeless toroidal torsional modes. We also find that the essential modal properties represented by the eigenfunctions are not strongly affected by introducing general relativity in the sense that the relativistic eigenfunctions of a mode have almost the same properties as the corresponding Newtonian eigenfunctions. An exception may be the two interfacial modes whose eigenfunctions differ from their Newtonian eigenfunctions, and we have shown that this can be explained as a result of the mode exchange between the two modes through an avoided crossing that occurs when we move from Newtonian dynamics to relativistic dynamics.

Needless to say, nonradial oscillations of neutron stars should be treated within the framework of general relativity because of their strong gravity. At present, however, we cannot investigate the pulsations of neutron stars in all the aspects by using fully general relativistic formalism because of the complexity general relativity brings about.222One of the problems that appears in general relativistic formalism is the existence of gravitational waves and a necessity to deal with a complex eigenvalue problem corresponding to quasi-normal modes. Some authors have treated this problem (e.g., Lindblom & Detweiler 1983; Leins, Nollert & Soffel 1993; Andersson, Kokkotas, & Schutz 1995). One of the simplest ways to study the oscillations of a neutron star is to treat the problems within the framework of Newtonian dynamics, that is, to solve the Newtonian oscillation equations for neutron star models constructed with the Newtonian hydrostatic equations. However, it is obvious that this Newtonian treatment cannot be fully approved for neutron stars since the Newtonian equilibrium models do not give us correct mass and radius for the stars. To reduce the distance between fully relativistic and Newtonian calculations of nonradial oscillations of the stars, we may use relativistic equilibrium models and solve the Newtonian oscillation equations, which was the strategy taken by McDermott, van Horn, & Hansen (1988). In this paper, as an extension of the study by McDermott, van Horn, & Hansen (1988), we solved relativistic oscillation equations derived in the relativistic Cowling approximation for relativistic neutron star models with a solid crust. Our calculation confirmed that the classification scheme McDermott, van Horn, & Hansen (1988) employed for the oscillation modes of neutron stars with a solid crust is valid even if we integrate relativistic oscillation equations.

###### Acknowledgements.

S.Y. acknowledges financial support from the Portuguese FCT through a Sapiens project, number 36280/99.


format: text

text:
# Signatures of interband scattering in spectroscopic experiments on 𝑀𝑔𝐵₂

## Abstract

Within a two-band, strong coupling model we analyze $`SIS`$-tunneling conductance, Raman scattering and optical conductivity measured on $`MgB_2`$ samples. We find that features observed in tunneling and Raman spectroscopy at intermediate energies $`\omega =\mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi 8÷10meV`$ can be consistently explained when incoherent scattering between $`\sigma `$ and $`\pi `$ band is considered. We calculate the optical conductivity and find that strong coupling effects are important for explaining the presently available data, and predict the location and magnitude of interband scattering features for this spectroscopic probe as well.

*Introduction.* The discovery of superconductivity in $`MgB_2`$ with a critical temperature $`T_C=39K`$ Nagamatsu  has initiated intense research efforts on elucidating the properties of this surprisingly simple compound. The main reason behind the relatively high transition temperature is still unknown. However, there is a growing consensus that strong electron-phonon interaction Kortus  is responsible for superconductivity in $`MgB_2`$. The large isotope effect measured by Bud’ko et al. Budko  for boron ($`\alpha =0.26\pm 0.03`$) clearly shows that phonons associated with boron vibrations play a significant role in superconductivity in this compound. Tunneling and point contact spectroscopy experiments Karapetrov ; Rubio ; Szabo ; Sharoni ; Kohen ; Schmidt  reveal a distribution of energy gap $`2\mathrm{\Delta }(0)/T_C`$ between 1.1 and 4.5. Recent experiments of Giubileo et al. Giubileo  pointed the presence of two gaps. Evidence of two gaps were found also in the specific heat Bouquet  and recently in Raman scattering experiments Chen .

According to Liu et al. Liu  the four Fermi surface sheets in $`MgB_2`$ can be grouped into quasi two-dimensional $`\sigma `$ bands and three dimensional $`\pi `$ bands, and the normal and superconducting properties of $`MgB_2`$ can be described by an effective two-band model. Within this model they calculated the coupling constants and energy gap ratios in the weak coupling regime. Later Golubov et al.Golubov  provided the two-band decomposition of the superconducting Eliashberg functions $`\alpha ^2F_{ij}(\omega )`$ (where $`i`$ and $`j`$ denote $`\sigma `$ and $`\pi `$ bands), which describe the electron-phonon coupling in $`MgB_2`$ as function of frequency. While electron-phonon interaction was clearly investigated in large spectrum of approximations, the role of scattering between the bands was not analyzed in detail.

We turn our attention to the theory of spectroscopic properties of $`MgB_2`$ because, similarly to the situation in cuprate superconductors, the most severe constrains on theoretical models are imposed by spectroscopic measurements. In fact, the aim of this paper is to point out that a simultaneous explanation of the main features in $`SIS`$ tunnel spectroscopySchmidt  and Raman scattering Chen  can only be done if one allows for the existence of incoherent (momentum non-conserving) single-particle scattering events between. We have also analyzed recent optical conductivity results Kaidl . We find, as shown in detail below, that agreement between theory and experiment is substantially improved over the weak-coupling analysis by considering strong coupling effects alone. We have also calculated the effect of incoherent interband scattering on such measurements. While the presently available data does not allow us to firmly conclude, that there is evidence for incoherent interband scattering in the optical conductivity data as well, it does not contradict the presence of the process either. We nevertheless present our predictions for optical conductivity in the presence of interband scattering, since we believe that our predictions for this probe will be useful in analyzing future, more accurate measurements.

*Theoretical Model.* We will adopt a two-band Eliashberg scheme to calculate the various spectroscopic quantities. The weak coupling version of the model for two-band superconductivity was originally proposed by Suhl, Matthias and WalkerSuhl . However, due to the lack of a genuine two-band superconductor at that time, interest in this model gradually decreased. More recently, the Suhl-Matthias-Walker model was revived and applied to cuprates by KresinKresin . We have used a similar modelMoca  to calculate the penetration depth in $`MgB_2`$, and we are now applying this framework to derive the spectroscopic properties as well. The starting point of our analysis is the two-band Hamiltonian $`=_0+_I+_{\mathrm{i}\mathrm{b}}`$. Here $`_0`$ is describing the free electrons,

$$_0=\underset{𝐤,\sigma }{}\epsilon _𝐤c_{𝐤\sigma }^+c_{𝐤\sigma }+\underset{𝐤,\sigma }{}\xi _𝐤d_{𝐤\sigma }^+d_{𝐤\sigma }+\underset{𝐪}{}\omega _𝐪b_𝐪^+b_𝐪.$$
(1)
In the above expression $`c_{𝐤\sigma }^+`$ $`\left(c_{𝐤\sigma }\right)`$ creates (annihilates) electrons in $`\sigma `$ band and $`d_{𝐤\sigma }^+`$ $`\left(d_{𝐤\sigma }\right)`$ are similar operators acting in $`\pi `$ band. The last term describes the non-interacting phonons with the dispersion given by $`\omega _𝐪`$. The second term in the Hamiltonian, $`_I`$ corresponds to the interaction of the conduction electrons with the phonons,

$`_I`$ $`=`$ $`{\displaystyle \underset{𝐤,𝐪,\sigma }{}}V_{\sigma \sigma }(𝐤,𝐪)c_{𝐤+𝐪\sigma }^+c_{𝐤\sigma }(b_𝐪+b_𝐪^+)`$ (3)
$`+{\displaystyle \underset{𝐤,𝐪,\sigma }{}}V_{\pi \pi }(𝐤,𝐪)d_{𝐤+𝐪\sigma }^+d_{𝐤\sigma }(b_𝐪+b_𝐪^+)`$
$`+{\displaystyle \underset{𝐤,𝐪,\sigma }{}}V_{\sigma \pi }(𝐤,𝐪)[c_{𝐤+𝐪\sigma }^+d_{𝐤\sigma }(b_𝐪+b_𝐪^+)+h.c.].`$
Finally we introduce a new term $`_{\mathrm{i}\mathrm{b}}`$ to account for the  *incoherent (momentum non-conserving) interband scattering*

$$_{\mathrm{i}\mathrm{b}}=\underset{𝐤,\sigma }{}(T_{\mathrm{𝐤𝐪}}c_{𝐤\sigma }^+d_{𝐪\sigma }+h.c.)$$
(4)
. The effect of this last term is the main focus of our study, and we will show below, that it plays an important role in explaining the observed features in spectroscopic measurements.

Our approach is based on the self-consistent solution of the Eliashberg equations corresponding to the Hamiltonian (1-4). Self-consistency is reached on the imaginary axis, and Padé approximation Vidberg  is then used to analytically continue to the real axis. This procedure gives accurately the frequency dependence of the gap functions $`\mathrm{\Delta }_\sigma (\omega )`$ and $`\mathrm{\Delta }_\pi (\omega )`$ and the renormalization factors $`Z_\sigma (\omega )`$ and $`Z_\pi (\omega )`$. The normal self-energy contribution to the $`i`$ band due to the direct scattering from and in the $`j`$ band can be written as $`1/(2\tau _{ij})g_j(i\omega _n)`$ where $`g_i(i\omega _n)=\omega _n/\sqrt{\omega _n^2+\mathrm{\Delta }_i^2(\omega _n)}`$ and the anomalous self-energy can be written as $`1/(2\tau _{ij})f_j(i\omega _n)`$ where $`f_i(i\omega _n)=\mathrm{\Delta }_i/\sqrt{\omega _n^2+\mathrm{\Delta }_i^2(\omega _n)}`$. The structure in $`\mathrm{\Delta }(\omega )`$ usually gives rise to observable variations in the density of states at a frequency comparable with the sum of phonon energy and gap value. However, in $`MgB_2`$ the phonon mode most relevant to superconductivity is $`E_{2g}67meV`$, and observation of density of states features at such high energies are difficult to be observed experimentally. For the present calculation we model the phonon density with a lorentzian: $`F(\omega )=A/((\omega \omega _1)^2+\omega _2^2)A/(\omega _3^2+\omega _2^2)`$ for $`\left|\omega \omega _1\right|<\omega _3`$ and zero otherwise. Here the frequency parameters are: $`\omega _1=67meV`$, $`\omega _2=3.5meV`$, and $`\omega _3=30meV`$. With these parameters, and taking a cut-off frequency at $`\omega _c=200meV`$, we fix $`A=1.31`$ by imposing $`\underset{0}{\overset{\omega _c}{}}F(\omega )𝑑\omega =1`$. The best fit to the experimental temperature dependence of the gap functionBuzea  and with the observed $`T_C=39K`$ was found for $`\alpha _{\sigma \sigma }^2=17.5`$, $`\alpha _{\pi \pi }^2=7.0`$, $`\alpha _{\sigma \pi }^2=1.0`$, $`\alpha _{\pi \sigma }^2=2.0`$. The effect of Coulomb pseudopotential have been considered through the renormalized parameters $`\mu _{ii}^{}=0.1`$ and $`\mu _{ij}^{}=0`$ for $`ij`$. The effect of cross-band interaction can be analyzed in terms of $`\mathrm{\Gamma }_{ij}=\tau _{ij}^1`$ which directly describe the scattering rate. This is similar with solving the Eliashberg equation for an alloy where similar terms appear due to scattering on impurities Carbotte . Throughout of our calculations we have considered that $`\mathrm{\Gamma }_{ii}=0`$ (no intra-band scattering) and that only $`\mathrm{\Gamma }_{ij}0`$ for $`ii`$.

*Break Junction Tunneling.* We start the comparison to experiment by analyzing the break junction $`SIS`$ tunneling experimental data of Schmidt et al. Schmidt . A careful analysis of the data reveals the presence a dip at a frequency $`\omega \mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi 10meV`$ and no measurable effects are seen near the frequency corresponding to the larger gap. This suggests that, due to its $`3D`$ nature, only the $`\pi `$ band gives a contribution to the tunneling current in break junctions, while the $`\sigma `$ band, due to its quasi $`2D`$ nature, has little overall contribution. The electron-phonon interaction alone does not give a satisfactory explanation for the above mentioned feature, and a new term is needed. This term describes incoherent single particle tunneling between bands, and is formally similar to the proximity term first proposed by McMillan McMillan . Within this extended Suhl-Matthias-Walker model, the superconductivity in $`\pi `$ band can be thought of as being induced by the $`\sigma `$-band. Scattering between bands diminishes the larger gap while the smaller one is enhanced. More importantly, however, self-energy contributions in lowest order in interband scattering will give rise to intermediate states of energies $`\mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi 810\mathrm{meV}`$, precisely at the voltage where the dip in the tunnel conductance is seen.

In order to study the effect of the direct inter-band interaction we have first calculated the density of states for each band by using the well known relation $`N_S^{(i)}(\omega )=N_N^{(i)}(\omega )Re\left(\left(\omega \mathrm{\Gamma }_T\right)/\sqrt{(\omega \mathrm{\Gamma }_T)^2\mathrm{\Delta }_i^2(\omega )}\right)`$, where $`N_S(\omega )`$ and $`N_N(\omega )`$ are the energy dependent density of states of the material in the superconducting and normal state. We have calculated first the superconducting density of states and then using these results we investigate the $`SIS`$ conductance by calculating first the tunneling current by a direct convolution:

$$I_{i,j}(eV)\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left[n_F(\omega )n_F(\omega +eV)\right]N_S^{(i)}(\omega )N_S^{(j)}(\omega +eV)$$
(5)
and then taking the derivative as function of bias. In Eq. (5) $`n_F(\omega )`$ is the Fermi-Dirac function and $`N_S^{(i)}(\omega )`$ with $`i=\pi ,\sigma `$ is the density of states corresponding to one of the bands. The conductance can be found by differentiating the tunneling current as function of bias and a direct quantitative comparison with the experimental data is possible. Fig. 1 represents the results obtained for the conductance for the $`\pi `$-band together with the experimental data of Schmidt et al. Schmidt . The interband self-energy leads to modification in the structure of the gap function at energies $`\omega =\mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi `$. Increasing the scattering rates between the bands lead to $`(i)`$ an enhance of the smaller gap corresponding to the $`\pi `$ band and to a decreasing of the $`\sigma `$-band gap; $`(ii)`$ increasing of the dip feature at the energy $`\mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi `$ in the $`SIS`$ conductance spectrum. The last effect can be verified experimentally by comparing samples of different impurity concentrations. The presence of impurities will increase such interband scattering processes, and an interesting and somewhat counter-intuitive increase in gap size is therefore expected with the increase of impurity concentration. This would be a qualitative consistency check on our theory.

*Raman Spectroscopy.* This is the second spectroscopic experiment that we analyze in the same theoretical framework and find an excellent agreement between the experiment and our model. Early Raman experiments performed on $`MgB_2`$ were somewhat inconsistent due to differences in sample quality: measurements on poly-crystaline samplesChen  reveal the presence of two gaps while single crystal measurements restricted to only $`ab`$-plane polarized spectra show only the presence of a one superconducting gap Quilty . Recently a more careful analysis resolved this inconsistency and the agreement is that the Raman spectrum does reveal the presence of two bandsQuilty1 .

We have considered the effect of scattering between the bands when calculating the Raman susceptibility Jiang  $`X(𝐪0,\omega )`$. The Raman spectrum was calculated separately for each band and then the spectrum was calculated as a weighted sum of contributions from each band. The inset of Fig.2 presents the normalized contributions coming from each band. There is a complete depletion of the Raman spectrum below $`2\mathrm{\Delta }_i`$ for each band and the spectrum has a vertical slope at $`2\mathrm{\Delta }_i`$. The weighting coefficients are 14.61 for $`\pi `$-band and 4.39 for the $`\sigma `$-band leading to a ratio of 3.32 for the contributions of each band (in contrast with the 2.97 ratio obtained in the Suhl-Matthias-Walker model by Chen et al.Chen . In our calculations the Raman vertex was calculated for simplicity in the limit $`\gamma (𝐩)=1`$. Even in this limit the effect of scattering between different bands is important, and leads to significant modifications of the spectrum. Our results for the Raman spectrum are presented in Fig.2 and compared to the experimental data of Chen et al. Chen . The dashed line is the weak coupling $`BCS`$ two-band model fit proposed by the same authors. The significant difference between our results and theirs is the presence of the hump at an energy $`\omega =\mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi `$, a feature that theoretically can be obtained by considering only the scattering between the bands, and does not appear in the simple Suhl-Matthias-Walker model. In our calculations the best fit of the experimental data was obtained when $`\mathrm{\Gamma }_{\sigma ,\pi }=0.5`$ and $`\mathrm{\Gamma }_{\pi ,\sigma }=3.0`$. The other parameters used in the calculations are r $`\alpha _{\sigma \sigma }^2=17.5`$, $`\alpha _{\pi \pi }^2=5.5`$, $`\alpha _{\sigma \pi }^2=1.0`$, $`\alpha _{\pi \sigma }^2=2.0`$ and the same values for the Coulomb pseudopotential as in the case of analysis of the tunneling experiment.

*Optical Conductivity.* The last spectroscopic experiment that we focus on is the optical conductivity, where the importance of the scattering between the bands was analyzed in the framework of Mattis-Bardeen theory Mattis . A quantitative comparison with the experimental data on optical conductivity was possible in this case as well, as shown in Fig. 3. One immediate conclusion is already apparent: the overall fit to the experimental data is much better in the Eliashberg theory than in the constant gap model.

Due to the small critical temperature obtained for the measured sample Kaidl  $`(T_C=30K)`$ we had to modify the phonon coupling constants in order to obtained the desired critical temperature. Therefore we consider $`\alpha _{\sigma \sigma }^2=15.5`$, $`\alpha _{\pi \pi }^2=5.0`$, $`\alpha _{\sigma \pi }^2=1.0`$, $`\alpha _{\pi \sigma }^2=2.0`$, and for the scattering rates the best fit with the experimental data was obtained for $`\mathrm{\Gamma }_{\sigma ,\pi }=0.5`$ and $`\mathrm{\Gamma }_{\pi ,\sigma }=2.5`$. The small feature appearing in the real part of the optical conductivity at $`\omega =\mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi 8meV`$ (in this case this energy is smaller due to the smaller critical temperature $`\mathrm{\Delta }_\pi (T=0K)=2.51\mathrm{meV}`$, $`\mathrm{\Delta }_\sigma (T=0K)=5.48\mathrm{meV}`$), is due to the modifications in the frequency dependence of the gap functions induced by the effect of inter-band scattering in the electronic self-energy. The small feature in the data near this frequency appears to be an extrinsic effect Carnahan , and is unlikely to correspond to interband scattering. However, we believe that future experiments will be able to detect the interband structure we predict here.

*Conclusions.* We have investigated the role of the incoherent inter-band scattering in the framework of two-band strong coupling model. We found results consistent with three experimental measurements. The main results of our work can therefore be summarized as follows: $`(i)`$ superconductivity in $`MgB_2`$ can be well explained in the framework of the two-band model; $`(ii)`$ The features observed in tunneling spectroscopy, and Raman scattering at energies $`\omega =\mathrm{\Delta }_\sigma +\mathrm{\Delta }_\pi 8÷10meV`$ can be consistently explained by considering the incoherent scattering between the bands. While there is no clear feature observed in optical conductivity at this frequency, we predict that such features will eventually be observed, $`(iii)`$ Overall, the fit to the experimental data is much better if a strong coupling theory is used, rather than a two-band, weak coupling, constant gap model. Schmidt .

*Acknowledgments.* We are especially indebted to Drs. H. Schmidt and M. Carnahan for supplying us with tunneling and optical spectroscopy data, respectively. We wish to thank them as well as Drs. G. Arnold, G. W. Crabtree, K. Gray, M. Iavarone, W.-K. Kwok, G. Karapetrov, A. Koshelev, J. Orenstein, and J. Zasadzinski for useful discussions.


text:
# Feshbach resonances in rubidium 87: Precision measurement and analysis

## Abstract

More than 40 Feshbach resonances in rubidium 87 are observed in the magnetic field range between 0.5 and 1260 gauss for various spin mixtures in the lower hyperfine ground state. The Feshbach resonances are observed by monitoring the atom loss, and their positions are determined with an accuracy of 30 mG. In a detailed analysis, the resonances are identified and an improved set of model parameters for the rubidium interatomic potential is deduced. The elastic width of the broadest resonance at 1007 G is predicted to be significantly larger than the magnetic field resolution of the apparatus. This demonstrates the potential for applications based on tuning the scattering length.

A Feshbach resonance is an exciting tool for controlling the atom-atom interaction in ultracold atomic gases. The elastic $`s`$-wave scattering length $`a`$ can be tuned over orders of magnitude simply by applying a magnetic field. Feshbach resonances have been observed in various alkali atoms Inouye et al. (1998); Courteille et al. (1998); Roberts et al. (1998); Vuletic et al. (1999); Loftus et al. (2002); Khaykovich et al. (2002); Strecker et al. (2002). They have been used to induce a controlled implosion of a Bose-Einstein condensate (BEC) Donley et al. (2001), to create a coherent superposition of an atomic BEC and a molecular state Donley et al. (2002), and to realize a bright soliton in a BEC Khaykovich et al. (2002); Strecker et al. (2002). Future applications could include experiments with the Mott-insulator phase transition Greiner et al. (2002), a Tonks gas Girardeau and Wright (2002), effects beyond the mean-field theory Dalfovo et al. (1999), and the creation of a molecular BEC kokkelmans:etal . In addition, the binding energies of ro-vibrational molecular states close to the dissociation threshold can be determined with high accuracy from the position of Feshbach resonances leading to a precise knowledge of the interatomic potential chin:etal . A Feshbach resonance showing up in elastic collisions is often accompanied by strong changes in the inelastic collision properties Stenger et al. (1999); Roberts et al. (2000). This offers a strategy for searching for new Feshbach resonances by monitoring the resulting atom loss.

Surprisingly, no Feshbach resonance has been observed in 87Rb, which is the isotope used in most of today’s BEC experiments. In 1995, a search with atoms in the $`|f,m_f=|1,1`$ state was carried out, but no resonances were found Newbury et al. (1995). Meanwhile various experiments Myatt et al. (1997); Wynar et al. (2000); Seto et al. (2000); Roberts et al. (2001) greatly improved the knowledge of the Rb interatomic potential. Recent models based on this Vogels et al. (1997); van Kempen et al. (2002) are consistent with the observations of Ref. Newbury et al. (1995). But these models predict four Feshbach resonances in the $`|1,1`$ state.

This letter reports the observation of more than 40 Feshbach resonances in 87Rb with most (but not all) atoms prepared in the $`|1,1`$ ground state. The theoretical model is extended to various spin mixtures and to bound states with rotational quantum number $`l3`$. Thus, all except one of the resonances can be clearly identified. Moreover, the measured position of one Feshbach resonance is used for an improved fit of the model parameters. The relative deviation between the predicted and observed positions is $`1.6\times 10^3`$ (rms). The observed loss might be due to two- or three-body inelastic collisions, but should be purely three-body for those resonances which involve only the $`|1,1`$ state. The broadest resonance at 1007 G offers the possibility to tune the scattering length and investigate its loss mechanism.

The experiment is performed with a new set-up similar to our previous one Ernst et al. (1998), but with all relevant components significantly improved. In particular, atoms are captured in a vapor-cell magneto-optical trap (MOT) at a loading rate of $`7\times 10^{10}\mathrm{s}^1`$. The atoms are transferred to a second MOT, in which $`6\times 10^9`$ atoms are accumulated by multiple transfer within 2 s. The atoms are then loaded into a Ioffe-Pritchard magnetic trap with a lifetime of 170 s. The measured bias-field drift of less than 1 mG/h illustrates the excellent stability of the magnetic trap. After 26 s of evaporative cooling, a BEC with up to $`3.6\times 10^6`$ atoms is formed in the $`|1,1`$ state.

The $`|1,1`$ state, in which the resonances are predicted, cannot be held in a magnetic trap. Therefore the atoms are now loaded into an optical dipole trap made of a single beam from an Yb:YAG laser at a wavelength of 1030 nm. A laser power of 45 mW is focused to a waist of 15 $`\mu `$m, resulting in measured trap frequencies of 930 Hz and 11 Hz, and an estimated trap depth of $`k_B\times 20\mu `$K.

Once the atoms are in the optical trap, a radio-frequency field is used to transfer the atoms to the desired $`|1,1`$ state. With a Stern-Gerlach method, the fraction of atoms that end up in the $`|1,1`$ state is determined to be roughly 90%. Almost all other atoms are in the $`|1,0`$ state. Next, a homogeneous magnetic field is applied in order to observe a Feshbach resonance. The field is created using the compensation coils of the magnetic trap which are in near-perfect Helmholtz configuration. Up to 1760 A of current are stabilized with a home-built servo to a few ppm. The magnetic field is held at a fixed value for typically 50 ms and then quickly switched off. Although many Feshbach resonances are rapidly crossed when the magnetic field is turned on or off, no significant loss of atoms is observed from these rapid crossings. After switching off the magnetic field, the atoms are released from the optical trap, and 14 ms later an absorption image of the expanded cloud is taken.

The search for Feshbach resonances was typically performed with a purely thermal cloud of $`4\times 10^6`$ atoms in the optical trap at a temperature of 2 $`\mu `$K, corresponding to a peak density of $`2\times 10^{14}`$ cm-3. The magnetic-field range between 0.5 and 1260 G was scanned and 43 Feshbach resonances were found. The magnetic-field values of these resonances $`B_{exp}`$ are listed in Tabs. 1 and 2. The fields were calibrated with 30 mG precision using microwave spectroscopy in the vicinity of each resonance.

In order to analyze these results theoretically, both continuum and bound-state calculations are carried out using an accurate description of the atomic interaction between two Rb atoms van Kempen et al. (2002). An introduction to the theory of Feshbach resonances can be found in Refs. Heinzen (1999); Burke (1999). The strong central part of the interaction, consisting of singlet and triplet potentials, conserves the orbital quantum numbers $`l,m_l`$ and the total-spin magnetic quantum number $`m_F`$ separately. (In bound-state spectroscopy, $`l`$ is often called $`N`$.) In the ultracold regime, with this part of the interaction included, Feshbach resonances are expected to occur only when an $`l=0`$ bound state with a certain value of $`m_F`$ crosses the dissociation threshold of an entrance channel with the same $`m_F`$. Figure 1 shows these crossings in the $`E`$-$`B`$ plane as solid dots along various $`(m_{f1}^{in},m_{f2}^{in})`$ thresholds. The resonance fields obtained from these bound-state calculations are listed in Tabs. 1 and 2 ($`l`$ = 0 quantum number). The same field values also follow from the $`B`$-dependence of the elastic $`S`$-matrix element, as obtained in a continuum calculation.

In Fig. 1, the bound states are labelled $`(f_1,f_2)v^{}`$, with the vibrational quantum number $`v^{}`$ = -1, -2, -3, … counting from the corresponding $`(f_1,f_2)`$ threshold. At $`B=0`$, the exchange interaction couples the atomic spins $`f_1`$, $`f_2`$ to a total molecular spin $`F`$ and causes a splitting between states with the same $`(f_1,f_2)v^{}`$ but different $`F`$. In the presence of a strong external magnetic field, however, $`m_{f1}`$ and $`m_{f2}`$ become good quantum numbers instead of $`F`$, while $`m_F`$ is always a good quantum number. What constitutes a strong field in this sense depends on the size of the $`F`$-splitting at $`B=0`$. For small $`|v^{}|`$, the $`F`$-splitting is hardly visible in Fig. 1.

The much weaker spin-spin interaction $`V_{ss}`$ consists of the magnetic dipole-dipole interaction of the valence electrons together with a second-order spin-orbit term Mizushima (1975). Due to its tensor form it breaks the spatial spherical symmetry and allows a redistribution between the angular momenta of the spin and spatial degrees of freedom so that the sum $`m_l+m_F`$ is the only conserved quantum number. $`V_{ss}`$ admixes an $`l`$ = 0 component in otherwise pure $`l`$ = 2 bound states and therefore induces additional $`s`$-wave Feshbach resonances in the ultracold regime. Resonances of this type have previously been observed in Cs chin:etal , but up to now never in 85Rb or 87Rb. The resulting resonance fields are again listed in the tables ($`l`$ = 2 resonances). For mixed species resonances, such as $`|1,0|1,1`$, the $`l=1`$ partial wave can be populated in the entrance channel, thus opening up the possibility to observe resonances due to $`l=1`$ or $`l=3`$ bound states.

The set of potential parameters used in this calculation is obtained as follows. The field value 911.74 G of one resonance is added to the set of eight experimental data already included in the analysis of Ref. van Kempen et al. (2002). This particular resonance was chosen, because the corresponding bound state is a pure triplet state and the previous experimental constraints on the singlet potential Seto et al. (2000) are much stronger than those on the triplet potential. Then, a least-squares analysis is applied with the column A parameters of Tab. I in Ref. van Kempen et al. (2002) as starting values. This leads to an adjusted set of parameter values differing from the starting values by less than $`1\sigma `$; in atomic units: $`C_6=4.707\times 10^3`$, $`C_8=5.73\times 10^5`$, $`C_{10}=7.665\times 10^7`$ (Ref. Marinescu et al. (1994)), $`J=0.486\times 10^2`$, $`a_S(^{87}\mathrm{Rb})=+90.6`$, $`a_T(^{87}\mathrm{Rb})=+98.96`$, $`v_{DS}(^{87}\mathrm{Rb})=0.454`$, and $`v_{DT}(^{87}\mathrm{Rb})=0.4215`$. Coupled-channel calculations Heinzen (1999); Burke (1999) based on these fine-tuned parameters then yield theoretical resonance fields $`B_{th}`$. The relative deviation of all these positions from the observations is $`1.6\times 10^3`$ (rms). Compared to the four resonance positions predicted in van Kempen et al. (2002), this is an improvement of a factor of 6. Since only one resonance position was included in the fit, the excellent agreement with all other positions demonstrates the accuracy of the model. Note that the observed position of maximum loss might deviate from the pole of the elastic scattering length Fedichev et al. (1996). However, such deviations are typically not larger than the elastic widths of the resonances; and they are small as discussed below.

Some resonances are so close together that they cannot be identified merely from their positions. In these cases, a Stern-Gerlach method was used to experimentally determine which $`m_f`$ states incurred the strongest atom loss. Thus the entrance channel could be identified, in particular when additionally varying the initial spin mixture. The entrance channel of the non-identified resonance at 1236.73 G was also determined with this method.

An interesting property of a Feshbach resonance – besides its position – is its strength. The strength of the elastic resonance is proportional to the field width $`\mathrm{\Delta }`$ over which the scattering length has opposite sign Burke (1999), as listed in Tab. 1. The strength of the inelastic scattering properties is quantified by the loss rate. The dominant loss for the resonances in Tab. 1 arises from inelastic three-body collisions. This is because for the parameters of the experiment, single-body loss is negligible; and since the $`|1,1`$ entrance channel is the absolute ground state of atomic 87Rb, inelastic two-body collisions cannot occur Stenger et al. (1999). (This is not the case for the resonances in Tab. 2.) The three-body loss is characterized by the coefficient $`K_3`$ in the rate equation $`\dot{N}=K_3n^2N`$, where $`N`$ is the atom number and $`n`$ the density. The depth $`d`$ listed in Tab. 1 is a non-linear, yet monotonic function of $`K_3`$. In Tab. 1, one finds a clear trend that stronger resonances with larger $`\mathrm{\Delta }`$ cause faster loss, i.e. larger $`d`$. This is plausible, although the theory of three-body losses is not yet fully understood, especially in the vicinity of Feshbach resonances Burke (1999); Fedichev et al. (1996).

Interestingly, the widths $`\mathrm{\Delta }`$ of the $`l=0`$ resonances in Tab. 1 are much smaller than for other alkali atoms Stenger et al. (1999); Vuletic et al. (1999); Roberts et al. (2001); Loftus et al. (2002). This is due to the approximate phase equality of the waves reflected from the short-range singlet and triplet potentials, that is also responsible for other remarkable 87Rb phenomena: the coexistence of condensates in different hyperfine states Myatt et al. (1997) and the smallness of the fountain-clock frequency shift van Kempen et al. (2002).

All resonances display a nearly symmetric loss feature. A Gaussian is fit to the atom loss at a 50-ms hold time in order to extract the depth and the width of the atom loss. The Gaussian was chosen for convenience and because it fits well to the data. The obtained rms-widths are between 20 mG and 100 mG for almost all resonances. For small $`d`$, this width is identical to the width of the resonance in $`K_3`$ (if two-body loss is absent). For large $`d`$, however, the non-linear dependence of $`d`$ on $`K_3`$ broadens the observed width at a given hold time, somewhat similar to saturation broadening of spectral lines. In addition, the finite temperature of the cloud gives rise to a broadening of typically 20 mG and leads to a shift of similar size. The observed rms-width of 24(4) mG at 965.96 G sets an experimental upper limit to thermal and technical broadening.

The broadest resonance is centered at 1007 G (see Fig. 2). The theoretical prediction for its elastic width, $`\mathrm{\Delta }=170(30)`$ mG, is large compared to the upper bound on the experimental broadening of 24(4) mG (see above). This demonstrates the potential for a controlled variation of the scattering length with the present set-up.

For various field values, the decay of the atom number was also measured as a function of time. If one assumes that the loss can be described by a rate equation, three-body loss dominates (see above). Values of $`K_3`$ determined from a fit are shown in Fig. 2. Note that the value $`K_3=3.2(1.6)\times 10^{29}`$ cm6s-1 obtained away from the resonance is consistent with the measured value $`K_3=4.3(1.8)\times 10^{29}`$ cm6s-1 for the $`|1,1`$ state Burt et al. (1997). The absence of two-body loss makes this resonance an ideal candidate for testing theories of three-body loss Fedichev et al. (1996). The values of $`K_3`$ in Fig. 2 are much smaller than those for 85Rb Roberts et al. (2000), where exciting experiments have been performed Donley et al. (2001, 2002). Three-body loss will therefore not be a substantial problem for applications.

All data obtained here fit well to a three-body decay, but the initial atom number extrapolated from the fit is a factor of up to 2 lower at the 1007 G resonance as compared to the off-resonance value. Since no data were taken for hold times shorter than 10 ms, this suggests the existence of an additional loss mechanism acting on much faster timescales. This could be molecule formation and dissociation during the first 10 ms of hold time Claussen et al. (2002). Note, however, that here loss from a purely thermal cloud is observed, while loss from a BEC is discussed in Ref. Claussen et al. (2002).

To summarize, more than 40 Feshbach resonances have been observed in 87Rb. All except one were identified by theory. After including one of the observed resonances in fitting the model parameters, theory and experiment are in excellent agreement. The magnetic field control is so accurate that it should be possible to resolve changes in the elastic scattering length on the broadest resonance at 1007 G.

The work in Garching is supported by the European-Union network “Cold Quantum Gases”, the DFG-Forschergruppe “Quantengase”, and GIF. The work at Eindhoven is part of the research program of the Stichting FOM, which is financially supported by NWO.