text:
# Depairing currents in the superconductor/ferromagnet proximity system Nb/Fe

## I Introduction

A still relatively little explored area of research in non-equilibrium superconductivity concerns phenomena involving spin-polarized quasiparticles. Pioneering work on spin-polarized tunneling in conventional s-wave superconductors was performed by Meservey and Tedrow , who studied different ferromagnets (F) in F/Al2O3Al tunnel junctions and found that the tunnel current can show varying degrees of spin polarization. More recently, experiments were performed by different groups in order to establish whether superconductivity can be suppressed by injecting spin-polarized quasiparticles . In these cases the combinations existed of a d-wave high-Tc superconductor (XBa2Cu3O7, with X = Y, Dy) and a fully spin-polarized ferromagnetic manganite (A0.67B0.33MnO3 with A = La, Nd and B = Ca, Sr), either with or without a barrier of a different oxide; measured was the change in the zero-field critical current density $`J_{c0}`$ of the superconducting films upon applying a current bias through the ferromagnet. The results are not fully conclusive, and certainly not quantitative. Although generally a suppression of $`J_{c0}`$ was observed, heating effects could not always be fully ruled out since the manganites are highly resistive metals (see the discussion in ref. ), and the geometry did not always allow to determine the area of the current injection, and therefore the injected current density. Moreover, since $`J_{c0}`$ in high Tc superconductors generally is not the depairing current $`J_{dp}`$ but involves flux motion, $`J_{c0}`$ is not a direct measure for the amount of depression of the superconducting order parameter. Similar experiments have to our knowledge not been performed with combinations of conventional metals, although that would have some clear advantages. The interpretation of results would not be complicated by e.g. inhomogeneous currents in the ferromagnet or anisotropic gaps in the superconductor; lithographic techniques could be brought to bear in order to have well defined superconducting bridges and injection contacts; and it should be possible to identify the effects of the spin-polarized quasiparticles on $`J_{dp}`$.

Still, two points deserve special interest. The first is that, in planning such an F/I/S experiment, there is the potential problem of insufficient knowledge of the tunneling process. This was already apparent in the work of Meservey and Tedrow cited above  since the experiments always showed a positive sign for the spin polarization, even in the cases of e.g. Co and Ni where a negative sign was expected. Recently, this was explained by demonstrating that the choice of barrier material can strongly influence and even reverse the spin polarization of the tunneling current , with obvious consequences for the interpretation of the injection experiments. It may be advantageous to also contemplate an (F or N)/I/F/S configuration; in this case the barrier is only used to increase the energy of the electrons coming from an N or F contact, while the polarization now takes place in a thin F layer between barrier and superconductor. The disadvantage here is that the F layer in connection with the superconductor will suppress the order parameter and therefore $`J_{dp}`$ in the S layer. Still, since the proximity effect (which gives rise to the depression) for S/F systems is understood reasonably well , the effect on $`J_{dp}`$ may also be quantifiable. The second point for consideration is that even in the case of conventional superconductors the determination of $`J_{dp}`$ need not be straightforward. The difficulty lies in the fact that the superconducting bridge must have a width of no more than both the superconducting penetration depth $`\lambda `$ and the superconducting coherence length $`\xi `$. The first is needed to avoid current pile-up near the edges (as a consequence of screening of the self-field), the second is required in order to avoid vortex nucleation and flow, which gives rise to dissipation before $`J_{dp}`$ is reached. These conditions can be met e.g. for Al, which has a Bardeen-Cooper-Schrieffer (BCS) coherence length $`\xi _0`$ of about 1.5 $`\mu m`$, while $`\lambda `$ can also be made of the order of 1 $`\mu m`$ by making the film thin enough. For Al-bridges of less than 1 $`\mu m`$ wide it was shown by Romijn et al.  that the measured $`J_{dp}`$ agreed very well with the theoretical calculations by Kuprianov and Lukichev  based on the Eilenberger equations and therefore valid in the whole temperature regime below $`T_c`$. For a material such as Nb, with $`\xi _0`$ and $`\lambda `$ of the order of 50 nm, such agreement need not be expected.

In this paper we show that, at least close to $`T_c`$, the values of the zero-field critical current $`J_{c0}`$ measured on bridge-structured Nb samples are essentially the values expected for the depairing current. Furthermore, we measure the depression of $`J_{c0}`$ in trilayers of Fe/Nb/Fe as function of the thickness $`d_{Nb}`$ of the Nb layer. We compare the behaviour of $`J_{c0}(d_{Nb})`$ with the behaviour of $`T_c(d_{Nb})`$, and also with calculations of the proximity effect and the pairbreaking velocity using a two-parameter formalism based on the Usadel equations. We find that $`J_{c0}(d_{Nb})`$ is well described by the same two parameters which describe the behaviour of $`T_c(d_{Nb})`$. The conclusion is that the suppression of the depairing current as a consequence of the depression of the order parameter in S/F structures can be well described by proximity effect theory, making (F,N)/I/F/S injection experiments a distinct possibility.

## II Depairing current; theory

Close to $`T_c`$, the classical Ginzburg-Landau (GL) result for the temperature dependence of the depairing current of a thin film, under the assumption of a homogeneous superconducting order parameter over the film thickness, is given by :

$$J_{dp}^{GL}(t)=J_{dp}^{GL}(0)(1t)^{3/2},$$
(1)
with $`t=T/T_c`$. The prefactor $`J_{dp}`$ is of the order of $`H_c/\lambda `$, with $`H_c`$ the thermodynamic critical field, and will be given more precisely below. For arbitrary temperatures, calculations were performed by Kupriyanov and Lukichev, who essentially solved the Eilenberger equations for a superconductor carrying a current, with the velocity of the condensate leading to a phase gradient . Their results recover the GL behaviour near $`T_c`$ :

$`J_{dp}^{GL}(t)`$ $`=`$ $`{\displaystyle \frac{16}{9\sqrt{7\zeta (3)}}}\sqrt{\chi (\rho _d)}[eN(0)v_Fk_BT_c]`$ (2)
$`\times (1t)^{3/2}.`$
Here, the constants have their usual meaning, $`N(0)`$ is the density of states at the Fermi level per spin direction, and $`\chi (\rho _d)`$ is the G’orkov function controlled by the ’dirt parameter’ $`\rho _d=(\mathrm{}v_F)/(2\pi k_BT_c\mathrm{}_e)`$, with $`\mathrm{}_e`$ the electronic mean free path. In the dirty limit, ( $`\rho _d\mathrm{}`$) $`\chi (\rho _d)1.33\mathrm{}_e/\xi _0`$, this becomes :

$$J_{dp}^{GL}(0)=\mathrm{\hspace{0.33em}1.26}[eN(0)v_F\mathrm{\Delta }(0)]\sqrt{\frac{\mathrm{}_e}{\xi _0}}$$
(3)
which is equivalent to the expression given by Romijn et al.  :

$$J_{dp}^{GL}(0)=\frac{16\pi ^2\sqrt{2\pi }}{63\zeta (3)}[eN(0)v_Fk_BT_c]\sqrt{\frac{k_BT_c\mathrm{}_e}{\mathrm{}v_F}},$$
(4)
which can also be written in terms of experimental parameters as

$$J_{dp}^{GL}(0)=7.84\left[\frac{(k_BT_c)^3}{e^2\mathrm{}v_F(\rho \mathrm{})}\frac{1}{\rho }\right]^{1/2}.$$
(5)
This way of writing also emphasizes the proportionality $`J_{dp}^{GL}(0)\sqrt{1/\rho }`$, since the product $`\rho \mathrm{}`$ is a materials constant. At low temperatures the value of $`J_{dp}`$ saturates, reaching a zero-temperature value of :

$`J_{dp}(0)`$ $`=`$ $`1.491eN(0)\sqrt{{\displaystyle \frac{D}{\mathrm{}}}}\mathrm{\Delta }^{3/2}(0)`$ (6)
$`=`$ $`0.486[eN(0)v_F\mathrm{\Delta }(0)]\sqrt{{\displaystyle \frac{\mathrm{}_e}{\xi _0}}},`$
with $`D=1/3v_F\mathrm{}_e`$ the diffusion constant. Comparison with eq. 3 shows that the ratio between the saturation value and the GL-extrapolated value equals $`J_{dp}(0)/J_{dp}^{GL}(0)`$ = 0.385.
In the case of F/S (or N/S) multilayers, the superconducting order parameter is depressed near the interfaces, and this has to be taken into account in calculating $`J_{dp}`$. For this we use the proximity effect model, based on the Usadel equations (dirty limit conditions), which was also used for calculating the depression of $`T_c`$ with decreasing thickness of the superconductor . Details will be given in Appendix A but here we briefly introduce the main parameters of the theory. In principle, the shape of the order parameter on both sides of the interface depends on the bulk transition temperature $`T_{c0}`$, on the coherence lengths $`\xi _{S,F}`$, on the normal state resistivities $`\rho _{S,F}`$, and on the transparency $`T`$ of the interface. From the boundary conditions for the order parameter (see eqs. A4) it follows that, apart from $`T_{c0}`$, only two independent parameters are needed, the proximity strength parameter $`\gamma `$ and the transparency parameter $`\gamma _b`$. The value of $`\gamma `$ = $`(\xi _S\rho _S)/(\xi _F\rho _F)`$ can be fully determined from experiment; the only free parameter is $`\gamma _b`$ ($`0\gamma _b\mathrm{}`$), which is approximately connected to the transparency $`T`$ (with $`0T1`$) by

