text:
# Inhomogeneous LOFF phase revisited for surface superconductivity

Also at \]L.D. Landau Institute for Theoretical Physics, Chernogolovka, 142432, Russia

## Abstract

We consider 2D surface superconductivity in high magnetic fields parallel to the surface. We demonstrate that the spin-orbit interaction at the surface changes the properties of the inhomogeneous superconducting Larkin-Ovchinnikov-Fulde-Ferrell state that develops above fields given by the paramagnetic criterion. Strong spin-orbit interaction significantly broadens the range of existence of the LOFF phase, which takes the form of periodic superconducting stripes running along the field direction on the surface, leading to the anisotropy of its properties. In connection with experiments by J.H. Schön et al. \[Nature 914, 434 (2001)\] on superconductivity of electrically doped films of the cuprate material CaCuO2, we also discuss this problem for the d-wave pairing to indicate the possibility of a re-orientation transition as the magnetic field direction is rotated in the plane parallel to the surface. Our results provide a tool for studying surface superconductivity as a function of doping.

There has been renewed experimental and theoretical interest in the properties of metallic states localized at a surface, and surface superconductivity (SSC). Surface states (a.k.a. Tamm’s levels) are well-known from the physics of semiconductors. Numerous ARPESLaShell ; Rotenberg  data and STM studies of Friedel oscillationsPetersen  have now proven the existence of 2D metallic bands even at surfaces of metals. The bands are well-separated from the bulk and possess clear-cut 2D Fermi surfaces. For heavy enough elements such Fermi surfaces are split further by strong spin-orbit (SO) interactions. For example, SO energy for electrons at the Fermi level for Au is estimated at 0.1 eV LaShell , while for Li-doped surfaces of Mo and W its value increases up to 0.13 eV and 0.5 eV, correspondingly Rotenberg .

Islands of a surface superconducting phase were also observed for the surface-doped tungsten bronzes, WO3:Na, at $`T_c=91.5K`$Reich  . The ARPES results mentioned above suggest that SSC may actually be a rather widespread phenomenon.

It is noteworthy that the bulk WO3 is an insulator at low dopingGoodenough . Recently a remarkable breakthrough has been accomplished by doping insulator surfaces electrically in the so-called Field-Effect Transistor (FET) geometrySchon . SSC with high Tc was induced by both electronic and hole doping of films of the prototype cuprate material, CaCuO2Schon1 .

All this makes us believe that the search for SSC emerges as a new and important development in studies of the properties of surfaces, especially their metallic properties. As for the SSC itself, its mechanisms are unknown, and may have nothing to do with the ones in the bulk. Of a special challenge is the possibility of superconductivity at surfaces of ordinary metals such as Cu Petersen  or lithium-doped Mo and WRotenberg , and the influence of substrate on superconductivity in thin films. Thus, it is necessary to return to a more careful investigation of the low temperature properties of surfaces with absorbed atoms and the role of adsorbed atoms as dopants of carriers into the metallic surface zone (e.g., see Ref.Rotenberg ) like it was discovered for WO3:Na Reich .

The major challenge to such a program lies in discerning SSC and studying its properties. Even for insulators where doping by FET, ideally, provides an effective control on surface properties, experimental tools to probe SSC are limited in numbers. The Meissner effect due to surface superconducting islands would probably never produce a bulk screening. Thermodynamical probes are also difficult because of the smallness of contributions from surface layers of atomic thicknesses. So far, e.g., in Ref. Schon , surface superconductivity has only been detected by measuring resistivity dependence on temperature in the perpendicular-to-plane magnetic fields.

In what follows we focus on destroying SSC by magnetic fields applied parallel to the surface. In high enough fields one should expect the appearance of a 2D version of inhomogeneous superconducting state known as Larkin-Ovchinnikov-Fulde-Ferrell phase (LOFF) LOFF . It was also shown that two-dimensionality broadens the region of the LOFF state on the B-T phase diagramBulaevskii . Motivated by experimental findings mentioned above, we investigate the peculiarities introduced into this phenomenon by the SO effect or non-s-wave pairing. Experiments by Schön et al. not only have proved FET to be an effective doping process. Different levels of doping result in different $`T_c`$-s. Our results for SSC in parallel magnetic fields are expressed in terms of this $`T_c`$, providing the tool for comparing theoretical predictions with experiments by controlling the doping level. Theoretically, there is no long range order in a 2D superconductor. However, correlations are destroyed on the exponentially large spatial scale, $`R\xi _0exp(E_F/T)`$, which would exceed the size of the film.

We employ below the weak-coupling BCS-like scheme by assuming that electrons interact via a weak short-ranged interaction, $`U(𝒓,𝒓^{})`$. Then $`T_cϵ_F`$, and only a narrow vicinity of the Fermi surface is involved. Thus, the interaction in the momentum representation can be taken in the form:

$$U(𝒑,𝒑^{})=\underset{l}{}U_l\chi _l(𝒑)\chi _l(𝒑^{\mathbf{}}),$$
(1)
where $`𝒑`$, $`𝒑^{\mathbf{}}`$ lie on the Fermi surface; the angular dependence is expressed through a complete set of basis functions $`\chi _l(𝒑)`$ (index $`l`$ enumerates different representations, as in expansions over the spherical functions in a 3D isotropic model). Superconducting order parameter, the ”gap”, $`\widehat{\mathrm{\Delta }}_{\alpha \beta }(𝒑)`$ is defined by the equation:

$$\widehat{\mathrm{\Delta }}_{\alpha \beta }(𝒑)=|U_l|\chi _l(𝒑)\frac{d^3p^{}}{(2\pi )^3}\chi _l(𝒑^{\mathbf{}})\left(T\underset{\omega _n}{}F_{\alpha \beta }(𝒑^{\mathbf{}};i\omega _n)\right),$$
(2)
where $`F_{\alpha \beta }(𝒑^{\mathbf{}};i\omega _n)`$ stands for the Fourier component of Gor’kov anomalous function:

$$\widehat{F}_{\alpha \beta }(𝒓𝒓^{},\tau \tau ^{})=T_\tau \left(\widehat{\mathrm{\Psi }}_\alpha (𝒓,\tau )\widehat{\mathrm{\Psi }}_\beta (𝒓^{},\tau ^{})\right),$$
(3)
and $`U_l<0`$ is a constant in Eq.(1) corresponding to the selected pairing channel. When the field operators for electrons are re-written in momentum space, $`\widehat{\mathrm{\Psi }}_\alpha (𝒓,\tau )=_𝒑\widehat{\mathrm{\Psi }}_\alpha (𝒑,\tau )e^{i𝒑𝒓},`$ Eq.(2) becomes:

$$\widehat{\mathrm{\Delta }}_{\alpha \beta }(𝒑)=|U_l|\chi _l(𝒑)\frac{d^3p^{}}{(2\pi )^3}\chi _l(𝒑^{\mathbf{}})\widehat{\mathrm{\Psi }}_\alpha (𝒑^{},\tau )\widehat{\mathrm{\Psi }}_\beta (𝒑^{},\tau ).$$
(4)
The two operators inside brackets in Eq.(4) anti-commute. In the presence of the center of inversion (CI) the behavior of $`\chi _l(𝒑^{\mathbf{}})`$ at $`𝒑^{\mathbf{}}𝒑^{\mathbf{}}`$ alone determines the symmetry of $`\widehat{\mathrm{\Delta }}_{\alpha \beta }(𝒑)`$ (even (singlet) vs odd (triplet) parity pairing). For an s- or d-pairing the order parameter below $`T_c`$ has the form:

$$\widehat{\mathrm{\Delta }}_{\alpha \beta }(𝒑;𝒒)=\mathrm{\Delta }(𝒒,T)(i\sigma _y)_{\alpha \beta }\chi _l(𝒑),$$
(5)
where the momentum $`𝒒`$ stands for the spatial dependence of the gap amplitude, $`\mathrm{\Delta }(𝒒,T)`$.

Surface always breaks the CI symmetry due to the difference between the ”top” and the ”bottom”. The direction bulk-to-surface determines $`𝒏`$, a unit vector normal to the surface. Qualitative changes in the surface electronic spectrum come about from the well-known Rashba termR60 :

$$\widehat{h}_{SO}=\alpha (𝝈\times 𝒑𝒏),$$
(6)
which specifies SO interactions at a surface. Eq.(6) lifts the two-fold spin degeneracy for band electrons. The electron spectrum now consists of two branches with two Fermi surfaces:

$$ϵ_\pm (𝒑)=v_F(pp_{F\pm });p_{F\pm }=p_F(1\pm \frac{\alpha }{v_F})$$
(7)
Even though SO splitting, $`2\alpha p_F`$, may be on a scale of tenths of eV LaShell ; Rotenberg , we assume that $`2\alpha p_Fϵ_F`$. Below we speak of a strong or weak SO meaning the relative values of $`2\alpha p_F`$ and Tc.

Some theory issues regarding superconductivity without CI due to the presence of the SO term Eq.(6) were first considered in Ref.Edelstein  and more recently in Ref. GR . The Gor’kov function (in the momentum representation), $`F_{\alpha \beta }(𝒑,+0)=\widehat{\mathrm{\Psi }}_\alpha (𝒑)\widehat{\mathrm{\Psi }}_\beta (𝒑)`$, that stands under integral in Eq.(4) represents the wave function for Cooper pairs in the condensate. In the presence of CI symmetry the latter can be classified according to the parity:

