text:
# Monte Carlo studies of the one-dimensional Ising spin glass with power-law interactions

## I Introduction

Two main theories of the spin-glass state have been proposed: the “droplet picture,”Fisher and Huse (1986, 1987, 1988); Bray and Moore (1986) and the replica symmetry-breaking (RSB) picture.Parisi (1979, 1980, 1983); Mézard et al. (1987) The droplet picture states that the minimum-energy excitation of linear dimension $`l`$ containing a given spin (a droplet) has an energy proportional to $`l^\theta `$, where $`\theta `$ is positive (assuming that the transition temperature $`T_c`$ is finite). It follows that in the thermodynamic limit, excitations that flip a finite fraction of the spins cost an infinite energy. In addition, the droplet picture states that these excitations are fractal with a fractal dimension $`d_s<d`$, where $`d`$ is the space dimension. By contrast, RSB follows Parisi’s exact solution of the Sherrington-Kirkpatrick (SK) infinite-range model and predicts that excitations involving a finite fraction of the spins cost a finite energy in the thermodynamic limit, and are space filling,Marinari and Parisi (2000) i.e., $`d_s=d`$.

There have been several numerical attemptsKrzakala and Martin (2000); Palassini and Young (2000); Marinari and Parisi (2000); Katzgraber et al. (2001); Katzgraber and Young (2001, 2002); Houdayer et al. (2000); Houdayer and Martin (2000) to better understand the nature of the spin-glass state for short-range spin glasses. The data, which are on small system sizes, are consistent with a picture in which the “stiffness exponent” for dropletlike excitations is different from that for domain-wall excitations, although it is still unclear if the difference is due to the limited range of system sizes,Moore (2002) or if it persists in the thermodynamic limit. Recently Krzakala and Martin,Krzakala and Martin (2000) as well as Palassini and Young,Palassini and Young (2000) find evidence for an intermediate picture (called “TNT” which stands for trivial-nontrivial) in which only a finite amount of energy is needed to generate system-size excitations in the thermodynamic limit, but their surface has a nontrivial fractal dimension less than $`d`$.

In addition to finite-$`T`$ Monte Carlo simulations, there has been a substantial amount of numerical work studying the energy of domain walls at $`T=0`$.Bray and Moore (1984); McMillan (1984a, b); Rieger et al. (1996); Hartmann (1999); Palassini and Young (1999); Hartmann and Young (2001); Carter et al. (2002) The energy of the domain wall is found to vary as $`L^\theta `$ where $`L`$ is the system size and $`\theta `$ is positive (for systems with $`T_c>0`$). The droplet theory assumes that the stiffness exponent for domain walls $`\theta _{\mathrm{domain}\mathrm{wall}}`$ is the same as the stiffness exponent for dropletlike excitations $`\theta _{\mathrm{droplet}}`$. In addition, $`\theta _{\mathrm{domain}\mathrm{wall}}>0`$. This is in contrast to the RSB and TNT pictures where the exponents are expected to be different and $`\theta _{\mathrm{droplet}}=0`$. To simplify the notation,com  from now on we will denote $`\theta _{\mathrm{domain}\mathrm{wall}}`$ by $`\theta `$ and $`\theta _{\mathrm{droplet}}`$ by $`\theta ^{}`$.

In this work we consider the question of whether there are one or two stiffness exponents for the case of Ising spin-glass models in one-dimension in which the interactions fall off with a power $`\sigma `$ of the distance. Depending on $`\sigma `$, the system can have either a finite transition temperature $`T_c`$ (with the long-range part of the interaction relevant), or $`T_c=0`$, in which case it can be either in the short-range or long-range universality class.Bray et al. (1986); Fisher and Huse (1988) The advantages of this system are the following: (i) a large range of system sizes $`L`$ can be studied, (ii) there are analytical predictions for the stiffness exponent within the droplet picture, and (iii) within a single model, simply by changing a parameter, one can cover the full range of behavior, from the absence of a spin-glass phase to the infinite-range SK model. We perform both Monte Carlo simulations at low $`T`$ (which indirectly gives $`\theta ^{}`$), for a value of $`\sigma `$ that has a finite $`T_c`$, and ground-state calculations for a range of values of $`\sigma `$ corresponding to all three regimes: $`T_c>0,T_c=0`$ (long range), and $`T_c=0`$ (short range). We use the parallel temperingHukushima and Nemoto (1996); Marinari (1998) approach for both the finite-$`T`$ and $`T=0`$ studies.

In Secs. II, III, and IV, we introduce the model, observables, and details of the Monte Carlo technique, respectively. In addition we describe the phase diagram of the model in detail. The critical region of the disordered long-range Ising spin glass is analyzed in Sec. V in order to determine the critical temperature of the system and to ensure that in Sec. VI, where we study the system at finite $`T`$, the temperature is sufficiently far below $`T_c`$ such that critical-point fluctuations do not affect the data. The zero-temperature results are described in Sec. VII and our conclusions are summarized in Sec. VIII.

## II Model

The Hamiltonian for the one-dimensional long-range Ising spin glass with power-law interactions is given by

$$=\underset{i,j}{}J_{ij}S_iS_j,$$
(1)
where the sites $`i`$ lie on a ring of length $`L`$, as shown in Fig. 1, and $`S_i=\pm 1`$ represent Ising spins. We place the spins on a ring in order to ensure periodic boundary conditions. The sum is over all spins on the chain and the couplings $`J_{ij}`$ are given by

$$J_{ij}=c(\sigma )\frac{ϵ_{ij}}{r_{ij}^\sigma },$$
(2)
where the $`ϵ_{ij}`$ are chosen according to a Gaussian distribution with zero mean and standard deviation unity:

$$𝒫(ϵ_{ij})=\frac{1}{\sqrt{2\pi }}\mathrm{exp}(ϵ_{ij}^2/2).$$
(3)
The constant $`c(\sigma )`$ in Eq. (2) is chosen to give a mean-field transition temperature $`T_c^{\mathrm{MF}}=1`$, where

$$\left(T_c^{\mathrm{MF}}\right)^2=\underset{ji}{}[J_{ij}^2]_{\mathrm{av}}=c(\sigma )^2\underset{ji}{}\frac{1}{r_{ij}^{2\sigma }}.$$
(4)
Here $`[\mathrm{}]_{\mathrm{av}}`$ denotes an average over disorder.

The distance between two spins on the chain $`r_{ij}`$ is determined by $`r_{ij}=2R\mathrm{sin}(\alpha /2)`$, where $`R`$ is the radius of the chain and $`\alpha `$ is the angle between the two sites on the circle (see Fig. 1). The previous expression can be rewritten in terms of $`L`$ and the positions of the spins to obtain

$$r_{ij}=\frac{L}{\pi }\mathrm{sin}\left(\frac{\pi |ij|}{L}\right).$$
(5)
The long-range Ising spin glass with power-law interactions has a very rich phase diagram in the $`d`$-$`\sigma `$ plane. This is summarized in Fig. 2, which is based on the work by Bray et alBray et al. (1986) and by Fisher and HuseFisher and Huse (1988) who present a detailed analysis of the role of long-range interactions within the droplet model. Spin-glass behavior is controlled by the long-range part of the interaction if $`\sigma `$ is sufficiently small, and by the short-range part if $`\sigma `$ is sufficiently large. More precisely, one has long-range behavior if the stiffness exponent of the long-range (LR) universality class, $`\theta _{\mathrm{LR}}`$, is greater than that of the short-range (SR) universality class, $`\theta _{\mathrm{SR}}`$, and vice versa. In addition, there is an exact result for $`\theta _{\mathrm{LR}}`$, namely,Bray et al. (1986); Fisher and Huse (1988)

$$\theta _{\mathrm{LR}}=d\sigma ,$$
(6)
so long-range behavior occurs if

$$\sigma <\sigma _c(d)=d\theta _{\mathrm{SR}}(d).$$
(7)
Equation (6) indicates that critical exponents depend continuously on $`\sigma `$ in the long-range region, even though they are independent of $`\sigma `$ in the region controlled by the short-range part of the interaction.

The condition for a finite-temperature transition is $`\theta >0`$, where $`\theta `$ refers here to the greater of $`\theta _{\mathrm{SR}}`$ and $`\theta _{\mathrm{LR}}`$. For the short-range model, there is a finite-temperature transition (i.e., $`\theta _{\mathrm{SR}}>0`$) for $`d`$ larger than the lower critical dimension $`d_l`$, which is found numerically to lie between 2 and 3.Bray and Moore (1984); McMillan (1984a, b); Kawashima and Young (1996); Ballesteros et al. (2000) In Figs. 3 and 4, we show the expected variation of $`\theta `$ with $`\sigma `$ both for $`d>d_l`$ and $`d<d_l`$.

For $`\sigma <d/2`$ the model does not have a thermodynamic limit ($`T_c`$ diverges) unless the interactions are scaled with an appropriate (inverse) power of $`L`$, i.e., $`c(\sigma )0`$ for $`L\mathrm{}`$. Note that $`\sigma =0`$ corresponds to the SK model.

For $`d=1`$, the case studied in our numerics, one hasBray and Moore (1984)

$$\theta _{\mathrm{SR}}(1)=1,$$
(8)
and so from Eq. (7), $`\sigma _c=2`$, in agreement with rigorous results by van Enter and van Hemmen.van Enter and van Hemmen (1985) This situation corresponds to Fig. 4 with $`\theta _{\mathrm{SR}}=1`$. To obtain a finite $`T_c`$ in $`d=1`$, one needs $`1/2<\sigma <1`$ and, in the finite-$`T`$ simulations, we choose $`\sigma =0.75`$ (see Fig. 2). For the zero-temperature calculations we use a range of $`\sigma `$ between $`0.10`$ and $`3.00`$. From Figs. 2 and 4, we see that this includes all the possible types of behavior: (i) $`T_c>0`$, (ii) $`T_c=0`$ (LR), (iii) $`T_c=0`$ (SR), and (iv) SK-like behavior for $`\sigma 0`$.

## III Observables

In our finite-$`T`$ simulations we primarily focus our attention on the spin overlap defined by

$$q=\frac{1}{L}\underset{i=1}{\overset{L}{}}S_i^\alpha S_i^\beta ,$$
(9)
where “$`\alpha `$” and “$`\beta `$” refer to two copies (replicas) of the system with the same disorder. In order to determine the critical temperature of the system, we study the Binder ratioBinder (1981) $`g`$ defined by

$$g=\frac{1}{2}\left(3\frac{[q^4]_{\mathrm{av}}}{[q^2]_{\mathrm{av}}^2}\right)=\stackrel{~}{g}[L^{1/\nu }(TT_c)].$$
(10)
Here, $`\mathrm{}`$ represents a thermal average, $`\nu `$ is the correlation length exponent, and $`T_c`$ is the critical temperature of the system. Because $`g`$ is dimensionless, curves for $`g(T,L)`$ cross at the same $`T_c`$ for all $`L`$. We also compute the spin-glass susceptibility $`\chi _{_{\mathrm{SG}}}`$, where

$$\chi _{_{\mathrm{SG}}}=L[q^2]_{\mathrm{av}}=L^{2\eta }\stackrel{~}{\chi }[L^{1/\nu }(TT_c)],$$
(11)
and $`\eta `$ describes the power-law falloff of the spin-glass correlations at $`T_c`$.

## IV Equilibration

For the finite-$`T`$ simulations, we use the parallel tempering Monte Carlo methodHukushima and Nemoto (1996); Marinari (1998) as it allows us to study larger systems at lower temperatures than is possible with conventional “single-spin-flip” Monte Carlo methods. We test equilibration by the traditional technique of requiring that different observables are independent of the number of Monte Carlo steps $`N_{\mathrm{sweep}}`$. Figure 5 shows data for the spin overlap as a function of Monte Carlo steps for different system sizes. We equilibrate by repeatedly doubling the number of Monte Carlo sweeps until the last three agree within error bars.

We also checked the equilibration of some of our data using a generalization of a test that was developed earlierKatzgraber et al. (2001) for nearest-neighbor models. For a Gaussian distribution of bonds it is straightforward to show, by integrating by parts with respect to the bonds, that the energy per spin, $`U=N^1_{i,j}[J_{ij}S_iS_j]_{\mathrm{av}},`$ can be related to a generalized “link overlap” $`q_l`$, defined by

$$q_l=\frac{2}{N}\underset{i,j}{}\frac{[J_{ij}^2]_{\mathrm{av}}}{(T_c^{MF})^2}[S_iS_j^2]_{\mathrm{av}},$$
(12)
as follows:

