text:
# Phase Separation close to the density-driven Mott transition in the Hubbard-Holstein model

## Abstract

The density driven Mott transition is studied by means of Dynamical Mean-Field Theory in the Hubbard-Holstein model, where the Hubbard term leading to the Mott transition is supplemented by an electron-phonon (e-ph) term. We show that an intermediate e-ph coupling leads to a first-order transition at $`T=0`$, which is accompanied by phase separation between a metal and an insulator. The compressibility in the metallic phase is substantially enhanced. At quite larger values of the coupling a polaronic phase emerges coexisting with a non-polaronic metal.

Since long time it is known that strong electron-electron (e-e) interactions can drive a metallic system insulating mott . The correlation caused Metal-Insulator transition (MIT), usually known as the Mott transition, occurs in several compounds like $`V_2O_3`$, $`CaTiO_3`$ and is often theoretically investigated in the framework of the Hubbard Hamiltonian. However, although this latter certainly captures the fundamental properties of correlated electronic systems, it lacks the lattice degrees of freedom, which can play a crucial role in many respects. Therefore, not only the interplay between the electronic correlation and the electron-phonon (e-ph) interaction is a relevant issue in general cg ; MILLIS , but interesting specific effects can arise from it, like in doped fullerenes, where a substantial e-ph coupling associated to Jahn-Teller distortion can be synergetic with the large electron-electron interaction to produce the observed rather large superconducting temperatures capone . Many experimental results also point toward an important role of phonons in the strongly correlated cuprates. In this light, it has also been proposed that the e-ph coupling could induce charge instabilities in the cuprates cg ; cdg  thereby relating the interplay between strong e-e and e-ph interactions to the issue of electronic phase separation (PS), stripe formation, and to the isotopic dependence of the pseudogap crossover temperatures in these materials andergassen . This connection between e-ph coupling and charge instabilities suggest that lattice degrees of freedom may greatly influence the thermodynamical stability of a strongly correlated electron system. The present work precisely focuses on the role of the e-ph coupling in modifying or emphasizing the physical properties of electrons close to their Mott insulating status. Specifically we study the role of phonons in changing the order of the density-driven MIT and in inducing PS close to the transition.

From the technical point of view we exploit the great progress made in this field with the dynamical mean-field theory (DMFT), which emerged as the first theoretical tool able to give a complete characterization of the Mott transitionrevdmft . This technique, which becomes exact in the infinite coordination limitmetzner , is in fact able to treat with the same accuracy both the metallic and the insulating phase, and it has been successfully applied to the Mott transition in the Hubbard model. Moreover, DMFT allows to treat on the same footing the e-e and e-ph interactions, which is crucial to our study. The main drawback of this approach is the momentum-independence of the self-energy, and the neglect of non-local interactions. This is not expected to introduce significant biases when the interactions are local, like in the model we study.

One of the most striking features of the DMFT treatment of the Mott transition is that both the interaction-driven transition at half-fillinggeorges ; kotliar  and the density-driven transitionfisher  are characterized by a coexistence of solution in some region of parameters. In the latter case, which we address in this paper, there is a region of chemical potentials, delimited by two curves $`\mu _{c1}(U)`$ and $`\mu _{c2}(U)`$, in which a metallic solution with $`n1`$ and an insulating one with $`n=1`$ coexistfisher . Naturally, the stable solution is the one that minimizes the grand-canonical free energy. At $`T=0`$ the metal is stable in the whole coexistence region, and is continuously connected to the insulator fisher . Moreover, the metal joins onto the insulator with a finite slope in the $`n`$-$`\mu `$ curve, i.e., it has a finite compressibility $`\kappa =n/\mu `$ fisher ; kajueter . At finite temperature the energetic balance gets more involved, and a first-order transition occurs along an intermediate line between $`\mu _{c1}`$ and $`\mu _{c2}`$, leading to PS between the two phasessahana .

In this work we study how this scenario is modified by the inclusion of a short-range interaction with phonons. We consider the e-ph interaction in its simplest model realization, the Holstein molecular crystal model, in which tight-binding electrons interact with local modes of constant frequencyholstein . The Hamiltonian is

$`H=`$ $`t{\displaystyle \underset{i,j}{}}c_i^{}c_j+H.c.+U{\displaystyle \underset{i}{}}n_in_i`$ (1)
$`g{\displaystyle \underset{i}{}}n_i(a_i+a_i^{})+\omega _0{\displaystyle \underset{i}{}}a_i^{}a_i,`$
where $`c_i`$ ($`c_i^{}`$) and $`a_i`$ ($`a_i^{}`$) are, respectively, destruction (creation) operators for fermions and for local vibrations of frequency $`\omega _0`$ on site $`i`$, $`t`$ is the hopping amplitude, $`g`$ is an e-ph coupling.

Previous studies have identified two different dimensionless e-ph couplings, $`\lambda =2g^2/\omega _0t`$, which measures the energetic convenience to form a bound state, and $`\alpha =g/\omega _0`$ which controls the number of excited phonons. The relevance of each coupling depends on the adiabatic ratio $`\gamma =\omega _0/t`$: If $`\gamma `$ is small $`\lambda `$ is the control parameter, while for large $`\gamma `$ the physics is mainly controlled by $`\alpha `$csg . Here we work with $`\gamma =0.2`$, and therefore we use $`\lambda `$ to measure the strength of the e-ph coupling.

In DMFT, the lattice model is mapped onto an impurity problem subject to a self-consistency condition, which contains all the information about the lattice. In the case of the Hubbard-Holstein model we have an Anderson impurity with a local phonon on the impurity site. If we work in the Bethe lattice of half-bandwidth $`t`$ the self-consistency enforcing the DMFT solution is given by

$$\frac{t^2}{4}G(i\omega _n)=\underset{k}{}\frac{V_k^2}{i\omega _nϵ_k},$$
(2)
where $`ϵ_k`$ and $`V_k`$ are the energies and the hybridization parameters of the Anderson impurity model. We use exact diagonalization to solve the impurity model caffarel . The technique consists in restricting the sum in Eq. (2) to a finite and small number of levels $`N_s1`$. The discretized model is then solved using the Lanczos method, which allows to compute the Green’s function at $`T=0`$. The method converges rapidly as a function of $`N_s`$, and just a few levels are enough to achieve convergence. Here we mainly present results for $`N_s=10`$, having checked in specific cases that no measurable change occurs in the physical quantities by increasing $`N_s`$note-ph .

As anticipated above, we discuss how the density-driven Mott transition is affected by the e-ph interaction. The main trend, which can be anticipated on intuitive grounds, is that the e-ph interaction favors the insulating phase with respect to the metal. We note that in Ref. sahana  the same effect is obtained by raising the temperature, because in that case the insulating state has a larger entropy. Here, when the e-ph interaction is turned on in the metallic state, two competing effects take place. On one side, the electron gains potential energy if the lattice is distorted, while, on the other side, the effective mass is enhanced, leading to a loss in kinetic energy. As a result of this competition, the electrons can not completely exploit the e-ph coupling to lower the energy. On the other hand, the electrons in the Mott insulator are already localized, so that they can gain potential energy by coupling with phonons without losing kinetic energy. This simple argument already tells us that the energetic balance between the two phases will be changed in favor of the insulator. For the same reasons we can also expect the region of $`\mu `$ where the insulator exists to be larger, i.e., that $`\mu _{c1}`$ decreases with increasing $`\lambda `$ notadoping .

We notice that simplified approaches, as the one based on Lang-Firsov and squeezing transformationsfeinberg  also suggest that the insulating behavior may be favored by the electron-phonon interaction. Within this approach, the Hubbard-Holstein model is transformed into a Hubbard model with reduced hopping, in which the insulating behavior is clearly favored with respect to the case where the electron-phonon interaction is absent.

These naive arguments are confirmed and based on solid ground by the DMFT calculations. We have studied in detail both the solutions as a function of $`\lambda `$ for values of $`U/t`$ larger than the critical value $`U/t3`$ for the Mott transition at half-filling revdmft . In Fig. 1 we plot the density as a function of the chemical potential for $`U/t=10`$ and various values of $`\lambda `$. For $`\lambda =0`$ (not shown) we recover the known results for the Hubbard model, with two solutions in a really small $`\mu `$ interval, and the metal is stable in the whole coexistence region, as it can be checked by comparing the grand canonical potentials $`\mathrm{\Omega }=E\mu N`$ of the two phase. In the insets of the various panels we plot the difference $`\mathrm{\Omega }_{met}\mathrm{\Omega }_{ins}`$. For small values of $`\lambda `$ we observe a slight increase of the coexistence region, but the metal stays at lower energy than the insulator (see upper panel of Fig. 1). For $`\lambda >1`$ the two curves start to cross and the grand canonical potential of the insulator becomes lower than the one of the metal for some range of chemical potential (in the middle panel we show the result for $`\lambda =2`$), leading to a jump of the occupation number by continuously varying the chemical potential. The transition becomes therefore of first order, and the chemical potential has a plateaux with respect to the density. As a result, the system is not stable in the interval of density between which the jump occurs, and it undergoes PS between the insulating solution with $`n=1`$ and a metallic solution with some density $`n_{ps}(\lambda )`$. $`n_{ps}`$ is found to be a decreasing function of $`\lambda `$. In Fig. 2 we show the PS region for $`U/t=10`$.