$$F_{\alpha \beta }(𝒑;+0)=\{\begin{array}{c}i(\sigma _y)_{\alpha \beta }f(𝒑);(f(𝒑)even)\\ i[(𝒅(𝒑)\widehat{𝝈})\widehat{\sigma }_y]_{\alpha \beta };(𝒅(𝒑)odd)\end{array}$$
(8)
With CI broken by non-zero SO term Eq.(6), the pairing wave function becomes a mixture of even and odd terms. It is important to realize that while this mixing changes physical properties of the SC phase, the gap order parameter, $`\widehat{\mathrm{\Delta }}_{\alpha \beta }(𝒑;𝒒)`$, preserves its singlet form Eq.(5). For instance, s-pairing indeed induces the non-zero triplet component Eq.(8), as shown in Refs.Edelstein ; GR . However, the latter does not automatically generate a ”triplet” gap, $`\widehat{\mathrm{\Delta }}_{\alpha \beta }^t(𝒑;𝒒)`$. Indeed, re-writing the integration over $`𝒑^{}`$ in Eq.(4) as $`d^3p^{}dS_F^{}d\xi `$, where $`\xi =v_F(pp_F)`$, we notice that the triplet $`F`$-component is odd in particle-hole transformation, $`\xi \xi `$, and, hence, the integrals of the form Eq.(4) would only give small terms of order $`(\alpha p_F/ϵ_F)1`$. In other words, while SO interaction may significantly change spin structure of the normal and anomalous Green functions, the ”gap”, $`\widehat{\mathrm{\Delta }}_{\alpha \beta }(𝒑,𝒒)`$, Eq.(2) preserves its usual form Eq.(5) with $`\chi _l(𝒑)=const`$ for isotropic pairing, and $`\chi _l(𝒑)(p_x^2p_y^2)`$ for the d-wave pairing.

We have calculated $`T_c`$ numerically for the 2D superconductor as a function of magnetic field and spin-orbit interaction. The result is shown in Fig.1 for the magnetic field strictly parallel to the surface (to exclude diamagnetic currents)Bulaevskii .

A 1st order phase transition was initially expected between superconducting and normal states, defined by comparing their free energies: $`F_s(T)=F_n\chi _N\frac{B^2}{2}`$, which would determine the so-called paramagnetic critical field, $`H_{par}`$ Clog  ($`\chi _N`$ \- the spin susceptibility in the normal phase). The transition from normal to superconducting state is actually (at lower temperatures) a second order transition into the LOFF state. The details of the phase diagram in the vicinity of the $`H_{par}`$ were studied numerically in Ref.br . The LOFF phase boundary is determined by Eq.(2) or Eq.(4), linearized in $`\mathrm{\Delta }_{\alpha \beta }(𝒒)`$ at an extremal $`𝒒`$. The corresponding expression for anomalous function linear in $`\mathrm{\Delta }_{\alpha \beta }(𝒒)`$ is obtained by solving the proper Gor’kov equations:

$$F_{\alpha \beta }(𝒑,𝒒;i\omega _n)=\widehat{G}_{\alpha \nu }^{(0)}(𝒑;i\omega _n)\widehat{\mathrm{\Delta }}_{\rho \nu }(𝒑,𝒒)\widehat{G}_{\beta \rho }^{(0)}(𝒑+𝒒;i\omega _n),$$
(9)
where $`\widehat{G}_{\alpha \beta }^{(0)}(𝒑;i\omega _n)`$ is the normal state Green function at non-zero $`\widehat{h}_{SO}(𝒑)`$ of Eq.(6), together with the Zeeman term, $`\mu _B(𝝈𝑩)`$:

$$[i\omega _n\xi \widehat{h}_{SO}(𝒑)\mu _B\widehat{𝝈}𝑩]\widehat{G}_{\alpha \beta }^{(0)}(𝒑;i\omega _n)=\widehat{1}.$$
(10)
The spin Hamiltionian on the left side, $`\widehat{H}=\widehat{h}_{SO}(𝒑)+\mu _B\widehat{𝝈}𝑩`$, may be easily diagonalized:

$$\stackrel{~}{ϵ}_\lambda (p)\xi =\lambda \sqrt{\alpha ^2p_F^2+2\alpha p_{Fy}\mu _BB+(\mu _BB)^2}\lambda \stackrel{~}{ϵ}(p)$$
(11)
($`\lambda =\pm 1`$ for the two branches, $`p_{Fy}`$ is the y-axis projection of the Fermi momentum, and $`Bx`$) with the eigenfunctions, spinors $`\eta ^\lambda (𝒑)`$ of the form:

$$\eta ^\lambda (𝒑)=\frac{1}{\sqrt{2}}\left\{\begin{array}{c}1\\ \frac{\mu _BBie^{i\phi (𝒑)}p_F\alpha }{\lambda \stackrel{~}{ϵ}(p)}\end{array}\right\}.$$
(12)
Substitution of Eq.(9) into Eq.(2) making use of Eq.(5) results in a rather cumbersome expression which generalizes the corresponding Eq.(7) of Ref.Bulaevskii . We sketch, therefore, only a few results for the low-T part of the phase diagram in Fig.1. Below we discuss the main changes in the shape of the phase diagram, as introduced by SO coupling or anisotropy. As for the structure of the LOFF phase itself, we assume that the numerical analysis done in Ref.br  remains applicable, i.e. the order parameter in the LOFF state has the structure of periodic stripes.

We state in more detail our results for the limiting cases of strong and weak SO interaction, which significantly simplify all calculations.

a) Strong SO: $`\alpha p_F\mathrm{\Delta }(0)`$; subsequent analysis, which we do not provide here, leads after short calculations to our final results:

$`𝒒𝑩,`$ $`|𝒒|={\displaystyle \frac{2\mu _BH_{c2}}{v_F}};`$
$`\mu _BH_{c2}`$ $`=`$ $`\sqrt{2\mathrm{\Delta }(0)\alpha p_F}\sqrt{{\displaystyle \frac{2\pi }{\gamma }}T_{c0}\alpha p_F}`$ (13)
One sees that SO interaction not only enhances the value of $`H_{c2}`$ in comparison with the LOFF for 2D model of Ref.Bulaevskii  , but it also fixes the direction of the structure vector. The resulting stripe structure is parallel to the magnetic field direction and has the perpendicular space periodicity $`L=\sqrt{\frac{\pi \gamma }{2T_{c0}\alpha p_F}}v_F`$.

According to GR , in the limit of strong SO interaction the spin susceptibility for parallel fields in superconducting state, $`\chi _S(T)`$ is non-zero and equal to $`\frac{1}{2}\chi _N`$. This increases the critical paramagnetic field Clog  only by a factor of $`\sqrt{2}`$, $`\mu _BH_{par}=\mathrm{\Delta }(0)`$. Comparing this with Eq.(13), one sees that strong SO significantly increases the area occupied by the LOFF state by a factor of $`\sqrt{\alpha p_F/T_{c0}}1`$. Strong SO scattering by defects also enhanses $`H_{c2}`$Klemm , but the LOFF state does not exist in presence of disorder. Analytical expressions can also be obtained at low temperatures. For the dependence of the transition temperature on magnetic field, $`T_c(B)`$, one obtains, at $`TT_{c0}\sqrt{T_{c0}/(\alpha p_F)}`$:

$$T_c(B)=5.784\left(\sqrt{2\pi T_{c0}\alpha p_F/\gamma }\mu _BB\right)^3/(\alpha p_F)^2,$$
(14)
and

$$T_c(B)=\pi T_{c0}^2/(2\gamma \mu _BB),T_{c0}\sqrt{T_{c0}/(\alpha p_F)}TT_{c0}$$
(15)
The expression for $`T_c(B)`$ also simplifies for small magnetic fields near $`T_{c0}`$:

$$\frac{T_{c0}T_c}{T_{c0}}=\frac{7\zeta (3)}{8\pi ^2}\frac{(\mu _BB)^2}{T_{c0}^2}$$
(16)
Note that while the spin-orbit interaction splits the Fermi surface, quasiparticles with the same band spin index Eq.(11) form Cooper pairs in the superconducting state, so that spin-orbit interaction alone does not change $`T_c`$.

b) Weak SO: $`\alpha p_F\mathrm{\Delta }(0)`$; unlike in case of strong SO, the Cooper pair is formed mainly by pairing of electrons from the FS’s with different spin indices. The LOFF phase in a 2D superconductor with no SO interaction was first analysed in Refs Bulaevskii ; br . According to Ref.Bulaevskii , $`\mu _BH_{c2}=\mathrm{\Delta }(0)=\sqrt{2}H_{par}`$, $`v_Fq=2\mu _BH_{c2}`$.

This result is zero order in spin-orbit interaction. The direction of $`𝒒`$ is not fixed with respect to $`B`$. Analysis to the second order in $`\alpha p_F`$ results in an anisotropy term which again fixes the vector $`𝒒`$, as in strong SO case Eq.(13), perpendicular to the direction of the magnetic field. Indeed, for the critical field as a function of the angle $`\beta `$ between $`𝒒`$ and $`𝑩`$ we find:

$$\mu _BH_{c2}=\mathrm{\Delta }(0)\frac{\alpha ^2p_F^2}{2\mathrm{\Delta }(0)}cos^2\beta ,$$
(17)
i.e. the maximum for $`H_{c2}`$ is reached for $`\beta =\pm \frac{\pi }{2}`$.

For small magnetic fields near $`T_c`$ we get:

$$\frac{T_{c0}T_c}{T_{c0}}=\frac{7\zeta (3)}{4\pi ^2}\frac{(\mu _BB)^2}{T_{c0}^2}$$
(18)
A quadratic dependence on $`B`$ remains valid for fields $`\mu _BB\alpha p_F`$. However, note the factor two difference between Eqs.(16) and (18): suppression of $`T_c`$ by magnetic field turns out to be slower in case of strong spin-orbit interaction than in case when the spin-orbit interaction is weak.

Although a more complicated LOFF periodic superstructure is possible, the energy considerations of Ref.LOFF  have shown in 3D that the stripe phase is energetically more favorable. A detailed study done numerically in 2D br  has shown a more complicated than just a sinusoidal shape of the order parameter. We assume that the results of Ref.br  remain valid in all cases considered above, so that the LOFF state preserves its striped order parameter form.

The major role of SO is in fixing of the LOFF superstructure. Anisotropy fixes the orientation of the LOFF stripes as well (see below). A 1st order re-orientation transition may be expected corresponding to an abrupt change in the direction of the superconducting stripes (similar to a spin-flop transition), as the magnetic field is rotated in the 2D plane.

In the above discussion so far we have neglected any anisotropy at all. Meanwhile, the anisotropy is, of course, important. We address the issue of the effect of pinning stripe direction to the particular crystal axis only for the d-wave order parameter, since the latter is intrinsically anisotropic. Stripes may orient themselves along the directions of the gap maximums. It can be easily shown that $`H_{c2}`$ for the $`d_{x^2y^2}`$ superconducting order parameter takes the form: $`\mu _BH_{c2}=\frac{\pi }{\gamma }e^{1/4}T_{c0}`$. The critical field for the LOFF phase of d-wave is somewhat higher than for s-wave (without the SO term Eq.(7)), as supeconducting stripes get pinned to the crystal axes by the form of the order parameter and the direction of the magnetic field. Reorientation transition/twinning is also expected for any other cause.

In summary, we have shown: i) that inhomogeneous state in parallel fields extends considerably the low temperature phase diagram of surface superconductivity with increased spin-orbit interaction; ii) all SC characteristics of the phase diagram in the $`(B,T)`$ plane can be expressed in terms of $`T_{c0}`$, the critical temperature in the absence of the field, which is controlled by doping, and, hence, not only allows a data comparison in the broad range of FET-doping, but may also serve as a method to extract the value of the SO interaction; iii) LOFF state properties are strongly anisotropic in the plane with respect to field direction - this, for example, can be seen by measuring the anisotropy of the AC susceptibility signal; iv) the indispensible feature of the LOFF state must be the re-orientation transitions at the field rotation in the plane caused by locking of the LOFF order parameter by anisotropy.

We would like to thank L. N. Bulaevskii for useful discussions and V. M. Edelstein for a few references, including Ref.Edelstein . This work was supported by NHMFL through the NSF Cooperative agreement No. DMR-9521035 and the State of Florida.


format: text

text:
# Complex X-ray spectral variability in Mkn 421 observed with XMM-Newton