$$T=\frac{1}{1+\gamma _b}$$
(7)
As was shown in ref , in F/S systems, $`T`$ can be quite low for a high magnetic moment in the F-layer, which is presumably due to the suppression of Andreev reflections by the exchange splitting of the spin subbands. Fig. 1 gives the results of some typical calculations, performed for the system Fe/Nb/Fe with the appropriate proximity effect parameters $`\gamma =`$ 34.6 and $`\gamma _b=`$ 42 (see section IV). Shown is $`J_{dp}(t)`$ for two different thicknesses ($`d_S/\xi _S=20,7.5`$), normalised on the bulk value $`J_{dp}^{bulk}(0)`$ as given by eq. 6. Note that this involves a factor $`(T_c/T_c^{bulk})^{3/2}`$.

The thickness dependence of $`T_c`$ and the normalised depairing current at T = 0 (see the insert of Fig. 1) are quite different, with a much stronger depression of the depairing current at relatively high thickness of the superconductor. This can be qualitatively understood by noting that $`T_c`$ is a measure for the maximum value of the superconducting order parameter in the layer, while the depairing current comes from an average over the layer thickness, which also involves lower values of the order parameter.

## III Experimental

Samples were grown on Si(100) substrates, by DC sputtering in a system with a base pressure of 10-9 mbar in an Ar-pressure of $`6\times 10^3`$ mbar. Sputtering rates were of the order 0.1 nm/s for Nb and 0.03 nm/s for Fe. One series of samples consisted of trilayers Nb/Fe/Nb with Nb thickness $`d_{Nb}`$ = 5 nm and the Fe thickness $`d_{Fe}`$ varying between 2 nm and 25 nm. These were used to determine the magnetization $`M_{Fe}`$ of the Fe layers in the presence of Fe/Nb interfaces with a commercial SQUID-based magnetometer. The behavior of $`M_{Fe}`$ versus $`d_{Fe}`$ could be well described with a straight line, yielding a magnetic moment per Fe atom of 2.36 $`\mu _B`$ ($`\mu _B`$ being the Bohr magneton), slightly above the bulk value of 2.2 $`\mu _B`$ and a magnetically dead layer per interface of 0.1 nm. This value is somewhat lower than reported for MBE-grown samples . Two other series of samples consisted of trilayers of Fe/Nb/Fe with $`d_{Fe}`$ = 5 nm and varying $`d_{Nb}`$. One set was structured by Ar-ion etching into strips with a width $`w`$ = 100 $`\mu m`$, the other into strips with a width $`w`$ = 6 $`\mu m`$, or sometimes 10 $`\mu m`$ or 20 $`\mu m`$. In both cases the length between the voltage contacts was 1 mm. The first set (deposited in two different runs) was used for measuring $`T_c(d_{Nb})`$, the second set for both $`T_c(d_{Nb})`$ and $`J_c(d_{Nb})`$. In all cases, the typical width of the resistive transitions to the superconducting state was 50 mK.

Also measured were single films of Fe and Nb with different strip widths in order to establish values for the specific resistivity $`\rho _{Fe,Nb}`$ (at 10 K), for $`T_c`$ and for the upper critical field $`B_{c2}(T)`$. On average, we find $`\rho _{Fe}`$ 7.5 $`\mu \mathrm{\Omega }cm`$, $`\rho _{Nb}`$ 3.7 $`\mu \mathrm{\Omega }cm`$, $`T_c`$ = 9 K and $`S`$ = $`dB_{c2}/dT`$ = 0.24 T/K, yielding $`\xi _{GL}(0)=`$ $`\sqrt{\mathrm{\Phi }_0/(2\pi ST_c)}`$ = 12.2 nm. This corresponds to $`\xi _S`$ = 7.8 nm. No special precautions were taken to shield residual magnetic fields. The zero-field critical current $`I_c`$ was determined at different temperatures $`T`$ by measuring current ($`I`$)- voltage ($`V`$) characteristics. For this, a DC current was switched on for the time of the order of 1 s and the voltage recorded, to prevent heating via the contacts. All samples showed a clear transition from the superconducting to the normal state, with a large and almost instantaneous increase in voltage at $`I_c`$. Upon detecting this rise, the current was also turned off since the sample then started to heat immediately. Most samples also showed a small rise in voltage prior to the major transition, probably due to vortex motion. We shall come back to this point in the discussion. Important for the theoretically expected value of $`J_{dp}(0)`$ is the value of the resistivity of the superconducting layer (see Eq. 4). This value, $`\rho _{Nb}`$, was extracted from the normal state resistance $`R_n`$ at 10 K of the patterned samples by assuming that the Nb layer and the 10 nm thick Fe layer ($`\rho _{Fe}`$ = 7.5 $`\mu \mathrm{\Omega }cm`$) contribute as parallel resistors.

The resulting values for $`\rho _{Nb}`$ are given in Table 1, together with the strip width $`w`$ and $`T_c`$. The values for the thinner films (around 50 nm) are somewhat larger than what we usually find for single Nb films, and approach that value for the thick films.

## IV Results

Fig. 2 shows the measured values for $`T_c(d_{Nb})`$ for both sample sets, with the the two types of open symbols denoting the two deposition runs for that set, and the solid symbols denoting the samples used for measuring $`J_{c0}`$.