$$q_l=1\frac{2T|U|}{(T_c^{\mathrm{MF}})^2}.$$
(13)
Note that this definition of $`q_l`$ is the natural generalization, to long-range interactions, of the link overlap discussed in earlier workPalassini and Young (2000); Katzgraber et al. (2001) on nearest-neighbor models, and that $`q_l`$ is normalized so that it tends to unity for $`T0`$ where $`S_iS_j^21`$.

As discussed in Ref. Katzgraber et al., 2001, the two sides of Eq. (13) are only equal to each other in equilibrium and approach the equilibrium value from opposite sides. This is illustrated in Fig. 6 which shows $`q_l`$ and the quantity on the right-hand side (RHS) of Eq. (13) for an increasing number of sweeps for $`L=128,T=0.1`$. The two quantities do indeed approach each other from opposite sides, and once they agree they do not change if the simulation is run for longer. We applied this test for several lattice sizes using the parameters in Table 1, and found that the data were well equilibrated in each case.

We also require that the acceptance ratios for the global moves which interchange the different temperatures in the parallel tempering scheme be greater than $`0.4`$ on average and roughly constant as a function of temperature. Table 1 shows the number of samples $`N_{\mathrm{samp}}`$ and the number of Monte Carlo sweeps $`N_{\mathrm{sweep}}`$ performed by each replica for each lattice size. The lowest temperature used is $`T=0.10`$, the highest $`T=1.40`$, well above the mean-field critical temperature $`T_c^{\mathrm{MF}}=1`$.

We also use parallel tempering as an optimization algorithm to determine the ground state in the $`T=0`$ calculations. This is described in Sec. VII.

## V Finite-Size Scaling in the Critical Region

In this section, we determine the critical temperature $`T_c`$ and the critical exponents for the disordered chain with $`\sigma =0.75`$. In Fig. 7, we show rescaled data for the Binder ratio $`g`$ as a function of temperature for different system sizes \[see Eq. (10)\]. The data scale well for $`T_c=0.62\pm 0.03`$ and $`1/\nu =0.30\pm 0.03`$, and agree with previous results by Leuzzi.Leuzzi (1999) The inset of the aforementioned figure shows the unscaled data for $`g`$ which cross at $`T_c0.62`$. The error bars are estimated by varying $`g`$ slightly until the data do not collapse well.

We use the same approach below for the spin-glass susceptibility, for which the data, rescaled according to Eq. (11), are presented in Fig. 8. The data scale well for $`\eta =1.34\pm 0.03`$ and the previously mentioned values of $`1/\nu `$ and $`T_c`$. These values can be compared with the results of Kotliar et al.:Kotliar et al. (1983) $`\eta =32\sigma (=1.5)`$, and $`1/\nu =(1/3)4ϵ`$ to first order in $`ϵ`$, where $`ϵ=\sigma 2/3`$. The latter gives $`1/\nu =0`$ suggesting that truncating the expansion to first order in $`ϵ`$ is not valid for $`\sigma =3/4`$.

Having studied the critical region of the model to estimate $`T_c`$ we are now able to simulate the system at low temperatures ($`T0.16T_c`$) in order to probe the spin-glass phase and test if $`\theta =\theta ^{}`$.

## VI Results at Low Temperatures

In this section we study the model with $`\sigma =0.75,d=1`$, marked by a cross in Fig. 2, at temperatures well below $`T_c0.62`$, in order to get information about large-scale, low-energy excitations in the spin-glass state.

Differences between the various models for the spin-glass state can be quantified by studyingMarinari and Parisi (2000); Katzgraber et al. (2001); Reger et al. (1990); Marinari et al. (1998); Marinari and Zuliani (1999) $`P(q)`$, the distribution of the spin-glass overlap $`q`$. The RSB picture predicts a nontrivial distribution with a finite weight in the tail down to $`q=0`$, independent of system size. On the contrary, the droplet picture predicts that $`P(q)`$ should be trivial with only two peaks at $`\pm q_{\mathrm{EA}}`$, where $`q_{\mathrm{EA}}`$ is the Edwards-Anderson order parameter. For finite systems there is also a tail down to $`q=0`$ but this vanishes in the thermodynamic limit likeFisher and Huse (1986); Bray and Moore (1986); Moore et al. (1998) $`L^\theta ^{}`$ with $`\theta ^{}=\theta `$.

Figures 9 and 10 show data for $`P(q)`$ at temperatures 0.10 and 0.23, respectively. There is clearly a peak for large $`q`$ and a tail down to $`q=0`$. At both temperatures one sees that the tail in the distribution is essentially independent of system size. A more precise determination of the size dependence of $`P(0)`$ is shown in Fig. 11 where, to improve statistics, we average over the $`q`$ values with $`|q|<q_{}`$, with $`q_{}=0.50`$. The expectedFisher and Huse (1986); Bray and Moore (1986); Moore et al. (1998) behavior in the droplet model is $`P(0)L^\theta ^{}`$, with $`\theta ^{}=\theta `$ whereBray et al. (1986); Fisher and Huse (1988) $`\theta =d\sigma `$. The dashed line in Fig. 11 has slope $`0.25`$, the expected value for $`\sigma =0.75`$ according to the droplet model. The size dependence is consistent with a constant $`P(0)`$, which implies that the energy to create a large excitation does not increase with size, at least for the range of sizes studied here.

We perform a two-parameter fit of the data in Fig. 11 of the form $`aL^\theta ^{}`$ to find for $`T=0.10`$, $`\theta ^{}=0.005\pm 0.009`$, for $`T=0.16`$, $`\theta ^{}=0.020\pm 0.007`$, and for $`T=0.23`$, $`\theta ^{}=0.016\pm 0.006`$. The fitting probabilitiesPress et al. (1995) $`Q`$ for the different temperatures are $`0.135`$, $`0.186`$, and $`0.08`$, respectively, which is reasonable. Fixing $`\theta ^{}=0`$, we obtain fitting probabilities of $`0.214`$, $`0.053`$, and $`0.076`$, respectively. On the contrary, fixing $`\theta ^{}=0.25`$, the predicted value for the droplet picture, we find $`Q10^{31}`$ for all three temperatures. Hence we are confident that, within error bars, $`\theta ^{}0`$ for the range of sizes studied.

Moore Moore (1986) has argued that $`P(q)`$ should be trivial in the thermodynamic limit (i.e., just two delta functions at $`q=\pm q_{EA}`$) to first order in an expansion about the “lower critical dimension” $`\sigma =1`$. This seems to disagree with our numerics, although possibly even larger sizes might be necessary to see the asymptotic behavior.

## VII Domain-wall energies at zero temperature

Since there are analytic predictions, Eqs. (6) and (8), for $`\theta `$ for the present model, it is useful to directly calculate $`\theta `$ from zero-temperature domain-wall calculations for different values of $`\sigma `$.

We use parallel tempering here as an optimization algorithm. We find that if we take the lowest temperature $`T_{\mathrm{min}}`$ to be $`0.05`$, then the minimum-energy state found at this temperature is the ground state. To test whether the true ground state has been reached, two criteria were adopted: (i) the same minimum-energy state has to be reached from two replicas at $`T_{\mathrm{min}}`$ for all samples, and (ii) this state has to be reached during the first $`1/4`$ of the sweeps in both copies. These conditions were satisfied for the parameters used in the simulations which are shown in Table 2.

Having found the ground state with periodic (P) boundary conditions, we then flip the boundary conditions to “antiperiodic” (AP) by the following prescription. Consider one of the nearest neighbor (nn) bonds, e.g., that between sites 1 and $`L`$ in Fig. 1. Then, for all $`i`$ and $`j`$, change the sign of the bond between $`i`$ and $`j`$ if the shorter path between those sites goes through the chosen nn bond. The new ground state is calculated and note is taken of the difference in energy

$$\delta E=E_{\mathrm{AP}}E_\mathrm{P}$$
(14)
for each sample. On average, there is no preference for AP or P and so $`[\delta E]_{\mathrm{av}}=0`$. Hence we take the average of the absolute value

$$\mathrm{\Delta }E=[|\delta E|]_{\mathrm{av}}$$
(15)
as a measure of the characteristic domain-wall energy.

Data for a range of sizes and values of $`\sigma `$ are shown in Fig. 12. Our program automatically scales the interactions to set $`T_c^{\mathrm{MF}}=1`$, so no modifications are needed for the case $`\sigma 0.5`$ since the fact that the scale factor $`c(\sigma )`$ \[in Eqs. (2\] and (4)) tends to zero for $`L\mathrm{}`$ in this range, is taken care of automatically.

A smaller range of sizes has been studied for larger $`\sigma `$ as parallel tempering is less efficient in finding the ground state in this limit. At first sight this seems surprising because the nearest-neighbor interactions dominate and so there is no frustration, apart possibly from that due to the boundary conditions. The spins become stuck when a frustrated nn bond is larger in magnitude than the neighboring bonds. Parallel tempering does not seem to help much in this situation, perhaps because the local minima do not have a hierarchical structure. It would be interesting to investigate in more detail the conditions under which parallel tempering is efficient.

One expects that $`\mathrm{\Delta }EL^\theta `$, where $`\theta =1\sigma `$ in the long-range region $`1/2<\sigma <2`$, and $`\theta =1`$ in the short-range region $`\sigma >2`$. Fits to the data obtained in Fig. 12, along with the analytic prediction, are shown in Fig. 13 and in Table 3. In particular, there is a distinctive positive slope for $`\sigma =0.75`$ in Fig. 12 (solid line), corresponding to a positive value of $`\theta `$. This is in clear contrast to the finite-$`T`$ data shown in Fig. 11.

We see that the overall trends in the analytic prediction are reproduced by the data, but there is some rounding where the analytical result has a change in slope. Note that our results for $`\theta `$ always lie well below $`1/2`$, the predictedBray et al. (1986); Fisher and Huse (1988) value for $`\sigma =1/2`$, and appear to saturate at about $`0.3`$ as $`\sigma 0`$, the SK limit.

In recent work Aspelmeier et al.Aspelmeier et al. (2002) have given analytic arguments, based on the RSB theory, that the defect energy should vary as $`N^\mu `$ for short-range models in high dimension $`d`$, where $`N=L^d`$ and $`\mu `$ is the exponent for the rms sample-to-sample fluctuations of the free energy of the SK model. A number of authorsCabasino et al. (1988); Bouchaud et al. (2002); Palassini (2003) have found that $`\mu 1/4`$, at least at $`T=0`$. We have therefore calculated $`\mu `$ for our one-dimensional model at $`T=0`$ for different values of $`\sigma `$ and show the results in Fig. 14. Indeed we find that $`\mu 1/4`$ for $`\sigma 0`$ (the SK limit). Also note that it crosses over to the expected “self-averaging” value of $`\mu =1/2`$ in the region $`\sigma >1/2`$ where we no longer need to scale the interactions with system size. Since we work in $`d=1`$ the expectation from Aspelmeier et al.Aspelmeier et al. (2002) is that $`\theta \mu `$ for $`\sigma 0`$. Figure 14 also shows $`\theta `$ and, indeed, it is quite striking how $`\theta `$ becomes close to $`\mu `$ at small $`\sigma `$.

## VIII Conclusions

We have performed finite-$`T`$ Monte Carlo simulations and ground-state calculations of the long-range disordered Ising chain in one dimension with power-law interactions. For the ground-state calculations we have studied a range of values of $`\sigma `$ that encompasses all the possible regimes, from large $`\sigma `$ where the system is controlled by the short-range part of the interactions and has $`T_c=0`$, all the way to near the SK model limit, $`\sigma =0`$. For the finite-temperature Monte Carlo runs we have concentrated on a single value of $`\sigma `$, namely $`0.75`$, which is in the range where there is a finite $`T_c`$.

For $`\sigma =0.75`$ we find the stiffness exponent for droplets, $`\theta ^{}`$, inferred from Monte Carlo simulations for sizes up to $`L=512`$, to be close to zero, whereas the stiffness exponent for domain walls, $`\theta `$, is positive. We find $`\theta =0.173\pm 0.005`$ (only statistical errors included) somewhat lower than the analytic predictionBray et al. (1986); Fisher and Huse (1988) that $`\theta =0.25`$. These results are consistent with the RSB and TNT scenarios according to which $`\theta ^{}=0`$ and $`\theta \theta ^{}`$. Similar results have been found in earlier work on short-range models in three and four dimensions.Krzakala and Martin (2000); Palassini and Young (2000); Katzgraber et al. (2001) According to the droplet picture, $`\theta ^{}=\theta `$ asymptotically, and there has been discussion as to whether the failure of the numerics to see this is an artifact of the small lattice sizes,Moore (2002) or whether it will persist in the thermodynamic limit. In the present work, the disagreement persists even though we are able to study a much larger range of $`L`$.