## 1  Introduction

Mkn 421 has a long history of observations at X–ray and other wavelengths. It is nearby ($`\mathrm{z}=0.031`$) and is generally the brightest BL Lacertae object in the X–ray sky. BL Lac objects are a sub–class of blazars; sources whose broad–band continuum emission are dominated by relatively featureless nonthermal radiation and which show rapid spectral variability with the characteristic timescale of variations typically decreasing with increasing energy.

The spectral energy distribution (SEDs) of blazars have two distinct peaks in $`\nu F_\nu `$ space; one in the IR/optical or UV/X–ray regime and one at higher energies in the X–ray/Gamma–ray regime. The lower energy component is probably incoherent synchrotron radiation and the high energy component Compton up–scattering of lower energy photons. It seems likely that the synchrotron photons provide the seeds for the Compton emission (the synchrotron self–Compton process), but this has not been unambiguously demonstrated. BL Lacs are often classified as LBLs (low–energy peaked BL Lacs) if the synchrotron component peaks in the IR/optical regime or HBLs if it peaks in the UV/X–ray, although the classes may be the extrema of a single population (e.g. Fossati et al. 1998). Mkn 421 is an HBL.

The rapid flux variability in blazars would require the luminosity to be highly super–Eddington if the emission was isotropic and it is generally assumed that the emission region moves at relativistic speeds towards the observer; an interpretation supported by observations of superliminal motion of radio knots in the jet–like structure of some radio galaxies.

Although the relativistic jet model provides a plausible explanation for the continuum properties of blazars, detailed physical constraints on the structure of the jets remain an elusive goal; snap–shot multi–wavelength spectra can be reproduced by a variety of jet models with different underlying assumptions. Studying the variability behaviour of both the spectrum and flux is therefore the most profitable way to constrain the various jet models. The relationship of the variations within different energy bands can, in principle, help to constrain the geometry of the emission region on the macroscopic scale and constrain the physical radiation and loss processes on the microscopic scale. For instance, X–ray variations in the HBL, PKS 2155–304, have been seen to lead those in the optical/UV (e.g., Edelson et al. 1995; Urry et al. 1997), supporting the view that the X-rays arise from the synchrotron component and not from Compton scattering from these lower–energy bands.

The X–ray spectral slope is especially sensitive to electron cooling timescales within the various models and observations of interband lags within this regime have often been interpreted in this context. Chiappetti et al. (1999) found that soft X–ray (0.1–1.5 keV) variations in PKS 2155–304 led hard X–ray (3.5–10 keV) variations by $``$0.3-4 hr and suggested it could occur through radiative cooling if the electron population has a sharp break or low–energy cut–off. Historically, however, the source has not exhibited a consistent behaviour in this respect. X–ray flares have also been observed with the hard band leading the soft (Sembay et al. 1993).

Mkn 421 has also been the subject of a number of monitoring campaigns and shows a similar variety of behaviour between observations. Takahashi et al. (1996) found that the hard X–ray band led the soft, whereas Fossati et al. (2000a) found the soft band led the hard. Takahashi et al. (2000) reported both types of behaviour seen within a single very long observation ($`11`$ days) of the source by BeppoSAX and ASCA.

Interpretation of interband lags can, however, be complicated by the fact that these observations have been made with low–Earth orbit satellites, resulting in periodic interruptions in the light curves every $`1.6`$ hr due to Earth–occultation. Edelson et al. (2001; hereafter E01) analysed uninterrupted XMM-Newton data on PKS 2155–304 and found no significant interband lag down to a reliable limit of $`0.3`$ hr, attributing the previous claims of lags on short time scales to a spurious signal induced into the CCF at the orbital period of the satellites. Takahashi et al. (2000) also noted the possibility of such an effect in their paper on Mkn 421. One advantage XMM-Newton has over previous X–ray missions for monitoring purposes is that it has a highly eccentric orbit with a 48 hour period, thus allowing much longer continuous observation periods. Nominally the only data gaps occur during passage through the radiation belts at perigee which takes around 7–8 hrs and, in the first year or so of the mission, a $`1`$ hour gap at apogee due to incomplete ground station coverage (Jansen et al. 2001).

XMM-Newton has observed both PKS 2155–304 and Mkn 421 on several occasions since launch. We have reported on the XMM-Newton observation of PKS 2155–304 in E01. In this paper we now present the results of several observations of Mkn 421. The observations and data reduction are discussed in the next section. The temporal analysis of these data are reported and compared with other data in § 3. The spectral analysis is discussed in § 4 and a concluding discussion is given in § 5.

## 2  Observations and Data Analysis

### 2.1  Observations

The XMM-Newton observations of Mkn 421 were performed as part of the performance–verification and calibration programme. The source was observed on–axis on four occasions at the epochs listed in Table 1. XMM-Newton has three co–aligned X–ray telescopes with focal plane cameras of two types; two EPIC–MOS imaging CCD cameras and one EPIC–pn imaging CCD camera (Turner et al. 2001; Strüder et al. 2001). Within the telescopes containing the MOS cameras, gratings divert approximately half the incident flux to two Reflection Grating Spectrometer (RGS) instruments (den Herder et al. 2001).

In this paper our analysis concentrates on data taken with the high–throughput EPIC CCD cameras. The EPIC cameras have superior sensitivity to short timescale flux variations and a broader energy response than the RGS instruments. During orbital revolutions 165, 171 and 259 (henceforth, O165, O171 and O259) both EPIC cameras were in their respective Small Window (SW) imaging modes throughout the observation. In orbit 84 (O84), however, the first part of the observation was performed in timing mode for the pn and MOS1 cameras, and SW mode for MOS2, followed by SW mode for the pn and MOS1 and timing mode for MOS2. There is, however, only a small overlap of about 5 ksec between the pn and MOS imaging data. Due to a radiation alert the MOS cameras were switched off shortly after the beginning of the pn SW mode observation (further details of the O84 observation can be found in Brinkmann et al. (2001) which gives the first analysis of the data from this observation). A variety of blocking filters were employed throughout the observations and these are listed in Table 1.

The data analysis software and calibration were generally in a more mature state for imaging mode than timing mode when the analysis was performed, so this paper deals with the imaging data only.

The imaging observations were performed in each instrument’s SW mode because they have a faster CCD frame readout time than their respective Full–Frame modes. This mode reduces the effects of pulse pile–up which causes spectral distortion in bright sources. Pulse pile–up occurs when two (or more) X–rays interact within a given frame close enough such that their deposited charge distributions overlap. This can lead to a loss of flux if the resultant charge distribution does not have a valid pattern within the instruments X–ray pattern library, or, if it still has a valid pattern, a new event whose energy is the sum of the contributing events. Hence, pile–up leads to an artifical hardening of the spectrum and can produce spurious correlations between soft and hard–band variations. The simplest way to correct for pile–up is to exclude the core of the point–spread function where pile–up is greatest. This is explained in more detail in the following Section.

The pn SW mode has a $`63\times 64`$ pixel window (equivalent to $`4.3\times 4.4`$ arcminutes on the sky) with a frame readout time of 5.7 ms. The common telescope boresight places the source off–centre within the window, typically less than an arcminute from the CCD boundary. For the MOS cameras the window is $`100\times 100`$ pixels ($`1.8\times 1.8`$ arcminutes), more centered with respect to the boresight, and the mode has a frame readout time of 0.3 s.

### 2.2  Data Reduction

A standard reduction of the raw pn and MOS event lists was performed using the Science Analysis System (SAS) Version 5.2. The SAS removes events from the location of defective pixels, corrects for charge transfer losses across the CCDs and applies a gain calibration to convert the recorded charge in ADC units to a pulse–invariant energy scale. Standard analysis packages are then used to create images, light–curves and spectra from the calibrated event lists.

Pulse pile–up in CCD detectors is nominally a purely statistical effect which can in principle be corrected for by modifying the instrument response function appropriately. However, without having a suitable tool available within the SAS, we have used the more direct approach of excluding events from the core of the point–spread function where, naturally, pile–up is the most extreme. Although this leads to a loss of flux, the high throughput of XMM means that we are still left with spectra and light curves with high statistical significance.

The size of the core to be excluded was determined by examining the distribution of mono–pixel to multi–pixel events as a function of energy. As the radius of the excluded region is increased, the distribution tends towards that consistent with the theoretical distribution for the non piled–up case. By this method we excluded core radii of diameter 8”, 6”, 10” and 10” respectively for the pn O84, O165, O171 and O259 observations and 8”, 12” and 12” for the MOS O165, O171 and O259 observations. The maximum radius of each annulus was 36.75” during O84 increasing to 44” for the remaining observations; the difference being due to a change in the relative boresight with respect to the pn window.

The point-spread function of each telescope is energy dependant so a correction needs to be applied to the effective area of the response file when fitting spectra extracted from an annulus with a gven inner and outer radius. To do this we used a description of the psf developed by the EPIC instrument teams using in-flight data from several on-axis point sources (Ghizzardi & Molendi 2002). The robustness of this procedure and other calibration related issues are discussed in more detail in Section 4.

In table 1 we list the mean observed $`0.110.0`$ keV source count rates (obtained after excluding the psf cores) during each observation. The source was most luminous during O171, but only $`1/3`$ as bright twelve days earlier during O165. O84 and O259 are at intermediate levels. The higher observed count rate in the soft band during O259 is primarily a consequence of the higher throughput of the Thin blocking filter compared with the Medium and Thick filters used earlier.

## 3  Timing Analysis

In this section, we compare gross variability in different energy bands for each of the observing epochs. Much of this analysis parallels that on PKS 2155–304 by E01. Although most of the analyses will be restricted to higher signal-to-noise pn data, the relevant MOS data are also included for consistency.

The light curves were constructed as follows: first, data were divided into soft (S; 0.1–0.75 keV), medium (M; 0.9–1.7 keV) and hard (H; 2–10 keV) bands. Then, the data were initially extracted at 1 sec intervals. They were then rebinned to 300 sec resolution and standard methods were used to measure the mean and standard errors from the spread of data in each bin. Due to data loss, there were two small gaps in the O165 light curves. The four points during 1.75–2.00 hrs and three points during 5.08–5.25 hrs after the beginning of the integration were instead determined by interpolation from the surrounding points.

In all bands, the source count rate exceeded the background by at least a factor of $`50100`$, so no background light curves were extracted. The resulting light curves are shown in Figure 1.

### 3.1 Variability Amplitudes

Variability amplitudes were measured for each light curve using the fractional variability amplitude, defined as