If we further increase the coupling, the effect of the e-ph interaction becomes even more dramatic. In the case of $`U/t=10`$, and $`\lambda 2.2`$, the insulating solution continuously evolves into a phase with density slightly smaller than 1. Also for this coupling the insulator is energetically favored in the coexistence region (see instet of lower panel in Fig. 1). An inspection to other physical quantities allows us to better characterize this strong-coupling phase as a polaronic state. In fact, if we compute the phonon displacement distribution function $`P(X)=X|\psi _0\psi _0|X`$, the evolution from the $`n=1`$ insulator to this state is accompanied by a change from a unimodal distribution of the phonon displacement to a bimodal distribution characteristic of the polaronic state. On the other hand, the metallic solution does not display polaronic features for these values of $`\lambda `$ and so the PS in this region, takes place between a polaronic insulator and a non-polaronic metal.

Our DMFT results show that the Hubbard-Holstein model displays PS close to half-filling, due to the first-order phase transition between the insulator and the metal. Even if PS occurs close to half-filling, just like in the large$`N`$ treatment of the same model, the present result is physically quite different. In Ref. becca  a divergent compressibility was indeed found within the metallic phase. Here this instability is prevented by the stabilization of the insulator. Nevertheless, we can follow the metallic solution in the regime where it is metastable, and check whether the e-ph interaction may lead to a divergence, or at least to a sizeable enhancement of the charge compressibility.

In the infinite-$`U`$ limit, the large$`N`$ approach has shown that the charge compressibility is well represented by a simple RPA-like formula,

$$\kappa (\lambda )=\frac{2N^{}}{1+(F_0^S)_e+(F_0^S)_{ph}},$$
(3)
where $`\kappa (\lambda )`$ is the compressibility for a given value of the e-ph coupling $`\lambda `$, $`N^{}`$ is the quasiparticle density of states per spin at the Fermi level, $`(F_0^S)_{e(ph)}`$ is the electronic (phononic) contribution to the symmetric Landau amplitude, and $`(F_0^S)_{ph}=4N^{}g^2/\omega _0`$, which differs from our $`\lambda `$ for the presence of the quasiparticle DOS instead of the real DOS. From (3) it is easy to derive $`\kappa (\lambda )/\kappa (0)=1/(1\lambda \kappa (0)t)`$, in which all the correlation effects are contained in the purely electronic compressibility calculated at $`\lambda =0`$. This implies that, if the electronic compressibility is large, a small or moderate e-ph coupling is sufficient to make the system unstable in the charge channel.

Coming back to the results shown in Fig. 1, the $`n`$-$`\mu `$ curves of the metallic solutions have a small finite slope when they approach $`n=1`$. However, if we consider the compressibility at fixed density as a function of $`\lambda `$, it turns out that the e-ph interaction is substantially increasing the compressibility. In Fig. 3 we present the inverse of the normalized compressibility as a function of $`\lambda `$ at $`U/t=5`$ for different values of the electron density, ranging from $`n=0.70`$ to $`n=0.95`$. An enhancement of $`\kappa (\lambda )/\kappa (0)`$ that varies from almost 6 in the $`n=0.70`$ case, to around 3 in the $`n=0.95`$ case corresponds to the decreasing linear behavior of all the curves seen in the figure. We also calculated the values of the renormalization factor $`Z`$ which is inversely proportional to the quasiparticle density of states. Since we found that $`Z`$ varies much less than $`\kappa `$ as a function of $`\lambda `$, such an enhancement of the compressibility cannot be due to mass renormalization alone. This also indicates that, within the Landau Fermi liquid approach, the total $`F_0^S`$ is expected to be renormalized by the e-ph interaction like in the large$`N`$ calculations. We extract informations on this renormalization of the Landau amplitudes by simply assuming that a relation similar to Eq. (3) holds also for the DMFT solution. We therefore fitted the DMFT data with the relation

$$\frac{\kappa (\lambda )}{\kappa (0)}=\frac{1}{1\beta \lambda \kappa (0)t},$$
(4)
where $`\beta `$ is the only fitting parameter. The results of this fit, together with the values of $`\kappa (\lambda )`$ in the inset, are shown in Fig. 3.

For almost the whole range of densities the fit is extremely accurate up to $`\lambda =2.2`$, with just a small correction to the ideal large-N result $`\beta =1`$. The closer we get to half-filling the more $`\beta `$ approaches 1, even if at $`n=0.95`$ the strong-coupling data display some fluctuations. The values of $`\beta `$ are reported in the figure.

In light of the ability of the simple Eq. (4) to describe the way the e-ph interaction modifies the scenario determined by electron correlation, the small values of $`\lambda `$ needed to obtain a divergent compressibility in the large$`N`$ approach cg ; becca  can be related to the smaller value of the electronic $`F_0^S`$ found in the large$`N`$ with respect to the one found in the present DMFT approach.

Our findings could be related to some experimental observations. First of all the presence of a (phonon-induced) PS close to a Mott insulating phase is suggestive of a chemical potential plateaux inferred from photoemission experiments in underdoped cuprates ino . The formation of doping-induced metallic states with a chemical potential slowly varying with doping inside a Mott gap was also recently deduced from tunneling experiments bozovic . In the same experiment, the small mixing of the metallic (superconducting) phase with the antiferromagnetic insulating one indicates that a spatial separation of the two phases can have small interfaces suggesting that the presence of long-range Coulomb forces can easily induce small-scale (presumably stripe-like) domain formation. We also remarkably found that for rather large values of the e-ph coupling ($`\lambda >2`$ in the $`U=10t`$ case) a PS can occur between an insulator away from half-filling with a finite density of polarons and a good metal with nearly free carriers. Also this situation of coexisting (but spatially separated) polarons and free carriers could find a realization in underdoped cuprates MUELLER .

In conclusion, we have studied the effect of the e-ph interaction on the density-driven Mott transition for large values of $`U`$. The e-ph interaction favors the insulating solution with respect to the metal. This effect, besides a quantitative modification of the coexistence region between the two solutions, determines a first-order transition between an insulator and a metal at some value of the chemical potential $`\mu _c`$, and to a PS between the Mott insulator with $`n=1`$ and a metal with a finite density $`n1`$. This PS occurs for $`\lambda 1`$. For substantially larger values of the e-ph coupling a polaron crossover occurs very close to half-filling and PS takes place between this polaronic state and a non-polaronic metal.

We acknowledge fruitful discussions with S. Ciuchi and support from MIUR Cofin 2001, prot. 2001023848, and from INFM, also through PA\_G0\_4 project.


format: text

text:
# NASCA LINES: A MYSTERY WRAPPED IN AN ENIGMA

## I Introduction

The well-known Nasca geoglyphs in the central Andes of Perú ($`14^o40^{}`$ to $`14^o55^{}`$ S lat., $`75^o00^{}`$ to $`75^o10^{}`$ W long.) have been one of the XX-th century rediscoveries that remain a mystery in many ways; our study of the lines that form these artistic figures shows that they may also be a riddle to modern mathematics. These geoglyphs were made by the pre-Inca inhabitants of the dry pampa by removing all the dark material (mainly stones) and revealing the clear granular substratum. Some indirect methods indicate that the geoglyphs were made before the year $`600`$ (pottery shreds on the surfaces Silver ), $`525\pm 80`$ (radiocarbon of a wooden post found nearby Reiche ), or from $`190`$ BC to $`660`$ (organic material “varnish” Dorn ); but due to their constitution, it is not possible to apply the usual direct dating methods and therefore, one is forced to look for alternative ways to estimate their age. The Nasca region has been subjected to immense perturbations (ignorant tourism mainly) and some of the geoglyphs are terribly damaged; although fortunately there are some good photographs that were taken when the figures had just been rediscovered, presumably by the mathematician Maria Reiche Reiche , who donated her lifelong work to humankind.

Research in various and diverse fields has shown that there is a tendency for some natural systems that evolve with time to increase their complexity Boyajian ; animal physiology and behavior Nat1 , and human cultural manifestations not being the exception Nat2 ; Nat3 . Having observed that there are different degrees of simplicity in the Nasca figures, even for the subset of biomorphs with deca- and hectometric dimensions, we have tried to establish their relative ages by measuring their fractal dimension. From our present day perspective, we conclude that simpler figures would have to be older than more complex ones, in agreement with the assumption of a tendency to greater complexity as their creators perfected their activity in time.

## II Method

Before measuring the fractal dimension of the available figures: three birds, two whales, a dog, a monkey, a spider, a pair of hands, and a tree (some shown in Fig. 1), we applied a careful digital process that involved intensity filtering, contrast magnification, and background-noise elimination, that resulted in a significantly increased contour-definition of the photographic images; two other methods for processing the figures were also employed to confirm consistent results.