Also noteworthy is our finding that $`\theta `$ tends to a value of about 0.3 for $`\sigma 0`$, the SK model limit. This is below the droplet theory prediction that $`\theta `$ should tend to 1/2 for $`\sigma 0`$, but is close to the prediction of Aspelmeier et al.Aspelmeier et al. (2002) that $`\theta 1/4`$ in this limit.

###### Acknowledgements.

We would like to that M. A. Moore for stimulating discussions and for bringing Ref. Moore, 1986 to our attention, and G. Blatter for carefully reading the manuscript. H.G.K. acknowledges support from the National Science Foundation under Grant No. DMR 9985978 and the Swiss National Science Foundation. A.P.Y. acknowledges support from the National Science Foundation under Grant No. DMR 0086287, and the EPSRC under Grant GR/R37869/01. He would also like to thank David Sherrington for hospitality during his stay at Oxford. This research was supported in part by NSF cooperative agreement ACI-9619020 through computing resources provided by the National Partnership for Advanced Computational Infrastructure at the San Diego Supercomputer Center. We would like to thank the University of New Mexico for access to their Albuquerque High Performance Computing Center. This work utilized the UNM-Alliance Los Lobos Supercluster. Part of the simulations were performed on the Asgard cluster at ETH Zürich.


format: text

text:
# A Model for Neutrino Masses and Dark Matter

## I Introduction

The discovery of very small, but non-zero neutrino masses via the observation of oscillations in atmospheric and solar neutrino experiments superk ; SNO  is one of the most exciting results in particle physics in recent years, and may provide important new information to help guide us in the search for new physics beyond the standard model. As we describe here, efforts to uncover a mechanism for generating such small mass differences may also shed light on the nature of non-baryonic dark matter.

The simplest way to introduce neutrino masses into the standard model is to postulate a new right-handed fermion, a singlet under the SM gauge group, that mixes with the SM fermions via a Dirac mass term. Unfortunately, in order to obtain neutrino masses below an eV, it is necessary to fine-tune the corresponding Yukawa coupling to one part in $`10^{11}`$, which is unacceptable. Thus, more elaborate models are needed.

Currently the most popular way to generate small neutrino masses is the seesaw mechanism seesaw . In this model one introduces a right handed neutrino $`N_R`$, and allows both Majorana and Dirac mass terms. Since the right handed neutrino is a singlet under the SM gauge group $`SU(2)_L\times U(1)_Y`$, its mass $`M_R`$ can be much larger than the electroweak scale without contributing large quantum corrections to the other couplings in the theory. One may integrate out such a heavy field leaving a Majorana neutrino with mass $`m_\nu \frac{m_D^2}{M_R}`$, where $`m_D`$ is the Dirac neutrino mass. The seesaw mechanism arises naturally in many extensions of the standard model gauge group, for example in left-right (L-R) models with group $`SU(3)_c\times SU(2)_L\times SU(2)_R\times U(1)_{BL}`$, the Pati-Salam model with group $`SU(4)_c\times SU(2)_L\times SU(2)_R`$ or in $`SO(10)`$ grand unified theories (GUTs).

Another possible way of achieving small neutrino masses is to extend the Higgs sector of the theory by adding a very heavy triplet scalar field $`\mathrm{\Delta }`$, with mass $`M_\mathrm{\Delta }`$, which couples to the SM fields via

$$_\mathrm{\Delta }=f_{\alpha \beta }L_\alpha ^TCi\tau _2\mathrm{\Delta }L_\beta +\mu \mathrm{\Phi }^Ti\tau _2\mathrm{\Delta }^+H\mathrm{\Phi }+\mathrm{h}.\mathrm{c}.,$$
(1)
where, $`L_\alpha `$ is the left-handed lepton doublet, $`\mathrm{\Phi }`$ is the standard model Higgs doublet and $`C`$ is the charge-conjugation matrix. Here the Yukawa couplings $`f_{\alpha \beta }`$ are antisymmetric in the generation indices $`\alpha `$ and $`\beta `$.

Integrating out our heavy triplet $`\mathrm{\Delta }`$ one obtains an effective theory whose Lagrangian density contains the term

$$_{LH}=\frac{1}{\mathrm{\Lambda }}LLHH,$$
(2)
where the scale $`\mathrm{\Lambda }M_\mathrm{\Delta }^2`$. Due to this term, after electroweak symmetry breaking the neutrino gets a mass $`m_\nu \frac{V_{EW}^2}{M_\mathrm{\Delta }}`$ where $`V_{EW}`$ is the vacuum expectation value (VEV) of the standard model higgs field.

Clearly, in both of these sets of models for generating neutrino masses, it is crucial that the masses of the right-handed neutrino and the scalar triplet respectively are much higher than the electroweak scale. As a consequence, these models are difficult to test in collider experiments. In addition, the existence of a large hierarchy between the electroweak scale and the new physics scale generates its own set of naturalness problems.

Here we reconsider the idea that SM neutrinos are massless at tree level, but acquire a nonzero mass at the one-loop level. A particularly simple example of this idea is provided by the Zee model zee  in which one augments the SM higgs doublet $`\mathrm{\Phi }_1`$ with a second higgs doublet $`\mathrm{\Phi }_2`$, and a charged field $`S`$ which transforms as a singlet under $`SU(2)`$. The Lagrangian density then contains the following relevant pieces

$$_{\mathrm{Zee}}=f_{\alpha \beta }L_\alpha ^TCi\tau _2L_\beta S^++\mu \mathrm{\Phi }_1^Ti\tau _2\mathrm{\Phi }_2S^{}+\mathrm{h}.\mathrm{c}..$$
(3)
Other extensions of the Zee model where small neutrino masses arise at loop levels have also been proposed Nasri:2001ax  that may have interesting implications for the generation of the baryon asymmetry of the universe Nasri:2001nb ; Luty:un .

In this letter we will demonstrate how a variant of the Zee model allows the generation of small phenomenologically acceptable neutrino masses and at the same time produces a new viable dark matter candidate near the electroweak scale.

## II The Model

We consider a model with sufficient symmetry that neutrino masses appear only at the three loop level. To achieve this we supplement the SM fields with two charged scalar singlet scalars $`S_1`$ and $`S_2`$, with masses $`M_{S_1}`$ and $`M_{S_2}`$ respectively and one right handed neutrino $`N_R`$, with mass $`M_R`$. We break lepton number explicitly by including a Majorana mass term for the right-handed neutrino, and impose a discrete $`Z_2`$ symmetry under which the SM fields and $`S_1`$ are singlets but $`S_2`$ and $`N_R`$ transform as

$$Z_2:\{S_2,N_R\}\{S_2,N_R\},$$
(4)
forbidding Dirac masses for the neutrinos. Given this symmetry, the most general renormalizable terms that may be added to the SM Lagrangian density

$`_{\mathrm{new}}=`$ $`f_{\alpha \beta }L_\alpha ^TCi\tau _2L_\beta S_1^++g_\alpha N_RS_2^+l_{\alpha _R}`$ (5)
$`+M_RN_R^TCN_R+V(S_1,S_2)+\mathrm{h}.\mathrm{c}.,`$
in which the potential $`V(S_1,S_2)`$ contains a $`(S_1S_2^{})^2`$ coupling $`\lambda _s`$. Let us assume a mild hierarchy of masses $`M_R<M_{S_1}<M_{S_2}`$ TeV and that the Yukawa couplings $`f_{\alpha \beta }`$, $`g_\alpha `$ are of order one. This implies that $`N_R`$ is stable if the above discrete symmetry is unbroken. Also, since lepton number is broken, we expect a left handed Majorana neutrino mass at the quantum level, and indeed this arises from the three loop diagram in figure 1.

Once again, we integrate out the heavy degrees of freedom $`N_R`$, $`S_1`$ and $`S_2`$ to obtain an effective theory containing the dimension five operator

$$_{LH}=\frac{1}{\mathrm{\Lambda }}(L_\alpha ^TC\overline{\tau }\mathrm{\Phi }_1)(\varphi _1^T\overline{\tau }L_\alpha ).$$
(6)
Here the scale $`\mathrm{\Lambda }`$ is given in terms of the heaviest singlet scalar mass $`M_{S_2}`$ and the Yukawa couplings $`\lambda _l=\frac{m_l}{V_{EW}}`$, $`f_{\alpha \beta }`$, $`g_\alpha `$ and $`\lambda _s`$ as

$$\mathrm{\Lambda }\frac{(4\pi ^2)^3}{f^2\lambda _s\lambda _{l}^{}{}_{}{}^{2}g_\alpha ^2}M_{S_2},$$
(7)
where a subscript $`l`$ denotes a leptonic quantity. For $`M_{S_2}`$ TeV, we obtain $`\mathrm{\Lambda }>10^9`$ GeV and so the model provides a hierarchy between the electron and neutrino masses and yields neutrino masses at the $`0.1`$ eV scale without involving mass scales significantly larger than a TeV.

## III Dark Matter

As we commented above, in our model the right-handed neutrino is stable, since its mass is smaller than $`M_{S_1}`$ and the discrete symmetry $`Z_2`$ is unbroken. It is natural to wonder whether this stable, neutral, weakly-interacting particle might be a reasonable cold dark matter (CDM) candidate. To determine this we need to calculate the relic abundance of right-handed neutrinos produced as the universe cools.

We begin by noting that the right-handed neutrino annihilates into right-handed charged leptons with cross-section given by

$$\sigma v\frac{g^4}{\pi }\frac{M_R^2}{M_{S_2}^4}.$$
(8)
We also need the standard result that the number density of right-handed neutrinos at the decoupling temperature $`T_D`$ is

$$n_{N_R+\overline{N}_R}=\frac{2T_D^3}{(2\pi )^{3/2}}\left(\frac{M_R}{T_D}\right)^{3/2}\mathrm{exp}\left(\frac{M_R}{T_D}\right).$$
(9)
Freeze-out occurs when

$$n_{N_R}\sigma \sqrt{g_{}}\frac{T_D^2}{M_{pl}},$$
(10)
where $`g_{}`$ is the effective number of massless degrees of freedom of at decoupling. Therefore, using our expressions for $`n_{N_R+\overline{N}_R}`$ and the cross-section we solve for $`\frac{M_R}{T_D}`$ to obtain

$`{\displaystyle \frac{M_R}{T_D}}=`$ $``$ $`\mathrm{ln}\left({\displaystyle \frac{1.66(2\pi )^{\frac{3}{2}}M_{S_2}^4}{M_{pl}g^4\mathrm{GeV}^3}}\right){\displaystyle \frac{1}{2}}\mathrm{ln}\left({\displaystyle \frac{M_R}{T_D}}\right)`$ (11)
$`+`$ $`3\mathrm{ln}\left({\displaystyle \frac{M_R}{\mathrm{GeV}}}\right)\mathrm{ln}\left({\displaystyle \frac{\sqrt{g_{}}}{g_R}}\right),`$
where the minus sign in front of the second term arises because the annihilation of $`N_R`$ is a p-wave process, and where we can ignore the masses of the annihilation products with masses less than $`M_R`$. For typical values $`g^20.1`$ and $`M_{S_2}`$ TeV we obtain

$$\frac{M_R}{T_D}20.$$
(12)
The ratio of the energy density in right-handed neutrinos to the critical density is approximately given by

$$\mathrm{\Omega }_{N_R}\frac{10^{37}\mathrm{cm}^2}{\sigma v}.$$
(13)
Note that if we take typical values $`g^20.1`$, $`M_R`$ TeV then $`\sigma v10^{36}\mathrm{cm}^2`$ and so one quite easily finds

$$\mathrm{\Omega }_{N_R}0.11,$$
(14)
as required for a viable CDM candidate.

Note also that because the right handed neutrino has no direct neutral current couplings to quarks, it would not be detected in existing WIMP detection experiments sensitive to elastic scattering off of nuclei.

## IV Current Experimental Constraints

We next consider whether what current phenomenological constraints on lepton number violation imply for this model.

First note that stringent constraint on the parameters $`f_{\alpha \beta }`$ can be derived from the new contribution to the muon decay effective lagrangian mediated by $`S_1`$. We find

$$\frac{f_{e\mu }}{M_{S_1}}10^4\mathrm{GeV}^2$$
(15)
The Lagrangian $`(6)`$ can also lead to the flavor violating process $`\mu e+\gamma `$. Using the experimental bound on its branching ratio we find:

$$\frac{f_{e\alpha }f_{\mu \alpha }^{}}{M_{S_1}^2}\mathrm{2.8.10}^9\mathrm{GeV}^2$$
(16)
The above constraints can be satisfied for the $`f_{e\alpha }0.1`$ and $`M_{S_1}`$ few TeV, and are also consistent with limits from neutrinoless double $`\beta `$-decay experiments, and are thus comparable with the requirements needed to generate neutrino masses and dark matter, as described above.