$$f_v=\frac{\sqrt{S^2\sigma _{err}^2}}{X},$$
(1)
where $`S^2`$ is the total measured variance, $`\sigma _{err}^2`$ is the mean squared errors, and $`X`$ is the mean count rate (see E01 for details). The fractional variability amplitude is also simply the square-root of the “excess variance” parameter used e.g., in Nandra et al. (1999); $`f_v=\sqrt{\sigma _{}^{2}{}_{xs}{}^{}}`$. These quantities are computed and shown in Table 2. Examination of this and Figure 1 indicate that stronger variations are clearly seen in the harder bands. This increasing variability with energy band is common in blazars (e.g., Ulrich, Maraschi & Urry 1997) and is a consequence of the fact that high energy electrons have shorter cooling timescales than low energy ones. Finally, different variability levels are seen in each of the epochs with O84 exhibiting the greatest overall degree of flux variations.

### 3.2 Comparison of the Hard and Soft Light Curves

Figure 2 shows a variety of quantities derived directly from the light curves in Figure 1. The top panel gives the geometric mean flux of the three bands, $`F_g=\sqrt[3]{F_H\times F_M\times F_S}`$. The second panel is the hardness ratio, $`HR=F_H/F_S`$. Note that the mean and hardness ratio light curves generally track each other very well, confirming that the source generally gets harder as it gets brighter, with the exception of the end of O165 where we observe a softening of the source as it continues to brighten.

Below that is an “overplot” diagram for the hard and soft bands. In these the light curves are subjected to the first two steps of cross-correlation function (CCF) analysis: at each point in the light curve, the mean count rate is subtracted and then the data are divided by the standard deviation. The resulting plot, with zero mean and unit variance for each light curve, provides a both a simple way to compare variations between bands as well as a convenient “sanity check” of the final CCFs by showing the intermediate step. Finally, the bottom panel is the difference between the two scaled light curves.

The temporal relation between hard and soft band variations shows a range of behaviors: O84: almost perfect agreement to within the errors; O165 and 259: good general agreement but significant divergances on short ($`>`$hour) time scales; O171: poor overall agreement and even a divergence between soft and hard bands. As seen below, the light curves that track the best (O84 and O259) are particularly useful because they allow the finest measurement/limits on the interband lag.

### 3.3 Cross-Correlation Analysis

Cross-correlation functions were measured using the discrete correlation function (DCF; Edelson & Krolik 1989) and the interpolated correlation function (ICF; White & Peterson 1994). These are shown in Figure 3, where a positive $`\tau `$ indicates that the softer band leads the harder band. They are also tabulated in Table 3. The errors, measured using the bootstrap method of Peterson et al. (1998), indicate a 2$`\sigma `$ upper limit of $`\tau 0.08`$ hr.

The best-tracking epoch (O84) also showed the highest correlation coefficients ($`r0.99`$) and the strongest limits on the interband lags ($`\overline{\tau }=0.01\pm 0.03`$). Such similar limits are seen in the next two strongest correlations (O165 S/M and O259 M/H, with $`r=0.96`$) and for the eight highest correlations (out of 12, those with $`r0.89,\overline{\tau }=0.01\pm 0.04`$). Recall that the bin size is 0.08 hr (300 sec); these indicate very secure ($`23\sigma `$) limits of $`|\tau |<0.08`$ hr).

However, it is also clear especially in O171 that the hard and soft light curves are not always tracking each other with high precision. In this case, the hard and soft actually diverge, and the peak correlation function is less than zero. This reveals a further complexity in the spectral variability that cannot be addressed with these short snippets of data, and indicates the need for much longer uninterrupted observations to fully describe the observational picture.

### 3.4 Comparison with Previous Results

Previous observations of Mkn 421 have claimed significant hard lags, in which the hard band lagged behind the soft by $``$0.6 hr (Fossati et al. 2000a), as well as soft lags, in which the soft band lagged behind the hard by $``$1 hr. By comparison, this paper finds that only during one of the four epochs are the data sufficiently well correlated to measure a meaningful CCF, and in that case, there is no measurable interband lag down to much tighter 2$`\sigma `$ limits: $`|\tau |0.08`$ hr (1 bin width).

In order to understand this discrepancy, it is important to keep in mind that the previous claims of interband lags all involved observations made with low-Earth orbit satellites (BeppoSAX and ASCA). As pointed out by E01, CCFs measured from such data are corrupted on time scales less than the orbital period (1.6 hr) and thus can yield spurious apparent lags. This can be seen in Figure 11 of Fossati et al. (2000a): there is a periodic feature running through the CCF with a period of 1.6 hr, causing excursions larger than any claimed local peaks. Indeed, the highest point in that CCF is at zero lag, and the next two highest points are immediately beside it. On the other hand, XMM-Newton produces uninterrupted light curves that do not suffer from this systematic errors, so the CCFs are intrinsically more reliable. Thus, we conclude that there is no reliable evidence for interband lags, hard or soft, in Mkn 421.

## 4  Spectral Analysis

The instrinsic X–ray spectra of BL Lacs are basically smooth and featureless. Over an energy band of roughly a decade the continuum can be well represented by a power law model whereas over a wider band the continuum exhibits a degree of curvature requiring fitting by broken–power law or continuously curved (CC) spectral models (e.g. Inoue & Takahara 1996; Fossati et al. 2000b). In this paper we have confined our analysis to the high energy 2.0 to 10 keV spectral band. There are unresolved cross–calibration issues between the EPIC instruments in the soft band which lead to somewhat contradictory results regarding the slope of the low–energy continuum. A full analysis of the entire spectral range will be left to a later paper. In the chosen band, the derived spectra can be fit a single power law with the column density fixed at the Galactic line–of–sight value of $`1.6\times 10^{20}\mathrm{cm}^2`$ (Lockman & Savage 1995) which is in agreement with previous spectral analyses of the broad–band X–ray spectrum (e.g. Fossati et al. 2000b).

Our spectral analysis on Mrk 421 proceeded as follows. Having removed the piled–up section of the point spread function as described earlier, we subdvided each observation into 5 ksec segments. This allowed us to extract spectra with good statistical significance and still track spectral variations throughout each observation. Between six and nine spectra were extracted per observation and per camera. For completeness, background spectra were extracted, but the background rate was found to be negligible with respect to the source and hence the background spectra were ignored. The pn data were analysed separately to the MOS, but the data from both MOS cameras were combined in the spectral fitting. All spectral fitting was done with Xspec Version 11.00.

Table 4 lists the derived spectral parameters from each fit. There is an excellent agreement between the pn and MOS parameters in general. The pn photon indeces for example are typically within $`\mathrm{\Delta }\alpha 0.05`$ of the MOS. The overal range in 2-10 keV source flux is similar to that observed in the BeppoSAX campaigns of 1997 and 1998 (Fossati et al. 2000b).

In Figure 4 we have plotted the spectral index from all observations, against the 2 to 10 keV (deabsorbed) flux. For clarity, the epic-pn parameters only are plotted. From this plot we can see that although the variation of slope with flux follows in general the classical trend of synchrotron emission, namely, that steeper slopes correspond with lower flux states, there is no strong correlation between these parameters globally. As expected from the previous discussion of the light curves, the variability behaviour is quite different between observations. In O84 the slope is stable but the flux varies. In O165, when the source was weakest and the spectrum steepest, neither parameter varied strongly. O171 and O259 show a similar behaviour in that a fairly linear correlation between slope and flux is evident.

## 5  Discussion

The temporal analysis found that in O84, the data were highly correlated with no evidence for an interband lag down to time scales of $`|\tau |0.08`$ hr, much tighter limits than previous claims of both hard and soft lags. Following on from similar results on PKS 2155–304 by E01 this provides further strong evidence that the detection of lags on short time scales by low–Earth orbit satellites are probably spurious artefacts caused by Earth occultation.

Measurements of a lag could potentially give the cooling time scale and therefore a limit on the magnetic field strength assuming that the lag is due to the electron population softening due to radiative losses. From Chiappetti et al. (1999),

$$B\delta ^{1/3}=300\left(\frac{1+z}{\nu _1}\right)^{1/3}\left(\frac{1(\nu _1/\nu _0)^{1/2}}{\tau }\right)^{2/3}G,$$
(2)
where $`B`$ is the magnetic field strength, $`\delta `$ is the Doppler factor, $`z`$ is the redshift, and $`\nu _1`$ and $`\nu _0`$ are the frequencies (in units of $`10^{17}`$ Hz) at which a lag, $`\tau `$, is measured. For the Orbit 84 measurement, we have an upper limit on the lag so we can derive a lower limit on the value of $`B\delta ^{1/3}`$. The most stringent limit is given by the hard and soft light curves with $`\nu _014.5`$ and $`\nu _11.0`$, $`z=0.031`$ and $`\tau \stackrel{<}{}300`$s giving $`B\delta ^{1/3}\stackrel{>}{}4.7`$ G.

This lower limit is rather higher (typically by a factor of $`510`$) than obtained from some estimates of $`B`$ and $`\delta `$ in the literature derived from fitting a variety of homogeneous SSC models to measurements of the broad–band SED of Mkn 421 observed at different epochs. For example, Ghisellini et al. (1998) obtain $`B=0.09`$ G and $`\delta =12`$ giving $`B\delta ^{1/3}=0.2`$ G for a pure SSC model. If the seed photons of the Compton emission arise from an external source (e.g. an accretion disk or the broad–line region) then the model allows for a somewhat higher magnetic field strength. In this case Ghisellini et al. (1998) obtain $`B=0.22`$ G and $`\delta =11`$ giving $`B\delta ^{1/3}=0.49`$ G. In comparison, Tavecchio, Maraschi and Ghisellini (1998), obtained $`B=0.25`$ G and $`\delta =25`$ giving $`B\delta ^{1/3}=0.7`$ G, also with a homogeneous SSC model, and Krawczynski et al. (2002) argue for $`B=0.22`$ G and $`\delta 50`$ which gives $`B\delta ^{1/3}=0.8`$ G.

Despite the differences in the underlying assumptions behind the specific homogeneous models employed a rather narrow range of values of $`B\delta ^{1/3}`$ are predicted. It is possible that source variability is the cause of the discrepancy with, perhaps, the magnetic field strength during O84 being unusually high compared with earlier epochs. Unfortunately we have no way to directly test this because we have no simultaneous observations in the gamma-ray regime with the O84 observation. The Krawczynski et al. (2002) paper does deal with simultaneous X-ray and TeV gamma-ray observations (obtained with RXTE, the Rossi X-ray Timing Explorer and HEGRA, the High Energy Gamma Ray Astronomy experiment) taken in February and 3-8th May 2000, which is close to, but not contemporaneous with the 25th May 2000 O84 observation. Appealing to source variability, however, seems contrived given the fact that observations from different epochs are returning a similar range of magnetic field strengths. The previous claims of lags on longer, $`1`$ hr, timescales lead to values of the $`B`$ and $`\delta `$ combination that are consistent with these models and this has been used as an argument in their favour on the grounds of self-consistancy. Our much tighter upper limit on the lag may, however, indicate a possible incompatibility with the homogeneous SSC model although the results are not conclusive.