The determination of the fractal dimension for each figure was carried out by applying the widely used box-counting method Yorke ; Nezadal . Briefly, the dimension of a black object over a white background or vice versa, and denoted by $`D`$, is defined as:

$$D=\underset{\epsilon 0}{lim}\frac{\mathrm{ln}N(\epsilon )}{\mathrm{ln}(1/\epsilon )},$$
(1)
where $`N(\epsilon )`$ is the number of boxes in a square grid of side-size $`\epsilon `$ required to cover the object in question. In our case, the side of the box was increased from one image pixel ($`0.32mm`$) to the whole image size ($`72dpi`$ density resolution), and all images were previously scaled to a box $`160mm`$ wide (the corresponding minimum height is $`220mm`$); this enabled us to avoid the loss of details when the small scales were analyzed. Fig. 2 shows the fitting of two straight lines to the measurements obtained at different scales. As a uniform rule for avoiding the uncertainties due to upper and lower scale cut-offs, we discarded the three lowest and three highest points in each graph before performing the fitting; the linear correlation coefficients for the ten analyzed figures varied from $`0.995`$ to $`0.999`$. Our results clearly show discernible differences in the values of the fractal dimension of distinct figures that imply an ordering of the biomorphs (Fig. 3). We will now use these differences and try to establish a chronological ordering of the figures in the following section.

## III Results and Conclusions

If one now adopts the point of view Tay1 ; Tay2  that complexity and richness increases steadily from dimension one to dimension two, and the hypothesis of an increasing complexity in drawing as time evolves is accepted, one then has to place lines before trapezoids in time, i. e., lines older than plane figures (Fig. 3). This result is in agreement with some recent archaeological findings Aveni ; Morrison , and indicates that the Nasca artists started their cultural activity with straight lines and, acquiring experience, moved on to the use of more complex geometries. The chronological order shown in Fig. 3 places the figures which essentially are combinations of straight, parallel lines before the ones that use mainly curves and nonparallel draws; the suggested relationship would then imply that the first figures are older and less sophisticated than the second set of figures, which look more complex in terms of the geometric elements just mentioned.

Although our study goes a different way when using the results –trying to establish a chronological ordering– it is similar to the analysis carried out for J. Pollock’s paintings which enabled the establishing of an initial link between the visual preference or aesthetic quality and the value of the fractal dimension of an object; it has also opened a new field of research on the physiological effects of visual perception (see TaylorSA  for a review).

We do not expect our premises to be accepted as unquestionable but await for them to be contested and challenged; all we can ascertain at present, is that the Nasca imagery thus reinflames the debate on the meaning of art, a problem which may well continue as an elusive question even for modern mathematical tools.

The authors sincerely acknowledge the invaluable suggestions made by R. P. Taylor.


text:
# Hole Doping Dependence of the Coherence Length in 𝐿𝑎_{2-𝑥}𝑆𝑟_𝑥𝐶𝑢𝑂₄ Thin Films

## Abstract

By measuring the field and temperature dependence of magnetization on systematically doped $`La_{2x}Sr_xCuO_4`$ thin films, the critical current density $`j_c(0)`$ and the collective pinning energy $`U_p(0)`$ are determined in single vortex creep regime. Together with the published data of superfluid density, condensation energy and anisotropy, for the first time we derive the doping dependence of the coherence length or vortex core size in wide doping regime directly from the low temperature data. It is found that the coherence length drops in the underdoped region and increases in the overdoped side with the increase of hole concentration. The result in underdoped region clearly deviates from what expected by the pre-formed pairing model if one simply associates the pseudogap with the upper-critical field.

It is well established that the superconductivity in the cuprates is originated from the insulating anti-ferromagnetic parent compound by doping holes into the system. The superconducting transition temperature $`T_c`$ increases in the underdoped region and decreases in the overdoped side, showing a parabolic shape of $`T_c`$ vs. doped hole concentration $`p`$, i.e., $`T_c/T_c^{max}=182.6(pp_c)^2`$ with $`p_c0.150.16`$, $`T_c^{max}`$ being the maximum $`T_c`$ at the optimal doping point $`p=p_c`$. At the underdoped side, an incomplete gap, namely the pseudogapPseudogap , has been observed near the Fermi surface at a temperature $`T^{}`$, which is much higher than $`T_c`$. It has been a core and longstanding issue to know the correlation between the pseudogap and the superconductivity. Based on the Uemura plotUemura  and other experimental data, Emery and KievelsonEmeryKievelson  put forward a picture to attribute the superconducting transition in the underdoped region to a Bose-Einstein condensation, while that in the overdoped side a BCS-type. This model suggests that the pseudogap corresponds to the energy scale of pre-formed pairing ( perhaps in the sense of resonating-valence-bond singlet pairingRVB  ) and the superconducting condensation occurs at $`T_c`$ when the long range phase coherence is established. This picture has been cited to explain the anomaly strong Nernst effect in normal stateXuZA ; WangYY1 , the time-domain optical conductivityOrenstein , the Andreev reflection experimentDeutscher , etc. although the upper temperature for these effects are still far below $`T^{}`$.

In the conventional BCS scenario, the coherence length $`\xi `$, or the upper critical field $`H_{c2}=\mathrm{\Phi }_0/2\pi \xi ^2`$, represent the pairing strength. Meanwhile, $`\xi `$ measures the size of Cooper pairs and vortex core, and near $`H_{c2}`$ the normal vortex core with diameter of 2$`\xi `$ overlap with each other. In the BEC-BCS picture mentioned aboveEmeryKievelson , the pairing and the condensation occur as two steps leading to a non-trivial issue of defining the superconducting coherence length. However, the vortex core size can be determined from experiment directly. In order to get a deeper insight to the underlying mechanism of cuprate superconductors, it is necessary to determine the doping dependence of the vortex core size or the superconducting coherence length $`\xi `$. This has turned out to be a very difficult task due to the very high $`H_{c2}(0)`$ in cuprates. In early 90s, by measuring the temperature dependence of diamagnetization or excess conductivity under different magnetic fields, many groupsHc2Early1 ; Hc2Early2  determined the slope of $`H_{c2}`$ near $`T_c`$ based on the mean-field fluctuation theoryScaling . Then according to the theory of Werthamer-Helfand-Hohenberg ( WHH )WHH , the value $`H_{c2}(0)`$ is derived by $`H_{c2}(0)=0.69|dH_{c2}/dT|_{T_c}`$. It was found that the coherence length obtained in this way decreases with the increase of p in the underdoped regionHc2Early1 . Because the true $`H_{c2}(0)`$ is very high, the method based on the WHH theory is indirect, one naturally argues whether the $`H_{c2}(0)`$ determined in this way is meaningful. In this Letter we present a new approach to derive the zero-temperature coherence length or vortex core size directly by analyzing the vortex collective pinning energy on systematically doped $`La_{2x}Sr_xCuO_4`$ thin films.

Seven $`La_{2x}Sr_xCuO_4`$ films made by pulsed laser ablation technique, with (00l) orientation and nominal composition x = 0.06, 0.10, 0.125, 0.15, 0.175, 0.20 and 0.25 have been measured. It is found that $`p`$ 0.07 for the sample $`x=0.06`$ due to easy oxygen absorption, and $`px`$ for other samples. The thickness of the films are typically 100nm and all samples have rectangular shape with the same size. The films are smooth and uniform without any traces of screw dislocations as shown by the atomic force microscope. It is understood that the vortex pinning in these films are from dense point-like disordersGriessen . Details about the fabrication were published previouslyFilm . A vibrating sample magnetometer(VSM 8T, Oxford 3001 )is used to measure the magnetization hysteresis loop( MHL ) of the films with the magnetic field perpendicular to the surface. In Fig.1 we show the $`T_c`$ vs. $`p`$ of these films, which can be roughly described by $`T_c/T_c^{max}=182.6(pp_c)^2`$ except for a small plateau near $`p=1/8`$ as widely known. Fig.2 shows the temperature dependence of the magnetic critical current density $`j_c`$ at 0.03 T determined from the MHLs via $`j_c=A\mathrm{\Delta }M`$, where $`\mathrm{\Delta }M`$ is the width of the MHL, and $`A`$ is a $`T`$ and $`H`$ independent pre-factor being the same for all films since they have the same shape and sizes. We intentionally choose the low field data for the present work because the single vortex creep condition is satisfied hereBlatter . In this semi-logarithmic plot, it is clear that $`j_c`$ drops slowly with increasing temperature in the low temperature region showing a linear-like behavior $`Logj_c=Logj_c(0)\alpha T`$ with $`\alpha `$ clearly depending on the hole concentration $`p`$.

In addition to slight quantum tunnelling behavior below about 5 K, the flux motion in HTS is dominated by the thermally activated flux creep ( TAFC ) with the dissipative electric field $`E=v_0Bexp(U/k_BT)`$, where $`v_0`$ is the attempting speed for flux creep, $`U`$ is the effective activation energy, $`B`$ the local magnetic inductance. A general relation $`U(T,j_c(T))=(U_p(T)/\mu )[(j_c^0(T)/j_c(T))^\mu 1]`$ is usually assumed, where $`\mu `$ is the glassy exponentVG , $`U_p(T)`$ is the collective pinning energy, $`j_c^0(T)`$ and $`j_c(T)`$ are the unrelaxed and the measured ( relaxed ) critical current density, respectively. In the low temperature region, since $`dlnj_c^0(T)/dT0`$ and $`U_p(T)U_p(0)`$, combining the TAFC equation leads to