Although the right-handed neutrino is a singlet under the SM gauge group, it may be produced in our model through the scatterings of SM particles and the subsequent decay of $`S_2`$.

If the center of mass energy of the next generation linear collider exceeds twice the mass of $`S_2`$ then these particles would be produced in, for example, $`e^+e^{}`$ annihilation. In this case, the dominant decay channel for the $`S_2`$ particle would be into a charged lepton and the lightest right-handed neutrino, providing a possible collider signature of the model.

## V Conclusions

It would be particularly exciting if a single mechanism might resolve two of the most important outstanding puzzles in particle physics and cosmology: the nature of dark matter, and the origin of neutrino masses. As we demonstrate here, it is possible to extend standard model of particle physics in a way so that both puzzles may find a common resolution. The symmetry structure of our model is such that neutrino masses occur only at three-loop level, and so are naturally small. As a consequence, the right-handed neutrino is not constrained to be have a large mass, as in other examples of neutrino mass generation, such as the seesaw model. Such a light right-handed neutrino is a possible weakly interacting massive particle (WIMP) dark matter candidate.

We have explored the experimental constraints on our model and it currently passes all standard model tests. More exciting perhaps is the fact that the TeV-scale scale energies required should be accessible at colliders in the near future. The model may also be tested by future neutrinoless double $`\beta `$-decay experiments. A new sort of direct detection WIMP technology would be required to directly detect the dark matter candidate we propose.

###### Acknowledgements.

LMK and SN are supported by the US Department of Energy (DOE) LMK acknowledges the hospitality of the Kavli Institute for Theoretical Physics during this work. The work of MT is supported by the National Science Foundation (NSF) under grant PHY-0094122.


text:
# Elastic Tensor of YNi2B2C

## Abstract

The complete elastic tensor of YNi2B2C was determined by application of the resonant ultrasound spectroscopy technique to a single-crystal sample. Elastic constants were found to be in good agreement with partial results obtained from ‘time-of-flight’ measurements performed on samples cut from the same ingot. From the measured constants, the bulk modulus and Debye temperature are calculated.

The discovery of superconductivity in the borocarbide quaternary intermetallic compoundsCava  has generated a great deal of interest due to the observed interplay between superconductivity and magnetic order. Compounds in this family have formulae of the form *R*Ni2B2C, where choice of the rare earth ion, *R*, strongly affects the physical properties, which range from superconducting with no magnetic ordering (*R*=Y) to non-superconducting with a Néel temperature of 19 K (*R*=Gd).Canfield

Recent heat transport measurements on non-magnetic borocarbides suggest that the superconductivity is very anomalous, with a gap that is either highly anisotropicBoaknin  or has point nodes.Matsuda  To elucidate the nature of the order parameter, a powerful tool is ultrasound attenuation (see for example Ref. Lupien\_PE, ). As a first step towards measuring the attenuation of sound it is important to establish the full elastic tensor of the material.