It interesting that even in these short XMM-Newton observations the spectral variability behaviour of Mkn 421 is quite complex. Of special interest is the lack of any obvious spectral hysteresis in the variation of intensity with slope as has been commonly observed in a number of blazars, including Mkn 421 (e.g. Sembay et al. 1993, Takahashi et al. 1996). When observed it most commonly follows a clockwise loop in the sense of Figure 4. This behaviour naturally follows whenever the slope is controlled by cooling processes which act faster at higher energies, as is the case when synchrotron cooling dominates. If the system is in a regime where the cooling and acceleration times are equal, however, the reverse behavour is expected and the slope will track the intensity in an anticlockwise loop (Kirk, Rieger & Mastichiadis 1998).

There are a number of possible reasons why our observations show no clean examples of hysteresis. First, it may simply be that the observations are too short relative to the characteristic flaring timescale and therefore we are observing only partial segments of a given loop. This may explain the difference in the 2-10 keV spectral behaviour between O84 and O171 which both show a decay in flux. In the latter case the observation may be starting from near the peak intensity where a rapid change in slope would be expected during the initial decay. In the former we may be well below the peak intensity in a segment where we expect only a minor change in slope (see Figure 3 from Kirk, Rieger & Mastichiadis 1998). Second, there may be contributions from multiple flares at different phases of their cycle; the lack of a correlation between the soft and hard bands in O171 appears to show an independent soft flaring component. Finally, the system overall may physically be in a transition zone between characteristically clockwise and anticlockwise spectral behaviour.

Attempts have been made to monitor Mkn 421 and other BL Lacs with long observations by BeppoSAX and ASCA however, these are limited on what they could achieve on analysing short time scales by being interrupted by Earth occultation. XMM-Newton can provide the uninterrupted light curves on the short time scales required to unambiguously test for the presence of X–ray lags, but a source such as Mkn 421 needs to be monitored for longer than the observations presented in this paper (which were designed primarily for instrument calibration purposes) in order to track the source behaviour over a wider range of time scales.

The authors thank the XMM-Newton team for the huge amount of effort required to build and operate this complex and powerful instrument. Sembay, Griffiths and Turner acknowledge financial support from the UK Particle Physics and Astronomy Research Council. Edelson and Markowitz are supported by NASA grants NAG 5-7317 and NAG 5-9023.


text:
# New Constraints on Neutron Star Cooling from Chandra Observations of 3C58

## 1 Introduction

Neutron stars are macroscopic manifestations of processes that otherwise occur only in individual atomic nuclei. Formed hot in the core collapse that terminates the life of a massive star, they are supported against gravitational implosion by neutron degeneracy pressure. However, details of the interior structure of neutron stars (NSs) remain poorly understood, largely due to our incomplete understanding of the strong interaction at ultrahigh densities. In the early stages of their lives, energy loss is dominated by neutrino emission. However, the neutrino production rate is highly dependent upon the structure of the interior. In the “normal” cooling scenario, neutrino production proceeds via the modified URCA process, neutrino bremsstrahlung in the stellar crust, and plasmon decay. The residual heat diffuses from the core to the surface, manifesting itself as blackbody-like emission (modified by effects of any residual atmosphere) which peaks in the soft X-ray band. The rate at which the surface temperature declines depends critically upon the neutrino emission rate, and thus provides constraints on hadronic physics at high densities.

In some models, the energy density of nucleons in the core is sufficiently high that the production of pion or kaon condensates, or quark matter, are predicted. The presence of such exotic particles enhances the neutrino rate immensely, resulting in much more rapid cooling. Consequently, measurements of the surface temperature for NSs of known age can provide important constraints on NS interiors, and thus on our knowledge of the strong interaction. Here we present observations of 3C58, a young Crab-like nebula believed to be associated with the supernova of 1181 AD. Given this age, 3C58 is presumably powered by one of the youngest neutron stars in the Galaxy, and one for which the residual cooling emission should be measurable if standard cooling models apply – a prospect that has motivated many searches for the pulsar and its cooling emission. The pulsar in 3C58 was recently discovered by Murray et al. (2002). To date, however, the strongest upper limit on cooling from the NS surface yielded a value well in excess of standard cooling curves (Helfand et al. 1995), while Torii et al. (2000) reported hot blackbody emission associated with a small emitting area that they associated with hot polar caps - a result inconsistent with XMM studies by Bocchino et al. (2001). Here we present Chandra observations and modeling of the emission from the NS which yield a surface temperature well below that predicted by standard cooling curves, providing strong support for the presence of enhanced neutrino production rates in the stellar interior.

## 2 Observations

3C 58 was observed with the ACIS detector onboard the Chandra Observatory on 04 September 2000. The aimpoint was placed on the S3 chip of the detector, and a 1/2-subarray mode was selected in order to minimize pileup from the compact central source while still imaging most of the pulsar wind nebula (PWN). Standard cleaning of the data to remove episodes of high background resulted in a final exposure of 33.8 ks. The image clearly reveals a compact central core, as reported by Murray et al. (2001) based on Chandra HRC observations, as well as other apparent structure in the inner nebula – most notably an arc of emission apparently emanating from the central region (Figure 1).

A total of eight relatively bright point sources are observed in the field defined by the S3 and S4 ACIS chips. Based on identification of stellar counterparts for four of these, we obtained a position reconstruction error of $`\mathrm{\Delta }\mathrm{RA}=0.17\pm 0.12`$ arcsec and $`\mathrm{\Delta }\mathrm{Dec}=+0.64\pm 0.08`$ arcsec, within typical accuracies for Chandra positions.

## 3 Data Analysis

### 3.1 Spatial Modeling

The central region of 3C58 is shown in Figure 2. The emission is clearly extended, with elongation in the N-S direction, perpendicular to the long axis of the main nebula. The pulsar PSR J0205+6449 (Murray et al. 2001) resides at the center of this compact nebula. We have carried out spatial modeling of the emission using a three-component model consisting of a 2-D gaussian combined with a point source, both superposed on a flat background. For the point source, we used a model of the Chandra-ACIS point spread function for the mean position of the central region. The free parameters in the fit are the central position, widths, rotation angle, and normalization of the 2-D gaussian, the position and normalization of the PSF model, and the surface brightness of the background. Our best-fit parameters for the 2-D gaussian yield a FWHM (along the long axis) of $`4.2^{\prime \prime }\pm 0.04^{\prime \prime }`$ and an ellipticity of $`0.49\pm 0.01`$, with the long axis running almost exactly north-south. In fact, Figure 2 shows that the full extent of the elongated core is $`25^{\prime \prime }`$ in the N-S direction. Murray et al. (2002) find similar results using HRC data, though a somewhat smaller extent which may be the result of the factor of $`2.5`$ fewer counts in the HRC image, as well as the fact that the 2-D gaussian model is only a rough approximation of the more extended structure. The best-fit point source location, when corrected for the offsets described above, is RA2000: 02h05m37.92s, Dec2000: $`+64^{}`$4942.8′′; the total number of counts from this component is $`880\pm 78`$ (90% confidence).

The western edge of the core nebula lies directly along a radio filament (Frail & Moffett 1993), shown as contours in Figure 2, and there appears to be a slight flattening of nebula along this side. Frail & Moffett (1993) suggested that this filament may represent the position of the termination shock where the pulsar wind is confined by the interior pressure of the PWN. Adopting the generally accepted distance $`d=3.2d_{3.2}`$ kpc (see Section 4.1), integration of the radio synchrotron spectrum yields a pressure of $`P_{\mathrm{neb}}=3.2\times 10^{10}d_{3.2}^1\mathrm{dyne}\mathrm{cm}^2`$ under the assumption of equipartition between the electron and magnetic energy densities. The ram pressure of the pulsar wind $`P_{\mathrm{wind}}=\dot{E}/4\pi \eta cr_w^2`$, where the wind covers a fraction $`\eta `$ of a sphere, $`\dot{E}`$ is the energy loss rate of the pulsar, and $`r_w`$ is the distance from the pulsar at which wind confinement occurs. The spin-down properties of the pulsar give $`\dot{E}=2.6\times 10^{37}I_{45}\mathrm{erg}\mathrm{s}^1`$ (Murray et al. 2001) where $`I_{45}`$ is the NS moment of inertia in units of $`10^{45}\mathrm{gm}\mathrm{cm}^2`$. Pressure balance should thus occur at $`r_w=5.5\times 10^{17}I_{45}^{1/2}\eta ^{1/2}d_{3.2}^{1/2}`$ cm, or at an angular distance $`\theta =11.4I_{45}^{1/2}\eta ^{1/2}d_{3.2}^{1/2}`$ arcsec. This is in good agreement with the the $`12`$ arcsec radial extent of the core X-ray emission in the NS direction, suggesting that the compact nebula surrounding the pulsar is bounded by the pulsar wind termination shock. It is possible that this is actually a toroidal structure, much like that seen in the Crab Nebula, and that the elongated surface brightness distribution is the result of the inclination angle. In this interpretation, the axis of the toroid, which presumably lies along the rotation axis of the pulsar, lies in the east-west direction when projected onto the sky. We note that the long axis of 3C58 itself, as well as an extended jet-like feature shown in Figure 2, are both aligned in this direction. Assuming that the radio filament lies along one side of the torus, its separation from the pulsar ($`4.5`$ arcsec) implies an inclination angle of $`70^{}`$. We defer detailed discussion of this picture, including the jet-like feature in Figure 2 and other possible outflow-related structures seen in Figure 1, to a future publication.

### 3.2 Spectral Modeling

Previous measurements of the column density for 3C58 yield $`N_H=(3.3\pm 0.4)\times 10^{21}\mathrm{cm}^2`$ (Torii et al. 2000). This is consistent with values reported by Bocchino et al. (2001). Using the entire nebula, and allowing the photon index to vary with radius111A variation in spectral index with radius is expected from synchrotron cooling. These results will be presented in a future publication. we obtain $`N_H=(3.75\pm 0.11)\times 10^{21}\mathrm{cm}^2`$. Using this value, the spectrum from the central $`3\times 3`$ pixel region222Each ACIS pixel is 0.492 arcsec on a side. in 3C58 is well fit ($`\chi _r^2=0.82`$) by an absorbed power law with a photon index $`\mathrm{\Gamma }=1.73\pm 0.07`$ (Figure 3); blackbody or thermal plasma models are ruled out at high significance. In particular, contrary to results from ASCA studies (Torii et al. 2000), we find no evidence of a soft blackbody component from the central source, which is consistent with the results of Bocchino et al. (2001). The majority of equations of state yield an effective NS radius larger than 12 km for any range of masses (Haensel 2001). To establish an upper limit to the effective blackbody temperature of any such a component, we have fixed the column density at the upper end of the 90% confidence interval quoted above, added a blackbody model with $`R_{\mathrm{}}=12`$ km to the power law, and varied the temperature of this component until the spectral fit probability fell to $`10^4`$. In re-fitting the data, we allowed the power law index to vary as well. The upper limit we determine is $`T_{\mathrm{}}<1.08\times 10^6`$ K. The model containing this thermal component is shown in Figure 3, and clearly exceeds what can be accommodated in the data. For comparison, we also plot the expected flux for a blackbody component whose temperature is that expected from standard neutron star cooling curves.