$$U_p(0)=ln(v_0B/E)/[\frac{dlnj_c(T)}{dT}]_{T0}$$
(1)
Here $`ln(v_0B/E)`$ depends on $`p`$ very weakly and can be taken as constant ($`10`$)Wen95 . In this case $`U_p(0)`$ can be determined by taking the slope of $`Logj_c(T)`$ vs. $`T`$ in low temperature region.Quantum . A linear fit to the low temperature data gives rise to $`U_p(0)`$ ( the slope ) and $`j_c(0)`$ ( the intercept ). The pinning energy $`U_p(0)`$ ( normalized at $`p=0.175`$ ) is presented in Fig.3 and the zero-temperature critical current density $`j_c^0(0)=j_c(0)`$ is shown in Fig.1. A clear peak of $`U_p(0)`$ is observed at $`p=0.175`$ rather than at $`p=0.15`$, as illustrated by a smaller slope of $`dLogj_c(T)/dT`$ for the sample p=0.175 than for 0.15.

In the single vortex pinning regime, if the pinning is induced by point-like disorders, like in our present thin films, the low temperature collective pinning energy $`U_p(0)`$ isBlatter

$$U_p(0)ϵ^{2/3}\xi ^3H_c^2\sqrt{j_c(0)/j_0}$$
(2)
where $`ϵ`$ is the anisotropy ratio $`ϵ=m_{ab}/m_c`$, and the depairing current $`j_0H_c/\lambda \sqrt{E_c\rho _s}`$ with $`H_c`$ being the thermodynamic field, $`\lambda `$ the penetration depth, $`E_c`$ the condensation energy, $`\rho _s`$ the superfluid density, and $`E_cH_c^2`$. All quantities in eq.(2) and the following discussions are valued at T = 0 K, and for clarity we neglect the notes ”0” for $`E_c`$ and $`\rho _s`$, i.e., $`E_c(0)=E_c`$, $`\rho _s=\rho _s(0)`$. It is important to note that eq.(2) is valid for both $`\delta T_c`$-pinning and $`\delta l`$-pinningBlatter ; Griessen , pinnings induced by the spatial fluctuation of $`T_c`$ and mean free path respectively. Thus the coherence length $`\xi `$ can be obtained via

$$\xi \frac{U_p(0)^{1/3}\rho _s^{1/12}}{E_c^{1/4}j_c(0)^{1/6}ϵ^{2/9}}$$
(3)
In above equation, $`j_c(0)`$ and $`U_p(0)`$ are determined in our experiment and shown in Fig.1 and Fig.3. The superfluid density $`\rho _s`$ and the condensation energy $`E_c`$ for many systems have been determined by Bernhard et al.Bernhard  and the normalized values show roughly similar doping dependence in different systems, thus we adopt these values in the following analysis. For $`p=0.07`$ no data of $`E_c`$ is available, and we follow the trend between $`p=0.1`$ and $`p=0.14`$, and extrapolate the data down to zero at $`p=0.05`$ where $`T_c`$ vanishes ( shown by the dotted line ). The doping dependence of the anisotropy ratio $`ϵ`$ is taken from the published values by Nakamura et al.Nakamura  and Kao et al.Kao . The computed results of $`\xi `$ are shown in Fig.4 together with the experimental data for YBCO-123 systemHc2Early1  determined by analyzing the diamagnetization data based on the critical fluctuation theoryScaling . Here the calculated value of $`\xi `$ has been normalized to $`20\AA `$ at $`p_c=0.15`$ according to the STM measurementPanSH . Clearly the coherence length $`\xi `$ drops down with the increase of $`p`$ in underdoped region. It is important to note that our approach is from the low temperature data, but what turns out is very close to that from the critical fluctuation analysis near $`T_c`$. The doping dependence of $`ϵ`$ does not seem to be very influential to the final result as shown by the full circles with $`ϵ(p)`$ considered, and open circles with $`ϵ`$ as a constant. Similar result is obtained by Ando et al.Ando1  in analyzing the high field magneto-conductivity in YBCO-123 system. In Fig.4 we also plot the doping dependence of $`\xi `$ based on the preformed pairing model, i.e., by assuming a simple linear relation $`H_{c2}(p)=H_{c2}(0)(1p/0.28)`$ and $`\xi =\sqrt{\mathrm{\Phi }_0/2H_{c2}\pi }`$, with $`\mu _0H_{c2}(0)=360T`$. It is clear that our result deviates from the expectation of the preformed pairing model when simply relating the pseudogap as the upper critical field. In the following we try to understand the data within the BEC-BCS picture.

If the $`H_{c2}(0)`$ is defined as the critical field that destroys the superconducting state at zero K, it should have very different meanings within different schemes for the superconductivity. In the conventional BCS picture, the energy scale for pairing, i.e. the superconducting gap $`\mathrm{\Delta }_s`$, is much smaller than the energy corresponding to the phase stiffness $`E_{ph}\rho _s\mathrm{}^2/m^{}`$, therefore the condensation occurs as soon as the pairs are formed. The $`H_{c2}`$ in this case is understood as the critical field for pair breaking. The condensation energy $`E_c=1/2N_F\mathrm{\Delta }_s^2`$, thus one can have a simple estimation that $`E_cN_F\mathrm{\Delta }_s^2\rho _s\mathrm{\Delta }_s`$, where $`N_F`$ is the electronic density of states near Fermi surface at $`T_c`$. Since $`\rho _s1/\lambda ^2`$, $`H_{c2}1/\xi ^2`$, $`E_cH_c^21/\xi ^2\lambda ^2`$, one has then $`H_{c2}\mathrm{\Delta }_s`$. Near $`T_c`$, the superconducting gap vanishes, thus $`H_{c2}`$ becomes zero and the coherence length diverges.

In the underdoped region, it is conjectured that the condensation is BEC type with the Cooper pairs formed above $`T_c`$. If we assume that $`E_c=\rho _s\mathrm{\Delta }_s`$, since now the energy scale $`\mathrm{\Delta }_s`$, the pair breaking energy is much larger than $`E_{ph}`$, $`H_{c2}`$ is now controlled by $`\rho _s`$. From the known correlationUemura  $`\rho _sT_c`$, we can expect that $`H_{c2}\rho _sT_c`$. This is very close to what we obtained in the unerdoped region ( inset of Fig.4 ). It is interesting to note that this simple relation $`H_{c2}T_c`$ has also been obtained by Muthukumar and WengWeng  within the bosonic resonating-valence-bond theory. Following the BEC picture for underdoped region, the vortex structure is also very different from that expected from the BCS picture. Firstly, the core state may be still gapped or partly gapped, reflecting the pseudogap feature. This may explain why the zero bias conductance peak ( ZBCP ) has not been observed within the vortex core by STMPanSH ; STMCore , NMRNMRCore  and specific heatWenCore . Secondly, the pairing amplitude keeps almost constant across the vortex core, while the superfluid density $`\rho _s`$ now drops to zero in the core center and gradually recovers its full value over a distance ( now the vortex core size $`\xi `$). In this sense the vortex core size, or the newly defined coherence length $`\xi `$ may be proportional to the London penetration depth $`\lambda `$ although $`\xi `$ is much smaller than $`\lambda `$. Thirdly, when $`T_c`$ is approached from below, both the vortex core size or the coherence length and the penetration depth will diverge due to the vanishing static superfluid density. At low temperatures when the magnetic field is increased, the vortex core tough each other at $`H_{c2}`$ in the conventional sense and the extremely diluted superfluid density will not be strong enough to sustain the long-life vortices. However the short-life vortices can still be generated above $`T_c(H)`$ which gives rise to the dynamic phase stiffness as observed by the ultra-fast optical conductivity measurementOrenstein . The anomalous strong Nernst signal in the normal state can also be interpretated in this wayXuZA ; WangYY1 ; WangYY2  although the upper magnetic field for observing the Nernst signal may be more close to the pair breaking field, which is proportional to the pseudogap. This interpretation is consistent with the result of high field measurement by Shibauchi et al.Shibauchi  who found that the pseudogap can be closed by the Zeeman splitting under a high magnetic field, while the field corresponding to the superconducting peak on the c-axis magneto-resistance $`H_{sc}`$ scales with $`T_c`$ in the entire doping regime. From our present data and analysis, we cannot evaluate whether there is a quantum critical transition of ground statesBernhard , or the transition concerning only the pairing symmetryDagan  near p=0.19. However in any case, in the overdoped side the pairing gap drops with $`p`$ and the BCS picture is restored. Our results on the doping dependence of the vortex core size call for future work on observing the vortex structure directly by scanning tunnelling microscope.