The 2nd order elastic constants relate stresses placed on a material along various directions to the resulting deformations of the crystal lattice. As such, knowledge of a compound’s elastic constants can be useful when investigating a number of solid state properties, such as interatomic potentials, equations of state, phonon spectra, the Debye temperature, and the bulk modulus. In the most general case, a highly anisotropic solid can have 21 independent second-order elastic constants, but the crystal structure for all compounds in the *R*Ni2B2C group is body-centered-tetragonal (space group *I*4/*mmm*), consisting of *R*-C planes separated by Ni2B2 layers stacked along the *c* axis,Lynn  reducing the number of independent constants to 6.Migliori  In Voigt notation, these constants are denoted *C*11, *C*33, *C*13, *C*12, *C*44, and *C*66; all other tensor entries can be retrieved from these six.

In this article, the elastic constants of YNi2B2C are determined at room temperature via two methods: resonant ultrasound spectroscopy and ultrasonic time-of-flight. Resonant ultrasound spectroscopy (RUS)Migliori  allows the determination of the complete elastic tensor of a material in one measurement, and was used as the primary experimental technique of this study. As a complement to the RUS measurements, elastic constants were also determined by measuring sound velocities using the more conventional ‘time-of-flight’ method. Note that this study is not simply a confirmation of existing measurements: it adds a missing elastic constant and provides a very accurate second measure to secure time-of-flight results in slight disagreement.

The RUS method utilizes the normal modes of vibration of a solid sample of known density, geometry and crystal structure to determine the elastic tensor. To perform an RUS measurement, the normal mode frequency spectrum of the sample in question is first obtained. This is done by holding the sample between two piezo-electric transducers, generally by its corners in order to reduce damping of the vibrational modes. One transducer is driven at a frequency which is swept across the region of the spectrum containing the normal modes, thus driving the sample resonances, while the other picks up the resulting amplitudes of crystal vibration (see Fig. 1).

Once a list of normal mode frequencies has been compiled, the elastic constants can be calculated. While it is possible to directly compute the normal mode frequencies of a sample, given the dimensions, mass and elastic constants, there is no known analytic method to do the reverse, so a computerized fitting algorithm has been developedMigliori  to numerically complete the calculation. This algorithm changes the elastic constants, calculates resulting normal mode frequencies, compares these calculated frequencies to the experimentally measured list and then iteratively repeats the process until the best fit to the data is obtained.

In the time-of-flight method, short ultrasound pulses are propagated through a sample carefully prepared with two flat, parallel faces, and later detected in order to obtain the sound velocity in the direction of pulse propagation. The diagonal entries of the elastic tensor can be easily calculated using the velocities measured along various principal crystal directions. The off-diagonal constants, however, require measurements of sound velocities in non-principal directions, which can be difficult. The equation (*C*11 \+ *C*12 \+ 2*C*66)/2 = $`\rho `$ $`v^2`$\[110/L\], where $`\rho `$ is the mass density of the material, and $`v_{[110/L]}`$ is the velocity of sound propagating in the  direction with a longitudinal polarization, is an example of how a sound velocity is related to the elastic constants for a tetragonal crystal.Truell  In this notation scheme, a transverse polarization, for example in the  direction, is indicated by $`T010`$.

With the complete elastic tensor, the bulk modulus of a tetragonal crystal, $`B_0=V\left(\frac{dP}{dV}\right)`$, is easily obtained via the equation:

$$B_0=\frac{C_{33}(C_{11}+C_{12})2C_{13}^2}{C_{11}+C_{12}+2C_{33}4C_{13}}$$
(1)
Also, an estimate of the Debye temperature can be made. For a crystal with tetragonal structure the following equationAlers  is employed:

$$\theta _D=\frac{h}{k_B}\left(\frac{9N}{4\pi V}\right)^{1/3}\rho ^{1/2}J^{1/2},$$
(2)
wherein *h* is Planck’s constant, *kB* is Boltzmann’s constant, *N* is the number of atoms per unit cell, *V* is the volume of the unit cell, $`\rho `$ is the mass density of the material, and *J* is the harmonic series expansion of an integration over sound velocities in all directions, written in terms of the elastic constants (Eqs. 32 and 33 in Ref. Alers, ).

A large high-quality single crystal specimen of YNi211B2C was grown using the floating zone method. From this specimen a sample was cut using spark erosion and then polished into the shape of a rectangular parallelepiped, with sides aligned along the principal crystal axes as determined by Laue x-ray backscattering measurements.

Dimensions of the polished sample were measured to be 0.0633(8) cm $`\times `$ 0.0873(9) cm $`\times `$ 0.238(2) cm. The sample mass was 0.00832(3) g, yielding a density of 6.3(1) g/cm3. Note that this is somewhat different from the density of 6.09(1) g/cm3 for YNi211B2C calculated using lattice constants obtained from x-ray/neutron measurements.Lynn ; Chakoumakos ; Yelon ; Buchgeister ; Jaenicke-Roessler ; Ignatov

A number of other measurements do indicate good sample quality. Specifically, a residual resistivity ratio (RRR) of 14.2, a superconducting width of 0.3 K, and a $`T_c`$ of 14.7 K, were observed via resistivity measurements on our specimen, which compare favourably with published values.Cava ; Chakoumakos ; Yelon ; Ignatov ; Rathnayaka ; Hong ; Godart ; Isida ; Andreone ; Movshovich

The experimental apparatus used to detect normal mode frequencies was a DRS Modulus I Resonant Ultrasound Spectrometer, with lead zirconate titanate transducers. This system has previously been used to successfully determine the elastic constants of Sr2RuO4.Paglione\_RUS

Once compiled, the list of normal mode frequencies was fitted, and appropriate elastic constants obtained. The measured dimensions of the sample were used as initial dimensions for the program input, but as the measuring errors in these values were the largest sources of error in this study, the dimensions were left as free parameters in the fit. Thus, 9 free parameters – 6 elastic constants and 3 dimensions – were fitted to the measured frequency data, while quantities such as the mass and density of the sample were held constant.

Two other samples were cut from the same ingot as the RUS sample and polished to facilitate time-of-flight measurements in the  and  directions. Sample thicknesses were measured to be 0.319(2) cm and 0.227(2) cm for the samples oriented along the  and  directions, respectively. In all cases, transducers were bonded to samples using phenyl salicylate and silicone for room and low temperature measurements respectively. Time-of-flight measurements were made at room temperature (300 K) and low temperature (2 K) using 20 and 30 MHz LiNbO3 transducers and a home-built spectrometer,Lupien\_PhD  which utilizes the phase comparison technique at frequencies between 20 and 500 Mhz. This apparatus has previously been used in a successful study of ultrasonic attenuation in Sr2RuO4.Lupien\_PE

Using the RUS technique, a total of 68 normal modes were measured between 0.5 MHz and 4.5 MHz. The frequencies of these measured modes, along with the frequencies of the calculated modes from the best fitting run and the percent difference between them, are shown in Table 1. The corresponding frequency spectrum is shown in Fig. 2, with a representative portion of the spectrum inset in more detail. In the full-range view, a large amount of noise was measured between approximately 0.7 MHz and 1.1 MHz, and is attributed to frequency characteristics of the measuring apparatus itself. A prototype RUS apparatus, with PVDF (polyvinylidene fluoride) film transducers,Paglione\_MSc  was used to verify that no normal mode frequencies of the sample lay in this range. In addition, the fits indicate no missing modes in this region.

The best fit to the data had a root-mean-square error between measured and calculated frequencies of 0.38%, indicative of a good fit. The fitted dimensions were 0.0636(6) cm $`\times `$ 0.0878(8) cm $`\times `$ 0.236(2) cm, which agree with measurements, also providing confirmation of a good fit. Note that only two modes were not detected, but were determined to be existent through the fitting procedure. Also, the second mode was not weighted in the fit, a common practice among the first 1-2 modes in RUS measurements.Migliori

The elastic constants determined via the best fit are listed in Table 2. Error values are obtained from the largest of either an estimation of the quality of the fit, as generated by the fitting program, or from the variation in *C*ij values over the range of densities allowed by the error in the measured dimensions. Specifically, the errors in *C*44 and all three fitted dimensions come from the error in density, while the rest come from quality-of-fit estimates.

Also included in Table 2 are the elastic constants determined from three different sources using the time-of-flight technique: room temperature measurements from this study, low-temperature measurements done by our groupPaglione\_2K  and low-temperature measurements by Isida *et al.*Isida  The elastic constants from our room temperature and low-temperature time-of-flight measurements were calculated from velocity data as discussed above.Truell  An example of the raw time-of-flight data used to obtain velocities is shown in Fig. 3. The time delay between each echo is equivalent to the total travel time of an ultrasound pulse going through the sample and back; when taken with the thickness of the sample this yields the velocity of sound for that particular direction of propagation/polarization combination.

The measured and calculated sound velocities with the appropriate polarization and direction are listed in Table 3. Note that whenever the density of YNi2B2C was required in a calculation, $`\rho =6.3(1)`$ g/cm3 was used for measurements performed on our specimen, and $`\rho =6.05`$ g/cm3 was used for the measurements by Isida *et al.*

Upon inspection of the data shown in Tables 2 and 3, it is immediately apparent that while the values obtained via both the RUS and time-of-flight methods of our group agree within error, the values published by Isida *et al.*Isida , with the exception of *C*66, are between 16% and 36% lower than the RUS results and both sets of time-of-flight results from our group. Reasons for this disagreement are not clear, but agreement of RUS and time-of-flight measurements on samples cut from the same specimen may suggest the significance of differences in sample composition between the two groups.

The RUS values of the elastic constants give a bulk modulus for YNi2B2C of 1.8(1)$`\times `$1012 dynes/cm2. This is in the range of the isothermal bulk modulus values of 1.2(1)$`\times `$1012 dynes/cm2 and 2.00$`\times `$1012 dynes/cm2 previously measured for this material by x-ray diffraction,Bud'ko ; Meenakshi  which gives us added confidence in our results. Furthermore, the Debye temperature was calculated via Eq. 2, yielding a value of 525(10) K, which falls within the wide range of published values of the Debye temperature for YNi2B2C: 310(20) KAndreone  from a resistivity measurement, and 415 K,Godart  489(5) KMovshovich  and 537 KHong  from specific heat measurements.

In conclusion, all six independent elastic constants of YNi2B2C have been determined using the resonant ultrasound spectroscopy method. Good agreement was found between these results and the elastic constants measured using the time-of-flight technique on a sample from the same ingot. A calculation of the bulk modulus and an estimation of the Debye temperature of YNi2B2C were performed with the elastic constants obtained in this study, and the results were found to be consistent with the range of values found in the literature for these quantities.

###### Acknowledgements.

This work was supported by the Canadian Institute for Advanced Research and funded by NSERC. J.P. also acknowledges the support of the Walter C. Sumner Foundation.


text:
# Diffusive Nuclear Burning in Neutron Star Envelopes

## 1 Introduction

The increasing number of observations of isolated neutron star (NS) atmospheres (see Pavlov, Zavlin & Sanwal 2002 for a recent review) and resulting constraints on their surface composition has highlighted the need to consider the role of nuclear physics during their cooling phase. The initial composition of the outer layer is neither known nor constrained by the theory of supernova explosions since the amount of matter needed to affect the outer envelope is small ($`10^{20}M_{}`$ for the photosphere). This much material can easily fallback and contaminate the outer atmosphere (Woosley & Weaver 1995). Even if the fallback material consists of heavy elements, spallation (see Bildsten, Salpeter & Wasserman 1992) will occur at some stage in the fallback, creating a plethora of lighter elements that rapidly gravitationally separate. Current evidence suggests that young pulsars ($`<10^{45}\mathrm{yrs}`$) possess photospheres of either hydrogen or helium, whereas older pulsars ($`>10^5\mathrm{yrs}`$) appear to have envelopes of heavier, more uniformly opaque elements (Yakovlev, Kaminker & Gnedin 2001). This may indicate an evolution of the outer envelope on the timescale of $`10^{45}\mathrm{yrs}`$, something we would like to understand better.

As part of our ongoing work on this problem, we present here a mechanism of nuclear processing of the outer envelope which we call diffusive nuclear burning (DNB). The basic picture is simple and was alluded to by Chiu & Salpeter (1964) as a mode of nuclear burning. Rosen (1968) carried out the initial calculations of the effect for surface temperatures in the range of $`(510)\times 10^6\mathrm{K}`$. Their work was in the limit where the diffusion time to larger depths was the rate-limiting step.

In this paper, we will consider much lower values for the surface temperature, where diffusion is not the limiting step. Consider a NS envelope that consists of hydrogen above a layer of carbon (or any other proton capturing material) in diffusive equilibrium (see illustrative Figure 1). For simplicity, we do not calculate the effect of an intervening helium layer in this paper, though we estimate its possible impact in the conclusions. As long as there has been adequate time to reach a diffusive equilibrium, the separation between the hydrogen and carbon is not strict and a diffusive tail of hydrogen penetrates deep into the carbon layer. The temperature rises with depth and can increase by 1-2 orders of magnitude 10 meters beneath the photosphere. At this location, the lifetime of a proton is very short. The convolution of the increasing nuclear burning rate and the decreasing hydrogen abundance creates a burning layer at this depth, where hydrogen burning peaks. As hydrogen burns, the diffusive tail is driven slightly away from diffusive equilibrium and a diffusive current is set up, $`𝑱_{\mathrm{DNB}}`$. We will always work in the limit where $`𝑱_{\mathrm{DNB}}`$ is far less than the maximum current that can be supported due to diffusive processes so that diffusive equilibrium remains an excellent approximation.

We find that, even when the photospheric temperatures and densities are so low that the local nuclear lifetime of a proton is in excess of a Hubble time, DNB can proceed on a much faster timescale. In this paper we make some simplifying assumptions that prohibits a direct application to radio pulsars, and thus the extension of the physics to radio pulsars and magnetars will be left to future work. For example, we require a magnetic field low enough so that there is a negligible effect on the temperature profile; $`B<10^9\mathrm{G}`$ (Ventura & Potekhin 2002; Potekhin & Yakovlev 2001). However, our work is applicable to some X-ray transients in quiescence, where the outer envelope consists of pure hydrogen, even though the depth of the hydrogen layer is unknown.

This paper is structured as follows. In § 2 we present the basic equations and microphysics governing the thermal and compositional structure of a NS envelope in diffusive equilibrium. The derivation of the electric field strength from the envelope structure equations is in § 3. Using the derived electric field strength, we numerically solve the thermal and compositional structure of the NS envelope. The analytic solution to the NS envelope structure is given in § 4. In § 5 we discuss the relevant hydrogen burning reactions that drive the NS envelope evolution, present numerical and analytic solutions for $`𝑱_{\mathrm{DNB}}`$, and apply our theory to the X-ray transients, specifically Cen X-4. We summarize our findings in § 6.

## 2 Atmospheric Structure and Microphysics

We only consider the top $`10^4\mathrm{cm}`$ of the envelope, where the density is $`<10^{10}\mathrm{g}\mathrm{cm}^3`$. Since the thickness of the envelope is much smaller than the NS radius ($`R10`$ km), we presume a plane parallel atmosphere with a constant downward gravitational acceleration, $`g=GM/R^2`$. We neglect all relativistic corrections. The structure of the NS envelope is determined by the equations of hydrostatic balance,

$`{\displaystyle \frac{dP_i}{dr}}`$ $`=`$ $`n_i\left(A_im_pgZ_ieE\right),`$ (1)
$`{\displaystyle \frac{dP_e}{dr}}`$ $`=`$ $`n_e\left(m_eg+eE\right),`$ (2)
where $`P_i`$, $`n_i`$, $`A_i`$, $`Z_i`$ are the pressure, number density, atomic number and charge of the $`i`$’th ion species and $`E`$ is the upward pointing electric field. The thermal structure is determined by the heat diffusion equation for a constant flux, $`acT_e^4/4`$, where $`T_e`$ is the effective temperature,

$$\frac{dT}{dr}=\frac{3\kappa \rho }{16T^3}T_e^4,$$
(3)
where $`\kappa `$ is the opacity. Finally, we demand charge neutrality,

$$n_e=\underset{i}{}Z_in_i.$$
(4)
Summing the pressure equations with the charge neutrality constraint recovers the familiar form of hydrostatic balance, $`dP/dr=\rho g`$.

In the outermost layers of a NS envelope, the nuclei are fully ionized and the radiative opacities are set by free-free absorption and Thomson scattering,

$$\kappa _{\mathrm{Th}}=\frac{n_e\sigma _{Th}}{\rho }\frac{0.4}{\mu _e}\mathrm{cm}^2/\mathrm{g},$$
(5)
where $`\mu _e=A/Z`$ is the electron mean molecular weight. We calculate the free-free opacity by adopting the formalism used in Schatz et al. (1999),

$$\kappa _{\mathrm{ff}}=7.53\times 10^6\frac{\rho _5}{\mu _eT_6^{7/2}}\frac{Z_i^2X_i}{A_i}g_{\mathrm{ff}}(Z_i,T,n_e)\mathrm{cm}^2/\mathrm{g},$$
(6)
where $`T_6=T/10^6\mathrm{K}`$, $`\rho _5=\rho /10^5\mathrm{g}\mathrm{cm}^3`$, $`X_i`$ is the mass fraction, and $`g_{\mathrm{ff}}`$ is the Gaunt factor. For our numerical calculations, we calculate the Gaunt factor using the Schatz et al. (1999) fitting formula \[given by their equation (A2)\]. For our analytic calculations we set $`g_{\mathrm{ff}}=1`$.

At very high densities, the dominant mode of heat transport is conduction by degenerate electrons. For strongly coupled Coulomb plasmas (SCCP), the electron thermal conductivity is determined by electron-ion scattering as given by the Wiedemann-Franz law (Ventura & Potekhin 2002; Brown, Bildsten & Chang 2002).

$$K_c=\frac{\pi ^2k_B^2Tn_e\tau _K}{3m_e^{}},$$
(7)
where $`m_e^{}=ϵ_\mathrm{F}/c^2`$ is the effective electron mass, $`ϵ_\mathrm{F}`$ is the electron Fermi energy, and $`\tau _K`$ is the effective relaxation time. For weak magnetic fields ($`B<10^9`$G), we use the thermal conductivity of electrons in SCCP given by Baiko et al. (1998) with analytic formulae given by Potekhin et al. (1999). In the liquid metal approximation, Brown et al. (2002) showed that this formalism agrees with the calculations of Urpin & Yakovlev (1980) and Itoh et al. (1983).

The transition of the envelope from primarily radiative to primarily conductive takes place over a narrow region called the sensitivity strip (Ventura & Potekhin 2002). Heat transport by radiation and conduction is of the same order in this region. The sensitivity strip also plays the dominant role in determining the temperature profile of the NS envelope (as discussed in § 4, also see Ventura & Potekhin 2002). In the sensitivity strip, the SCCP is always in the liquid metal regime, whose conductivity is fairly well modeled (Brown et al. 2002; Ventura & Potekhin 2002). Thus, we are confident in our results for the thermal profile.

The electron equation of state transitions from an ideal gas EOS to a degenerate EOS. To model this transition in our numerical solutions, we have used the Paczynski (1983) equation of state. In our analytic calculations, we solve our envelope for the two different limits and introduce a free parameter $`\mathrm{\Psi }`$ to join them at a point defined by

$$\frac{P_{\mathrm{e},\mathrm{deg}}}{P_{\mathrm{e},\mathrm{nondeg}}}\mathrm{\Psi },$$
(8)
where $`P_{\mathrm{e},\mathrm{deg}}`$ and $`P_{\mathrm{e},\mathrm{nondeg}}`$ are the nonrelativistic degenerate and nondegenerate electron equations of state respectively. For this calculation we take $`\mathrm{\Psi }=1.5`$, which gives a reasonable approximation to the full numerical solution in the region of the envelope where DNB occurs.

## 3 Elemental Distribution and Electric Field Strength

The composition of the envelope in diffusive equilibrium is found by simultaneously solving the equations of hydrostatic balance for each species and the flux equation. This approximation is excellent once enough time has passed for the diffusive equilibrium to be established (Brown et al. 2002).

We begin by calculating the local electric field. In early studies of white dwarfs and NSs (Rosen 1968; Fontaine & Michaud 1979), the electric field was presumed to be that for an ideal gas of ions and electrons,

$$eE=\frac{A}{Z+1}m_pg.$$
(9)
However, in the degenerate case where $`P_eP`$, the field is

$$eE=\frac{A}{Z}m_pg.$$
(10)
These results are only valid in certain limits, where the plasma consists of one dominant species and where the gas is either very degenerate or very non-degenerate. For multi-ion plasmas and for partial degeneracy, these simple approximations break down. For a more general solution, we calculate the electric field from the ion and electron equations of hydrostatic balance (eq. ), and the flux equation (eq. ). We presume an ideal gas equation of state for the ions and an arbitary electron equation of state. Finally we drop the $`m_eg`$ term.

Let us first consider the isothermal case. We first expand $`P_e/r`$,

$$\frac{n_e}{r}=\left(\frac{P_e}{n_e}\right)^1n_eeE,$$
(11)
and then use the ideal gas equation of state for the ions:

$$\frac{n_i}{r}=\frac{n_i}{kT}(A_im_pgeZ_iE).$$
(12)
Multiplying each ion equation (eq. ) by $`Z_i`$ and subtracting their sum from the electron equation (eq. ) and presuming charge neutrality (eq. ) gives

$$eE=\frac{m_pgn_iA_iZ_i}{n_iZ_i^2+n_ek_BT\left(P_e/n_e\right)^1}.$$
(13)
For the nonisothermal case, we get an additional term in the expansion of $`P_e/r`$ and $`P_i/r`$, e.g.

$$\frac{n_e}{r}=\left(\frac{P_e}{n_e}\right)^1\left(n_eeE\frac{P_e}{T}\frac{T}{r}\right),$$
(14)
and the electric field becomes

$$eE=\frac{n_iZ_i\left(A_im_pg+k_B\left(dT/dr\right)\right)\left(P_e/T\right)\left(P_e/n_e\right)^1\left(dT/dr\right)}{n_iZ_i^2+n_ek_BT\left(P_e/n_e\right)^1}.$$
(15)
This result agrees with the results of Macdonald, Hernanz & Jose (1998) and Althaus & Benvenuto (2000), but we have decoupled our solution of the electric field from the drift velocity and diffusion coefficient. Hence, the electric field is simply a function of local parameters. For the purposes of numerical calculation we use the Paczynski (1983) electron equation of state. For our analytic calculations, we use either the ideal gas or the nonrelativistic degenerate electron equation of state depending on the local conditions of the plasma.

At this point, it is helpful to illustrate a solution for the particular case of a $`\mathrm{H}`$/$`{}_{}{}^{12}\mathrm{C}`$ envelope. In Figure 2, we plot the electric field of the layer as a function of column depth, $`y=P/g`$. We integrate from the photosphere with a background of $`\mathrm{H}`$. At a certain depth we introduce a low carbon abundance and continue integrating inward. As expected, the electric field approximates that of the value given by equation (9) at the extreme limits where one ion dominates over the other. The slight discrepancy between the value give by equation (9) at the extreme limits is due to degeneracy effects from the Paczynski (1983) fit of the electron equation of state. The flux and structure equations can now be solved self consistently to model any arbitrary envelope in diffusive equilibrium. In Figure 3, we plot the solution for a $`\mathrm{H}`$/$`{}_{}{}^{12}\mathrm{C}`$ envelope. We plot the composition profile in terms of the concentration of hydrogen, $`f_H=n_H/n_t`$, where $`n_t`$ is the total number of ions.

## 4 Approximate Solutions for Temperature and Composition

The temperature and composition profiles can be represented well with approximate analytic solutions that illuminate the physics. Hernquist & Applegate (1984) and Ventura & Potekhin (2002) have presented analytic thermal calculations of the NS envelope, but their solutions are too complex for our analytic estimates. Hence, we present simplified approximate solutions in § 4.1. For an envelope in diffusive equilibrium, Fontaine & Michaud (1979) have presented an excellent approximate solution for the composition profile. In § 4.2, we extend their technique to the degenerate regime.

### 4.1 Analytic Temperature Profile

Roughly speaking, the NS envelope can be divided into three zones, a radiative outer zone, a sensitivity strip (where opacity changes from primarily radiative to conductive), and a conductive isothermal interior (also see Ventura & Potekhin 2002). In the radiative, nondegenerate outer zone, the opacity is mostly determined by free-free absorption. With an ideal gas equation of state, $`yg=\rho k_BT/\mu m_p`$, where $`\mu =A/(Z+1)`$ is the mean molecular weight of the plasma, the solution to the heat diffusion equation and hydrostatic balance is

$$T(y)=1.58\times 10^6y^{4/17}\left(\mu Zg_{14}T_{e6}^4\mu _e^2\right)^{2/17}\mathrm{K},$$
(16)
where for brevity we have written $`g_{14}=g/(10^{14}\mathrm{cm}\mathrm{s}^2)`$ and $`T_{e6}=T_e/(10^6\mathrm{K})`$. For a radiative nondegenerate envelope consisting of multiple layers of ions, the solution is modified via the introduction of constants of integration at the boundaries. If the boundaries are sufficiently far away in pressure from the sensitivity strip, the temperature profile of the underlying layers will approximate that of a pure component.

In the sensitivity strip, which is defined by $`\kappa _{\mathrm{rad}}(y_{\mathrm{ss}})=\kappa _{\mathrm{cond}}(y_{\mathrm{ss}})`$, the opacity transitions from primarily radiative to primarily conductive. Here we follow the formalism of Ventura & Potekhin (2002). To find the sensitivity strip, we take the approximate expression for conductivity (Ventura & Potekhin 2002),

$$K_c=2.3\times 10^{15}\frac{T_6}{\mathrm{\Lambda }Z}\frac{\chi ^3}{1+\chi ^2}\frac{\mathrm{ergs}}{\mathrm{cm}\mathrm{s}\mathrm{K}},$$
(17)
where $`\chi =p_F/m_ec(\rho _6/\mu _e)^{1/3}`$ is the relativity parameter, $`\mu _e=A/Z`$, and $`\mathrm{\Lambda }`$ is the Coulomb logarithm. We set $`\mathrm{\Lambda }=1`$ for simplicity’s sake and $`1+\chi ^21`$ since the relativity parameter is small in the nonrelativistic regime. Equating the two opacities, equations (6) and (17), we find that the sensitivity strip is defined by the condition,

$$y=11.2\frac{A^{1/3}T_6^{17/6}\mu _e^{2/3}}{Z^{1/3}g_{14}\mu }\mathrm{g}/\mathrm{cm}^2.$$
(18)
Taking this condition and inserting the solution for the radiative zone, we obtain expressions for the column depth of the sensitivity strip, $`y_{\mathrm{ss}}`$,

$$y_{\mathrm{ss}}=6.86\times 10^4\frac{ZT_{e6}^4\mu _e}{g_{14}^2\mu ^2}\mathrm{g}/\mathrm{cm}^2,$$
(19)
and temperature,

$$T_{\mathrm{ss}}=2.2\times 10^7\left(\frac{Z^3T_{e6}^{12}}{g_{14}^3\mu ^3}\right)^{2/17}\mathrm{K},$$
(20)
as a function of effective temperature, gravity, and composition. Beyond the sensitivity strip ($`y>y_{\mathrm{ss}}`$), we solve the constant flux equation using the conductive opacity,

$$T\frac{dT}{dy}=2.2\times 10^{16}\frac{ZT_{e6}^4\mu _e}{\rho ^2},$$
(21)
when $`1+\chi ^21`$. Integration of this equation is relatively straightforward given an equation of state. There are two cases that we need to consider. For the nondegenerate regime, the solution for $`y>y_{\mathrm{ss}}`$ is

$$T=T_{\mathrm{ss}}\mathrm{exp}\left[1.62\times 10^4\frac{Z\mathrm{\Lambda }T_{e6}^4\mu _e}{g_{14}^2\mu ^2}\left(y^1y_{\mathrm{ss}}^1\right)\right]\mathrm{K},$$
(22)
where we have joined it with the radiative solution at the sensitivity strip.

The electrons become more degenerate as the column increases. We define the degenerate boundary with our analytic degeneracy condition (eq. ), which gives the condition

$$y=22\frac{T_6^{5/2}\mathrm{\Psi }^{5/2}\mu _e^{5/3}}{g_{14}\mu ^{5/2}}\mathrm{g}/\mathrm{cm}^2,$$
(23)
where we have assumed that $`y_{\mathrm{deg}}y_{\mathrm{ss}}`$. Using our solution for $`T`$ beyond the sensitivity strip (eq. ), we find a system of equations for $`T_{\mathrm{deg}}`$ and $`y_{\mathrm{deg}}`$, which are solved numerically.

The solution beyond the degenerate boundary is

$$T^2T_{\mathrm{deg}}^2=1.37\times 10^{16}T_{e6}^4\frac{Z}{g_{14}^{6/5}\mu _e}\left(y_{\mathrm{deg}}^{1/5}y^{1/5}\right)\mathrm{K}^2.$$
(24)
Taking the limit of equation (24) as $`y\mathrm{}`$ gives an analytic expression for the core temperature, $`T_c`$, as a function of $`T_e`$, which reproduces previous scalings (Gudmundsson, Pethick & Epstein 1983; Potekhin, Chabrier & Yakovlev 1997). Ventura & Potekhin (2002) present a more generalized solution in terms of the relativity parameter $`\chi `$.

In Figure 3 we present the approximate and numerical solutions for the temperature and composition which we discuss in the next section. For the thermal profile, the agreement is excellent (to within 10%). It also agrees well with the analytic solution given by Ventura & Potekhin (2002).

### 4.2 Analytic Composition Profiles

We find analytic solutions for the composition when one ion species is dominant and the electrons are either nondegenerate or very degenerate. For the nondegenerate case, we first derive the result of Fontaine & Michaud (1979). We start with the hydrostatic balance equations for the ions, where “1” is the background and “2” is the trace,

$`{\displaystyle \frac{dP_1}{dr}}=n_1\left(A_1m_pgZ_1eE\right),`$ (25)
$`{\displaystyle \frac{dP_2}{dr}}=n_2\left(A_2m_pgZ_2eE\right).`$ (26)
For the ions, $`P_i=n_ik_BT`$. Dividing each pressure equation above by their respective $`P_i`$ and subtracting “1” from “2”, we have:

$$k_BT\frac{d\mathrm{ln}\left(P_2/P_1\right)}{dr}=\left(A_2A_1\right)m_pg+\left(Z_2Z_1\right)eE,$$
(27)
which in the trace limit ($`n_2n_1`$) gives

$$\mathrm{ln}\left(\frac{P_2}{P_1}\right)=\mathrm{ln}\left(\frac{n_2}{n_1}\right)\mathrm{ln}\left(\frac{n_2}{n_1+n_2}\right)=\mathrm{ln}f_2,$$
(28)
where $`f_2`$ is number fraction of the trace.

In the nondegenerate regime, the electric field is given by equation (9). Using this result and using the trace approximation so that $`AA_1`$ and $`ZZ_1`$, equation (27) becomes

$$\frac{d\mathrm{ln}f_2}{dr}=\frac{A_1m_pg}{(Z_1+1)k_BT}\left[\frac{A_2}{A_1}(Z_1+1)Z_21\right]$$
(29)
Using hydrostatic balance $`dP/dr=\rho gA_1m_pgn_1`$ and $`P(Z_1+1)n_1kT`$, we have the elegant result,

$$\frac{d\mathrm{ln}f_2}{d\mathrm{ln}P}=\left[\frac{A_2}{A_1}(Z_1+1)Z_21\right].$$
(30)
The concentration of a trace species is a power law of pressure with an exponent, $`\delta =A_2(Z_1+1)/A_1Z_21`$. It is widely applicable because it gives the concentration in terms of a pressure contrast,

$$f_2(P)=f_2(P_\mathrm{b})\left(\frac{P}{P_\mathrm{b}}\right)^\delta ,$$
(31)
where $`f_2(P_\mathrm{b})`$ is the concentration at an arbitrary boundary $`P_\mathrm{b}`$.

We now extend this result to the highly degenerate regime, where the electron pressure dominates over ion pressure, and the electric field is

$$eE=\frac{m_pg}{\mu _e}=\frac{A}{Z}m_pg.$$
(32)
Placing this into equation (27) gives

$$\frac{d\mathrm{ln}f_2}{dP}=\frac{A_1m_p}{\rho (P,T)k_BT(P)}\left[\frac{A_2}{A_1}\frac{Z_2}{Z_1}\right].$$
(33)
We solve for $`\rho (P,T)`$ via the electron equation of state and $`T(P)`$ from the flux equation. This solution is best when the electrons are highly degenerate. However, this is also the regime where the temperature begins to be isothermal (Ventura & Potekhin 2002), allowing us to replace $`T(P)`$ by $`T_{\mathrm{deg}}`$, where $`T_{\mathrm{deg}}`$ is the temperature at the degenerate boundary, giving

$$\frac{d\mathrm{ln}f_2}{dP}=\frac{A_1m_p}{\rho _{\mathrm{deg}}\left(P/P_{\mathrm{deg}}\right)^{1/\gamma }k_BT_{\mathrm{deg}}}\left[\frac{A_2}{A_1}\frac{Z_2}{Z_1}\right],$$
(34)
where $`\gamma `$ is the polytropic index of the degenerate electron equation of state, $`P\rho ^\gamma `$. The condition of degeneracy (eq. ) determines $`P_{\mathrm{deg}}`$ and $`\rho _{\mathrm{deg}}`$.

The electrons are nonrelativistic at the degenerate boundary so $`\gamma =5/3`$ and we integrate (eq. ) to get

$$f_2=f_{2,\mathrm{deg},0}\mathrm{exp}\left\{\eta \frac{A_1m_pP_{\mathrm{deg}}}{\rho _{\mathrm{deg}}k_BT_{\mathrm{deg}}}\left[\frac{Z_2}{Z_1}\frac{A_2}{A_1}\right]\left[\left(\frac{P}{P_{\mathrm{deg}}}\right)^\eta 1\right]\right\}$$
(35)
where $`f_{2,\mathrm{deg},0}`$ is the concentration where the envelope becomes degenerate and $`\eta =1\frac{1}{\gamma }=0.4`$. Note that $`\eta `$ is positive for all degenerate equations of state: the concentration always decreases exponentially with increasing pressure. Our solution reproduces the Boltzmann solution for a particle experiencing a upward pointing force $`F`$ in the isothermal limit, $`f_2\mathrm{exp}\left(Fr/kT\right)`$ since $`rP/\rho P^\eta `$. Thus, in the degenerate regime the concentration falls faster than the power law relation in the nondegenerate regime. In Figure 3, we compare the numerical solution with our analytic calculation. The agreement is poor for higher densities due to the increasing importance of relativistic equation of state, but gives reasonable agreement in the burning layer $`y_{\mathrm{burn}}10^6\mathrm{g}\mathrm{cm}^2`$.

## 5 Hydrogen Burning in Diffusive Equilibrium

The relevant nuclear processes that consume hydrogen depend on the temperature and composition of the underlying matter. We consider the case of a NS which has a pure hydrogen layer on top of a layer of proton capturing material (i.e. carbon, nitrogen, oxygen, etc) in diffusive equilibrium. We have not yet calculated the impact of an intervening He layer, but we think it will be small as long as it does not penetrate to a depth greater than the burning layer.

For simplicity, we first consider carbon. For hydrogen on carbon, there are only two possible nuclear reactions: $`{}_{}{}^{12}\mathrm{C}(\mathrm{p},\gamma ){}_{}{}^{13}\mathrm{N}`$ and $`\mathrm{p}(\mathrm{p},e^++\nu _e)\mathrm{D}`$. The D produced in the p-p capture undergoes $`\mathrm{D}(\mathrm{p},\gamma ){}_{}{}^{3}\mathrm{He}`$ which in general is the endpoint of this process. Since we always have a region where T $`>10^7\mathrm{K}`$, the reaction rate for p-p capture is very slow compared to the proton capture rate onto carbon. The $`{}_{}{}^{13}\mathrm{N}`$ decays to $`{}_{}{}^{13}\mathrm{C}`$ which has $`A/Z=13/62.167`$, greater than the local $`A/Z=2`$ for pure carbon. The $`{}_{}{}^{13}\mathrm{C}`$ will thus sink through the outer layers of the NS, and will not reside long enough in the burning layer to facilitate a catalytic cycle. Hence, the complete exhaustion of $`\mathrm{H}`$ by this process requires an excess of $`{}_{}{}^{12}\mathrm{C}`$ compared to $`\mathrm{H}`$.

### 5.1 Diffusive Nuclear Burning

We calculate the burning rate presuming that the burning time is slow compared to the time for hydrogen to diffuse down through the carbon to the burning layer. We express this hierarchy of diffusion time to nuclear burning time by studying the $`\mathrm{H}`$ equation of continuity,

$$\frac{n_H}{t}+\mathbf{}\mathbf{}𝑱_H=\frac{n_H}{\tau _H},$$
(36)
where $`\tau _H=\sigma v^1n_C^1`$ is the local lifetime of $`\mathrm{H}`$ to $`{}_{}{}^{12}\mathrm{C}`$ capture. The condition of diffusive equilibrium is $`n_H/t=0`$ or more specifically that $`n_H/t`$ changes only on the timescale associated with depletion of the hydrogen column, which is long compared to both the local nuclear burning time $`\tau _H`$ and the diffusion time and therefore is dropped. With the condition of steady state $`𝑱_H`$ and $`𝑱_{\mathrm{DNB}}`$ are equivalent. Writing the result in one dimension for $`𝑱_H=𝑱_{\mathrm{DNB}}=𝒟n_H/zn_Hv_{\mathrm{dr}}`$, we have

$$𝒟\frac{^2n_H}{z^2}v_{\mathrm{dr}}\frac{n_H}{z}=\frac{n_H}{\tau _H},$$
(37)
, where $`𝒟`$ is the diffusion coefficient and $`v_{\mathrm{dr}}`$ is the drift velocity. We presume that $`𝒟`$ and $`v_{\mathrm{dr}}`$ change slowly with $`z`$. Our presumption that nuclear burning does not affect diffusive equilibrium implies that locally the timescale associated with nuclear burning is much longer than the timescale associated with diffusion to the burning layer, or $`\tau _H\tau _{\mathrm{ion}}=l^2/𝒟`$, where $`l`$ is the ion scale height for $`\mathrm{H}`$in the $`{}_{}{}^{12}\mathrm{C}`$ layer. Therefore, taking $`𝒟`$ from Brown et al. (2002), we have

$$\tau _{\mathrm{ion}}=800\frac{T_6^{0.7}Z_i^{1.3}\rho _5^{0.6}}{A_i^{0.1}g_{14}^2\left(A_iZ_i\right)^2}\mathrm{s},$$
(38)
where $`A_i`$ and $`Z_i`$ are the atomic number and charge of the background proton capturing element.

For proton capturing reactions, the proton lifetime (Clayton 1983) is

$$\tau _\mathrm{H}^1=2.45\times 10^{16}\left(\frac{S_0}{1\mathrm{keV}\mathrm{barns}}\right)\frac{\rho X_i}{A_i}T_6^{2/3}\mathrm{exp}\left(BT_6^{1/3}\right)\mathrm{yrs}^1$$
(39)
where $`X_i`$ and $`A_i`$ is the mass fraction and atomic number of the proton capturing element and $`S_0`$ and $`B`$ are parameters determined by the proton capturing element. We derive an analytic expression for the condition of diffusive equilibrium by expanding equation (39) about $`T_6=40`$,

$$\tau _H=1.32\times 10^5\rho _5^1\left(\frac{T_6}{40}\right)^{14}\mathrm{s}.$$
(40)
Comparing equations (38) and (40), our condition for carbon is

$$T_6<45.5\rho _5^{0.11},$$
(41)
for $`g_{14}=2`$.

If the conditions for DNB are fulfilled for a sufficiently large pressure range, the total burning rate of hydrogen, defined as $`\zeta _H=𝑱_{\mathrm{DNB}}(\mathrm{H}/{}_{}{}^{12}\mathrm{C}\mathrm{boundary})`$ is

$$\zeta _H=\frac{y_H}{\tau _{\mathrm{col}}}=\frac{n_Hm_p}{\tau _H(n_H,n_C,T)}𝑑z,$$
(42)
where $`y_H`$ is the integrated column of hydrogen, $`\tau _{\mathrm{col}}`$ is the characteristic time for that column to be consumed and $`n_H`$ and $`n_C`$ are the local number density of hydrogen and carbon respectively. This equation is incorporated into our flux and structure equations and solved simultaneously, yielding $`\tau _{\mathrm{col}}`$ as a function of our stellar properties and $`y_H`$. We plot $`d\zeta _H/d\mathrm{lg}_{10}y`$ for an envelope with $`y_H=100\mathrm{g}\mathrm{cm}^2`$ in Figure 4. In the bottom graph we show the local concentration and temperature. Due to the inverse dependence of the local hydrogen burning rate on temperature and concentration, the burning layer is concentrated in a narrow pressure range at a depth where the $`\mathrm{H}`$ abundance is small. The diffusion timescale to this depth is $`\tau _{\mathrm{ion}}=10^4\mathrm{s}`$, which is much shorter than the nuclear timescale of $`1.5\times 10^7\mathrm{s}`$. *This observation is the essence of DNB, namely, that practically all the burning occurs in the exponentially suppressed diffusive hydrogen tail.*

In Figure 5, we plot the characteristic time, $`\tau _{\mathrm{col}}`$, and characteristic mass burning rate, $`\dot{M}_{\mathrm{DNB}}=4\pi R^2y_H/\tau _{\mathrm{col}}`$, for $`g_{14}=2`$ and fixed core temperature, $`T_c`$. For a central temperature of $`5\times 10^7\mathrm{K}`$, our calculations are valid only up to $`3\times 10^7\mathrm{g}\mathrm{cm}^2`$. Above this column, the assumption of constant flux breaks down since the energy released from DNB is comparable to the flux.

Two results are obvious from this figure. First the relation between lifetime (and hence mass burning rate) and column is related via a simple power law between $`\tau _{\mathrm{col}}`$ and $`y_H`$. We derive this power law in § 5.2. Secondly, the lifetime of the entire envelope is set by the lifetime of the photosphere.

The location of the burning layer must satisfy our derived diffusive equilibrium condition, (eq. ) for our calculation to be valid. In Figure 6, we plot various models against our condition for DNB in diffusive equilibrium (eq. ). The burning layers is the region between the two vertical lines in each model. For models with $`T_e>10^6\mathrm{K}`$ (or $`T_c>6\times 10^7\mathrm{K}`$), our present calculation is not valid as DNB will occur in the slow diffusion limit. We will expand on this point and study the slow diffusion limit in a future paper.

### 5.2 Approximate Solution

The competition between the falling $`\mathrm{H}`$ abundance and rising temperature with depth forces much of the hydrogen burning to occur in a narrow zone. As a result, we can solve the burning integral (eq. ) over the burning layer by approximating it with the method of steepest descents.

Since DNB occurs deep below the hydrogen layer, $`X_i1`$, equation (42) then becomes

$$\zeta _H=\frac{y_H}{\tau _{\mathrm{col}}}=2.45\times 10^{16}A_i^2\left(\frac{S_0}{1\mathrm{keV}\mathrm{barns}}\right)\rho yf_HT_6^{2/3}\mathrm{exp}\left(BT_6^{1/3}\right)d\mathrm{ln}(y)\mathrm{yrs}^1$$
(43)
where $`f_H`$ is the local number concentration of hydrogen, which decreases with increasing $`y`$. Note that we have made a change of variables of $`dz=y\rho ^1d\mathrm{ln}(y)`$.

This integral should be evaluated in both the nondegenerate, radiative regime and the degenerate, conductive regime. However, since the majority of the burning occurs in the degenerate, conductive regime, we will evaluate the integral there by the method of steepest descents (Clayton 1983),

$$g(x)\mathrm{exp}(h(x))𝑑xg(x_0)\mathrm{exp}(h(x_0))\sqrt{\frac{2\pi }{h^{\prime \prime }(x_0)}},$$
(44)
where $`g(x)`$ is a slowly varying function function of $`x`$, $`h(x)`$ is a peaked function, and $`x_0`$ is the solution of $`h^{}(x_0)=0`$, the location of the peak of the burning rate. In order to carry out this approximation for equation (43), we make the following identifications,

$`\mathrm{exp}(h(x))`$ $``$ $`y\rho f_H\mathrm{exp}\left(BT_6^{1/3}\right),`$
$`g(x)`$ $``$ $`T_6^{2/3}.`$ (45)
For the degenerate, conductive regime, the peaked function is

$$h\left[\mathrm{ln}(y)\right]=\frac{8}{5}\mathrm{ln}(y)BT_6^{1/3}\eta \frac{A_im_pgy_{\mathrm{deg}}}{\rho _{\mathrm{deg}}k_BT_{\mathrm{deg}}}\left(\frac{1}{Z_i}\frac{1}{A_i}\right)\left[\left(\frac{y}{y_{\mathrm{deg}}}\right)^\eta 1\right].$$
(46)
The last term comes from the analytic solution for the composition profile (eq. ), where the background is the proton capturing element and the trace is hydrogen. The temperature, $`T`$, in the degenerate regime is given by equation (24). Since the burning peak $`y_0`$ is defined by a transcendental equation, $`h^{}(\mathrm{ln}(y_0))=0`$, we solve it numerically and then calculate $`h^{\prime \prime }(\mathrm{ln}(y_0))`$. For a total column of hydrogen of $`y_H=100\mathrm{g}\mathrm{cm}^2`$, $`T_e=8\times 10^5\mathrm{K}`$, and $`g_{14}=2`$, we find by the method of steepest descents, $`\tau _{\mathrm{col}}=468\mathrm{yrs}`$, which compares well with the numerical answer of $`\tau _{\mathrm{col}}=428\mathrm{yrs}`$. The burning layer is centered around $`y_0=7.2\times 10^5\mathrm{g}\mathrm{cm}^2`$ with a local temperature of $`T(y=y_0)=2.6\times 10^7\mathrm{K}`$ and central core temperature $`T_c=3.3\times 10^7\mathrm{K}`$.

Power law scalings of $`\tau _{\mathrm{col}}`$ with $`y_H`$, $`T_e`$ and $`g`$ can also be determined from this integral. The simplest power law scaling involves $`y_H`$, which we determine as follows. The rate, $`y_H/\tau _{\mathrm{col}}`$, is directly proportional to $`f_H`$ as in equation (43). From equation (31) and (35), $`f_Hf_{\mathrm{H},\mathrm{deg},0}=(y_{\mathrm{deg}}/y_H)^\delta `$, where $`f_{\mathrm{H},\mathrm{deg},0}`$ is the hydrogen number fraction at the degenerate boundary and $`\delta =A_2(Z_1+1)/A_1Z_21=17/12`$ for H on $`{}_{}{}^{12}\mathrm{C}`$. Therefore, $`\tau _{\mathrm{col}}y_H^{1+\delta }`$, which is confirmed from comparisons to the numerical results. This scaling is accurate as long as $`y_H<y_{\mathrm{deg}}`$.

The scaling of $`\tau _{\mathrm{col}}`$ against $`T_{e6}`$ and $`g_{14}`$ can be derived from a power-law expansion of burning integral (eq. ) in terms of $`T`$. The expansion of the burning integral (eq. ) gives

$$\frac{y_H}{\tau _{\mathrm{col}}}f_HT_{0,6}^{2/3}\mathrm{exp}\left(BT_{0,6}^{1/3}\right),$$
(47)
where $`T_{0,6}=T_0/10^6`$ and $`T_0`$ is the temperature of the burning layer. Since we are in the degenerate conductive regime, $`T_0T_c`$. The location of the burning region, $`y_0`$ and $`\rho _0`$, scale weakly with $`T`$ compared to the temperature of the burning region, $`T_0`$. Taking the analytic solution to the temperature profile (eq. ) and degeneracy point (eq. ), the power law dependence between $`T_c`$, $`T_{e6}`$ and $`g_{14}`$ is $`T_0T_cT_{e6}^2g_{14}^{1/2}`$. Therefore, the scaling of $`f_H`$ in $`{}_{}{}^{12}\mathrm{C}`$ goes as,

$$f_H\mathrm{exp}\left\{\eta \frac{m_pP_{\mathrm{deg}}}{\rho _{\mathrm{deg}}k_BT_0}\left[\left(\frac{P}{P_{\mathrm{deg}}}\right)^\eta 1\right]\right\}.$$
(48)
Expanding this equation at $`T_0=10^8K`$ gives $`f_H(T_0/10^8)^{0.475}g_{14}^{0.285}T_{e6}^{0.95}g_{14}^{0.0475}`$. Expanding the exponential in equation (47) for carbon ($`B=136.96`$) gives $`\mathrm{exp}\left(BT_c^{1/3}\right)T_c^{9.8}T_{e6}^{19.7}g_{14}^{4.9}`$. Thus, the scaling for the rate is $`y_H/\tau _{\mathrm{col}}T_{e6}^{19.3}g_{14}^{4.52}`$. Hence, the dominant scaling comes from the exponential dependence on temperature of the nuclear burning rate with small corrections from the other factors. Putting all the scalings together with our analytic calculation, we find

$$\tau _{\mathrm{col}}6.5\left(\frac{y_H}{100\mathrm{g}\mathrm{cm}^2}\right)^{5/12}T_{e6}^{19.3}\left(\frac{g_{14}}{2}\right)^{4.52}\mathrm{yrs}.$$
(49)
In Figure 7, we compare the scaling found numerically with this analytic result and find reasonable agreement. The agreements for the scalings of $`y_H`$ and $`g_{14}`$ and the numerical results are good. However, the agreement for the scaling of $`T_e`$ and numerical results diverge for low $`T_e`$. This is due to the changing power law with respect to expanding around different values of $`T_c`$. Expanding around a value of $`T_c`$ appropriate for these lower values of $`T_e`$ resolves this discrepancy.

These scalings can be calculated for other proton capturing elements. The scaling of $`\tau _{\mathrm{col}}`$ and $`y_H`$ remains the same with the exponent $`1\delta _i`$, where $`\delta _i`$ takes on different values for the differing backgrounds $`(Z_i,A_i)`$. We can use the same strategy as earlier, expanding our results in $`T_0`$ and inserting the scalings for $`T_{e6}`$ and $`g_{14}`$ afterwards. Hence, $`f_H`$ has the same scaling as before, $`f_H(T_0/10^8)^{0.475}g_{14}^{0.285}`$. However, the dominant scaling from the exponential dependence of temperature on the burning rate changes. Hence we choose to ignore the scalings associate with $`f_H`$. In terms of the exponential factor $`B_i`$, this dependence becomes $`\mathrm{exp}\left(B_iT_{0,6}^{1/3}\right)\left(T_0/10^8\mathrm{K}\right)^{B_i/14}T_{e6}^{B_i/7}g_{14}^{B_i/28}`$. Since $`f_H`$ does not scale strongly with $`T_0`$ compared to the nuclear reaction rate (especially with heavier elements where $`B_i`$ gets progressively larger), we ignore the scaling of $`f_H`$ with $`T_0`$ and $`g_{14}`$. To put this another way, the value of the number density, $`f_H`$, at the burning layer does not change as drastically as the value of the nuclear timescale at the burning layer with $`T_0`$. Putting these scalings all together, we have

$$\tau _{\mathrm{col}_\mathrm{i}}=\tau _{\mathrm{col}_\mathrm{i},0}\left(\frac{y_H}{100\mathrm{g}\mathrm{cm}^2}\right)^{1+\delta _i}T_{e6}^{B_i/7}g_{14}^{B_i/28},$$
(50)
Table 1 lists the values of the scalings, the prefactor, calculated both numerically and analytically, and the largest effective temperature where DNB in diffusive equilibrium is valid. Above $`T_{e,6\mathrm{DNB}}`$, equation (50) is not valid since DNB no longer occurs in diffusive equilibrium.

### 5.3 Application to Cen X-4 in Quiescence

The large number of transiently accreting NS also provide an astrophysical site for our work. The work of Rutledge et al. (1999) (see Bildsten and Rutledge 2001 for an overview) first showed that much of the quiescent emission from these objects was thermal emission from the surface, allowing for many $`T_e`$ measurements (Rutledge et al. 2001a, b, 2000, 1999). Of all these LMXBs, only Cen X-4 is safely in the DNB regime. Its effective temperature is low enough ($`T_e=8.8\times 10^5`$ Rutledge et al. 2001a) so that DNB in diffusive equilibrium is a good approximation for a hydrogen column $`y_H<10^8\mathrm{g}\mathrm{cm}^2`$. The column depth of H on Cen X-4 after an outburst is expected to be $`y_H<10^8\mathrm{g}\mathrm{cm}^2`$; set by ignition conditions of type I X-ray bursts (Brown et al. 2002).

Accretion outbursts were observed in 1969 and 1979 from Cen X-4, but the latter outburst was particularly weak. Presuming no outburst has occurred in the intervening years, this gives an age of 23 years. For a hydrogen column of $`y_H=10^8\mathrm{g}\mathrm{cm}^2`$ after the 1979 outburst, we find that the present parameters of Cen X-4 would be $`y_H2\times 10^6\mathrm{g}\mathrm{cm}^2`$ and $`T_c4\times 10^7`$ K if there has been no accretion in quiescence. A much larger H column ($`y_H2\times 10^7\mathrm{g}\mathrm{cm}^2`$) is needed to burn matter via DNB at the supplied rate if the NS is accreting in quiescence at a rate $`\dot{M}_q=10^{13}M_{}\mathrm{yr}^1`$.

Presuming little accretion in quiescence and fixed $`T_c`$, we evolve $`y_H`$ and the resulting flux as a function of the time since the last outburst. The age of the envelope or time since the last outburst, $`t_{\mathrm{age}}`$ is the age derived from the scaling of $`\tau _{\mathrm{col}}`$ with $`y_H`$ in the analytic solution,

$$t_{\mathrm{age}}=\left(\frac{dy_H}{dt}\right)^1𝑑y_H=\frac{12}{5}\tau _{\mathrm{col}}.$$
(51)
In Figure 8, we plot the resultant flux from a total integrated column of $`\mathrm{H}`$, $`y_H`$, against the lifetime of that column for fixed core temperature, $`T_c`$. Evolution via DNB predicts that the outgoing flux will drop by $`12\%`$ over a twenty year period roughly 100 years after the outburst. This drop is due to the consumption of $`\mathrm{H}`$ that makes up the sensitivity strip of the *$`\mathrm{H}`$ layer*. At present (twenty years after the outburst) the flux will vary by only 3% over a ten year timescale. The associated $`\dot{M}_{\mathrm{DNB}}=2\times 10^{15}M_{}\mathrm{yr}^1`$ for this time after the outburst. Observations of Cen X-4 have placed constraints on the variation of $`T_e`$ over a five year timescale of $`<10\%`$ (Rutledge et al. 2001a). Observational confirmation of DNB on Cen X-4 is therefore unlikely.

## 6 Conclusions and Discussions

We have shown that diffusive nuclear burning (DNB) is an effective mechanism that burns surface hydrogen on a NS on astrophysically relevant timescales. Our numerical and analytic solutions for the rates of DNB are in the limit where the H is in diffusive equilibrium with the underlying proton capturing elements. Our work is similar to that originally carried out by Michaud and Fontaine (Michaud et al. 1984; Michaud & Fontaine 1984) for diffusion-induced hydrogen burning in white dwarfs and also by Iben & MacDonald (1985a, b). However, there are several crucial differences: the NS temperatures are much higher and the local gravity much stronger. Hence, we were able to make simplifying assumptions and avoid the active diffusion calculation with nuclear burning performed by Iben & MacDonald (1985a). We have also solved for the concentration profile in the degenerate regime self-consistently, whereas the effects of degeneracy were neglected in Michaud and Fontaine’s paper.

The largest remaining uncertainty in applying our work is the unknown composition of the material underlying the hydrogen. Our assumptions of low magnetic fields and surface temperatures ($`T_e<10^6\mathrm{K}`$ for DNB to be nuclear-rate limited) limit the DNB applications to quiescent NSs in transient low-mass X-ray binaries and millisecond radio pulsars. Little is known about the surface composition of millisecond radio pulsars, though in the absence of our work one would definitely expect H to be dominant on the surface since these objects underwent mass accretion prior to becoming a pulsar. If there is enough carbon in the underlying composition from previous H/He burning during type I X-ray bursts (Schatz et al. 1999), then there is clearly adequate time to burn off the H. A clear detection of H at the photosphere of a millisecond radio pulsar would thus constrain the underlying abundances.

The impact of an intervening helium layer depends on its thickness. A thin layer of helium, which does not penetrate into the burning layer and remains in the nondegenerate regime, enhances the burning rate. This is because the electric field in a nondegenerate helium plasma ($`E=4m_pg/3e`$) is smaller than the electric field in a nondegenerate carbon plasma ($`E=12m_pg/7e`$). Hence the diffusive tail of hydrogen can penetrate more easily into the helium layer onto the underlying layer of proton capturing elements and the number fraction of hydrogen is enhanced at the burning layer. For a thick helium layer which goes beyond the burning layer, DNB would be effectively shut off. However, a thick helium layer would not likely sit on top of a thick proton capturing layer, if the proton capturing elements have the same A/Z as helium and both layers are degenerate. This is because the electric field is no longer a differentiating factor allowing the helium to mix downward with the proton capturing material. This would effectively dilute the abundance of proton capturing elements, thereby reducing the rate, but not shutting off DNB. The timescale for this mixing and its exact impact on the overall burning rate still needs more careful attention.

The recent observations of a pair of X-ray spectral lines on 1E1207.4-5209 (Sanwal et al. 2002; Mereghetti et al. 2002) show that hydrogen is not present on the surface of this young NS (Sanwal et al. 2002; Hailey & Mori 2002). These spectral lines appear to be mid-atomic elements like oxygen or neon (Hailey & Mori 2002) or helium if $`B10^{14}\mathrm{G}`$ (Sanwal et al. 2002). Since the age of the associated SN remnant is 7 kyrs, the mechanism of hydrogen removal must be extremely fast. However, our calculation is not directly applicable to this system. The spindown of 1E1207.4-5209 implies a dipole B-field of $`B3\times 10^{12}\mathrm{G}`$ (Pavlov et al. 2002). This object also has a high temperature of $`T1.41.9\mathrm{MK}`$ at which our assumption of diffusive equilibrium breaks down. However, the physics remains the same and we expect this mechanism to be active in this system.

Since we have explored one special limit of DNB, it is not surprising that we have not found a physical system in which DNB is directly observable. In a future paper we will consider the problem of diffusion limited DNB, which is applicable to systems where the assumption of diffusive equilibrium breaks down ($`T_e>10^6\mathrm{K}`$ for $`\mathrm{H}`$/$`{}_{}{}^{12}\mathrm{C}`$). We will also consider the problem of DNB in highly magnetized sources such as young radio pulsars ($`B10^{1213}\mathrm{G}`$) and magnetars ($`B10^{1415}\mathrm{G}`$).

We thank D. Uzdensky for showing us a nice mathematical trick, and R. Sunyaev for reminding us of the importance of the magnetic field on the thermal structure in radio pulsars. We would also like to thank S. W. Davis, C. J. Deloye, A. Socrates and D. Townsley for discussions and the referee for important clarifications. P.C. would like to thank the Department of Physics and Department of Astronomy at Columbia University, where the early and latter parts of this work was done, for their hospitality. This research was supported by NASA via grant NAG 5-8658 and by the NSF under Grants PHY99-07949 and AST01-96422. L. B. is a Cottrell Scholar of the Research Corporation.