As with all stars, the emission from the surface of a NS is not a blackbody; rather, it is modified by the presence of whatever atmosphere might exist. One expects the surface of the NS to be covered with Fe, but an atmosphere consisting of H, He, and/or intermediate-mass elements acquired either from ejecta fallback following the neutron star’s formation, or from material accreted from the ISM, is also a possibility. From models of nonmagnetic atmospheres, the primary effect of H or He atmospheres is a considerable deviation of the high energy end of the spectrum relative to the Wien tail of a pure blackbody. The result is that the best-fit blackbody model overestimates the effective temperature – typically by as much as a factor of two. For atmospheres dominated by heavier elements the blackbody fit gives a good approximation to the temperature (Lloyd, Hernquist, & Heyl 2002). The temperature upper limit derived above is thus conservative.

### 3.3 Spatial/Temporal Modeling

The temporal information from PSR J0205+6449 provides an additional constraint on any thermal emission originating at the surface of the NS. The spectral results summarized above clearly establish that the bulk of the flux from the pulsar is nonthermal. The narrow pulse profile measured by Murray et al. (2001) is consistent with this picture. While temperature nonuniformities on the NS surface can, in principle, lead to a pulsed X-ray signal, light bending effects from the strong gravitational field result in considerable broadening of any such pulse profile (Yancopoulos, Hamilton & Helfand 1994; Zavlin, Shibanov & Pavlov 1995) The very sharp pulse profile of PSR J0205+6449 is inconsistent with such an origin. We can thus use the strength of the pulsed signal along with the modeled brightness of the pointlike component to obtain an independent limit on any remaining emission arising from surface cooling.

In the 33 ks HRC-S observation of 3C58, the pulse profile for PSR J0205+6449 yields a total of $`127.2\pm 17.2`$ counts (90% confidence) above the mean flux level in the narrow phase bins containing the bulk of the pulsed power. The count rate conversion between ACIS-S and HRC-S based on the measured power law spectrum reported above is $`R_{\mathrm{HRC}\mathrm{S}}=0.26R_{\mathrm{ACIS}\mathrm{S}}`$. The lower limit to the number of pulsed counts in the $`34`$ ks ACIS-S observation is thus $`435`$. Our spatial modeling yields a total of $`880\pm 78`$ counts from the point source in ACIS-S. The 90% confidence upper limit to the number of unpulsed counts is thus $`523`$, for an ACIS-S count rate of $`1.54\times 10^2\mathrm{ct}\mathrm{s}^1`$. For a radius at infinity of 12 km, the upper limit to the blackbody temperature is $`1.13\times 10^6`$ K, consistent with the limit obtained through spectral modeling. As we discuss below (and show in Figure 4), this limit falls well below values predicted for standard cooling models assuming the historical age of 3C58. We note that the limit depends on uncertainties in the distance and column density. However, even at the most extreme values of acceptable column density based on fits to the entire nebula, the distance required to match the standard cooling predictions is $`6`$ kpc, a value inconsistent with that inferred from H i measurements.

## 4 Discussion

Our measurement of an upper limit to the temperature of the young NS in 3C58 provides the strongest constraint on NS cooling yet derived, and poses significant problems for standard cooling models. We explore here the implications of our result for the nuclear equation of state and the structure and evolution of NS.

### 4.1 Neutron Star Cooling

The cooling rate of isolated NSs has been a subject of considerable theoretical work predating even the discovery of the first pulsars (e.g. Bahcall & Wolf 1965a,b). The poorly understood properties of the strong nuclear potential at the densities found in NS interiors make these calculations difficult, and lead to a wide range of predictions based on different assumptions for the equation of state, composition, and details of superconductivity (see, e.g., the review by Tsuruta 1998, and references therein). While there is a clear consensus that the early cooling proceeds via neutrino emission from the NS core, the time scale over which this dominates depends critically on the neutrino production rate which, in turn, can vary by orders of magnitude depending upon the state of matter in the interior.

Broadly speaking, models can be divided into “standard” and “non-standard” cooling scenarios. The former class is based upon neutrino emission via the modified URCA process, n–n and p–p neutrino bremsstrahlung, crust neutrino bremsstrahlung, and plasmon neutrino processes. Assumptions must still be made in the cooling calculations, such as the form of the the nucleon-nucleon force and the many-body technique used to determine the equation of state, but the differences between stiff and soft equations of state are not large for this class of models. In Figure 4, we plot a variety of cooling curves associated with different models for the NS interior and its properties (see Page 1998 and references therein)333Tabulated values are available at http://www.astroscu.unam.mx/neutrones/NS-Cooler/NS-Cooler.html. The solid curve corresponds to “standard” cooling using an equation of state of moderate stiffness.

Non-standard cooling models incorporate neutrino emissivities associated with other processes that may operate in NS interiors. These include the presence of pion condensates which may form at sufficiently high densities. The resulting pion-induced beta decay leads to very a high neutrino emissivity and a correspondingly shorter cooling time for the NS interior (Bahcall & Wolf 1965a,b, Maxwell et al. 1977). Similar processes involving kaon condensates or quark matter may operate as well. Alternatively, strong magnetic field effects in the crust may allow the normal URCA process to proceed, which also leads to extremely high neutrino production rates. These exotic cooling mechanisms modify the NS cooling curves substantially. The dashed curves in Figure 4 represent approximations for several nonstandard cooling models to illustrate the associated rapid cooling (Page 1998). The effects of superfluidity can substantially moderate the rapid cooling because the significantly reduced heat capacity of the superfluid particles reduces the neutrino rate considerably.

We have indicated on Figure 4 the upper limit for surface cooling emission from PSR J0205+6449 derived above. It is evident that the limit falls considerably below predictions from standard cooling models, suggesting the presence of some exotic cooling contribution in the interior. We note that the placement of this point on the cooling curves depends upon the age assumption for 3C58. Here we have used the historical age based on association with SN 1181. Four independent written records – two from Japan, and one each from the northern and southern Chinese empires – record the presence of a guest star that appeared in the late summer of AD 1181 and was visible to the naked eye for nearly 6 months. Detailed analysis of the star’s celestial coordinates are summarized by Stephenson and Green (1999) who identify a faint star near $`ϵ`$ Cas as the closest permanent star to the location of the event. This position lies less than $`2^{}`$ from 3C58. The only other young supernova remnant within $`10^{}`$ is Tycho’s SNR, the unambiguously identified remnant of SN 1572. The association of 3C58 with SN 1181 thus appears extremely secure. However, even if we increase the age estimate to 5000 yr, as has been suggested based on optical and radio expansion data (Fesen 1985, Bietenholz et al. 2001), the measurement still falls below standard cooling models (see Figure 4). Deeper high-resolution X-ray measurements of 3C58 are of considerable interest in order to further constrain this important diagnostic of the NS interior.

## 5 Summary

The size of the extended emission region immediately surrounding the pulsar in 3C58 is consistent with that expected for the termination shock where the pulsar wind pressure is balanced by the interior pressure of the nebula. The elongated morphology, along with the position and size of an apparently associated radio filament, suggest a toroidal morphology similar to that observed in the Crab Nebula. The orientation of the torus suggests a pulsar rotation axis aligned in the east-west direction, along the primary axis of extension of 3C58 itself. The inclination angle of the torus is large, and the resulting projected X-ray morphology does not permit us to resolve the toroidal structure with the current X-ray data.

Our Chandra observations of 3C58 and its associated pulsar provide important constraints on models for NS cooling. The upper limit to the thermal emission from the NS surface falls well below predictions from standard cooling models and appears to require the presence of some rapid cooling mechanism. Additional time-resolved observations of the pulsar could further constrain this limit, providing important new information on the equation of state of nuclear matter and the structure of NS interiors.

The authors would like to thank John Bahcall for helpful discussions on NS cooling and Dale Frail for supplying the VLA image of the central region in 3C58. We also thank Bryan Gaensler for his helpful suggestions and careful reading of the text. This work was supported in part by the National Aeronautics and Space Administration through contract NAS8-39073 and grant GO0-1117A (POS) and contract NAS8-38248 (SSM).


text:
# Specific heat and thermal conductivity in the mixed state of MgB2

## Abstract

The specific heat $`C`$ and the electronic and phononic thermal conductivities $`\kappa _e`$ and $`\kappa _{ph}`$ are calculated in the mixed state for magnetic fields $`H`$ near $`H_{c2}`$. The effects of supercurrent flow and Andreev scattering of the Abrikosov vortex lattice on the quasiparticles are taken into account. The resulting function $`C(H)`$ is nearly linear while $`\kappa _e(H)`$ exhibits an upward curvature near $`H_{c2}`$. The slopes decrease with impurity scattering which improves the agreement with the data on MgB$`_2`$. The ratio of phonon relaxation times $`\tau _n/\tau _s=g(\omega _0,H)`$ for phonon energy $`\omega _0`$, which is nearly a step function at $`\omega _0=2\mathrm{\Delta }`$ for the BCS state, is smeared out and tends to one for increasing $`H`$. This leads to a rapid reduction of $`\kappa _{ph}(H)`$ in MgB2 for relatively small fields due to the rapid suppression of the smaller energy gap.

Numerous experiments indicate that the superconducting state of MgB2 ($`T_c=40`$K)  is a conventional s-wave pairing state mediated by the electron-phonon interaction. Nevertheless, the theoretical explanation of physical quantities is demanding because two different energy gaps are formed where the larger gap $`\mathrm{\Delta }_1`$ is associated with the nearly cylindrical $`\sigma `$-sheets of the Fermi surface and the smaller gap $`\mathrm{\Delta }_2`$ with the three-dimensional $`\pi `$-sheets. Analysis of the temperature dependence of the measured thermal conductivity in the basal plane of MgB$`_2`$  gives supporting evidence for two different gaps on different Fermi surface sheets. The situation is complicated because the thermal conductivity $`\kappa `$ is a sum of an electronic part $`\kappa _e`$ and the phononic or lattice conductivity $`\kappa _{ph}`$. In Ref.  $`\kappa _e`$ and $`\kappa _{ph}`$ were calculated with help of the BCS theory  where $`\kappa _e`$ is limited by elastic scattering (scattering by phonons is omitted) and $`\kappa _{ph}`$ is limited by electron scattering. This theory for $`\kappa _{ph}`$ is supplemented by adding to the phonon relaxation rate due to scattering by electrons, the relaxation rates due to scattering by point defects  and other defects. We remark that the peak in $`\kappa _{ph}`$ below $`T_c`$ in the cuprate superconductors has been fitted by including scattering by sample boundaries, sheet-like faults, and dislocations .