In conclusion, the hole doping dependence of the coherence length or vortex core size has been determined directly from the low temperature data in wide doping regime. It is found that the coherence length drops in the underdoped region and increases in the overdoped side with the increase of hole concentration. This behavior clearly prevents the idea of relating the pseudogap with the upper-critical field in underdoped region. The BEC-BCS picture can be applied to interpret the data leading to consequences about a new structure and many novel features of the vortex line.

###### Acknowledgements.

This work is supported by the National Science Foundation of China (NSFC ), the Ministry of Science and Technology of China, and the Knowledge Innovation Project of Chinese Academy of Sciences. The work at Penn State is supported in part by NSF under Grant No. 9623889 and ONR under Grant No. N00014-00-1-0294. We thank J. L. Tallon for kind permission of using their published data ( Ref.20 ). We are grateful for fruitful discussions with Z. Y. Weng, T. Xiang and Q. H. Wang.


text:
# Dense quark matter in compact stars

## 1 Introduction

It has become clear over the last few years that the phase diagram of QCD is much richer than originally believed. In addition to the hadronic phase with which we are familiar and the quark gluon plasma (QGP) that is predicted to lie at temperatures above $`170\mathrm{MeV}`$, there is a whole family of “colour superconducting” phases, that are expected to occur at high density and low temperature . The essence of colour superconductivity is quark pairing, driven by the BCS mechanism, which operates when there exists an attractive interaction between fermions at a Fermi surface. The QCD quark-quark interaction is strong, and is attractive in many channels, so we expect cold dense quark matter to generically exhibit colour superconductivity. Moreover, quarks, unlike electrons, have colour and flavour as well as spin degrees of freedom, so many different patterns of pairing are possible. This leads us to expect a rich phase structure in matter beyond nuclear density.

Colour superconducting quark matter may occur naturally in the universe, in the cold dense cores of neutron stars. (It may also be possible to create it in low-energy heavy ion colliders, such as the proposed Compressed Baryonic Matter facility at GSI Darmstadt.) There densities are above nuclear density, and temperatures are of the order of tens of keV. Thus most work on signatures has focussed on properties of colour superconducting quark matter that would affect observable features of compact stars, and I will discuss some of these below.

## 2 Phase diagram of QCD

The strange quark plays a crucial role in the phases of QCD. In Fig. 1 I show two conjectured phase diagrams for QCD, one for a light strange quark and one for a heavy strange quark. In both cases, along the horizontal axis the temperature is zero, and the density rises from the onset of nuclear matter through the transition to quark matter. Compact stars are in this region of the phase diagram, although it is not known whether their cores are dense enough to reach the quark matter phase. Along the vertical axis the temperature rises, taking us through the crossover from a hadronic gas to the quark gluon plasma. This is the regime explored by high-energy heavy-ion colliders (see other contributions to these proceedings).

If the effective strange quark mass is low enough in quark matter at a few times nuclear density, then there is a direct transition from nuclear matter to colour-flavour-locked (CFL) quark matter . In the CFL phase the strange quark participates symmetrically with the up and down quarks in Cooper pairing—this is described in more detail below. If the strange quark is too heavy to pair symmetrically with the light quarks at these densities, then there will be an interval of some other phase. This may be crystalline colour superconductivity  or some form of single-flavour pairing .

## 3 Colour superconductivity

To understand the essence of colour superconductivity, we can borrow some intuition from condensed matter physics . QCD is asymptotically free—the interaction becomes weaker as the momentum grows—so at sufficiently high density and low temperature, there is a Fermi surface of almost free quarks. The interactions between quarks near the Fermi surface are certainly attractive in some channels (quarks bind together to form baryons) and it was shown by Bardeen, Cooper, and Schrieffer (BCS)  that if there is any channel in which the interaction is attractive, then there is a state of lower free energy than a simple Fermi surface. That state arises from a complicated coherent superposition of pairs of particles (and holes)—“Cooper pairs”.

We can understand the BCS mechanism in an intuitive way as follows. The Helmholtz free energy is $`F=E\mu N`$, where $`E`$ is the total energy of the system, $`\mu `$ is the chemical potential, and $`N`$ is the number of fermions. The Fermi surface is defined by a Fermi energy $`E_F=\mu `$, at which the free energy is minimized, so adding or subtracting a single particle costs zero free energy. Now switch on a weak attractive interaction. It costs no free energy to add a pair of particles (or holes), and the attractive interaction between them then lowers the free energy of the system. Many such pairs will therefore be created in the modes near the Fermi surface, and these pairs, being bosonic, will form a condensate. The ground state will be a superposition of states with all numbers of pairs, breaking the fermion number symmetry.

Since pairs of quarks cannot be colour singlets, the resulting condensate will break the local colour symmetry $`SU(3)_{\mathrm{colour}}`$. We call this “colour superconductivity”. Note that the quark pairs play the same role here as the Higgs particle does in the standard model: the colour-superconducting phase can be thought of as the Higgs phase of QCD. For more detailed discussions of calculations of the favored pairing patterns, see the reviews .

### 3.1 Three flavours: Colour-flavour locking (CFL)

The favored pairing pattern at high densities, where the strange quark Fermi momentum is close to the up and down quark Fermi momenta, is “colour-flavour locking” (CFL). In QCD with three flavours of massless quarks the Cooper pairs cannot be flavour singlets, and both colour and flavour symmetries are necessarily broken . Both NJL  and gluon-mediated interaction calculations  agree that the attractive channel exhibits a pattern called colour-flavour locking (CFL),

$$\begin{array}{c}q_i^\alpha q_j^\beta _{1PI}^{}C\gamma _5\left((\kappa +1)\delta _i^\alpha \delta _j^\beta +(\kappa 1)\delta _j^\alpha \delta _i^\beta \right)\\ [SU(3)_{\mathrm{colour}}]\times \underset{\left[U\left(1\right)_Q\right]}{\underset{}{SU(3)_L\times SU(3)_R}}\times U(1)_B\underset{\left[U\left(1\right)_{\stackrel{~}{Q}}\right]}{\underset{}{SU(3)_{C+L+R}}}\times _2\end{array}$$
(1)
Colour indices $`\alpha ,\beta `$ and flavour indices $`i,j`$ run from 1 to 3, Dirac indices are suppressed, and $`C`$ is the Dirac charge-conjugation matrix. The term multiplied by $`\kappa `$ corresponds to pairing in the $`(\mathrm{𝟔}_S,\mathrm{𝟔}_S)`$, which although not energetically favored breaks no additional symmetries and so $`\kappa `$ is in general small but not zero . The Kronecker deltas connect colour indices with flavour indices, so that the condensate is not invariant under colour rotations, nor under flavour rotations, but only under simultaneous, equal and opposite, colour and flavour rotations. Since colour is only a vector symmetry, this condensate is only invariant under vector flavour+colour rotations, and breaks chiral symmetry. The effect of including a strange quark mass is to distort the CFL pairing , and a sufficiently large strange quark mass will unpair the strange quark completely , yielding a two-flavour colour superconductor (2SC) , and at lower values it may induce a flavour rotation of the condensate known as “kaon condensation” . The features of the CFL pattern of condensation are

* The colour gauge group is completely broken. All eight gluons become massive. This ensures that there are no infrared divergences associated with gluon propagators.
* All the quark modes are gapped. The nine quasiquarks (three colours times three flavours) fall into an $`\mathrm{𝟖}\mathrm{𝟏}`$ of the unbroken global $`SU(3)`$, so there are two gap parameters. The singlet has a larger gap than the octet.
* A rotated electromagnetism (“$`\stackrel{~}{Q}`$”) survives unbroken. It is a combination of the original photon and one of the gluons.
* Two global symmetries are broken, the chiral symmetry and baryon number, so there are two gauge-invariant order parameters that distinguish the CFL phase from the QGP, and corresponding Goldstone bosons which are long-wavelength disturbances of the order parameter. When the light quark mass is non-zero it explicitly breaks the chiral symmetry and gives a mass to the chiral Goldstone octet, but the CFL phase is still a superfluid, distinguished by its baryon number breaking.
* The symmetries of the 3-flavour CFL phase are the same as those one might expect for 3-flavour hypernuclear matter , so it is possible that there is no phase transition between them.

In a real neutron star we must require electromagnetic and colour neutrality  (ignoring charge-separated phases). It turns out that this penalizes the 2SC phase relative to the CFL phase . The reason is that the CFL phase has already paid most of the cost of neutrality, since it brings the $`u`$, $`d`$, and $`s`$ Fermi surfaces close together , and the 2SC phase only pairs 4 of the 9 quark colours and flavours, so it has much less pair binding energy than the CFL phase.

The arguments made in Ref.  are model-independent, based on simplified assumptions about the dependence of the constituent strange quark mass $`M_s`$ on $`\mu `$ and expanding the free energy in powers of $`M_s/\mu `$. The NJL calculation of Ref.  handles $`M_s\mu `$ and includes the coupling between the chiral condensate and quark condensate gap equations. The net result is that once neutrality is imposed, there is no (or very little) density range in which 2SC is the phase with the lowest free energy.

## 4 Mass-radius relationship for compact stars