The overall data spread is small, and the data can be well described by the proximity effect theory for S/F systems we used for analyzing the behavior of V/(FexV$`{}_{1x}{}^{})`$ in ref. , with the two parameters $`\gamma `$ and $`\gamma _b`$ defined above. We use the same value for $`\xi _F`$ as in the case of V/Fe, $`\xi _{Fe}`$ = 0.14 nm and values for $`\xi _s`$, $`\rho _F`$ and $`\rho _S`$ as given in section III, yielding $`\gamma `$ = 34.6. The best description for $`T_c(d_{Nb})`$ then is for $`\gamma _b`$ = 42, as shown by the drawn line in Fig. 2. The critical thickness for the S-layer for onset of superconductivity $`d_{cr}^S`$ can be taken either from the lowest measured value for $`T_c`$ or from the extrapolated value of the calculated curve, $`d_{cr}^S`$ = 29 nm, corresponding to a ratio $`d_{cr}^S/\xi _S`$ = 3.7, which is somewhat higher than in the case of V/Fe where we found 3.2. Apparently, the effect of ferromagnet on superconductor is slightly stronger in the Nb/Fe case, but this is not the issue of the current paper.

In Fig. 3 $`J_{c0}=I_{c0}/(wd)`$ is plotted versus reduced temperature $`t=T/T_c`$ for $`d_{Nb}=42`$, 60 and 75 nm. All curves show a clear upturn with decreasing temperature in the region close to $`T_c`$, above $`t0.9`$. Plotting $`J_{c0}^{2/3}(t)`$ versus $`t`$ results in a straight line in this temperature regime, which can be extrapolated to $`t=0`$. The ensuing values for $`J_{c0}^{GL}(0)`$ are given in Table 1 for all samples, and comprise some of the main experimental results. They can also be used to normalize the data. Fig. 4 shows $`[J_{c0}(t)/J_{c0}^{GL}(0)]^{2/3}`$ versus $`t`$ together with the line $`1t`$ (the GL-behavior) and the result of the full theoretical calculation, which is now independent of the parameters. All data collapse on the universal curve above $`t=0.9`$. At lower temperatures, the thinnest films ($`d_s=`$ 36, 40, 42. 53 nm) follow the full calculation quite closely, even down to $`t`$ 0.6. The difference between the data of 36 nm and 40 nm is mainly due to the choice of the normalization value, and reflects the accuracy of that determination. For thicker films the first deviation progressively shifts to higher $`t`$.

## V Discussion

The first point to be discussed is whether the measured values of $`J_{c0}`$ agree with the theoretical estimates for $`J_{dp}`$. The absolute value of $`J_{dp}^{GL}(0)`$ can be calculated with Eq. 5. The materials constants for Nb are well documented ; we use the values $`v_F`$ = $`5.6\times 10^5`$ m/s and $`\rho \mathrm{}`$ = $`3.75\times 10^{16}\mathrm{\Omega }m^2`$. Eq. 5 then yields for the Nb film $`J_{dp,Nb}^{GL}(0)=1.70\times 10^{12}`$ A/m2, which is quite close to the experimentally determined value of $`J_{c0,Nb}^{GL}(0))=1.5\times 10^{12}`$ A/m2 (see Table 1). It is also in good correspondence with the data presented by Ando et al.  on films with a thickness of 100 nm and different strip widths between 0.1 $`\mu m`$ and 10 $`\mu m`$, who found a fitted value $`J_{dp,Nb}^{GL}(0)=1.26\times 10^{12}`$ A/m2. It appears that the depairing current is directly probed by the measurement of $`J_{c0}`$.
Next we consider the dependence of $`J_{c0}^{GL}(0)`$ on the superconducting Nb layer thickness $`d_{Nb}`$. As Eq. 5 shows, $`J_{dp}(0)`$ is proportional to $`\sqrt{1/\rho _{Nb}}`$. Since $`\rho _{Nb}`$ of the samples differs, this leads to some variation in the expected value for $`J_{dp}(0)`$ which can be taken into account by multiplying $`J_{c0}^{GL}(0)`$ by $`\rho _{Nb}^{1/2}`$. Normalizing this value to the single Nb film yields the dependence on $`d_{Nb}`$ as shown in Fig. 5. $`J_{c0}^{GL}(0)`$ in the trilayers is clearly reduced with respect to the bulk Nb value and increases with increasing $`d_{Nb}`$, but much more slowly than $`T_c`$ does. The correspondence with theory is good at low $`d_{Nb}`$, with some deviations above $`d_{Nb}75`$ nm. This coincides with the findings on the temperature dependence of $`J_{c0}(t)`$, shown in Fig. 4 : for small $`d_{Nb}`$ there are only small deviations in the whole measured temperature regime, for large $`d_{Nb}`$ the deviations are large below $`t=0.9`$. This suggests that at high $`d_{Nb}`$ the extrapolation for $`J_{c0}^{GL}(0)`$ leads to somewhat underestimated values. In essence, we conclude that the model used to describe the depression of $`T_c`$ in F/S/S trilayers also adequately describes the behaviour of $`J_{dp}`$.

This leaves one final point to be addressed. In order to determine the depairing current it is usually understood that two conditions have to be fulfilled  : the current has to be distributed uniformly over the strip, and the width $`w`$ should be small enough to preclude vortex formation and motion. In terms of penetration depth $`\lambda (t)`$, strip thickness $`d_s`$ and Ginzburg-Landau coherence length $`\xi (t)`$ this means :

$`w<`$ $`\lambda _{eff}(t)`$ $`=\lambda (t),d_s<\lambda _{eff}(t)`$
$`=\lambda ^2(t)/d_s,d_s>\lambda _{eff}(t)`$
$`w<`$ $`\mathrm{\hspace{0.17em}4.4}\xi (t)`$ $`=4.4\xi (0)/(1t)^{1/2}.`$ (8)
Estimating $`\lambda (0)`$ from $`\lambda (0)=1.05\times 10^3\sqrt{\rho _0/T_c}`$ we find it ranges between 67 and 113 nm. Both conditions mean for all samples $`1t<10^4`$, much smaller than the region where $`J_{c0}(t)(1t)^{3/2}`$, (Fig. 4), and the question is valid whether the current is uniform, as has implicitly been assumed in the analysis.
Qualitatively, current is expected to pile up at the edges of the strip in order to minimize the self-field inside. The edge current will then sooner reach the value of $`I_{dp}`$. By using $`J_{dp}=I_{dp}/(wd)`$, this would lead to underestimating the real value of $`J_{dp}`$. From the close agreement between the experimental and theoretical values this does not appear to be the case. Quantitatively, the situation can be assessed that the self-field of the sample is completely screened ($`B_z=0`$ in the sample). The current distribution is then given by  :

$$J(x)=\frac{I_T}{\pi d\sqrt{W^2x^2}}W<x<W,$$
(9)
where $`I_T`$ is the transport current through the sample, $`x`$ is in the direction of the width $`w`$ of the film, $`x=0`$ in the middle of the film and $`2W=w`$. According to this formula, the current diverges at the edges of the film. It can be assumed, however, that the field penetrates over a distance $`d/2`$ from the edges, but is kept out of the rest of the sample by the screening current. Then, the current within $`d/2`$ from the edges can be set equal to $`J_{dp}^{GL}(0)`$ and beyond $`d/2`$ it decreases according to Eq. 9. The following calculation can be done for the Nb film. The transport current $`I_T`$ in the screened part of the strip can be calculated from Eq. 9 by using

$$J(x=Wd/2)=J_{dp}^{GL}(0).$$
(10)
The total current $`I`$ including the edges is given by :

$`I/d`$ $`=`$ $`2d/2J_{dp}^{GL}(0)+{\displaystyle _{W+d/2}^{Wd/2}}{\displaystyle \frac{I_T}{\pi d\sqrt{W^2x^2}}}`$ (11)
$`=`$ $`dJ_{dp}^{GL}(0)\times `$
$`(1+\sqrt{{\displaystyle \frac{w}{2d}}}\left[\mathrm{sin}^1{\displaystyle \frac{x}{|W|}}\right]_{W+d/2}^{Wd/2})`$
The ratio $`I/(wd)/J_{dp}^{GL}(0)`$, which can be calculated from Eq. 11, gives the fraction of $`J_{dp}^{GL}(0)`$ which would be actually measured as the critical current under the given current distribution, where the depairing current is reached at the edges. It can be easily seen that it equals 1 when the current is uniform. For the Nb film with $`w`$ and $`d`$ as given in Table 1, Eq. 11 yields a fraction of 0.11, an order of magnitude below what is actually measured. The conclusion is that $`J(x)`$ is much more uniformly distributed than might be expected. The reason is probably that a field and moving vortices exist in the film, indicated by a voltage onset below the jump to the normal state. This breaks up the Meissner state, causes a much more uniform current distribution and allows the correct determination of the depairing current over a much larger region than expected on the basis of the condition $`w<\lambda ,\xi `$. Still, the deviations of $`J_{c0}(t)`$ compared to the theoretical behaviour at higher $`d_{Nb}`$ in the Fe/Nb/Fe trilayers may be due to the low value of $`\lambda `$ and inhomogeneities in the current distribution at these high thicknesses. At low thicknesses, there are two effects which increase $`\lambda `$ above the bulk value. Firstly, for $`d_s<`$ $`\lambda `$ (around $`d_{Nb}75`$ nm) the effective penetration depth increases according to $`\lambda _{eff}=\lambda (0)^2/d`$, and can become significantly higher than $`\lambda (0)`$. Secondly, the suppression of the order parameter as measured by the decrease of $`T_c/T_{c0}`$ results in a higher value for $`\lambda (0)`$. From that point of view the full agreement between the measured and calculated values of $`J_{c0}(t)`$ at the lowest thicknesses is not surprising.

## VI Summary

In this paper we have addressed the question of the value of the superconducting depairing current in F/S/F trilayers with varying $`d_s`$, where the superconducting order parameter is inhomogeneously suppressed by the pair breaking in the F-layers. The same model which is adequately describes the suppression of $`T_c`$ with decreasing $`d_s`$ with two parameters (proximity strength $`\gamma `$ and interface transparency $`\gamma _b`$ or $`T`$) can also be used to compute the suppression of the depairing current. Measurements of the zero-field critical current $`J_{c0}`$ (as defined by the current where the resistance jumps to the normal state value) in thin strips of Fe/Nb/Fe show that the temperature dependence near $`T_c`$ is as expected for the depairing current. Also the absolute value of $`J_{c0}`$ of single Nb films is close to the theoretically expected value and the measured suppression of $`J_{c0}`$ in the trilayers follows the calculated behaviour. We conclude that the current distribution is homogeneous and that the depairing current is measured, even though the strip widths are larger that the superconducting penetration depth and coherence length. Also, the proximity effect model correctly describes the shape of the order parameter, at least in the superconducting layer. These findings can be of use in experiments on the effect of injecting polarized quasiparticles.

## Acknowledgements

This work is part of the research program of the ’Stichting voor Fundamenteel Onderzoek der Materie (FOM)’, which is financially supported by NWO. A. A. G. acknowledges support by NWO in the framework of the Russian-Dutch collaboration programme, grant 047-005-01. We would like to thank P. H. Kes for helpful discussions.

## Appendix A : calculation of $`J_{dp}`$

We assume that the dirty limit conditions are fulfilled in both S and F layers, so that the F/S bilayer can be described by the Usadel equations. In the absence of a depairing current in the S layer, and in the regime of large exchange energy in the ferromagnet ($`E_{ex}k_BT_c`$) these equations were discussed extensively by Buzdin et al. (see also Demler et al. ). Here we rewrite these equation in $`\theta `$-parametrization ($`F=\mathrm{sin}\theta `$, $`G=\mathrm{cos}\theta `$) and include the pair-breaking effects by current along the S film :

$$\xi _S^2\frac{d^2}{dz^2}\theta _S(z)\stackrel{~}{\omega }\mathrm{sin}\theta _S(z)+\mathrm{\Delta }(z)\mathrm{cos}\theta _S(z)=0,$$
(A1)
$$\xi _F^2\frac{d^2}{dz^2}\theta _F(z)i\mathrm{sin}\theta _F(z)=0,$$
(A2)
$$\mathrm{\Delta }\mathrm{ln}(T/T_c)+\pi T\underset{\omega _n}{}(\frac{\mathrm{\Delta }}{\left|\omega _n\right|}\mathrm{sin}\theta _S)=0,$$
(A3)
where $`\omega _n=\pi (2n+1)T/T_c`$ is the normalized Matsubara frequency, $`\stackrel{~}{\omega }=\left|\omega _n\right|+Q^2\mathrm{cos}\theta (z)`$, $`\mathrm{\Delta }`$ is the pair potential in a superconductor normalized to $`\pi T_c`$, $`\xi _S=(\mathrm{}D_S/2\pi T_c)^{1/2}`$, $`\xi _F=(\mathrm{}D_F/2E_{ex})^{1/2}`$ and $`D_{F,S}`$ are the coherence lengths and the electronic diffusion coefficients in F and S metals. Moreover, $`Q=\xi _s\chi /x`$ is the normalized gradient-invariant superfluid velocity in the $`x`$-direction, with $`\chi `$ the phase of the pair potential $`\mathrm{\Delta }`$. There are two sources of pair-breaking in the problem, the volume one by the current and the surface one by the ferromagnet. The latter is described by the boundary conditions at the FS interface ($`z=0`$)

$`\xi _S{\displaystyle \frac{d}{dz}}\theta _S`$ $`=`$ $`\gamma \xi _F{\displaystyle \frac{d}{dz}}\theta _F`$ (A4)
$`\gamma _b\xi _F{\displaystyle \frac{d}{dz}}\theta _F`$ $`=`$ $`\mathrm{sin}(\theta _S\theta _F)`$ (A5)
where the parameter $`\gamma =\rho _S\xi _S/\rho _F\xi _F`$ describes the strength of the suppression of superconductivity in S by the ferromagnet.

The parameter $`\gamma _b`$ describes the effect of boundary transparency (coupling strength) between the layers. In the NS case, when the decoupling is due to the presence of an additional potential barrier at the interface, $`\gamma _b=R_B/\rho _F\xi _F`$, with $`R_B`$ the normal state resistance of the N/S interface . In the F/S bilayer there is no general microscopic derivation for $`\gamma _b`$, combining the effect of exchange splitting and an additional interface barrier. A simple estimate is still possible, when the exchange splitting is the main cause for intransparency . Then $`\gamma _b=(2/3)(l_F/\xi _F)<(1T_A)/T_A>`$ where $`T_A`$ is the transmission probability of scattering between the majority and minority spin subbands, i.e. the probability of Andreev reflection. This process is implicitly described by the boundary condition $`\gamma _b\xi _F\frac{d}{dz}\theta _F=\mathrm{sin}(\theta _S\theta _F)`$ since $`\theta _F`$ is off-diagonal in spin indices. Here the brackets $`<\mathrm{}>`$ denote the Fermi surface averaging, which is generally proportional to the overlap area of the projections of different spin subbands onto the contact plane . As a result, $`T_A`$ drops roughly linearly (for spherical Fermi surfaces) as a function of $`E_{ex}`$, both for ballistic and diffusive interfaces . The supercurrent density is given by

$$J_s(z,Q)=\frac{2\pi \sigma _s}{e}QT\underset{\omega _n}{}\mathrm{sin}^2\theta _s.$$
(A6)
Since the superconducting pair potential $`\mathrm{\Delta }`$ and the Green’s function $`\theta _s`$ are suppressed by the superflow $`Q`$, the dependence $`J_s(Q)`$ must be found selfconsistently. In the well known spatially homogeneous case  the function $`J_s(Q)`$ behaves non-monotonously : the supercurrent $`J_s`$ increases with $`Q`$ at small $`Q`$, then reaches a maximum and finally drops to zero, when $`\mathrm{\Delta }`$ is fully suppressed by current. The depairing current is defined as the maximum of $`J_s(Q)`$. A similar situation holds in the spatially inhomogeneous case considered here, with the difference that the solutions for $`\theta (z)`$ and $`\mathrm{\Delta }(z)`$ of the proximity effect problem (eqs. A1-A3) should be calculated selfconsistently for a given $`Q`$ using the boundary conditions at the FS interface (eqs. A4, A5). This problem is solved numerically by the method applied previously to NS bilayers and described in detail in . Then the local $`z`$-dependent supercurrent density $`J_s(z,Q)`$ is calculated from eq.A6 by summing the solutions $`\mathrm{sin}^2\theta _s`$ over $`\omega _n`$. Finally the density is averaged over film thickness $`J_s(Q)=`$ $`d_s^1_0^{d_s}J_s(z,Q)𝑑z`$ and the depairing current is found from the maximum of the dependence of $`J_s`$ on $`Q`$.


format: text

text:
# Magnetic Field Induced Gap and Kink Behavior of Thermal Conductivity

It is now well known from the study of field theoretical models that an external magnetic field can be a strong catalyst for a symmetry breaking leading to the generation of a fermion dynamical mass even at the weakest attractive interactions . This phenomenon of magnetic catalysis (MC) of symmetry breaking has been shown to be essentially model independent and to have wide applications in several physical systems , in particular, in condensed matter physics . The essence of the MC effect lies in the dimensional reduction of the fermion pairing dynamics when the pairing energy is much less than the Landau gap $`\sqrt{eB}`$ ($`B`$ is the magnitude of the magnetic field induction). In this case, the fermions are mostly confined to the lowest Landau level (LL) and their attractive interaction, regardless how small it may be, leads to mass generation.

In recent years, several authors  suggested that the magnetic catalysis may play a relevant role in the physics of high-$`T_c`$ superconductors in external magnetic fields. This idea was inspired in part by the outcome of various experiments  initiated by Krishana’s et al. on the thermal conductivity of high-$`T_c`$ compounds in the presence of an external magnetic field perpendicularly applied to the cuprate plane. According to these experiments, at temperatures significantly lower than the critical temperature of superconductivity, the profile of the thermal conductivity $`\kappa `$ versus the magnetic field displays a sharp break in its slope at a transition field $`B_c`$, followed then by a plateau region in which it ceases to change with increasing field up to the highest attainable fields $`14T`$.

To understand how the magnetic catalysis was related  to these experimental observations, let us recall that high-$`T_c`$ superconducting cuprates are d-wave superconductors which implies that their Fermi surface is characterized by nodal points where the gap function vanishes. Then, the low-energy quasiparticle (QP) dynamics is concentrated mainly around these nodal points because these states can be occupied even at very low temperatures. The QP excitations in the vicinity of each node can be described by a massless Dirac Lagrangian . Hence, a natural candidate to modelling the QP interactions in a d-wave superconductor is a $`3`$-dimensional relativistic four-fermion interaction (Nambu-Jona-Lasinio (NJL) type) model. It is well known  that the gap equations of such a kind of models in the presence of an external magnetic field lead to the generation of a dynamical fermion mass (QP gap), which scales (at zero temperature) with the field as $`m\sqrt{eB}`$. The critical temperature at which the dynamical mass vanishes is determined by the dynamical mass at zero temperature, and therefore, it scales with the magnetic field with the same law as the zero-temperature mass. One of the outcomes of the above mentioned experiments  was that the critical temperature for the appearance of the kink-like behavior of the thermal conductivity scales with the magnetic field as $`T_c\sqrt{B}`$. Then, the square root field dependence of the critical temperature in $`(2+1)`$dimensional NJL models was considered as a hint that the MC phenomenon could be important to explain the experimental results of Krishana et al. .

By now, several possible mechanisms for the arising plateau have been proposed . All of them are based on the adoption of the QP picture as the low-energy excitations at the nodes of the $`d`$wave symmetric order parameter. In the approach of Refs., following the ideas of Ref. , it was assumed that the QP spectrum in the magnetic field was characterized by Landau levels. Furthermore, in Ref.  it was supposed that a QP gap opens at the critical field $`B_c`$ leading to the exponential vanishing of the quasiparticle contribution in the thermal conductivity. Since the total conductivity is assumed to be the sum of the QP term $`\kappa _{\mathrm{QP}}`$ and the phonon term $`\kappa _{\mathrm{ph}}`$ with the phononic part independent of the field, this would lead to a plateau for high fields. Laughlin  relates the QP gap to a weakly first order phase transition leading to the development of an additional small parity-violating superconducting order parameter. In another line of reasoning , the drop in the QP’s conductivity, hence the plateau in the total conductivity, is ascribed to the opening of a field-induced energy gap in $`d`$-wave superconductors due to a second order phase transition.

On the contrary, the approach of Franz  does not appeal to the generation of any QP gap. Albeit the phononic part of the thermal conductivity is also considered to be field-independent, in Franz’ scenario the QPs conductivity itself becomes field-independent above a crossover field $`B_c`$. The effect is primarily due to the compensation of the enhancement of the QP density of states, arising from the Doppler shift of QP’s spectrum in the superfluid velocity field around vortices (Volovik effect ), and the growth of QP’s width caused by the scattering on vortices. The Volovik effect certainly plays a decisive role at weak magnetic fields and low temperatures (less than $`1`$ K), where the increase in the thermal conductivity was predicted theoretically , and later confirmed in the experiment , to follow the $`\sqrt{B}`$ dependence of the density of states. However, Volovik’s arguments are based on semiclassical calculations valid in case one has well isolated vortices with supercurrents circulating around them. At higher fields, when vortices begin to overlap, this semiclassical picture should be replaced by a more adequate quantum mechanical treatment.

It should be stressed that none of the above works address the question of the kink itself in the thermal conductivity. Moreover, we should mention that the applicability of the Landau level quantization scheme used in the QP picture of the d-wave superconductor in a constant magnetic field has been recently subjected to intense criticism by several authors . The Landau level-like spectrum proposed by Gorkov and Schiffer , and Anderson  was based on the assumption that the spatially dependent superfluid velocity can be neglected, but according to Melnikov , the superfluidity velocity strongly mixes individual Landau levels. In order to take into account the effects of spatially varying supercurrents, Franz and Tesanovich  introduced a singular gauge transformation to trade the phase of the order parameter for a ficticious $`U(1)`$ gauge field. As a result, the two fluxes, due to the external magnetic field and the $`U(1)`$ field, cancel each other on average, washing out the whole picture of Landau level quantization. In this framework the low-energy quasiparticles are in fact Bloch waves described by massless Dirac fermions moving in a vector potential of physical supercurrents but zero average magnetic field. While the arising scenario might seem promising, the approach of Ref.  itself is not free of problems: there is no good foundation from a physical point of view to use singular gauge transformations and to introduce fictitious gauge fields, both of which are needed in this scheme. Moreover, because of the singular transformations, the physical picture itself becomes strongly dependent on the way such transformations are implemented (i.e., different singular gauge transformations can give rise to different spectra, depending on the regularization scheme ). Finally, as far as we know, it has not been proved up to now that such a formulation can be applicable in the regime of weak-to-moderate magnetic fields, where the individual vortices overlap and the magnetic field can be considered to be fairly uniform.

The aim of the present paper is to discuss a mechanism that generates the kink effect within the framework of the MC phenomenon. We emphasize that although, as shown in the derivations that follow, our numerical curves are in close qualitative agreement with the experimental curves of Ref. , we do not pretend to claim that our results explain the experimental findings of Krishana et. al. , as our calculations are based on the Landau level quantization whose applicability to d-wave superconductivity, as mentioned above, is not well established. Our main result in this paper is then of a more general and theoretical character. That is, we show that the MC phenomenon can be responsible for a kink effect in the thermal conductivity of (2+1)-dimensional relativistic fermion systems, and we point out the very important finding that the kink-like effect close to the critical field is essentially model independent, since it is determined by the critical behavior of the dynamically generated mass near the phase transition point. This fact makes our result relevant beyond the particular model under consideration, linking it to the universality class of theories with such a critical behavior.

In connection with this, we would like to call the attention of the reader to a new class of recently discovered materials , which could be a playground for applications of the phenomenon of magnetic catalysis. The so-called highly oriented pyrolitic graphites (HOPG) have layered structure with two isolated points in the Brillouin zone where the dispersion is linear, so that the electronic states can be described here too in terms of relativistic Dirac particles . It was already suggested  that the magnetic catalysis in the system of Coulomb-interacting Dirac fermions can provide an explanation for the semimetal-insulator phase transition observed in HOPG in the presence of a magnetic field perpendicular to the layers. Since the quasiparticles in a graphite do not suffer from the Bloch waves versus Landau levels dilemma, our present results should have full strength there and we anticipate that the thermal conductivity of these systems should exhibit a behavior similar to the one we are reporting here.

In this paper we present for the first time a consistent calculation of the thermal conductivity in the presence of an external magnetic field in a model with the simplest four-fermion interaction for quasiparticles. We use the same constant magnetic field approximation that was already explored in Ref. to calculate the thermal conductivity. However, our calculation deviates considerably from what was done in . Not only we take into account the contribution of all Landau levels, but the definition of the heat current itself is different. In the limit of narrow width of quasiparticles ($`\mathrm{\Gamma }T,\sqrt{B}`$), after a gap is opened, the thermal conductivity exhibits a new term proportional to $`\sigma ^2`$ ($`\sigma `$ is the gap). This term originates from the compensation of the matrix elements of transitions between the zeroth and the first LLs (behaving as $`1/B`$) and the LL density of states which in turn is proportional to $`B`$. This is one more manifestation of the MC phenomenon: not only a gap is induced by the magnetic field, but the transitions between the zeroth and the first LLs are enhanced. In mean-field approximation and near the phase transition point the gap behaves like $`\sigma 0.523\sqrt{eBeB_c}`$ what yields a positive contribution into the slope of the thermal conductivity leading to a jump (kink effect) of $`\kappa (B)`$ at $`B=B_c`$.

We start from the following NJL model in an external magnetic field in $`(2+1)`$ dimensions,

$$=\overline{\psi }_a[i\gamma ^0\mathrm{}\frac{}{t}+v_D\gamma ^i(i\mathrm{}\frac{}{x^i}\frac{e}{c}A_i)]\psi _a+\frac{g}{2N}(\overline{\psi }_a\psi _a)^2,$$
(1)
In (1) the vector potential is taken in the symmetric gauge $`A_\mu =(0,\frac{B}{2}x_2,\frac{B}{2}x_1)`$, and $`v_D`$ is a characteristic velocity . In what follows we take $`\mathrm{}=v_D=k_B=1`$ and absorb $`c`$ into the charge $`e`$. We will restore these constants when needed. We assume also that the fermions carry an additional flavor index $`a=1,\mathrm{},N`$ ($`N=2`$ for realistic $`d`$-wave superconductors). The Dirac $`\gamma `$ matrices are taken in a reducible four-component representation.

The Lagrangian density (1) is invariant under the discrete (chiral) symmetry $`\psi \gamma _5\psi ,\overline{\psi }\overline{\psi }\gamma _5,`$ which forbids the fermion mass generation in perturbation theory. The mass generation can be studied introducing an auxiliary field $`\sigma =(g/N)\overline{\psi }_i\psi _i`$ by means of the Hubbard-Stratanovich trick which permits one to integrate over fermion fields in the path integral representation of the partition function. The field $`\sigma `$ has no dynamics at the tree level but it acquires the kinetic term due to fermion loops. Likewise, a nontrivial minimum of the effective potential (the expectation value of $`\sigma `$) gives mass to fermions and spontaneously breaks the discrete symmetry leading to a neutral condensate of fermion-antifermion pairs. Studying the minimum of the effective potential we find that, at a fixed temperature $`T`$, there is a critical value of the magnetic field $`\sqrt{eB_c}/T4.148`$ such that for subcritical fields $`eBeB_c`$ the gap is zero, while for $`eB>eB_c`$ a nontrivial gap appears. The critical curve equation has the form $`(v_D/c)^210^{10}B=21.5T^2`$ for $`B`$ measured in Tesla. Using the approximated value of the characteristic velocity $`v_D10^7cm/s`$ , we obtain the critical curve $`B=0.014T^2`$. Notice that the obtained critical curve fits well the experimental curve of Ref..

To derive an expression for the static thermal conductivity in an isotropic system, we follow the familiar linear response method and apply Kubo’s formula

$$\kappa =\frac{1}{TV}\mathrm{Im}\underset{0}{\overset{\mathrm{}}{}}𝑑ttd^2x_1d^2x_2u_i(x_1,0)u_i(x_2,t),$$
(2)
where $`V`$ is the volume of the system and $`u_i(x,t)`$ is the heat-current density operator. The brackets denote averaging in the canonical ensemble. Physically, the thermal conductivity $`\kappa `$ appears as a coefficient in the equation relating the heat current to the temperature gradient $`\stackrel{}{u}=\kappa \stackrel{}{}T`$ under the condition of absence of particle flow. If we neglect the chemical potential, the heat density coincides with the energy density $`ϵ`$, hence the quantity $`\stackrel{}{u}`$ that satisfies the continuity equation $`\dot{ϵ}(x)+\stackrel{}{}.\stackrel{}{u}(x)=0`$ can be interpreted as the heat current. From the Lagrangian density (1) we find $`u_i=\frac{i}{2}\left(\overline{\psi }\gamma _i_0\psi _0\overline{\psi }\gamma _i\psi \right)`$.

Neglecting vertex corrections  the calculation of the thermal conductivity reduces to the evaluation of the bubble diagram . In this approximation, making use of the spectral representation for the fermion Green’s function, we arrive at the following expression for the thermal conductivity in the Matsubara formalism

$$\kappa =\frac{1}{32\pi T^2}\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{d\omega \omega ^2}{\mathrm{cosh}^2\frac{\omega }{2T}}d^2k\mathrm{tr}\left[\gamma ^iA(\omega ,\stackrel{}{k})\gamma ^iA(\omega ,\stackrel{}{k})\right].$$
(3)
The spectral function $`A(\omega ,\stackrel{}{k})=\pi ^1\mathrm{Im}S^R(\omega +i\epsilon ,\stackrel{}{k})`$ is derived from the fermion propagator in a magnetic field decomposed over Landau level poles ,

$`A(\omega ,\stackrel{}{k})=e^{\frac{\stackrel{}{k}^2}{eB}}{\displaystyle \frac{\mathrm{\Gamma }}{2\pi }}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{M_n}}[{\displaystyle \frac{(\gamma ^0M_n+\sigma )f_1(\stackrel{}{k})+f_2(\stackrel{}{k})}{(\omega M_n)^2+\mathrm{\Gamma }^2}}`$ (4)
$`+{\displaystyle \frac{(\gamma ^0M_n\sigma )f_1(\stackrel{}{k})f_2(\stackrel{}{k})}{(\omega +M_n)^2+\mathrm{\Gamma }^2}}],M_n=\sqrt{\sigma ^2+2eBn},`$ (5)
$`f_1(\stackrel{}{k})=2\left[P_{}L_n(s)P_+L_{n1}(s)\right],f_2(\stackrel{}{k})=4(\stackrel{}{k}\stackrel{}{\gamma })L_{n1}^1(s).`$ (6)
Here $`P_\pm =(1\pm i\gamma ^1\gamma ^2)/2`$ are projector operators, $`L_n,L_n^1`$ are Laguerre’s polynomials ($`L_1^1=0`$), $`s=2\stackrel{}{k}^2/eB,`$ and $`\sigma `$ is the dynamical fermion mass obtained from the finite temperature gap equation in the magnetic field. We introduced also the quasiparticles width replacing $`\epsilon `$ in the definition of $`A(\omega ,\stackrel{}{k})`$ by finite $`\mathrm{\Gamma }`$, which is due to interaction processes, in particular, scatterings on impurities. Straightforward calculation of the trace and further integration over momenta in Eq. (3) produce Kronecker’s delta symbols $`\delta _{n,m1}+\delta _{m,n1}`$ due to orthogonality of Laguerre’s polynomials. This implies that only those transitions between neighbor Landau levels contribute to the heat transfer. After performing one of the sums in (3) we obtain

$`\kappa `$ $`=`$ $`{\displaystyle \frac{eB\mathrm{\Gamma }^2N}{\pi ^2T^2}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\omega \omega ^2}{\mathrm{cosh}^2\frac{\omega }{2T}}}{\displaystyle \frac{1}{(\omega ^2+M_n^2+\mathrm{\Gamma }^2)^24\omega ^2M_n^2}}`$ (7)
$`\times `$ $`{\displaystyle \frac{(\omega ^2+M_n^2+\mathrm{\Gamma }^2)(\omega ^2+M_{n+1}^2+\mathrm{\Gamma }^2)4\omega ^2\sigma ^2}{[(\omega ^2+M_{n+1}^2+\mathrm{\Gamma }^2)^24\omega ^2M_{n+1}^2]}}.`$ (8)
Further summation over $`n`$ in Eq. (8) can be carried out expanding the integrand in terms of partial fractions. The resulting sums are expressed through digamma functions and the final expression for $`\kappa `$ is written as follows

$`\kappa ={\displaystyle \frac{N\mathrm{\Gamma }^2}{2\pi ^2T^2}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\omega \omega ^2}{\mathrm{cosh}^2\frac{\omega }{2T}}}{\displaystyle \frac{1}{(eB)^2+(2\omega \mathrm{\Gamma })^2}}\{2\omega ^2+`$ (9)
$`{\displaystyle \frac{(\omega ^2+\sigma ^2+\mathrm{\Gamma }^2)(eB)^22\omega ^2(\omega ^2\sigma ^2+\mathrm{\Gamma }^2)eB}{(\omega ^2\sigma ^2\mathrm{\Gamma }^2)^2+4\omega ^2\mathrm{\Gamma }^2}}`$ (10)
$`{\displaystyle \frac{\omega (\omega ^2\sigma ^2+\mathrm{\Gamma }^2)}{\mathrm{\Gamma }}}\mathrm{Im}\psi \left({\displaystyle \frac{\sigma ^2+\mathrm{\Gamma }^2\omega ^22i\omega \mathrm{\Gamma }}{2eB}}\right)\}.`$ (11)
This formula is the main result of our paper and we can use it now to study the different asymptotic regimes. It is easy to show that in the limit of zero field $`\kappa `$ reduces to

$$\kappa _0=\frac{N}{4\pi ^2T^2}\underset{0}{\overset{\mathrm{}}{}}\frac{d\omega \omega ^2}{\mathrm{cosh}^2\frac{\omega }{2T}}\left[1+\frac{\omega ^2+\mathrm{\Gamma }^2}{\omega \mathrm{\Gamma }}\mathrm{arctan}\frac{\omega }{\mathrm{\Gamma }}\right],$$
(12)
where we put $`\sigma =0,`$ since the mass is not generated in the zero-field weak coupling NJL model. Eq. (12), up to an overall factor $`k_B^2(v_F^2+v_\mathrm{\Delta }^2)/\mathrm{}v_Fv_\mathrm{\Delta }`$, coincides with the corresponding expression obtained in Ref.. With the overall factor replacing $`N`$ ($`=2`$ in real $`d`$wave superconductor), Eq. (12) reproduces the universal (or residual) thermal conductivity at low $`T`$ in the so-called “dirty” limit, $`T\mathrm{\Gamma }`$. The residual conductivity was recently observed in experiments  giving explicit confirmation of the existence of gapless quasiparticles in $`d`$-wave cuprates at $`T<T_c`$.

Let us consider the case $`\mathrm{\Gamma }T`$ with $`B0`$. In this approximation we obtain

$$\kappa \frac{N\mathrm{\Gamma }}{4\pi T^2}\left\{\frac{\sigma ^2}{\mathrm{cosh}^2\frac{\sigma }{2T}}+4\underset{n=1}{\overset{\mathrm{}}{}}\frac{n(\sigma ^2+2eBn)}{\mathrm{cosh}^2\frac{\sqrt{\sigma ^2+2eBn}}{2T}}\right\}.$$
(13)
Asymptotically, at $`\sqrt{eB}\sqrt{eB_c}T`$, the dynamical mass behaves as $`\sigma \sqrt{eB}`$ what leads to an exponential decrease as in the case of gapless fermions. However, the term $`\mathrm{cosh}^2\sigma /2T`$ in Eq. (13) is of order one when $`\sigma 2T`$, thus, there is no suppression of this term for a certain range of fields where it is almost field independent (plateau region).

Near the phase transition point $`\sigma a\sqrt{eBeB_c}`$ and the first term in Eq. (13) gives positive contribution to the slope of the thermal conductivity at $`eBeB_c`$ leading to the jump in the slope of $`\kappa `$ (kink effect) at $`eB=eB_c`$. The parameter $`a`$ is model-dependent, but the scaling of $`\sigma `$ is determined by the universality class of the theory. For the NJL model (1) we find, in mean-field approximation, $`a0.523`$.

The explicit appearance of the kink has been corroborated by numerical calculations as shown in Fig.1. Notice the break in the slope of $`\kappa `$ (kink effect) at the critical value $`B_c`$ in the presence of $`\sigma `$. For $`B>B_c`$ the kink is followed by a region where $`\kappa `$ is only weakly dependent on the field. While decreasing the temperature, the position of the kink moves to the left in accordance with the critical line $`B_c=0.014T^2`$. Notice the similarity of our results (Fig. 1) with the experimental behavior reported in  for high-Tc superconductors. Moreover, Fig. 1 also reproduces the experimentally observed crossing of the curves that occurs with decreasing $`T`$ in such a way that the lower $`T`$ curve reaches the higher value at large fields.

Because of the important role played by the first term in (13) to get the kink effect, it is instructive to clarify its origin. From Eq. (8) one can see that while the contribution of transitions between Landau levels with $`n1`$ in the integrand behaves as $`1/(eB)^2`$ for large fields, the contribution of transitions between the zeroth and the first LLs decreases only as the first power of the field ($`1/eB`$). Since the density of LLs is proportional to $`eB`$, this implies that the transitions between the zeroth and the first LLs are not suppressed, even though the gap between the levels grows with the field. Note that for gapless QPs the transitions with $`n=0`$ are completely suppressed in the regime $`\mathrm{\Gamma }0`$: indeed, their contribution is proportional to $`\delta (\omega )`$ and the integrand in (8) contains an $`\omega ^2`$ factor. For gapped QPs, $`\delta (\omega )`$ is replaced by $`\delta (\omega \pm \sigma )`$ and the transitions with $`n=0`$ survive.

In summary, in the present paper the thermal conductivity at $`B0`$ in the framework of a relativistic four-fermion model was calculated. Assuming a uniform magnetic field approximation, we showed that the thermal conductivity exhibits a kink behavior when the field reaches the critical value $`B_c=0.014T^2`$. Two main features determine this effect: the generation of a QP gap in a magnetic field (MC phenomenon) and the lack of suppression of zeroth-first LLs transitions. We stress that the appearance of a kink in the thermal conductivity is model independent, being only determined by the critical behavior of the induced dynamical mass near the phase transition. The universality character of this result leaves an open window for possible future applications.

We would like to acknowledge V.A. Miransky for stimulating discussions and careful reading of the manuscript, and A. Hams and I. Shovkovy for help in numerical calculation. This research has been supported in part by NSF under Grants PHY-0070986 (E.J.F. and V.I.), PHY-9722059 (E.J.F., V.P.G. and V.I.) and POWRE-PHY-9973708 (V.I.). The work of V.P.G. is supported also by the SCOPES-project 7UKPJ062150.00/1 of the Swiss NSF and Grant-in-Aid of Japan Society for the Promotion of Science (JSPS) #11695030. V.P.G. wishes to acknowledge JSPS for financial support.


text:
# Small, Dense Quark Stars from Perturbative QCD

## Abstract

As a model for nonideal behavior in the equation of state of QCD at high density, we consider cold quark matter in perturbation theory. To second order in the strong coupling constant, $`\alpha _s`$, the results depend sensitively on the choice of the renormalization mass scale. Certain choices of this scale correspond to a strongly first order chiral transition, and generate quark stars with maximum masses and radii approximately half that of ordinary neutron stars. At the center of these stars, quarks are essentially massless.

Strongly interacting matter under extreme conditions can reveal new phenomena in Quantum Chromodynamics (QCD). Compact stars serve as an excellent observatory to probe QCD at large density, as their interior might be dense enough to allow for the presence of chirally symmetric quark matter, i.e., quark stars .

The usual model used for quark stars is a bag model, with at most a correction $`\alpha _s`$ from perturbative QCD . In the massless case, the first order correction cancels out in the equation of state, so that one ends up finally with a free gas of quarks modified only by a bag constant. If the bag constant is fit from hadronic phenomenology, then the gross features of quark stars are very similar to those expected for neutron stars: the maximum mass is $`2.M_{}`$, with a radius $`10`$ km.

In this article we consider quark stars, using the equation of state for cold, dense QCD in perturbation theory to $`\alpha _s^2`$ . These results are well known, and our only contribution is to use modern determinations of the running of the QCD coupling constant . At the outset, we stress that we do not suggest that the perturbative equation of state is a good approximation for the densities of interest in quark stars. Rather, we use it merely as a model for the equation of state of QCD.

To $`\alpha _s^2`$, there is significant sensitivity to the choice of the renormalization mass scale. Under our assumptions, we find that this choice is tightly constrained by the physics. We consider two illustrative values of this parameter. One choice corresponds to a weakly first order chiral transition (or no true phase transition), and gives maximum masses and radii very similar to that of neutron stars. The second choice corresponds to a strongly first order chiral transition , and generates two types of stars. One type has densities a few times that of nuclear matter, and looks like the stars of a weakly first order chiral transition. In addition, however, there is a new class of star , with densities much higher than that of nuclear matter . For this new class, the maximum mass is $`1.M_{}`$, with a radius $`5`$ km. Other models with nonideal behavior also generate small, dense quark stars .

Assume that the chiral phase transition occurs at a chemical potential $`\mu _\chi `$ . Our perturbative equation of state is applicable only in the chirally symmetric phase, when the quark chemical potential $`\mu >\mu _\chi `$. In this phase, the effects of a strange quark mass, $`m_s100`$ MeV , are small relative to the quark chemical potentials, $`\mu >300`$ MeV. Thus we take three flavors of massless quarks with equal chemical potentials .

The thermodynamic potential of a plasma of massless quarks and gluons was calculated perturbatively to $`\alpha _s^2`$ by Freedman and McLerran  and by Baluni , using the momentum-space subtraction (MOM) scheme. The MOM coupling is related to that in the modified mimimal substraction scheme, $`\overline{\mathrm{MS}}`$, as :

$$\frac{\alpha _s^{\mathrm{MOM}}}{\pi }=\frac{\alpha _s^{\overline{\mathrm{MS}}}}{\pi }\left[1+𝒜\frac{\alpha _s^{\overline{\mathrm{MS}}}}{\pi }\right];$$
(1)
$`\alpha _s=g^2/(4\pi )`$, with $`g`$ the QCD coupling constant, and $`𝒜=151/48(5/18)N_f`$, with $`N_f`$ the number of massless flavors. (While we take $`N_f=3`$, we give formulas for arbitary $`N_f`$.) In the $`\overline{\mathrm{MS}}`$ scheme, to $`\alpha _s^2`$ the thermodynamic potential is then

$`\mathrm{\Omega }(\mu )={\displaystyle \frac{N_f\mu ^4}{4\pi ^2}}\{12\left({\displaystyle \frac{\alpha _s}{\pi }}\right)`$ (2)
$``$ $`[G+N_f\mathrm{ln}{\displaystyle \frac{\alpha _s}{\pi }}+(11{\displaystyle \frac{2}{3}}N_f)\mathrm{ln}{\displaystyle \frac{\overline{\mathrm{\Lambda }}}{\mu }}]\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^2\},`$ (3)
where $`G=G_00.536N_f+N_f\mathrm{ln}N_f`$, $`G_0=10.374\pm .13`$ , and $`\overline{\mathrm{\Lambda }}`$ is the renormalization subtraction point. In $`\overline{\mathrm{MS}}`$ scheme, the thermodynamic potential is manifestly gauge invariant. We take the scale dependence of the strong coupling constant, $`\alpha _s\alpha _s(\overline{\mathrm{\Lambda }})`$ as :

$`\alpha _s(\overline{\mathrm{\Lambda }})`$ $`=`$ $`{\displaystyle \frac{4\pi }{\beta _0u}}[1{\displaystyle \frac{2\beta _1}{\beta _0^2}}{\displaystyle \frac{\mathrm{ln}(u)}{u}}`$ (4)
$`+`$ $`{\displaystyle \frac{4\beta _1^2}{\beta _0^4u^2}}((\mathrm{ln}(u){\displaystyle \frac{1}{2}})^2+{\displaystyle \frac{\beta _2\beta _0}{8\beta _1^2}}{\displaystyle \frac{5}{4}})];`$ (5)
$`u=\mathrm{ln}(\overline{\mathrm{\Lambda }}^2/\mathrm{\Lambda }_{\overline{\mathrm{MS}}}^2)`$, $`\beta _0=112N_f/3`$, $`\beta _1=5119N_f/3`$, and $`\beta _2=28575033N_f/9+325N_f^2/27`$. The scale $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ is fixed by requiring that $`\alpha _s=0.3089`$ at $`\overline{\mathrm{\Lambda }}=2`$ GeV ; for $`N_f=3`$, (4) gives $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}=365`$ MeV.

All thermodynamic quantities follow consistently from $`\mathrm{\Omega }(\mu )`$. The pressure is given by $`p(\mu )=\mathrm{\Omega }(\mu )`$, the quark number density by $`n(\mu )=(p/\mu )`$, and the energy density by $`ϵ=p+\mu n`$. Given our stated assumptions, the only freedom we have in the model is the choice of the ratio $`\overline{\mathrm{\Lambda }}/\mu `$. To illustrate this, we take the values $`\overline{\mathrm{\Lambda }}/\mu =1,2,3`$.

For reasons which will become clear later, we find the choice $`\overline{\mathrm{\Lambda }}=2\mu `$ especially interesting. Start with a very large chemical potential, such as $`\mu =100`$ GeV, for which $`\alpha _s.095`$ (for the purposes of discussion, assume $`N_f=3`$ at this scale). The first order term decreases the ideal gas pressure by $`6\%`$; the sum of the first and second order terms decrease the pressure by $`7\%`$ of the ideal gas value. Because the strong coupling constant runs relatively slowly with $`\mu `$ at large $`\mu `$, even at $`\mu =1`$ GeV, where $`\alpha _s.31`$, the first order term decreases the ideal gas pressure only by $`20\%`$; the first and second order terms, by $`30\%`$.

As can be seen from (2), the perturbative expansion of the thermodynamic potential is an expansion in a power series not just in $`\alpha _s`$, but in $`\alpha _s\mathrm{log}(\alpha _s)`$. The logarithm of $`\alpha _s`$ arises from the plasmon effect, where the Debye mass squared $`m_D^2\alpha _s\mu ^2`$. Because gluons at $`T=0`$ have four-dimensional phase space in loop integrals, however, the plasmon effect is relatively innocuous, and only produces logarithms, $`\mathrm{log}(m_D/\mu )\mathrm{log}(\alpha _s)`$.

This is in stark contrast to the perturbative expansion of the free energy at nonzero temperature, $`T0`$. While there is again a plasmon effect, $`m_D^2\alpha _sT^2`$, because in loop integrals static gluons at $`T0`$ have a three dimensional phase space, the perturbative expansion is not in $`\alpha _s`$, but in $`\sqrt{\alpha _s}`$. The series in $`\sqrt{\alpha _s}`$ is much worse behaved than that at $`\mu 0`$, $`T=0`$, and does not converge until very high temperatures . The convergence appears to improve after resummation , or by using Padé methods .

Consequently, the perturbative series for the thermodynamic potential may be much better behaved at $`\mu 0`$ (and $`T=0`$) than at $`T0`$ . This does not imply that a given value of $`\alpha _s`$, which is adequate to compute the thermodynamic potential, works equally well for all other quantities. In particular, the gaps for color superconductivity are nonperturbative, $`\varphi \mathrm{exp}(1/\sqrt{\alpha _s})`$ , and much smaller values of $`\alpha _s`$ appear to be required to reliably compute them . In QCD, effective models find that even when $`\mu 400`$ MeV, these gaps are at most $`100`$ MeV . As the relative change in the thermodynamic potential is only $`(\varphi /\mu )^2`$, then, for the equation of state in QCD, color superconductivity is never a large effect.

To truly know how well perturbation theory converges at $`\mu 0`$, it is imperative to compute the thermodynamic potential to $`\alpha _s^3`$. Unlike the case of $`T0`$, which is sensitive to nonperturbative effects from static magnetic gluons from $`\alpha _s^3`$ on, at $`\mu 0`$ (and $`T=0`$), the entire power series in $`\alpha _s`$ is well defined .

We now use the perturbative calculation of the thermodynamic potential for $`\mu <1`$ GeV. Since both terms $`\alpha _s`$ and $`\alpha _s^2`$ have negative coefficients, as $`\alpha _s(\mu )`$ increases with decreasing $`\mu `$, eventually the pressure vanishes. While it is clearly invalid using perturbation theory when $`p=0`$, it at least provides a well defined model of dense QCD. In Fig. 1 we show the pressure for $`\overline{\mathrm{\Lambda }}=2\mu `$; it vanishes at $`\mu _c=425`$ MeV, where $`\alpha _s.65`$. This corresponds to a quark density $`4.35\rho _0`$, where $`\rho _0`$ is the density of quarks in nuclear matter, $`3\times .16`$/fm3.

A weakness in our model is how to match the equation of state for massless quarks, (2), onto that for massive quarks and hadrons. The quark chemical potential must be larger than one third of the nucleon rest mass, minus one third the binding energy of nuclear matter, $`\mu >\mu _{\mathrm{min}}3135`$ MeV. While the pressure vanishes at $`\mu _{\mathrm{min}}`$, hadronic (or quark) matter certainly exists, with nonzero pressure, for all $`\mu >\mu _{\mathrm{min}}`$. Thus we imagine that a very “soft” equation of state for massive quarks (and hadrons) matches onto the equation of state in (2) at some $`\mu >\mu _c`$ (see also below). Consequently, $`\mu _c`$ cannot be much higher than $`\mu _{\mathrm{min}}`$.

It is this which limits the choice of $`\overline{\mathrm{\Lambda }}/\mu `$ in our model. For $`\overline{\mathrm{\Lambda }}=\mu `$, $`\mu _c=767`$ MeV. It is absurd to think that the pressure of massive quarks could be small to densities $`33\rho _0`$. Thus we do not consider this case further.

For $`\overline{\mathrm{\Lambda }}=3\mu `$, $`\mu _c=300`$ MeV when $`\alpha _s.6`$. By the Hugenholtz-van Hove theorem, when the pressure vanishes, the ratio of the total energy to the baryon number is $`E/A=3\mu _c`$. For iron, $`E/A=930`$ MeV. Thus for $`\overline{\mathrm{\Lambda }}=3\mu `$, $`E/A=3\mu _c=900`$ MeV, and, as suggested by Bodmer and Witten , strange quark matter is absolutely stable relative to hadronic matter.

While possible, we prefer an alternate view. Our perturbative equation of state is valid only in the chirally symmetric phase, for $`\mu >\mu _\chi `$. Perhaps when $`\mu <\mu _\chi `$, the true equation of state is close to our perturbative model, but vanishes smoothly as $`\mu \mu _{\mathrm{min}}`$. As discussed later, this is a model for a weakly first order (or no) chiral phase transition.

The structure of a quark star is determined by the solution to the Tolman-Oppenheimer-Volkov (TOV) equations . For the TOV equations, all that matters is the relationship between pressure and energy density. This is shown in Fig. 2 for $`\overline{\mathrm{\Lambda }}/\mu =2,3`$. The numerical solution of the TOV equations gives the mass radius relationships of Fig. 3. For any solution to the TOV equation, the chemical potential reaches its maximum value at the center of the star; as one goes out in radius, the chemical potential decreases, and equals $`\mu _c`$ (where the pressure vanishes) at the edge of the star. For $`\overline{\mathrm{\Lambda }}=3\mu `$, the maximum mass is $`M_{\mathrm{max}}2.14M_{}`$. At this mass, the radius is $`R_{\mathrm{max}}12`$ km; the chemical potential at the center of the star is $`\mu 456`$ MeV, which corresponds to a quark density of $`\rho _{\mathrm{max}}4.87\rho _0`$. When $`\overline{\mathrm{\Lambda }}=2\mu `$, $`M_{\mathrm{max}}1.05M_{}`$. At this mass, the radius is $`R_{\mathrm{max}}5.81`$ km; the chemical potential at the center is $`\mu 649`$ MeV, corresponding to a quark density of $`\rho _{\mathrm{max}}14\rho _0`$.

To help understand these results, it is useful to compare to the equation of state of a nonideal bag model:

$$\mathrm{\Omega }(\mu )=\frac{N_f}{4\pi ^2}a_{\mathrm{eff}}\mu ^4+B_{\mathrm{eff}};$$
(6)
$`B_{\mathrm{eff}}`$ is an effective bag constant, and the parameter $`a_{\mathrm{eff}}`$ measures deviations from ideality. A common choice is to take $`a_{\mathrm{eff}}`$ from the thermodynamic potential to one loop order, with a fixed value of the coupling constant: $`a_{\mathrm{eff}}=12\alpha _s/\pi `$ . In a bag model, the relationship between pressure and energy density is linear, $`p=(ϵ4B)/3`$, irrespective of the value of $`a_{\mathrm{eff}}`$. Thus we can uniformly scale $`p`$, $`ϵ`$, and $`B`$ together, so the maximum mass and radius satisfy a simple scaling relation, $`M_{\mathrm{max}}R_{\mathrm{max}}1/B^{1/2}`$ .

For the chemical potentials of relevance to a quark star, somewhat surprisingly we find numerically that the pressure in (2) can be well approximated by the effective bag model, (6). This can be seen from Fig. 2, where the relationship between pressure and energy is very nearly linear. When $`\overline{\mathrm{\Lambda }}=3\mu `$, the pressure agrees with a bag model with $`B_{\mathrm{eff}}^{1/4}=157`$ MeV and $`a_{\mathrm{eff}}=.626`$ to within $`2\%`$ for $`\mu :300470`$ MeV. This is close to the usual value in the MIT bag model, $`B^{1/4}=145`$ MeV . When $`\overline{\mathrm{\Lambda }}=2\mu `$, the pressure agrees with a bag model with $`B_{\mathrm{eff}}^{1/4}=223`$ MeV and $`a_{\mathrm{eff}}=.628`$ to within $`4\%`$ for $`\mu :425650`$ MeV.

Consequently, the mass radius relationships for our quark stars agree well with a bag equation of state. To $`5\%`$, the maximum masses and radii scale according to $`1/B_{\mathrm{eff}}^{1/2}`$. The shape of the mass radius curve is also the same as for a bag model. Notably, light quark stars have small radii. This is because for light stars, $`MM_{}`$, the chemical potential at the center of the star is near $`\mu _c`$, and the equation of state is controlled by that of massless fields, minus a bag constant.

Our results for $`\overline{\mathrm{\Lambda }}=2\mu `$ can be compared to other equations of state for dense QCD . All of these can be viewed as models in which there is nonideality at a scale significantly higher than nuclear matter densities. Ref.  uses a Schwinger-Dyson model, and finds $`M_{\mathrm{max}}0.7M_{}`$ and $`R_{\mathrm{max}}9.`$ km. Refs.  and  use models with massive quasiparticles, so that the masses act as a type of nonideality. Ref.  finds $`M_{\mathrm{max}}0.8M_{}`$ and $`R_{\mathrm{max}}4.`$ km; ref.  finds $`M_{\mathrm{max}}1.35M_{}`$ and $`R_{\mathrm{max}}10.`$ km. We note that in Nambu-Jona-Lasino models, stars with a quark core do not arise, even at the maximum mass .

What about the manifestly nonperturbative phase in which chiral symmetry is spontaneously broken, $`\mu <\mu _\chi `$? To understand this, consider an expansion in a large number of colors . The usual large $`N_c`$ limit is to let $`N_c\mathrm{}`$ at fixed $`N_f`$. Since quark loops are suppressed in this limit, however, gluons are only affected by quarks when $`\mu N_c^{1/4}`$ or larger. This is in contrast to the transition at a nonzero temperature, which occurs at a temperature $`N_c^0`$ . We then consider a generalized large $`N_c`$ limit, in which $`N_f\mathrm{}`$ at fixed $`N_f/N_c`$ . The quark thermodynamic potential is $`\mathrm{\Omega }N_fN_c\mu ^4`$, and the quark number density $`nN_fN_c\mu ^3`$. In this limit, $`\mu _\chi N_c^0`$, as at nonzero temperature.

For the purposes of our discussion, all that matters is that baryon masses $`m_BN_c`$. The baryon chemical potential is related to the quark chemical potential as $`\mu _B=N_c\mu N_c`$. When the baryons are nonrelativistic, so their Fermi momenta $`k_f1`$, the baryon number density is $`n_Bd_Bk_f^3`$. The degeneracy of baryons is at least $`N_f`$, and could easily be larger, $`N_f^2`$. The baryon thermodynamic potential is (naively) $`\mathrm{\Omega }_Bd_Bk_f^5/m_B`$. For $`k_f1`$, however, everything is fine: the baryon density is $`N_c`$ that of quarks, and $`\mathrm{\Omega }_B(1/N_c)\mathrm{\Omega }`$. In terms of the quark chemical potential, however, $`\mu =\mu _B/N_c`$; with $`m_q=m_B/N_c`$, $`\mu m_q+k_f^2/(2m_BN_c)+\mathrm{}`$. That is, for $`k_f1`$, the region in $`\mu `$ over which hadrons are a reasonable description is small, $`1/N_c^2`$. Also, the binding energy of nuclear matter is automatically $`1/N_c^2`$.

From this large $`N_c`$ argument, we conclude that a hadronic description is applicable only in a very narrow region of $`\mu `$; for larger $`\mu `$, still of order one, a quark description is appropriate. This need not be a true phase transition; rather, simply that the thermodynamic potential may be very difficult to compute in terms of hadrons, but relatively simple in terms of quarks. For example, when $`k_fN_c^{1/2}`$, so $`\mu m_q1/N_c`$, naively $`\mathrm{\Omega }_Bd_BN_c^{3/2}N_c^{1/2}\mathrm{\Omega }`$. This cannot be right — the thermodynamic potential of baryons cannot dominate that of quarks. The only resolution is that there are cancellations — analogous to those which occur for baryon-meson couplings — which greatly reduce the baryon thermodynamic potential, so that it is comparable to that of quarks. In other words, the hadronic thermodynamic potential must “soften” whenever $`\mu m_q>1/N_c^2`$.

Once one is away from this (narrow) hadronic window in $`\mu `$, the appropriate equation of state for $`\mu <\mu _\chi `$ is that for massive quarks. There are then two possibilities.

The first is that the thermodynamic potential for massless quarks matches, more or less smoothly, onto that of massive quarks. This requires either a weakly first order chiral transition, or perhaps just crossover. Below $`\mu _\chi `$, as $`\mu \mu _{\mathrm{min}}`$ the quark thermodynamic potential vanishes in a fashion typical of massive particles; within $`1/N_c^2`$ of $`\mu _{\mathrm{min}}`$, a hadronic description is applicable.

This is illustrated by the choice $`\overline{\mathrm{\Lambda }}=3\mu `$. For a star at its maximum mass, at the center $`\mu 456`$ MeV; as the radius increases, $`\mu `$ decreases, until one enters a phase in which chiral symmetry is broken (assuming $`\mu _\chi <456`$ MeV). Thus realistically, all of our “quark” stars have mantles with massive quarks, and then hadrons. For stars near the maximum mass, we assume this mantle is thin, and does not greatly alter its properties. As the mass decreases, however, the portion of the star in the chirally symmetric phase becomes small, until the entire star is composed entirely of massive quarks and hadrons. At this point, the relationship between the stars mass and radius is no longer like that of Fig. 3. Instead, it looks like that of nonrelativistic matter, for which the radius increases as the mass decreases. It is necessary for the chiral transition to be weakly first order, or a smooth crossover, for the mass-radius curve to be continuous.

The second possibility is that the equation of state for massless quarks does not match smoothly onto that for massive quarks, with a strongly first order chiral transition . As the thermodynamic potential approaches the ideal gas limit at large $`\mu `$, and vanishes at $`\mu _{\mathrm{min}}`$, this requires that the pressure is small at a value of $`\mu _\chi \mu _{min}`$. Below $`\mu _\chi `$, by construction the pressure of massive quarks is small, with a “soft” equation of state .

This occurs if $`\overline{\mathrm{\Lambda }}=2\mu `$. At the maximum mass, $`\mu 649`$ MeV at the center; if $`\mu _\chi 649`$ MeV, most of the star is composed of massless quarks. As the mass of the star decreases, so does the amount in the chirally symmetric phase. If the chiral phase transition is strongly first order, eventually one jumps to a second branch, in which the chemical potentials are always $`<\mu _\chi `$ throughout the star. Stars on this second branch are composed only of massive quarks and hadrons, with a maximum mass and radius like that of “ordinary” neutron stars. A strong first order chiral transition is necessary to ensure that there are two, distinct branches. Using toy models for the thermodynamic potential, numerically we obtained solutions to the TOV equations which display two branches: we patched a thermodynamic potential for massive quarks, for $`\mu <\mu _c\mu _\chi `$, onto that for massless quarks, for $`\mu >\mu _c`$. In this case, our stars of massless quarks constitute a third class of compact stars, after white dwarfs and “ordinary” neutron stars .

Most pulsars have masses $`1.5M_{}`$ . The MACHO project has also reported micro-lensing events for the Large Magellanic Cloud with masses $`M=0.15`$–0.9$`M_{}`$ . For a weakly first order chiral phase transition, if MACHO events are hadronic stars, they must have large radii. For a strongly first order chiral phase transition , MACHO events could be quark stars, with small radii, and pulsars might represent the second branch. We stress that our numbers for the maximum mass and radius are meant only to be suggestive. Even so, we believe that our conclusions are qualitatively correct; a weak (or no) chiral phase transition leads to one type of compact objects, a strongly first order chiral transition, to two.

We thank R. Harlander, J. Lenaghan, A. Peshier, K. Rajagopal, A. Rebhan, K. Redlich, D. Rischke, and especially L. McLerran for fruitful discussions. We thank the U.S. Department of Energy for their support under Contract No. DE-AC02-98CH10886; E.S.F., for the support of CNPq (Brazil); and J.S.B., for the support of RIKEN and BNL.