Recently, the ab-plane thermal conductivity of MgB2 has been measured as a function of magnetic field $`H`$ with orientations both parallel and perpendicular to the c-axis . At low temperatures $`\kappa (H)`$ drops steeply for increasing $`H`$ up to a relatively low field which we denote by $`H_{c2}^{(2)}`$ ($``$1 kOe), and then it rises continuously up to $`H_{c2}`$ ($``$ 30 kOe for $`𝐇𝐜`$). The first drop of $`\kappa `$ is interpreted as due to the behavior of $`\kappa _{ph}(H)`$ which is caused by a strong suppression of the smaller gap $`\mathrm{\Delta }_2`$ by relatively small fields. This leads to a rapid increase in the number of those quasiparticles that dominate the scattering of phonons . The rapid suppression of the smaller gap can be explained in terms of a two-band model with different energy gaps where the smaller gap is induced by Cooper pair tunneling from the band with the larger gap . The observed increase of $`\kappa (H)`$ above $`H_{c2}^{(2)}`$ is attributed to the increase of $`\kappa _e(H)`$ because $`\kappa _{ph}`$ has already reached its normal state value . Since this field dependence of $`\kappa _e(H)`$ at low temperatures is qualitatively similar to that of the specific heat coefficient $`\gamma (H)`$ , it is concluded that the growth of $`\kappa _e(H)`$ is due to the rapid increase of quasiparticles in the vortex cores associated with the larger energy gap.

In this Letter we present theories for the specific heat and the electronic and phononic thermal conductivities in the vortex state for applied fields near the upper critical field $`H_{c2}`$. For simplicity we consider only one isotropic s-wave pairing gap either on a cylindrical Fermi surface with $`𝐇𝐜`$, or on a sperical Fermi surface. Our goal is to explain the measured quantities $`\gamma (H)`$ and $`\kappa _e(H)`$ near $`H_{c2}`$ ($``$ 30 kOe), and $`\kappa _{ph}(H)`$ below the effective upper critical field $`H_{c2}^{(2)}`$ ($``$ 1 kOe) for the vortex lattice associated with the smaller energy gap. Our theories are based on the normal and anomalous Green’s functions $`G`$ and $`F`$ obtained from the Gorkov integral equations with kernels given by the product of Abrikosov vortex lattice order parameters $`\mathrm{\Delta }(𝐫_\mathrm{𝟏})\mathrm{\Delta }^{}(𝐫_\mathrm{𝟐})`$ and the phase factor due to the magnetic field . These Green’s functions depend sensitively on $`\mathrm{sin}\theta `$ where $`\theta `$ is the angle between the quasiparticle momentum $`𝐩`$ and the magnetic field $`𝐇`$. For $`\theta 0`$, $`G`$ and $`F`$ tend to the Green’s functions for a BCS superconductor while, for $`\theta \pi /2`$, they take into account the effects of supercurrent flow and Andreev scattering due to the vortex lattice. This theory has been applied previously to calculate $`\kappa _e(T,H)`$ for energy gaps with line nodes like those occuring in the high-$`T_c`$ cuprates and Sr2RuO4 .

A simplified version of the theory in Ref. , which has been derived from the Eilenberger equations , yields the following expression for the spatial average of the density of states:

$$N(\omega ,\theta )/N_0A(\mathrm{\Omega };\stackrel{~}{\mathrm{\Delta }},\theta )=\text{Re}\left[1+\frac{8\stackrel{~}{\mathrm{\Delta }}^2}{\mathrm{sin}^2\theta }\left[1+i\sqrt{\pi }zw(z)\right]\right]^{1/2}.$$
(1)
The quantities appearing in the equation are defined as follows: $`z=2[\mathrm{\Omega }\stackrel{~}{\mathrm{\Delta }}+i(\mathrm{\Lambda }/v)\mathrm{\Gamma }]/\mathrm{sin}\theta `$ with $`\theta =\mathrm{}(𝐩,𝐇)`$, $`\mathrm{\Lambda }=(2eH)^{1/2}`$, $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }\mathrm{\Lambda }/v`$, $`\mathrm{\Omega }=\omega /\mathrm{\Delta }`$, $`\stackrel{~}{\mathrm{\Delta }}^2=(H_{c2}H)/6\beta _AH`$, $`\mathrm{\Delta }^2=\mathrm{\Delta }_0^2(T)\left[1(H/H_{c2})\right]`$, and $`w(z)=\mathrm{exp}(z^2)\mathrm{erfc}(iz)`$. Here $`v`$ is the Fermi velocity, $`\beta _A`$ the Abrikosov parameter which we take as 1.2, $`\mathrm{\Delta }_0(T)`$ is the BCS gap, and $`\mathrm{\Gamma }`$ is the normal-state impurity scattering rate. In the Born and unitary scattering limits, $`\mathrm{\Gamma }`$ has to be replaced by $`\mathrm{\Gamma }A`$ or $`\mathrm{\Gamma }/A`$, respectively, where $`A(\mathrm{\Omega })`$ is calculated self consistently. The specific heat C is given by $`C=N_0T_0^{\mathrm{}}𝑑xx^2\text{sech}^2(x/2)A(Tx/\mathrm{\Delta })`$. At low temperatures ($`T\mathrm{\Delta }`$), $`C=\gamma _s(H)T`$, where the specific heat coefficient $`\gamma _s(H)`$ is proportional to $`A(\mathrm{\Omega }=0)`$. In Fig. 1 we have plotted $`A(\mathrm{\Omega }=0)`$ versus $`H/H_{c2}`$ for $`\mathrm{sin}\theta =1`$ corresponding to a cylindrical Fermi surface and $`𝐇𝐜`$, and impurity scattering rates $`\delta \mathrm{\Gamma }/\mathrm{\Delta }_0=`$0.1, 0.2, and 0.5. A is seen to be nearly linear near $`H_{c2}`$ with a slope at $`H_{c2}`$ that decreases for increasing $`\delta `$. In Fig. 1 we also show our results for the angular average $`\overline{A}=_0^{\pi /2}𝑑\theta \mathrm{sin}\theta A(\theta )`$ which corresponds to a three dimensional Fermi surface. Comparison with the solid curves in Fig. 1 shows that the slopes at $`H_{c2}`$ are decreased by the angular average. It should be pointed out that our results are only strictly valid in the vicinity of $`H_{c2}`$. However, solution of the Eilenberger equations for a vortex lattice shows that the spatial average of the resulting density of states is well approximated by $`A(\mathrm{\Omega })`$ for fields down to about $`(1/2)H_{c2}`$ . The solid curve of $`A`$ vs $`H/H_{c2}`$ for $`\delta =0.5`$ in Fig. 1 agrees qualitatively with the measured field dependence of the specific heat coefficient $`\gamma _s(H)`$ at low temperatures and fields near $`H_{c2}`$ . In the unitary impurity scattering limit $`A`$ increases slightly for decreasing field, however at lower fields our theory does not apply anyway.

We turn now to the theory of the electronic thermal conductivity $`\kappa _e`$ in the vortex state near $`H_{c2}`$ which has been developed in Ref. . These expressions are easily modified to apply to an isotropic s-wave pairing state and a field along the c-axis. Employing the expression valid at low temperatures we obtain for the ratio $`\kappa _{es}/\kappa _{en}`$ as a function of $`H/H_{c2}`$ the plots shown in Fig. 2 for constant $`\mathrm{sin}\theta =1`$ and impurity scattering rates $`\delta =`$0.1, 0.2, and 0.5. One sees that these plots exhibit upward curvatures with decreasing slopes at $`H_{c2}`$ for increasing values of $`\delta `$. In Fig. 2 we have also plotted the corresponding results for the angular averages over $`\theta `$ of $`\kappa _{es}/\kappa _{en}`$. Comparison with the measured $`\kappa (H)`$ near $`H_{c2}`$, which is presumably dominated by the field dependence of $`\kappa _e(H)`$ , shows that the measured upward curvature towards $`H_{c2}`$ is qualitatively best described by the upper solid curve in Fig. 2 for the relatively large impurity scattering rate $`\delta =0.5`$. It should be noted that the renormalization of the impurity scattering rate in the unitary limit, $`\mathrm{\Gamma }\mathrm{\Gamma }/A`$, slightly increases the upward curvature of $`\kappa _{es}/\kappa _{en}`$ in comparison to that obtained without renormalization or in the Born scattering limit. This is in somewhat better agreement with the data. Of course our theory cannot explain the measured field dependence of $`\kappa _e(H)`$ at lower fields which first rises steeply and then saturates until the upward curvature sets in at about $`(1/2)H_{c2}`$. Our theory also cannot describe the measured anomalous temperature dependence $`\kappa _eT^{1/2}`$ for fields below approximately $`(1/2)H_{c2}`$: $`\kappa _{en}`$ and thus $`\kappa _{es}`$ is proportional to $`T/\mathrm{\Gamma }`$.

It was argued in Ref.  that the similar field dependencies of the specific heat and electronic thermal conductivity can be explained in terms of the relationship $`\kappa _e=C_ev_F\mathrm{}/3`$ where $`\mathrm{}`$ is the mean free path. In fact, our theoretical expression for $`\kappa _e`$ has a similar form if one sets $`\mathrm{}=v_F\tau _e`$ where $`\tau _e=1/\text{Im}\xi _0`$ with $`\text{Im}\xi _0=\gamma +\gamma _A`$ . Here $`\gamma `$ is the impurity scattering rate and $`\gamma _A`$ is the scattering rate due to Andreev scattering of the quasiparticles by the vortices. $`\gamma _A`$ is shown as a function of angle $`\theta `$ for different values of $`\mathrm{\Omega }`$ and $`\stackrel{~}{\mathrm{\Delta }}`$ in Ref. . These plots show that, for extended states ($`\mathrm{\Omega }1`$), $`\gamma _A`$ increases from zero as $`\theta `$ increases from 0 to $`\pi /2`$. This means that the Andreev scattering rate has a maximum for quasiparticles moving perpendicular to the vortex axis. For increasing field (decreasing $`\stackrel{~}{\mathrm{\Delta }}`$), $`\gamma _A`$ decreases which has the effect that the curve for $`\kappa _{es}/\kappa _{en}`$ exhibits an upward curvature in comparison to $`C_e`$ which has a slight downward curvature for increasing field.

Until now there exists only a phenomenological theory for the phonon heat transfer in the mixed state where the phonons are scattered by vortices consisting of normal state cylinders . We develop here the microscopic theory of $`\kappa _{ph}`$ (limited by electron scattering) in the vortex state in close analogy to the theory for the BCS state . In that case the sum of probabilities for absorption and emission of phonons by quasiparticles yields a relaxtion rate $`1/\tau _s`$ for a phonon of energy $`\omega _0`$ which is proportional to the integral over quasiparticle energy $`E`$ of the expression $`|EE^{}/\epsilon \epsilon ^{}|[1\mathrm{\Delta }^2/EE^{}][f(E/T)f(E^{}/T)]`$. Here $`E/\epsilon =E/(E^2\mathrm{\Delta }^2)^{1/2}`$ is the density of states, $`E^{}=E+\omega _0`$, and $`f(E/T)`$ is the Fermi function. In the mixed state near $`H_{c2}`$, the density of states $`E/\epsilon `$ is replaced by $`A(\mathrm{\Omega })`$ as given in Eq. (1). The coherence term $`(\mathrm{\Delta }/\epsilon )(\mathrm{\Delta }/\epsilon ^{})`$, which arises from the matrix elements for quasiparticle scattering, is replaced by $`FF^{}`$ where the anomalous Green’s function $`F`$ is given in Ref. . The spectral function of $`F`$ yields the following analog of $`\mathrm{\Delta }/\epsilon `$:

$$B(\mathrm{\Omega };\stackrel{~}{\mathrm{\Delta }},\theta )=\text{Re}\left[\frac{i\sqrt{\pi }\mathrm{\hspace{0.17em}2}(\stackrel{~}{\mathrm{\Delta }}/\mathrm{sin}\theta )w(z)}{\left\{1+(8\stackrel{~}{\mathrm{\Delta }}^2/\mathrm{sin}^2\theta )[1+i\sqrt{\pi }zw(z)]\right\}^{1/2}}\right].$$
(2)
The argument $`z`$ was defined below Eq.(1). The functions $`A(\mathrm{\Omega })`$ and $`B(\mathrm{\Omega })`$ are even and odd in $`\mathrm{\Omega }`$. For $`\theta 0`$ they tend to the BCS functions $`\omega /(\omega ^2\mathrm{\Delta }^2)^{1/2}`$ and $`\mathrm{\Delta }/(\omega ^2\mathrm{\Delta }^2)^{1/2}`$, respectively. In this way we obtain for the ratio of phonon relaxation times in the normal and superconducting states (denoted by g in Refs. , , and ):

$`\tau _n/\tau _s=g(\mathrm{\Omega }_0)`$ $`=`$ $`\left[1\mathrm{exp}(\mathrm{\Omega }_0\mathrm{\Delta }/T)\right](2/\mathrm{\Omega }_0){\displaystyle _0^{\mathrm{}}}𝑑\mathrm{\Omega }f[(\mathrm{\Omega }\mathrm{\Omega }_0/2)\mathrm{\Delta }/T]f[(\mathrm{\Omega }+\mathrm{\Omega }_0/2)\mathrm{\Delta }/T]`$ (4)
$`\times \left[A(\mathrm{\Omega }\mathrm{\Omega }_0/2)A(\mathrm{\Omega }+\mathrm{\Omega }_0/2)B(\mathrm{\Omega }\mathrm{\Omega }_0/2)B(\mathrm{\Omega }+\mathrm{\Omega }_0/2)\right].`$
The quantity $`\mathrm{\Omega }_0=\omega _0/\mathrm{\Delta }`$ is the phonon energy $`\omega _0`$ divided by the effective gap $`\mathrm{\Delta }=\mathrm{\Delta }(H)`$. For $`T\mathrm{\Delta }`$, the function $`ff`$ in the integrand of Eq.(4) is approximately one in the range from $`\mathrm{\Omega }=0`$ to $`\mathrm{\Omega }_0/2`$ and zero above $`\mathrm{\Omega }_0/2`$. In Fig. 3 we show some examples of the function $`g(\mathrm{\Omega }_0)`$ for parameter values $`\theta =\pi /2`$, $`\delta =0.1`$, and $`\stackrel{~}{\mathrm{\Delta }}=`$0.1, 0.2, 0.3, and 0.6 corresponding to $`H/H_{c2}=`$0.93, 0.78, 0.61, and 0.28. For increasing $`\stackrel{~}{\mathrm{\Delta }}`$, or decreasing field, the function $`g`$ tends to the BCS step function which is zero in the range from $`\mathrm{\Omega }_0=0`$ to 2, and $`\pi /2`$ for $`\mathrm{\Omega }_0>2`$ (see Fig. 1 of Ref. ). The physical meaning of $`g`$ at low temperatures for the BCS superconductor is that the minimum energy of a phonon for creating a pair of quasiparticles is $`\omega _0=2\mathrm{\Delta }`$. Fig. 3 shows that, for increasing field, the relaxation rate is more and more smeared out and the ratio $`\tau _n/\tau _s`$ tends to one. A similar effect occurs for increasing impurity scattering rate $`\delta `$.

The expression for the phonon thermal conductivity limited by electron scattering and several other scattering processes is given by

$$\kappa _{ph}=At^3_0^{\mathrm{}}𝑑x\frac{x^4\text{e}^x}{(\text{e}^x1)^2}\left[1+\alpha t^4x^4+\beta t^2x^2+\delta tx+\gamma txg(xT/\mathrm{\Delta })\right]^1.$$
(5)
Here $`x=\omega _0/T`$, $`t=T/T_c`$, and the coefficients $`A`$, $`\alpha `$, $`\beta `$, $`\delta `$, and $`\gamma `$ refer to scattering by sample boundaries, point defects, sheet-like faults, dislocations, and quasiparticles, respectively. For $`t1`$ and $`T\mathrm{\Delta }`$ the argument of $`g(\mathrm{\Omega }_0)`$ in Eq.(5) is approximately equal to zero. Taking the limit $`\mathrm{\Omega }_00`$ in Eq.(4) and noting that the coherence term B vanishes at $`\mathrm{\Omega }_0=0`$, we obtain $`g(\mathrm{\Omega }_0=0)=[A(\mathrm{\Omega }_0=0,\stackrel{~}{\mathrm{\Delta }},\theta )]^2`$. For constant $`\theta =\pi /2`$ this expression yields the field dependence of $`\tau _n/\tau _s`$ for the case $`𝐇𝐜`$ and a cylindrical Fermi surface. For a three dimensional Fermi surface, which is presumably more appropriate for the quasiparticles associated with the smaller gap in MgB2, we have to multiply $`[A(0,\theta )]^2`$ by $`\mathrm{sin}\theta `$ and integrate over $`\theta `$ from 0 to $`\pi /2`$. In Fig. 4 we have plotted our results for $`g`$ at $`\theta =\pi /2`$ and for the angular average of $`g`$ vs $`H/H_{c2}`$ for impurity scattering rates $`\delta `$=0.1 and 0.5.

The measured $`\kappa _{ph}`$ at $`H=0`$  has been analyzed in terms of an expression which is equivalent to our Eq.(5) , apart from the term proportional to $`\beta `$. By inserting two different contributions for the electron scattering term $`\gamma `$ corresponding to different gaps $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$, the fitting procedure indicated that the dominant part of the scattering of phonons by electrons is provided by that part of the electronic excitation spectrum experiencing the smaller gap, $`\mathrm{\Delta }_2`$. According to the two-band theory of Ref. , this smaller gap is suppressed by a relatively small effective upper critical field $`H_{c2}^{(2)}`$ which leads to the measured fast drop of $`\kappa _{ph}(H)`$ as $`H`$ increases from 0 to $`H_{c2}^{(2)}`$ . We conclude that the measured field dependence of $`\kappa _{ph}(H)`$ below $`H_{c2}^{(2)}`$ at low temperatures can be explained by inserting in Eq.(5) the values of the various constants obtained for $`H=0`$ and, for $`\tau _n/\tau _s=g(H)`$, the angular average of $`g(\mathrm{\Omega }_0=0)=[A(\mathrm{\Omega }_0=0,\stackrel{~}{\mathrm{\Delta }},\theta )]^2`$. The behavior of $`g(H)`$, which determines the form of the fast reduction of $`\kappa _{ph}(H)`$, depends on the impurity scattering rate $`\delta `$ (see Fig. 4). We note that the model of phonon scattering by vortices consisting of normal state cores  yields a linear dependence of $`g`$ on $`H/H_{c2}`$.

Our results have been derived for an isotropic s-wave gap. Recently an anisotropic s-wave gap $`\mathrm{\Delta }(𝐩)/\mathrm{\Delta }=(1+a\mathrm{cos}^2\theta )^{1/2}`$ has been proposed which can explain the anisotropy of $`H_{c2}`$ in MgB$`_2`$ . For this gap with reasonable parameter values of a = 10 and 20 we find that the effect of this gap anisotropy on the field dependence of $`C`$ and $`\kappa _{es}/\kappa _{en}`$ is similar to that of impurity scattering. For instance, the curves of $`C`$ and $`\kappa _{es}/\kappa _{en}`$ versus $`H/H_{c2}`$ for $`\delta =0.1`$ and $`a=20`$ appear to yield better fits to the experimental points in Figs. 1 and 2 than the dashed curves for $`a=0`$ and $`\delta =0.5`$.

The other important effect on the upward curvature of $`\kappa _{es}/\kappa _{en}`$ arises from the decrease of the Andreev scatering rate $`\gamma _A`$ with increasing field. This decrease is given approximately by the expression $`\gamma _A=\sqrt{\pi }(\stackrel{~}{\mathrm{\Delta }}^2/\mathrm{sin}\theta )\text{exp}[(2\mathrm{\Omega }\stackrel{~}{\mathrm{\Delta }}/\mathrm{sin}\theta )^2]`$ where $`\stackrel{~}{\mathrm{\Delta }}^2(H_{c2}H)/H`$. It is interesting that the exponent in this expression, in particular the term $`(\mathrm{sin}\theta )^2`$, is similar to the exponent in the scattering rate $`1/\tau _v`$ which has been derived previously . In that model a quasiparticle is converted into a hole due to Andreev reflection by vortex screening currents. This process corresponds just to the imaginary part of the self energy calculated in Ref. .

Our theory is strictly valid only for fields near the upper critical field $`H_{c2}`$. However, recent numerical solutions of the quasiclassical Eilenberger equations for the Abrikosov vortex lattice  have shown that the analytical expressions for the density of states  (here the function $`A(\mathrm{\Omega })`$) yield at low energies impressively good results over the whole range of fields down to $`H_{c1}`$. These results for low energies are important for the calculation of thermodynamic quantities. We expect that this is also true for the function $`B(\mathrm{\Omega })`$, the spectral function of the anomalous Green’s function.

In conclusion, we have employed microscopic theories for the electronic specific heat $`C(H)`$ and the thermal conductivity $`\kappa _e(H)`$ and developed a microscopic theory for the phononic thermal conductivity $`\kappa _{ph}(H)`$ to explain the measured field dependence of these quantities in the mixed state of MgB$`_2`$. We find that impurity and Andreev scattering, as well as Fermi surface and gap anisotropy, give rise to large effects. For sufficiently large impurity scattering rate, or gap anisotropy, one can approximately fit the data for $`C(H)`$ and $`\kappa _e(H)`$ for MgB2 over a broad field range below $`H_{c2}`$. Our new theory for $`\kappa _{ph}(H)`$ is capable of describing the observed rapid reduction of the total $`\kappa `$ of MgB2 for small fields if this is actually caused by suppression of a second smaller energy gap.

We thank T. Dahm for helpful discussions.