The high density and relatively low temperature required to produce colour superconducting quark matter may be attained in compact stars, also known as “neutron stars”, since it is often assumed that they are made primarily of neutrons. Typical compact stars have masses close to $`1.4M_{}`$, and are believed to have radii of order 10 km. (For reviews see Ref. ).

Colour superconductivity affects the equation of state at order $`(\mathrm{\Delta }/\mu )^2`$. It also gives mass to excitations around the ground state: it opens up a gap at the quark Fermi surface, and makes the gluons massive. One would therefore expect it to have a profound effect on transport properties, such as mean free paths, conductivities and viscosities.

Although the effects of colour superconductivity on the quark matter equation of state are subdominant, they may have a large effect on the mass-radius relationship. The reason for this is that the pressure of quark matter relative to the hadronic vacuum contains a constant (the “bag constant” $`B`$) that represents the cost of dismantling the chirally broken and confining hadronic vacuum,

$$p\mu ^4+\mathrm{\Delta }^2\mu ^2B.$$
(2)
If the bag constant is large enough so that nuclear matter is favored (or almost favored) over quark matter at $`\mu 320\mathrm{MeV}`$, then the bag constant and $`\mu ^4`$ terms almost cancel, and the superconducting gap $`\mathrm{\Delta }`$ may have a large effect on the equation of state and hence on the mass-radius relationship of a compact star .

In Ref.  Sanjay Reddy and I explored the effect of quark pairing on the $`M`$-$`R`$ relationship at values of the bag constant that are consistent with nuclear phenomenology. Fig. 2 shows the mass-radius curve for a plausible model of dense matter: the Walecka nuclear equation of state, and quark matter with physically reasonable values of the bag constant $`B^{1/4}=180`$ MeV ($`B=137\mathrm{MeV}/\mathrm{fm}^3`$) and strange quark mass $`m_s=200`$ MeV . Curves for unpaired ($`\mathrm{\Delta }=0`$) and colour-superconducting ($`\mathrm{\Delta }=100\mathrm{MeV}`$) quark matter are shown. At these values the stars are typically “hybrid”, containing both quark matter and nuclear matter. The solid lines in Fig. 2 correspond to stars that either have no QM at all, or a sharp transition between NM and QM: the core is made of quark matter, which is the favored phase at high pressure, and at some radius there is a transition to nuclear matter, which is favored at low pressure. The transition pressure is sensitive to $`\mathrm{\Delta }`$, for reasons discussed earlier. The dashed lines are for stars that contain a mixed NM-QM phase. In all cases we see that light, large stars consist entirely of nuclear matter. When the star becomes heavy enough, the central pressure rises to a level where QM, either in a mixed phase or in its pure form, occurs in the core. As can be seen from the figure the transition density is very sensitive to $`\mathrm{\Delta }`$.

In Ref.  we kept the bag constant fixed, assuming that it could be fixed by other observations, and we treated the quark matter as free quarks with a pairing energy. It is interesting to see what happens when we relax these assumptions. To see how closely quark matter can mimic nuclear matter, we tune the bag constant to keep constant the density of nuclear matter at its transition to quark matter. To allow for effects of quark interactions beyond Cooper pairing, we follow the parameterization of Fraga et. al. , who find that the $`𝒪(\alpha _s^2)`$ pressure for three flavours over a finite interval in $`\mu `$ can be mimicked by a bag model inspired form given by

$$\begin{array}{ccc}\hfill P_{\alpha _s^2}(\mu _{\mathrm{low}}\mu \mu _{\mathrm{high}})& =& \frac{3}{4\pi ^2}a_{\mathrm{eff}}\mu ^4B_{\mathrm{eff}},\hfill \\ \hfill a_{\mathrm{eff}}& & 1c.\hfill \end{array}$$
(3)
Choosing the renormalization scale $`\mathrm{\Lambda }=2\mu `$ they find that $`a_{\mathrm{eff}}=0.628`$ ($`c=0.372`$) and $`B_{\mathrm{eff}}^{1/4}=223`$ MeV in the range $`\mu =425`$ to $`650`$ MeV and for $`\mathrm{\Lambda }=3\mu `$ they find that $`a_{\mathrm{eff}}=0.626`$ ($`c=0.374`$) and $`B_{\mathrm{eff}}^{1/4}=157`$ MeV in the range $`\mu =300`$ to $`470`$ MeV.

First, we explore the effect of a colour superconducting gap $`\mathrm{\Delta }`$ and perturbative correction $`c`$ on the mass-radius relationship. We fix the bag constant by requiring that that the nuclear to quark matter phase transition occur at nuclear matter baryon density $`\rho =1.5n_{\mathrm{sat}}`$. The resultant $`M(R)`$ curves are shown in Fig. 3.

1. The stars resulting from quark matter equations of state without perturbative correction ($`c=0`$,blue lines) are smaller and lighter. This is because they are inherently highly favored, and so requiring the phase transition to occur at $`\rho =1.5n_{\mathrm{sat}}`$ leads to a large bag constant. By contrast, the equations of state with perturbative correction ($`c=0.4`$, red lines) are inherently less favored, and require a low bag constant to bring the phase transition down to $`\rho =1.5n_{\mathrm{sat}}`$.
2. The stars with $`c=0.3`$ have mass-radius relationships that are very similar to the pure nuclear APR98 matter. In fact, for the case where there are perturbative corrections but no colour superconductivity the equations of state ($`p(\mu )`$) are so similar that our program found a series of phase transitions back and forth between CFL and APR98 up to $`\mu =546\mathrm{MeV}`$ (baryon density $`\rho =5.4n_{\mathrm{sat}}`$). This is why the $`c=0.3`$ red dotted curve lies almost exactly on top of the solid black (APR98) curve, even though there was a phase transition from APR98 to CFL at $`\rho =1.5n_{\mathrm{sat}}`$ (which is first attained when the APR98 star reaches a mass of $`0.315M_{}`$, $`R=13.3\mathrm{km}`$).
3. In Ref.  we showed that at fixed bag constant, colour superconductivity has a strong effect on the mass-radius relationship of compact stars. Here, by comparing the dashed lines with the dotted lines in Fig. 3, we see that it is difficult to distinguish the effect of colour superconductivity from a change in the bag constant. In Fig. 3, as we vary parameters $`c`$ and $`\mathrm{\Delta }`$ of the quark matter equation of state, the bag constant is tuned to maintain a constant value of the nuclear density at the transition to quark matter, and in this situation colour superconductivity only makes a small difference to the mass-radius relationship.

## 5 Other phenomenology

Mass and size are the most obvious observable properties of a compact star, but various groups are working on other possible signatures of colour-superconducting quark matter.

Interfaces and mixed phases. These were studied in Ref. , and it was found that a mixed phase only occurs if the surface tension of the interface is less than about $`40\mathrm{MeV}/\mathrm{fm}^2=0.2\times (200\mathrm{MeV})^3`$, a fairly small value compared to the relevant scales $`\mathrm{\Lambda }_{\mathrm{QCD}}200\mathrm{MeV}`$, $`\mu 400\mathrm{MeV}`$. A sharp nuclear-quark interface will have an energy-density discontinuity across it, which will affect gravitational waves emitted in mergers, and also the $`r`$-mode spectrum and the damping forces to which $`r`$-modes are subject.

Crystalline pairing (the “LOFF” phase). This is expected to occur when two different types of quark have sufficiently different Fermi momenta that BCS pairing cannot occur . This is a candidate for the intermediate phase of Fig. 1, where the strange quark mass, combined with requirements of weak equilibrium and charge neutrality, gives each quark flavour a different Fermi momentum. The phenomenology of the crystalline phase has not yet been worked out, but recent calculations using Landau-Ginzburg effective theory indicate that the favored phase may be a face-centered cubic crystal , with a reasonably large binding energy. This raises the interesting possibility of glitches in quark matter stars.

Cooling by neutrino emission. The cooling rate is determined by the heat capacity and emissivity, both of which are sensitive to the spectrum of low-energy excitations, and hence to colour superconductivity . CFL quark matter, where all modes are gapped, has a much smaller neutrino emissivity and heat capacity than nuclear matter, and hence the cooling of a compact star is likely to be dominated by the nuclear mantle rather than the CFL core . Other phases such as 2SC or LOFF give large gaps to only some of the quarks. Their cooling would proceed quickly, then slow down suddenly when the temperature fell below the smallest of the small weak-channel gaps. This behavior should be observable.

$`r`$-mode spindown. The $`r`$-mode is a bulk flow in a rotating star that, if the viscosity is low enough, radiates away energy and angular momentum in the form of gravitational waves. One can rule out certain models for compact stars on the grounds that they have such low damping that they could not support the high rotation rates observed in pulsars: $`r`$-mode spindown would have slowed them down. Madsen  has shown that for a compact star made entirely of quark matter in the CFL phase, even a gap as small as $`\mathrm{\Delta }=1`$ MeV is ruled out by observations of millisecond pulsars. It remains to extend this calculation to the more generic picture of a quark matter core surrounded by a nuclear mantle.

Acknowledgments

I thank the organizers of SQM 2003. The work reported in Section 4 was performed in collaboration with Sanjay Reddy, and was supported by the UK PPARC and by the U.S. Department of Energy under grant number DE-FG02-91ER40628.


text:
# The First Supernova Explosions in the Universe

## 1. INTRODUCTION

One of the most important challenges in modern cosmology is to understand how the cosmic dark ages ended (e.g., Barkana & Loeb 2001; Bromm & Loeb 2003a). Recent numerical simulations have presented evidence that the first (Pop III) stars, formed out of metal-free gas in low-mass halos at redshifts $`z20`$, were predominantly very massive with $`M_{}100M_{\mathrm{}}`$ (Bromm, Coppi, & Larson 1999, 2002; Abel, Bryan, & Norman 2000, 2002). As the first stars mark the crucial transition from a smooth, homogeneous universe to an increasingly complex and structured one, the important question arises: How did the first stars die?

The answer to this question sensitively depends on the precise mass of a Pop III star. In particular, if the star has a mass in the narrow interval $`140M_{}260M_{\mathrm{}}`$, it will explode as a pair-instability supernova (PISN), leading to the complete disruption of the progenitor (Fryer, Woosley, & Heger 2001; Heger et al. 2003). Pop III stars with masses below or above the PISN range are predicted to form black holes. This latter fate is not accompanied by a significant dispersal of heavy elements into the intergalactic medium (IGM), since most of the newly synthesized metals will be locked up in the black hole. The PISN, however, will contribute all its heavy element production to the surrounding gas.

How did the pristine IGM get enriched with the first heavy elements? To account for the widespread presence of metals in the Ly$`\alpha `$ forest at $`z5`$ (e.g., Songaila 2001), star formation in low-mass systems at $`z10`$ has been proposed as a likely source (Madau, Ferrara, & Rees 2001). The metals produced in these low-mass halos can more easily escape from their shallow potential wells than those at lower redshift (e.g., Aguirre et al. 2001a, b). In addition, the enriched gas has to travel much shorter distances between neighboring halos at these early times, and it might therefore have been easier to establish a uniform metal distribution in the IGM.

In this Letter, we present numerical simulations of the supernova explosions that end the lives of massive Pop III stars. In contrast to earlier work, both analytical (e.g., Larson 1974; Dekel & Silk 1986; Scannapieco, Ferrrara, & Madau 2002; Furlanetto & Loeb 2003) and numerical (e.g., Mori, Ferrara, & Madau 2002; Thacker, Scannapieco, & Davis 2002; Wada & Venkatesan 2003), we here focus on the minihalos that are the sites for the formation of the very first stars. The gas in these halos, with masses of $`M10^510^6M_{\mathrm{}}`$ and virializing at $`z20`$, can only form stars in the presence of a sufficient amount of H2, the only viable coolant in these systems (e.g., Haiman, Thoul, & Loeb 1996; Tegmark et al. 1997; Yoshida et al. 2003a).

An important motivation for this numerical study is to substantiate our earlier semi-analytical discussion of a fundamental transition in the character of star formation at high redshifts, from very massive Pop III stars to lower-mass Pop II stars at $`z15`$ (Mackey, Bromm, & Hernquist 2003; see also Schneider et al. 2002). This transition redshift is an important ingredient in determining the reionization history of the universe, which has very recently been constrained by the WMAP satellite (e.g., Cen 2003; Ciardi, Ferrara, White 2003; Haiman & Holder 2003; Wyithe & Loeb 2003a,b; Sokasian et al. 2003; Yoshida et al. 2003b).

## 2. Numerical Methodology

The evolution of the dark matter and gas components is calculated with the parallel version of GADGET (Springel, Yoshida, & White 2001) in its “conservative entropy” formulation (Springel & Hernquist 2002). This code combines the smoothed particle hydrodynamics (SPH) method with a hierarchical (tree) gravity solver. Here, we briefly describe the additions to GADGET which are necessary for the investigation of explosions in zero-metallicity gas.

We include all the relevant cooling processes for gas temperatures $`10^2T10^9`$K (Cen 1992; Bromm et al. 2002), as well as heating due to the photoionization of H, He, and He+. The chemical reaction network comprises the 9 species H, H+, H-, H2, H$`{}_{}{}^{+}{}_{2}{}^{}`$, $`e^{}`$, He, He+, and He++, including the reactions given in Haiman et al. (1996). Our implicit, backwards-differencing method to solve the coupled set of rate equations is fast and accurate (see Bromm et al. 2002).

To achieve the dynamical range required to follow the evolution of the gas from cosmological initial conditions down to the formation of high density clumps, which are the immediate progenitors of Pop III stars, and of the subseqent supernova (SN) blast waves, we have implemented a resampling technique within SPH (see Bromm & Loeb 2003b for details). This method of refining a coarser, parent simulation, and following the further collapse with increased resolution by ‘splitting’ the parent SPH particles into a number, $`N_{\mathrm{ref}}`$, of child particles has already been successfully applied to problems in Galactic star formation (Kitsionas & Whitworth 2002).

## 3. Simulations: Pre-Explosion Phase

### 3.1. Initial Conditions and Clump Formation

We initialize our simulations at $`z_i=100`$ according to a standard $`\mathrm{\Lambda }`$CDM model with parameters $`\mathrm{\Omega }_m=1\mathrm{\Omega }_\mathrm{\Lambda }=0.3`$, $`\mathrm{\Omega }_B=0.045`$, $`h=0.7`$, and $`\sigma _8=0.9`$, close to the values measured by WMAP (Spergel et al. 2003). Our computational box is periodic, has a comoving length of $`100h^1`$kpc, and contains $`128^3`$ particles each in dark matter and gas. For the initial gas temperature we adopt the value of 200 K. The fractional free-electron abundance is initialized as $`x_e=4.6\times 10^4`$, and the hydrogen molecule abundance as $`f_{\mathrm{H}_2}=2\times 10^6`$ (Anninos & Norman 1996).

We follow the evolution until $`z20`$, at which point the gas in one of the halos has collapsed to a density $`n10^4`$cm-3, thus violating the requirement that the local Jeans mass be resolved everywhere: $`M_J>M_{\mathrm{res}}500M_{\mathrm{}}`$ (e.g., Bate, Bonnell, & Bromm 2003). This first high-density clump to form in the simulation comprises a mass of $`10^3M_{\mathrm{}}`$ within a radius of $`1`$ pc. The clump lies at the center of a virialized dark matter minihalo of total mass $`10^6M_{\mathrm{}}`$ and virial radius $`150`$ pc. If allowed to collapse further, the gas in the clump would eventually form one very massive star (VMS), or a small multiple of them (see Abel et al. 2002; Bromm et al. 2002). Here, we assume that the clump will form a single VMS, producing ionizing radiation during its short lifetime before finally exploding as a PISN. Next, we describe the radiative feedback that precedes the explosion.

### 3.2. Radiative Feedback

To be able to simulate the propagation of the SN remnant into the surrounding gas with sufficient resolution, we have selected a Cartesian volume around the location of the high-density clump. The SPH particles within this volume of physical length 1 kpc serve as the parent particles for the resampling procedure, such that each parent spawns $`N_{\mathrm{ref}}=16`$ child particles (see Bromm & Loeb 2003b for details). After refinement, there are $`5\times 10^5`$ SPH particles in the dark matter minihalo that hosts the PISN explosion, and the improved mass resolution is $`M_{\mathrm{res}}20M_{\mathrm{}}`$. The dark matter component, on the other hand, is not refined.

It is crucial to take into account the radiative feedback from the VMS on the density structure in the minihalo. Otherwise, we would severely overestimate the radiative losses in the early stages of the explosion. Recently, Bromm & Loeb (2003c) have implemented an approximate scheme to address this radiative feedback that allows one to estimate the evolution of the central density field without carrying out a fully self-consistent radiative transfer calculation. We here briefly summarize this approach.

We assume that the gas within the minihalo is rapidly ionized on a timescale that is short compared to the dynamical time in the halo. Specifically, we initialize our feedback simulation, aimed at determining the hydrodynamic response to the internal ionizing radiation, with a fully-developed H II region already in place. We choose a Strömgren radius of $`R_\mathrm{S}200`$ pc, and assume that the ionizing radiation is unattenuated within $`R_\mathrm{S}`$ (Osterbrock 1989). We now include photoheating terms in the entropy equation, such that (in units of ergs s-1 cm-3): $`\mathrm{\Gamma }r^2d\nu f_\nu (R_{})\sigma _\nu (\nu \nu _{\mathrm{th}})/\nu _{\mathrm{th}}`$ for $`r<R_\mathrm{S}`$, and $`\mathrm{\Gamma }=0`$ otherwise. Here, $`f_\nu (R_{})`$ is the flux at the surface of a very massive Pop III star (Bromm, Kudritzki, & Loeb 2001b). We have also included the relevant photoreactions in the chemical reaction network (see Haiman et al. 1996). We allow the photoheating due to the central point source to go on for $`\tau _{\mathrm{MS}}2\times 10^6`$yr, the approximate main-sequence lifetime of a VMS. Under the influence of the Pop III radiation field, the gas inside the H II region relaxes to a temperature of $`T3\times 10^4`$K, and the density within a central core of $`R_{\mathrm{core}}20`$ pc reaches $`n1`$ cm-3.

At the end of the radiative feedback phase, we switch off the central point source, and initialize the SN explosion inside this photoheated density field. The detailed density structure very close to the VMS, at radii $`r1`$ pc, is not reliably calculated by our approximate procedure. As we are here primarily interested in the evolution of the SN shell at radii $`10`$ pc, this caveat is not expected to compromise the main results in this study.

## 4. Simulations: Post-Explosion Phase

We carry out two explosion simulations. In each case, we assume that the Pop III star in the center of the minihalo explodes as a PISN. To span the possible energy range, we consider explosion energies of $`E_{\mathrm{SN}}=10^{51}`$ and $`10^{53}`$ergs, corresponding to stars with masses $`M_{}150`$ and $`250M_{\mathrm{}}`$ (Fryer et al. 2001). We focus in particular on the latter case, which marks the upper mass limit of the PISN regime. We insert the explosion energy as thermal energy, distributing it amongst $`500`$ SPH particles. These particles lie within an initial radius of $`2`$ pc, such that the explosion is initialized with conditions appropriate for the adiabatic Sedov-Taylor (ST) phase. Since the PISN is completely disrupted without leaving a remnant behind, we can trace the subsequent fate of the metals by using the SPH particles that represent the stellar ejecta as markers.

The expansion of the SN blast wave proceeds as shown in Figure 1, where we compare our simulations with simple analytical calculations. Initially, the blast wave evolves into a roughly uniform medium at radii $`rR_{\mathrm{core}}20`$ pc, resulting in the usual ST scaling: $`R_{\mathrm{sh}}t^{0.4}`$. Once the blast wave reaches beyond the core, however, it encounters an isothermal density profile in the remainder of the halo, leading to the scaling (e.g., Ostriker & McKee 1988): $`R_{\mathrm{sh}}30\mathrm{pc}(E_{\mathrm{SN}}/10^{53}\mathrm{ergs})^{1/3}(t/10^5\mathrm{yr})^{2/3}`$. A few $`10^6`$yr after the explosion, radiative (inverse Compton) losses become important, and the SN remnant enters its final, snowplow phase. The significance of inverse Compton cooling distinguishes very high-$`z`$ SN explosions from those occuring in the present-day universe. From Fig. 1, it is evident that the numerical simulations nicely reproduce the analytical expectation.

After approximately $`10^5`$yr, when $`t_{\mathrm{cool}}R_{\mathrm{sh}}/u_{\mathrm{sh}}`$, a dense shell begins to form at a radius $`50`$ pc. The thermodynamic properties of the dense shell, the interior bubble, and the surrounding medium are summarized in Figure 2. The shocked, swept-up gas exhibits three distinct thermal phases, as can clearly be discerned in Fig. 2: the first at $`T10^5`$K, the second at $`T10^4`$K, and the third at $`T10^3`$K, with the latter two corresponding to the dense shell. The cooling from very high temperatures to $`10^4`$K, and subsequently to $`10^2`$K, proceeds in non-equilibrium, such that a fraction of free electrons persisits below $`8,000`$ K, allowing the re-formation of H2 (e.g., Shapiro & Kang 1987; Oh & Haiman 2002). The molecular fraction in the densest regions asymptotically assumes a ‘plateau’ value of $`f_{\mathrm{H}_2}3\times 10^3`$. Inside the expanding SN remnant lies a hot pressurized bubble (the red symbols in Fig. 2) that corresponds to the original stellar ejecta. This bubble drives the outward motion, and its gas is slowly cooled by adiabatic expansion and inverse Compton losses. Initially, cooling due to adiabatic expansion dominates (see the distribution of the red symbols in Fig. 2 as compared to the dotted line).

In Figure 3, we show the projected gas density $`10^6`$yr after the explosion for the $`E_{\mathrm{SN}}=10^{53}`$ergs case. Two features are apparent: Firstly, the blast wave has succeeded in completely disrupting the minihalo, dispersing most of the halo gas into the IGM. Comparing $`E_{\mathrm{SN}}`$ with the binding energy of the minihalo, $`E_{\mathrm{grav}}10^{51}`$ergs, such an outcome is clearly expected. It depends, however, on the inefficient radiative processes that cool the hot bubble gas in the initial stages of expansion. The evolution of the $`E_{\mathrm{SN}}=10^{51}`$ergs case indeed proceeds very differently. The blast wave effectively stalls at the virial radius, and the minihalo does not suffer complete disruption.

The second interesting feature in Fig. 3 is the vigorous fragmentation of the dense shell into a large number of small cloudlets. The emergence of these cloudlets, corresponding to gas at $`T10^3`$K, is not driven by gravity (e.g., Elmegreen 1994 and references therein), but is instead due to a thermal instability triggered by the onset of atomic cooling. The explosions simulated here do not result in the formation of lower-mass stars, or Pop II.5 stars in the terminology of Mackey et al. (2003). This is consistent with the prediction in Mackey et al. (2003) that Pop II.5 stars can only form in dark matter halos massive enough to be able to cool via atomic hydrogen (see also Salvaterra, Ferrara, & Schneider 2003). To further test the possibility of SN triggered star formation in the high-redshift universe, we will explore a range of cosmological environments, corresponding to different halo masses and collapse redshifts, in future work.

What is the fate of the metals that are produced in the PISN progenitor? PISNe are predicted to have substantial metal yields, of the order of $`y=M_Z/M_{}0.5`$, and for the largest stellar masses in the PISN range, most of $`M_Z`$ is made up of iron (Heger & Woosley 2002). We find that in our $`E_{\mathrm{SN}}=10^{53}`$ergs simulation, $`90`$% of the metals have escaped from the DM minihalo $`4\times 10^6`$yr after the explosion. The metals have filled most of the interior hot bubble (see inset of Fig. 3), and there has not yet been sufficient time for them to significantly mix into the dense surrounding shell. The eventual metal distribution is rather flattened, with most metals ending up in the low-density voids. A rough estimate of the resulting metallicity in the surrounding, contaminated IGM can be obtained as follows: $`Z=M_Z/M_{\mathrm{gas}}100M_{\mathrm{}}/10^5M_{\mathrm{}}10^2Z_{\mathrm{}}`$. This value is well above the critical metallicity threshold, $`Z_{\mathrm{crit}}10^{3.5}Z_{\mathrm{}}`$, such that very massive stars can only form out of gas with $`Z<Z_{\mathrm{crit}}`$ (Omukai 2000; Bromm et al. 2001a; Schneider et al. 2002). We conclude that Pop III stars can continue to form in the pristine gas within the neighboring halos, as the metals that are dispersed in the explosion do not reach them.

## 5. Summary and Conclusions

We have simulated the death of the first stars in minihalos at $`z20`$, assuming that they end their lives as PISNe. For the highest possible explosion energies, $`E_{\mathrm{SN}}=10^{53}`$ergs, we find that a large fraction ($`90`$%) of the stellar metals are dispersed into the surrounding IGM. The extent of the metal enriched region, with a physical size of $`1`$ kpc, is comparable to the radius of the relic H II region around the minihalo. In a companion paper (Yoshida, Bromm, & Hernquist 2003), we argue that if Pop III stars have led to an early partial reionization of the universe, as may be required by the recent WMAP results (e.g., Wyithe & Loeb 2003b), this will have resulted in a nearly uniform enrichment of the universe to a level $`Z_{\mathrm{min}}10^4Z_{\mathrm{}}`$ already at $`z15`$.

For the most energetic PISNe, these metals are predominantly in the form of iron, and the systems investigated in this Letter could have been responsible for an early burst of Fe enrichment. On nucleosynthetic grounds, such an initial episode of star formation has been inferred by Qian & Wasserburg (2002). Indeed, the ‘prompt inventory’ postulated by these authors could have arisen in the PISN explosions studied here. Our simulations may also be relevant for explaining the recent observations of quasars at $`z6`$, exhibiting surprisingly strong Fe II lines (Freudling, Corbin, & Korista 2003). These authors have interpreted their results as hinting at an early nucleosynthetic contribution from Type Ia SNe, but it might be challenging to accomodate the required long evolutionary timescales, first to form a white dwarf, and then to sufficiently accrete from a binary companion to trigger the explosion. In fact, if our picture is correct that normal Pop II stars can only form at $`z15`$ (Mackey et al. 2003), the available time between $`z15`$ and $`6`$ is only $``$0.7 Gyr, rendering the Type Ia scenario problematic. A PISN origin and the short timescales connected to it, however, could naturally account for the high Fe/$`\alpha `$ ratios observed (see Heger & Woosley 2002).

It would indeed be a remarkable feature of the universe, if the first stars had endowed the IGM with such a widespread, near-universal level of pre-enrichment.

We are grateful to Volker Springel for making available to us a version of GADGET. This work has been supported in part by NSF grant AST 00-71019. The simulations were performed at the Center for Parallel Astrophysical Computing at the Harvard-Smithsonian Center for Astrophysics.