text:
# Monte Carlo study of the superfluid weight in doped antiferromagnets

## I Introduction

The quest on the nature of the superconducting transition observed in doped antiferromagnets, such as the lamellar high-$`T_c`$ copper-oxides, has entered a new era. Indeed, a consensus is emerging on the importance of the zero-temperature superfluid weight $`D_s(0)`$ as the fundamental energy scale determining the corresponding transition temperature $`T_c`$. An early clue for this result came from the empirical Uemura relation, valid for underdoped cuprates, which shows roughly a proportionality between the independently measurable quantities $`T_c`$ and $`D_s(0)`$,

$$k_BT_c=AD_s(0),$$
(1)
where $`A`$ is a dimensionless parameter of order unity. Given that the superfluid weight $`D_s(T)`$, at a general temperature $`T`$, measures the “phase stiffness” of the condensate, i.e., the energy cost to produce spatial variations of its phase, Emery and Kivelson argued that (1) provides strong evidence for a phase-fluctuation driven superconducting transition. Of course, in strictly two-dimensional systems with continuous U(1) gauge (phase) symmetry, melting of long-range phase coherence with increasing temperature, and hence loss of superfluidity of the charge carriers, proceeds by thermal generation and subsequent unbinding of vortex-antivortex pairs. The detailed description of this transition is given by the Kosterlitz-Thouless (KT) theory in the context of the classical XY model. A characteristic feature of the KT transition is the discontinuous drop of $`D_s(T)`$ to zero at $`T=T_c`$, as $`T_c`$ is approached from below. Despite the inevitable rounding of this discontinuity by the weak coupling between the copper-oxide layers, Corson et al., in a remarkable experiment, have provided direct evidence of the KT nature of the superconducting transition by observing dynamic effects of thermally generated vortices in the frequency-dependent conductivity of underdoped cuprate thin films.

A common feature of theoretical models assuming a phase-fluctuation mechanism for the superconducting transition, is the prediction of a $`T`$-linear decrease in $`D_s(T)`$ at temperatures well below $`T_c`$, in accordance with the experimental observations. However, all the aforementioned models, do not incorporate explicitly the doping dependence of $`D_s(T)`$ and therefore provide no framework for a proper treatment of the high-$`T_c`$ cuprate superconductors as doped antiferromagnets. This is a major drawback since the proximity of the doped materials to the Mott insulating antiferromagnetic state at half-filling ($`n_e=1`$) leads to a vanishing $`D_s(0)`$ and $`T_c`$, as the hole concentration $`(1n_e)`$ tends to zero. The importance of the antiferromagnetic fluctuations, even well within the underdoped regime, has been recently confirmed in a dramatic way. Specifically, the elastic neutron-scattering measurements in pristine La2-xSrxCuO4, with $`x=0.10,0.12`$, have established the presence of static long-range antiferromagnetic order in the superconducting ground state, suggested by earlier work. This conclusion is supported also by the subsequent observation of conventional two-magnon Raman scattering in the same material.

The above discussion underlines the need for models of the high-$`T_c`$ cuprate superconductors that can account, on the same footing, for the intimately coupled phase and spin fluctuations of the condensate. Such a model has been suggested some time ago and shown, by one of the authors, to display flux quantization and a finite superfluid weight, i.e., superconductivity, in its ground state. Our model consists of a $`t`$-$`t^{}`$-$`J`$ Hamiltonian and a suitable $`1/N`$ expansion that provide a framework for the study of: (i) the ground-state properties, using standard “spin-wave” techniques which allow the incorporation of leading quantum-fluctuation effects, and (ii) the finite-temperature properties, using an associated classical (large-$`N`$ limit) energy functional and corresponding partition function in terms of which important physical quantities can be readily expressed and calculated by Monte Carlo simulation. Point (i) was the subject of earlier works and the results, especially for the optical properties, were found to provide support for our effective model when compared with experiment. In the present paper we consider point (ii), focusing on the study of the superfluid weight as a function of temperature and doping, a subject of current experimental and theoretical interest.

In Sec. II we give a brief description of our effective $`t`$-$`t^{}`$-$`J`$ Hamiltonian and the associated classical energy functional and partition function emerging in the large-$`N`$ limit. In Sec. III we derive an explicit expression for the superfluid weight $`D_s(T)`$ which we study by Monte Carlo simulation using the standard Metropolis algorithm. We present results for the shape of the scaled curve, $`D_s(T)/D_s(0)`$ vs $`k_BT/D_s(0)`$, the value of the dimensionless ratios $`k_BT_c/D_s(0)`$ and $`D_s(T_c)/k_BT_c`$, as well as the value of the zero-temperature slope $`D_s^{}(0)`$, clarifying their doping dependence in the regime of interest, i.e., close to half-filling. The KT nature of the transition to the superconducting state is supported by finite-size scaling analysis using the relevant Weber-Minnhagen scaling formula. In particular, our numerical results for the temperature dependence of $`D_s(T)`$, at 10% hole doping, compare reasonably well with recent measurements in underdoped La2-xSrxCuO4, for the corresponding concentration value $`x=0.10`$. Furthermore, the zero-temperature slope of the superfluid weight is predicted to approach, close to half-filling ($`n_e1`$), the universal value: $`D_s^{}(0)=k_B/2`$, consistent with the available experimental data. The latter limiting value is shown to be a hallmark of the strong antiferromagnetic correlations, present in this temperature and doping regime. Our concluding remarks are summarized in Sec. IV.

## II Effective model

The effective model under consideration is described by a $`t`$-$`t^{}`$-$`J`$ Hamiltonian expressed in terms of Hubbard operators $`\chi ^{ab}=|ab|`$ as

$$H=\underset{i,j}{}t_{ij}\chi _i^{0\mu }\chi _j^{\mu 0}+\frac{1}{2}J\underset{i,j}{}(\chi _i^{\mu \nu }\chi _j^{\nu \mu }\chi _i^{\mu \mu }\chi _j^{\nu \nu }),$$
(2)
where the index 0 corresponds to a hole, the Greek indices $`\mu ,\nu ,\mathrm{}`$ assume two distinct values, for a spin-up and a spin-down electron, and the summation convention is invoked. Here $`J`$ is the antiferromagnetic spin-exchange interaction between nearest-neighbor sites $`i,j`$ on a square lattice endorsed with periodic boundary conditions and a total number of sites $`\mathrm{\Lambda }=\mathrm{\Lambda }_x\times \mathrm{\Lambda }_y`$, where $`\mathrm{\Lambda }_x=\mathrm{\Lambda }_y`$. For the hopping matrix elements $`t_{ij}`$ we assume

$$t_{ij}=\{\begin{array}{cc}t& \text{if }i,j\text{ are nearest neighbors}\hfill \\ t^{}& \text{if }i,j\text{ are next nearest neighbors}\hfill \\ 0& \text{otherwise}.\hfill \end{array}$$
(3)
The conventions in (3) incorporate opposite signs for $`t`$ and $`t^{}`$ as it is appropriate for the hole-doped cuprates. In Ref.  we generalized the local constraint associated with (2) to $`\chi _i^{00}+\chi _i^{\mu \mu }=N`$, where $`N`$ is an arbitrary integer, and considered the commutation properties of the $`\chi ^{ab}`$ operators to be those of the generators of the U(3) algebra. A Holstein-Primakoff realization for the latter algebra in terms of hard-core bosons resolves then explicitly the local constraint and can be used to develop a perturbation theory based on the $`1/N`$ expansion, restoring the physical value $`N=1`$ at the end of the calculation.

In the presence of an external magnetic flux $`\mathrm{\Phi }`$, threading the two-dimensional lattice in an Aharonov-Bohm torus geometry, the hopping matrix elements $`t_{ij}`$ are modified by the well known Peierls phase factor and should be substituted in (2) according to

$$t_{ij}t_{ij}e^{iA_{ij}},\text{with}A_{ij}=\frac{2\pi \mathrm{\Phi }}{\mathrm{\Lambda }_x\mathrm{\Phi }_0}(𝐑_i𝐑_j)𝐞_x.$$
(4)
Here $`𝐑_i`$ is the position vector for lattice site $`i`$, $`𝐞_x`$ is the unit vector along the $`x`$-axis encircling the flux lines and $`\mathrm{\Phi }_0=2\pi \mathrm{}c/q`$ is the so-called flux quantum. As argued in Ref. , in the context of the present effective model, whereby carriers are treated as hard-core bosons, the charge $`q`$ entering $`\mathrm{\Phi }_0`$ should be set equal to $`q=2e`$, where $`e`$ is the electronic charge.

In the large-$`N`$ limit “condensation” occurs, i.e., the Bose operators become classical commuting fields. Considering only uniform density configurations, the corresponding classical energy functional resulting from (2)–(4) takes the form

$`H(\mathrm{\Phi })`$ $`=`$ $`n_e(1n_e){\displaystyle \underset{i,j}{}}t_{ij}`$ (8)
$`[\mathrm{cos}{\displaystyle \frac{\theta _i}{2}}\mathrm{cos}{\displaystyle \frac{\theta _j}{2}}\mathrm{cos}(A_{ij}+{\displaystyle \frac{\psi _i\psi _j\varphi _i+\varphi _j}{2}})`$
$`+\mathrm{sin}{\displaystyle \frac{\theta _i}{2}}\mathrm{sin}{\displaystyle \frac{\theta _j}{2}}\mathrm{cos}(A_{ij}+{\displaystyle \frac{\psi _i\psi _j+\varphi _i\varphi _j}{2}})]`$
$`+{\displaystyle \frac{n_e^2}{4}}J{\displaystyle \underset{i,j}{}}[\mathrm{cos}\theta _i\mathrm{cos}\theta _j+\mathrm{sin}\theta _i\mathrm{sin}\theta _j\mathrm{cos}(\varphi _i\varphi _j)1].`$
where $`n_e`$ is the average electronic density, the angles $`(\theta _i,\varphi _i)`$ determine the local spin direction, and the remaining parameter $`\psi _i`$ determines the local phase of the condensate. The above functional form makes apparent the coupling between the phase and spin variables of the condensate through the kinetic energy term, proportional to $`t_{ij}`$.

The description of the finite-temperature classical theory is now completed using the energy functional (LABEL:eq:H-classical) to construct the corresponding partition function $`Z(\mathrm{\Phi })`$ and free energy per lattice site $`F(\mathrm{\Phi })`$,

$$Z(\mathrm{\Phi })=e^{\beta \mathrm{\Lambda }F(\mathrm{\Phi })}=\left(\underset{i}{}\mathrm{sin}\theta _id\theta _id\varphi _id\psi _i\right)e^{\beta H(\mathrm{\Phi })},$$
(10)
where $`\beta =1/(k_BT)`$ and the integrations at each lattice site $`i`$ extend over the intervals: $`0\theta _i\pi `$, $`0\varphi _i2\pi `$, and $`0\psi _i4\pi `$. Invoking standard thermodynamic identities, important physical quantities can be readily expressed in terms of the partition function (10) and studied by Monte Carlo simulation. In the following we focus our discussion on the study of the superfluid weight $`D_s(T)`$.

## III Superfluid weight

At a finite temperature $`T`$, the superfluid weight (or helicity modulus) $`D_s(T)`$ is given by the curvature of the infinite lattice limit of the free energy $`\mathrm{\Lambda }F(\mathrm{\Phi })`$ at $`\mathrm{\Phi }=0`$,

$$D_s(T)=\mathrm{\Lambda }\left(\frac{\mathrm{\Phi }_0}{2\pi }\right)^2\left[\frac{^2F(\mathrm{\Phi })}{\mathrm{\Phi }^2}\right]_{\mathrm{\Phi }=0}.$$
(11)
Physically, $`D_s(T)`$ determines the ratio of the density of the superfluid charge carriers to their mass and hence can be related to the experimentally measurable in-plane magnetic (London) penetration depth, as noted later on in this section. Carrying out the second derivative with respect to $`\mathrm{\Phi }`$ in (11) we have more explicitly that

$`D_s(T)`$ $`=`$ $`n_e(1n_e){\displaystyle \frac{2}{z\mathrm{\Lambda }}}<{\displaystyle \underset{i,j}{}}t_{ij}|𝐑_i𝐑_j|^2`$ (17)
$`[\mathrm{cos}{\displaystyle \frac{\theta _i}{2}}\mathrm{cos}{\displaystyle \frac{\theta _j}{2}}\mathrm{cos}\left({\displaystyle \frac{\psi _i\psi _j\varphi _i+\varphi _j}{2}}\right)`$
$`+\mathrm{sin}{\displaystyle \frac{\theta _i}{2}}\mathrm{sin}{\displaystyle \frac{\theta _j}{2}}\mathrm{cos}\left({\displaystyle \frac{\psi _i\psi _j+\varphi _i\varphi _j}{2}}\right)]>`$
$`[n_e(1n_e)]^2{\displaystyle \frac{\beta }{\mathrm{\Lambda }}}<\{{\displaystyle \underset{i,j}{}}t_{ij}[(𝐑_i𝐑_j)𝐞_x]`$
$`[\mathrm{cos}{\displaystyle \frac{\theta _i}{2}}\mathrm{cos}{\displaystyle \frac{\theta _j}{2}}\mathrm{sin}\left({\displaystyle \frac{\psi _i\psi _j\varphi _i+\varphi _j}{2}}\right)`$
$`+\mathrm{sin}{\displaystyle \frac{\theta _i}{2}}\mathrm{sin}{\displaystyle \frac{\theta _j}{2}}\mathrm{sin}\left({\displaystyle \frac{\psi _i\psi _j+\varphi _i\varphi _j}{2}}\right)]\}^2>,`$
$`z=4`$ being the coordination number of the square lattice. As shown in Ref. , close to half-filling ($`n_e1`$) and for a sufficiently large $`t^{}`$, the ground state of (LABEL:eq:H-classical), in the absence of magnetic flux ($`\mathrm{\Phi }=0`$), is described by a planar spin configuration ($`\theta _i=\pi /2`$) in which the local twist angles and phases exhibit long-range order according to: $`\varphi _i=𝐐𝐑_i`$, $`\psi _i=𝐐^{}𝐑_i`$, where $`𝐐=(\pi ,\pi )`$ is the usual spin-modulating antiferromagnetic wavevector and $`𝐐^{}=(\pi ,\pi )`$ is a phase-modulating wavevector. The zero-temperature value of the superfluid weight follows then easily from (LABEL:eq:D(T)) as

$$D_s(0)=4t^{}n_e(1n_e).$$
(19)
For the typical two-dimensional model with continuous symmetry under consideration, we expect that at low but finite temperatures, the long-range order will be destroyed by the proliferation of excited Goldstone modes, leading to a $`T`$-linear decrease of $`D_s(T)`$. At higher temperatures we expect that the thermal generation and subsequent unbinding of vortex-antivortex pairs will lead eventually to a discontinuous drop of $`D_s(T)`$ to zero, at a critical point $`T=T_c`$, in a KT type of transition. In order to study numerically $`D_s(T)`$ in the whole temperature range and affirm the aforementioned physical picture, we performed a Monte Carlo simulation using the standard Metropolis algorithm. Our calculations were carried out on small lattices, with typical sizes $`\mathrm{\Lambda }_x=8,16,32,64`$. For a given temperature we performed of the order of $`10^4`$ thermalization steps and of the order of $`10^6`$ measurements. We considered values for the dimensionless ratios $`\epsilon =t^{}/t=0.45`$ and $`t/J=1.0`$, which are thought to be relevant for the copper-oxide layers, and restricted our study to the underdoped regime, i.e., to small $`(1n_e)`$ values up to 10% hole doping. The latter restriction is dictated by the fact that models of the $`t`$-$`t^{}`$-$`J`$ kind, being rather simple extensions of the Mott-Heisenberg antiferromagnetic insulator, cannot properly account for the nontrivial evolution of the electronic structure of the cuprates that occurs at higher doping values, namely, the closing of the pseudogap in the optimally doped and overdoped (Fermi liquid) regime.

Typical Monte Carlo results for the superfluid weight vs temperature are shown in Fig. 1(a), for $`(1n_e)=0.01`$, and Fig. 1(b), for $`(1n_e)=0.10`$. At low temperatures the superfluid weight has a weak finite-size dependence and displays the expected $`T`$-linear decrease. In particular, for $`T0`$ and $`n_e1`$, we have established the asymptotic form

$$\frac{D_s(T)}{D_s(0)}=1\frac{k_BT}{2D_s(0)}\text{for}T0\text{and}n_e1.$$
(20)
The analytic expression (20) is shown in Fig. 1 by a dotted line. Evidently, this asymptotic line is approached very closely from below by the low temperature numerical data already in the case of the 1% hole doping; see Fig. 1(a). From (20) follows that the zero-temperature slope of the superfluid weight approaches, close to half-filling, the parameter-free universal value:

$$D_s^{}(0)=k_B/2=0.043\text{meV}\text{K}^1\text{for}n_e1.$$
(21)
The upper limiting value (21) is a rather stringent prediction of our theory and seems, indeed, to be consistent with the available experimental data in the high-$`T_c`$ cuprate superconductors. A comparison of Fig. 1(a) with Fig. 1(b) reveals an increase in the magnitude of the slope $`D_s^{}(0)`$ upon doping, a trend also consistent with experiment.

We emphasize that the asymptotic form (20) is a physical consequence of the fact that close to half-filling the antiferromagnetic exchange energy, $`Jn_e^2`$, overwhelms the hole kinetic energy, $`t_{ij}n_e(1n_e)`$, and in particular the term $`D_s(0)`$ given by (19). Hence, in the limit $`n_e1`$ and in the relevant temperature range $`k_BTD_s(0)Jn_e^2`$ (so that $`T0`$), the thermal average (LABEL:eq:D(T)) may be simplified by freezing the spin variables $`(\theta _i,\varphi _i)`$ to their zero-temperature antiferromagnetic configuration and allowing fluctuations only in the phase variables $`\psi _i`$. In this case, the vanishing overlap between the opposite sublattice spin states leaves the direct hopping $`t^{}`$ between next-nearest-neighbor (n.n.n) sites as the only relevant process of charge transport. One can then easily show that the expressions for the energy functional (LABEL:eq:H-classical) and superfluid weight (LABEL:eq:D(T)) reduce to those of a classical XY model for the $`\psi _i`$ variables, but with only a n.n.n interaction $`I`$, where $`I=D_s(0)/2`$. The reduction of the structure of the phase fluctuations, close to half-filling, to that of the n.n.n XY model and not to the commonly assumed structure of the nearest-neighbor (n.n) XY model, is an important prediction and a direct consequence, in the context of our theory, of the presence of strong antiferromagnetic correlations in this regime. Numerically, the validity of our argument becomes apparent in Fig. 2 displaying almost coinciding $`D_s(T)`$ Monte Carlo data for the n.n.n XY model (opaque diamonds) and the $`t`$-$`t^{}`$-$`J`$ model with a very small hole concentration $`(1n_e)=0.01`$ (filled diamonds). Therefore, we may exploit the detailed results of Appendix A for the n.n.n XY model, see Eq. (A7), to conclude the asymptotic form (20) and hence the limiting value (21). The latter value, $`D_s^{}(0)=k_B/2`$, being twice that of the n.n XY model (see Appendix A), serves as a distinct hallmark of the sublattice structure of the strong antiferromagnetic correlations in the limit $`n_e1`$. Our observations here affirm also, by analogy to the well-known physics of the XY model, the presence of a KT transition for the superfluid weight of the $`t`$-$`t^{}`$-$`J`$ model, when $`n_e1`$, and allow to transcribe relevant results for the former model, established in Appendix A, to the latter, e.g.,

$`{\displaystyle \frac{D_s(T_c)}{k_BT_c}}={\displaystyle \frac{4}{\pi }}\text{for}n_e1,`$ (22)
(23)
$`k_BT_c/D_s(0)=0.4435(5)\text{for}n_e1.`$ (24)
In Fig. 1 and Fig. 2 we present $`(4/\pi )k_BT/D_s(0)`$ by a short dashed line. According to (22), in the limit $`n_e1`$ the latter line should intersect the corresponding Monte Carlo data curve of $`D_s(T)/D_s(0)`$ vs $`k_BT/D_s(0)`$ precisely at the $`k_BT_c/D_s(0)`$ value given by (24); see opaque (filled) diamonds in Fig. 2.

It should be noted that in the actual simulations of the $`t`$-$`t^{}`$-$`J`$ model we can only use finite, though possibly small, hole concentrations which inevitably lead to deviations from the limiting values (22)–(24). In order to obtain rather accurately the corresponding transition temperature $`T_c`$ we have used the finite-size scaling analysis of Weber and Minnhagen, which is appropriate for KT type of transitions. In this analysis one measures the chi-square values $`\chi ^2(T)`$ of the fitting of the Monte Carlo data for the superfluid weight, at each given temperature and for a sequence of small lattice sizes, to a certain scaling formula, derived by the latter authors from the Kosterlitz renormalization group equations. Specifically, one assumes at each temperature $`T`$ the following $`\mathrm{\Lambda }_x`$-dependence of the superfluid weight $`D_s(T,\mathrm{\Lambda }_x)`$,

$$\frac{\pi D_s(T,\mathrm{\Lambda }_x)}{2k_BT}=R(T)\left(1+\frac{1}{2\mathrm{ln}[\mathrm{\Lambda }_x/\mathrm{\Lambda }_0(T)]}\right),$$
(25)
where $`\mathrm{\Lambda }_0(T)`$ is some characteristic length of the order of the lattice constant and has no singularity at $`T=T_c`$. The logarithmic lattice size dependence involved in (25) is characteristic of the presence of vortices in the KT transition. Strictly speaking (25) is correct only at the critical point $`T=T_c`$. Given now a value for $`R(T_c)`$, the critical temperature can be determined by two alternative procedures: (a) we fix $`R(T)`$ to be $`R(T_c)`$ and use $`\mathrm{\Lambda }_0(T)`$ as the only adjustable parameter in (25) to measure the chi-square values $`\chi ^2(T)`$ of the fitting, in which case $`T_c`$ corresponds to the minimum of the $`\chi ^2(T)`$ curve, or (b) we use both $`R(T)`$ and $`\mathrm{\Lambda }_0(T)`$ as adjustable parameters in (25) and determine $`T_c`$ from the point where the $`R(T)`$ curve crosses the line $`R(T)=R(T_c)`$. The correct value for $`R(T_c)`$ should lead uniquely to the same value for $`T_c`$ in both procedures. The application of this finite-size scaling to the n.n.n XY model is summarized pictorially in Fig. 4 of Appendix A and justifies, as far as the $`t`$-$`t^{}`$-$`J`$ model is concerned, the limiting values (22)–(24) with $`R(T_c)2`$, for $`n_e1`$.

For the finite doping value $`(1n_e)=0.01`$ or $`(1n_e)=0.10`$ (and fixed $`\epsilon =t^{}/t=0.45`$, $`t/J=1.0`$), the application of the aforementioned finite-size scaling analysis, using the lattice size sequence $`\mathrm{\Lambda }_x=4,6,8,10,12`$, leads to: $`k_BT_c/D_s(0)=0.4346(9)`$ with $`R(T_c)=2.04`$, or $`k_BT_c/D_s(0)=0.3502(8)`$ with $`R(T_c)=2.51`$, respectively. To be sure, the presence in Fig. 1 of nonzero $`D_s(T)`$ values above the corresponding $`T_c`$, instead of a discontinuous drop to zero, is due to finite-size effects which grow rapidly above the estimated critical temperature; a typical behavior for a KT transition. Indeed, following an original argument by Weber and Minnhagen, we note that the success of the scaling formula (25) in the present model provides, ipso facto, strong evidence that the relevant phase transition is of the KT type. Our results show a modest increase of the jump ratio $`D_s(T_c)/(k_BT_c)=2R(T_c)/\pi `$ upon doping. The latter behavior seems physically similar to that known in the literature of the frustrated XY model given that, in the context of the $`t`$-$`t^{}`$-$`J`$ model, doping induces a form of dynamic frustration for the phase variables, via their inevitable coupling to the fluctuating spin variables.

Furthermore, our results show that the dimensionless parameter $`A=k_BT_c/D_s(0)`$, introduced in context of the empirical Uemura relation (1), is not doping independent but decreases modestly upon doping, while away from half-filling it also depends on the couplings $`\epsilon `$ and $`t/J`$. Nevertheless, for rough theoretical estimates of the KT transition temperature $`T_c`$ in terms of $`D_s(0)`$, in the underdoped regime, one may always use the universal limiting value $`A=0.4435(5)`$, for $`n_e1`$, quoted in (24). We remind that that the latter value is characteristic of the n.n.n XY model and equals to half the corresponding value of the n.n XY model (see Appendix A) commonly employed to this end.. Note that the use of $`A=0.4435(5)`$ in conjunction with (1) brings earlier theoretical overestimations of the KT transition temperature for the copper-oxides layers, derived with $`A0.9`$, down to more reasonable values. In all cases, the present analysis confirms $`D_s(0)`$ as the fundamental energy scale determining $`T_c`$ in the underdoped regime.

Having discussed the generic trends of the superfluid weight as a function of temperature and doping, in the underdoped regime, it is instructive to provide a more direct comparison of our theory with experiment. To this end we note that in the lamellar high-$`T_c`$ superconductors, the experimental value for the superfluid weight per copper-oxide plane, $`D_s^{(\mathrm{exp})}(T)`$, can be extracted from the directly measurable in-plane magnetic (London) penetration depth, $`\lambda _{ab}(T)`$, using the relation:

$$D_s^{(\mathrm{exp})}(T)=\frac{(\mathrm{}c)^2d}{4\pi q^2\lambda _{ab}^2(T)},$$
(26)
where $`d`$ is the average distance between planes and we remind that $`q=2e`$. Using the experimental data of Panagopoulos et al. for $`\lambda _{ab}(T)`$ on the underdoped La2-xSrxCuO4, with $`x=0.10`$, and the structural parameter $`d=6.64\text{Å}`$, we depict in Fig. 2 by triangles the corresponding experimental values (26) for the superfluid weight vs temperature. As a result of the weak coupling between the copper-oxide layers, the experimental data display no KT discontinuity but rather a continuous drop of the superfluid weight to zero, at a specific temperature value, that is not simply related to the ideal KT transition temperature of a copper-oxide monolayer. Corresponding Monte Carlo results for 10% hole doping are depicted in Fig. 2 by circles and calculated for a $`64\times 64`$ lattice, with $`\epsilon =0.45`$ and $`t/J=1.0`$. As noted earlier in this section, the theoretical KT transition temperature for the latter set of parameters is $`k_BT_c/D_s(0)=0.3502(8)`$, while nonzero Monte Carlo points above this value are due to finite-size effects. Evidently, our theoretical results (circles) in Fig. 2 compare reasonably well with the experimental data (triangles) throughout their common relevant temperature range, i.e., up to $`T_c`$.

In Fig. 2 we also depict Monte Carlo results for 1% hole doping (filled diamonds), as well as results for the n.n.n XY model (opaque diamonds), thus providing the theoretical lineshape of the superfluid weight vs temperature, in the limit $`n_e1`$. Clearly, it will be very interesting to have measurements of the in-plane magnetic penetration depth on La2-xSrxCuO4, with hole concentration $`x`$ as small as it is experimentally possible, to compare with the present definite theoretical prediction.

## IV Concluding remarks

In this paper we have presented a study of the temperature and doping dependence of the superfluid weight $`D_s(T)`$ in doped antiferromagnets described by the $`t`$-$`t^{}`$-$`J`$ model (2)–(3). Using Monte Carlo simulations and finite-size scaling analysis we have demonstrated that the phase fluctuations of the condensate, emerging in an appropriate classical (large-$`N`$) limit, drive superconductivity via a Kosterlitz-Thouless type of transition. Our theoretical results reproduce important generic experimental trends of $`D_s(T)`$, observed in the underdoped high-$`T_c`$ cuprate superconductors. This includes the $`T`$-linear decrease of $`D_s(T)`$ at low temperatures and the increase of the magnitude of the slope $`D_s^{}(0)`$ upon doping.

In particular, the sublattice structure of the strong antiferromagnetic correlations in the half-filled-band limit was shown to dictate the lineshape of $`D_s(T)/D_s(0)`$ vs $`k_BT/D_s(0)`$, for $`n_e1`$, to be identical to that of the n.n.n XY model. In order to check this definite theoretical prediction we have suggested measurements of the in-plane magnetic penetration depth in very lightly doped cuprates. Here we should add that higher order $`1/N`$ corrections are expected to renormalize downwards the fundamental energy scale $`D_s(0)`$ but, nevertheless, leave the lineshape of the scaled curve $`D_s(T)/D_s(0)`$ vs $`k_BT/D_s(0)`$ essentially intact.

The present study shares some common features with earlier works invoking a phase-fluctuation mechanism for the high-$`T_c`$ superconductivity. On the other hand, all these works including ours are radically different from phenomenological approaches that implicate the thermally excited nodal quasiparticles in a $`d`$-wave BCS superconducting state for the reduction of $`D_s(T)`$ with increasing temperature. The weak-coupling BCS type of approaches, however, are undermined by the absence of normal state quasiparticle peaks near the Brillouin zone points $`(0,\pm \pi )`$ and $`(\pm \pi ,0)`$ where superconductivity is supposed to originate. At present it is still difficult to discern experimentally whether the temperature and doping dependence of $`D_s(T)`$ is dominated by phase fluctuations or by nodal BCS-like quasiparticle excitations. However, recent experiments in cuprate thin films have provided strong evidence for the inherent two-dimensional character of superconductivity and for the KT nature of the superconducting transition.

## Acknowledgments

I would like to thank C. Panagopoulos for providing the penetration depth data and X. Zotos and G. Varelogiannis for stimulating discussions.

## A Generalized XY model

In the main body of the paper we noted that the structure of the phase fluctuations of the $`t`$-$`t^{}`$-$`J`$ model under study reduces, in the half-filled-band limit, to that of a classical XY model with only next-nearest-neighbor (n.n.n) interactions. In order to clarify the properties of the latter model, in juxtaposition to those of the more conventional nearest-neighbor (n.n) XY model, we consider briefly in this appendix the following Hamiltonian

$$H_{\mathrm{XY}}=\frac{1}{2}\underset{i,j}{}I_{ij}\mathrm{cos}(\psi _i\psi _j),$$
(A1)
assuming

$$I_{ij}=\{\begin{array}{cc}(1\alpha )I& \text{if }i,j\text{ are nearest neighbors}\hfill \\ \alpha I& \text{if }i,j\text{ are next nearest neighbors}\hfill \\ 0& \text{otherwise},\hfill \end{array}$$
(A2)
where $`\alpha `$ is a free parameter, with $`0\alpha 1`$, and $`I>0`$. At each lattice site $`i`$ the angle $`\psi _i`$ varies in the interval: $`0\psi _i2\pi `$. Evidently, $`\alpha =0`$ corresponds to the n.n XY model, while $`\alpha =1`$ corresponds to the n.n.n XY model.

The superfluid weight (or helicity modulus) for the generalized XY model (A1) reads

$`D_s(T)`$ $`=`$ $`{\displaystyle \frac{2}{z\mathrm{\Lambda }}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}I_{ij}|𝐑_i𝐑_j|^2\mathrm{cos}(\psi _i\psi _j)`$ (A4)
$`{\displaystyle \frac{\beta }{\mathrm{\Lambda }}}\left\{{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}I_{ij}[(𝐑_i𝐑_j)𝐞_x]\mathrm{sin}(\psi _i\psi _j)\right\}^2,`$
in agreement with corresponding early results. In view of (A2), the ground state configuration of (A1) is simply given by $`\psi _i=0`$, while the zero-temperature value of the superfluid weight follows immediately from (LABEL:eq:xy-D(T)) as

$$D_s(0)=(1+\alpha )I.$$
(A6)
Integrating the quadratic (Gaussian) fluctuations around the ground state configuration we obtain, after some lengthy algebra, the following low-temperature asymptotic expansion for the superfluid weight

$$\frac{D_s(T)}{D_s(0)}=1[1+G(\alpha )]\frac{k_BT}{zD_s(0)}\text{for}T0.$$
(A7)
Here $`G(\alpha )`$ is a dimensionless geometric factor given by

$$G(\alpha )=\frac{1}{\mathrm{\Lambda }}\underset{𝐪}{}\frac{\alpha (1\delta _𝐪)}{(1\alpha )(1\gamma _𝐪)+\alpha (1\delta _𝐪)},$$
(A8)
with

$$\gamma _𝐪=\frac{1}{2}(\mathrm{cos}q_x+\mathrm{cos}q_y),\delta _𝐪=\mathrm{cos}q_x\mathrm{cos}q_y.$$
(A9)
We emphasize that $`G(\alpha )`$ is an increasing function of $`\alpha `$ with end-point values: $`G(0)=0`$ and $`G(1)=1`$. Hence, from (A7) follows that the zero-temperature slope $`D_s^{}(0)`$ evolves monotonically from $`k_B/4`$ to $`k_B/2`$, as the parameter $`\alpha `$ varies from 0 (n.n XY model) to 1 (n.n.n XY model); see dash-dotted and dotted lines in Fig. 3.

Noting now that $`|𝐑_i𝐑_j|^2=2`$, for n.n.n sites $`i,j`$, whereas $`|𝐑_i𝐑_j|^2=1`$, for n.n sites $`i,j`$, a cautious inspection of (LABEL:eq:xy-D(T)) reveals that, in the thermodynamic limit and at each given temperature $`T`$, the superfluid weight of the n.n.n XY model should be twice as large the superfluid weight of the n.n XY model: $`D_s^{(\alpha =1)}(T)=2D_s^{(\alpha =0)}(T)`$. The latter property is explicit in the low-temperature analytic results (A6)–(A7) and we have confirmed its validity in the whole temperature range by Monte Carlo simulation using the standard Metropolis algorithm with the parameters quoted in Sec. III. Indeed, as shown in Fig. 3, the Monte Carlo data curves of $`D_s(T)/D_s(0)`$ vs $`k_BT/D_s(0)`$ for the n.n.n XY model coincide, within numerical error, with those for the n.n XY model, when the $`k_BT/D_s(0)`$ axis is scaled by a factor of 2. This agreement becomes better with increasing lattice size.

For the conventional n.n XY model it is well-known that in the thermodynamic limit we have the jump ratio: $`D_s(T_c)/(k_BT_c)=2/\pi `$, or equivalently in the notation of (25), the value $`R(T_c)=1`$. Hence, in view of the simple relation $`D_s^{(\alpha =1)}(T)=2D_s^{(\alpha =0)}(T)`$, we anticipate for the n.n.n XY model the jump ratio: $`D_s(T_c)/(k_BT_c)=4/\pi `$, or equivalently, $`R(T_c)=2`$. In order to demonstrate numerically the latter property we have carried out a finite-size scaling analysis for the n.n.n XY model based on the scaling formula (25), as described in detail Sec. III. The results are summarized in Fig. 4. We note that the minimum of the $`\chi ^2(T)`$ curve in Fig. 4(a) occurs, within numerical error, at the same point were the $`R(T)`$ curve in Fig. 4(b) crosses the line $`R(T)=2`$. Hence the assignment $`R(T_c)=2`$ leads, indeed, to a uniquely determined value for $`T_c`$ which from the crossing point in Fig. 4(b) is estimated to be $`k_BT_c/D_s(0)=0.4435(5)`$. The success of the finite-size scaling analysis validates then the assignment $`R(T_c)=2`$ for the n.n.n XY model.

It is worth emphasizing that in terms of the dimensionless parameter $`A=k_BT_c/D_s(0)`$, the Monte Carlo results of the present appendix imply the value $`A^{(\alpha =1)}=0.4435(5)`$ for the n.n.n XY model and, of course, twice as large corresponding value for the n.n XY model, i.e., $`A^{(\alpha =0)}=0.887(1)`$. Note that the latter value for the n.n XY model agrees with the original estimate of Weber and Minnhagen derived with the same finite-size scaling procedure. Pictorially, our results are manifest in Fig. 3 were the crossings of the short (long) dashed line with the Monte Carlo data curves for the n.n.n (n.n) XY model provide estimates for the value of $`A`$ of increasing accuracy, as the lattice size increases.


format: text

text:
# Event Rates for Binary Inspiral

## Introduction

Compact objects, neutron stars (NS) or black holes (BH), formed from relatively massive stars can spiral in and coalesce when found in tight binaries, the orbital evolution of which is driven by gravitational radiation. As angular momentum losses dominate, the orbit shrinks and the two compact objects can eventually merge as they revolve in orbit around each other. The prototype progenitor system of such inspiral events is the binary pulsar PSR B1913+16 (the “Hulse–Taylor” pulsar HT75 ). Sensitive pulsar timing measurements have revealed that the orbital period decreases at a rate comparable (to better than 1%) to that predicted by general relativity for the emission of gravitational waves TW82 , TW89 . The ultimate coalescence of the two neutron stars seems inevitable.

Although PSR B1913+16 will not reach coalescence for another 300 Myr, similar inspiraling systems in the Milky Way and nearby galaxies are thought to be primary sources of gravitational radiation for ground–based interferometric gravitational–wave detectors, currently under construction or commissioning (e.g., LIGO, VIRGO, GEO600). In addition to NS–NS close binaries, NS–BH and BH–BH binaries are also expected to form through the evolution of massive binaries and to contribute to the detection of inspiral events.

The expected detection rate of inspiral events depends on (i) the strength of the expected gravitational–wave signal, (ii) the gravitational–wave detector sensitivity, and (iii) the coalescence rate of each binary population. The first two considerations define a maximum distance $`D_{\mathrm{max}}`$, out to which different types of inspiral events and mergers could be detected. The coalescence rate for each population is estimated in two steps: first, the Galactic rate, and then its extrapolation out to the maximum distance of interest. Based on the current understanding of the LIGO sensitivities, the maximum distances out to which inspiral events could be detected by LIGO II (and LIGO I) are (approximately), 350 Mpc (20 Mpc) for NS–NS binaries, 700 Mpc (40 Mpc) for NS–BH binaries, and 1500 Mpc (100 Mpc) for BH–BH binaries (assuming 1.4 M NS and 10 M BH; Sam Finn, private communication). Given our current best knowledge (based on recent redshift surveys) of galaxy distributions out to those distances K00 , it can be estimated that, for a LIGO II detection rate of 1 event per year, the following Galactic coalescence rates are required: $`5\times 10^7`$ yr-1 for NS–NS binaries, $`5\times 10^8`$ yr-1 for NS–BH binaries, and $`5\times 10^9`$ yr-1 for BH–BH binaries.

Formation rates of coalescing compact binaries (systems with tight enough orbits that merge within a Hubble time $`10^{10}`$ yr) have been calculated so far using two very different methods: either entirely theoretically, based on binary evolution models, or, for NS–NS binaries, empirically, based on the observed NS–NS sample. In what follows we present an up–to–date review of current rate estimates, addressing in detail the most important uncertainties associated with them. We also discuss independent ways of obtaining upper limits to the coalescence rate of NS–NS binaries and possible implications of the formation of systems without recycled pulsars.

## Theoretical Rate Estimates

The formation rate of coalescing binary compact objects can be calculated, given a sequence of evolutionary stages leading to binary compact object formation. Over the years, a relatively standard picture has been formed describing the birth of such systems based on considerations of NS–NS binaries BvdH91 . More recently, variations of the standard evolutionary channel have also been discussed and suggested B95 , mainly based on worries about the fate of neutron stars in situations of hypercritical accretion (not limited to the photon Eddington rate), and their possible collapse into black holes. In all versions, however, the main picture remains the same: the initial binary progenitor consists of two binary members massive enough to eventually collapse into a NS or a BH. The evolutionary path involves multiple phases of stable or unstable mass transfer, common–envelope phases (where one or possibly two stellar cores spiral in the envelopes of evolved stars and eventually lead to the ejection of these envelopes), and accretion onto compact objects, as well as two core collapse events. The final outcome of interest is the formation of binary compact objects in close binary orbits.

Such theoretical modeling has been undertaken by a number of different groups by means of population syntheses. This provides us with ab initio predictions of coalescence rates. Monte Carlo numerical techniques are employed in following the evolution of a large ensemble of primordial binaries with certain assumed initial properties through a multitude of channels until compact object binaries are formed. The changes in the properties of the binaries at the end of each stage are calculated based on our current understanding of the various evolutionary processes involved: wind mass loss from massive hydrogen– and helium–rich stars, mass and angular–momentum losses during mass transfer phases, dynamically unstable mass transfer and common–envelope evolution, effects of highly super–Eddington accretion onto NS, and supernova explosions with kicks imparted to newborn NS or even BH. Given our limited understanding of some of these phases, the results of population synthesis are expected to depend on the assumptions made in the treatment of the various processes. Therefore, exhaustive parameter studies are required by the nature of the problem.

Recent studies of the formation of compact objects and calculations of their Galactic coalescence rates (LPP97 , FBB98 , PZ98 , BB98 , FWH99 , BBZ99  ) have explored the input parameter space and the robustness of the results at different levels of (in)completeness. Almost all of these groups have studied the sensitivity of the predicted coalescence rates to the average magnitude of the kicks imparted to compact objects at birth. The range of predicted NS–NS Galactic rates obtained by varying the kick magnitude alone is found in the range $`<10^75\times 10^4`$ yr-1. This large range indicates the importance of supernovae (two in this case) in the evolution of massive binaries. Variations in the assumed mass–ratio distribution for the primordial binaries can further change the predicted rate by about a factor of $`10`$, while assumptions of the common–envelope phase add another factor of about $`10100`$. Variation in other parameters typically affects the results by factors of two or less. Predicted rates for BH–NS and BH–BH binaries lie in the ranges $`<10^710^4`$ yr-1 and $`<10^710^5`$ yr-1, respectively when the kick magnitude to both NS and BH is varied. Other uncertain factors such as the critical progenitor mass for NS and BH formation lead to variations of the rates by factors of $`1050`$.

It is evident that recent theoretical predictions for coalescence rates cover a wide range of values (typically 3–4 orders of magnitude), because the various input parameters and assumptions affect strongly the absolute normalization (birth rate) of the modeled populations. Given these results, it seems fair to say that, at least at present, population synthesis calculations have a rather limited predictive power and provide fairly loose constraints on coalescence rates. One way to improve the reliability of such predictions is to study a number of different binary populations (with or without compact objects) and incorporate a number of independent observational constraints, such as star formation rate, supernova rates of different types, binarity of Wolf–Rayet stars, and others. A number of constraints on the population synthesis models could help restricting the predicted coalescence rates in narrower ranges BKB00 .

## Empirical Rate Estimates

The large range of theoretically predicted Galactic coalescence rates of double compact objects motivates us to examine other ways of obtaining rate estimates. The observed sample of coalescing NS–NS binaries found in the Galactic field (PSR B1913+16 and PSR B1534+12) provides us with alternative estimates of their coalescence rate. “Empirical” estimates can be obtained using the observed pulsar and binary properties along with models of selection effects in radio pulsar surveys P91 , N91 . For each observed object, a scale factor can be calculated based on the fraction of the Galactic volume within which pulsars with properties identical to those of the observed pulsar could be detected by any of the radio pulsar surveys, given their detection thresholds. This scale factor is a measure of how many more pulsars like the ones detected in the coalescing NS–NS systems exist in our galaxy. The coalescence rate can then be calculated based on the scale factors and estimates of detection lifetimes summed up for all the observed systems. Based on this method the first two studies concluded that the NS–NS Galactic coalescence rate is $`10^6`$ yr-1.

Since then, estimates of the NS–NS coalescence rate have known a significant downward revision primarily because of (i) the increase of the Galactic volume covered by radio pulsar surveys with no additional coalescing NS–NS being discovered CL95 , (ii) the increase of the distance estimate for PSR B1534+12 based on measurements of post-Newtonian parameters S98  (iii) revisions of the lifetime estimates vdHL96 , A98 . Recent estimates place the NS–NS rate for our Galaxy in the range $`13\times 10^7`$ yr-1. Further, it has been realized that a number of upward correction factors must be included, most importantly to account (i) for the beamed nature of pulsar emission and correct for all the binary pulsars with beams that our line of sight does not intersect, and (ii) for the faint end of the pulsar luminosity function and correct for those systems that are too faint to be detected. These two correction (multiplication) factors have so far typically been assumed to be $`3`$ and $`10`$, respectively.

In a study just recently completed K00 , we especially focused on all the uncertainties associated with these empirical estimates. We found that the upward correction factor for the faint end of the pulsar luminosity is the most important source of uncertainty. However, it is highly sensitive to the number of observed objects and its distribution function widens dramatically for small–number samples. For a sample of two objects (as the observed one) the faint–pulsar correction factor can vary from very small (close to unity) to as high as $`200`$ (see following subsection). Beyond the issue of faint pulsars, we considered a number of uncertainties and correction factors. Based on recent observational data for both PSR B1913+16 and PSR B1534+12, we found that the beaming correction factor is higher than previously thought ($`6`$) but with a rather small uncertainty ($`10`$%). Other factors, such as pulsar ages and lifetimes, and spatial distribution, lead to an uncertainty factor of about 2. We estimate the Galactic NS–NS coalescence rate in the range $`10^65\times 10^4`$ yr-1, which is still narrow compared to the range covered by the theoretical estimates.

### Small Number Sample and Pulsar Luminosity Function

One important limitation of empirical estimates of the coalescence rates is that they are derived based on only two observed NS–NS systems, under the assumption that the observed sample is representative of the true population, particularly in terms of their radio luminosity. Assuming that the recycled pulsars in NS–NS binaries follow the radio luminosity function of young pulsars and that therefore their true Galactic population is dominated in number by low–luminosity pulsars, it can be shown that the current empirical estimates most probably underestimate the true coalescence rate. If a small–number sample is drawn from a parent population dominated by low–luminosity (hence hard to detect) objects, it is statistically more probable that the sample will actually be dominated by objects from the high–luminosity end of the population. The result is that the population overall is thought to be brighter than it really is, and therefore, detectable over a larger Galactic volume. Consequently, the empirical estimates based on such a sample will tend to overestimate the detection volume for each observed system, and therefore underestimate the scale factors and the resulting coalescence rate.

This effect can be clearly demonstrated with a Monte Carlo experiment K00  using simple models for the pulsar luminosity function and the survey selection effects. As a first step, the average observed number of pulsars is calculated given a known “true” total number of pulsars in the Galaxy (thick-solid line in Figure 1). As a second step, a large number of sets consisting of “observed” (simulated) pulsars are realized using Monte Carlo methods. These pulsars are drawn from a Poisson distribution of a given mean number ($`<N_{\mathrm{obs}}>`$) and have luminosities assigned according to the assumed luminosity function. Based on each of these sets, one can estimate the total number of pulsars in the Galaxy using empirical scale factors, as is done for the real observed sample. The many (simulated) “observed” samples can then be used to obtain the distribution of the estimated total Galactic numbers ($`N_{\mathrm{est}}`$) of pulsars. We find that these $`N_{\mathrm{est}}`$ distributions are very strongly skewed and lead to possible correction factors for the faint pulsars in a wide range of values (covering typically a couple of orders of magnitude). The median and 25% and 75% percentiles of this distribution are plotted as a function of the assumed number of systems in the (fake) “observed” samples in Figure 1 (thin–solid and dashed lines, respectively).

It is evident that, in the case of small–number observed samples (less than $`10`$ objects), the estimated total number, and hence the estimated coalescence rate, can be underestimated by a significant factor. For observed samples with an expected number of objects equal to two, for example, the true rate may be much higher by more than a factor as high as $`200`$. This underestimation factor represents an upward correction factor that must be applied to the rate estimated using the observed sample of coalescing NS–NS binaries. However, we note that distribution of this correction factor covers a wide range and becomes highly skewed for small number samples (less than about 10 objects), and therefore it is currently quite uncertain. We conclude that correcting for the undetected, faint pulsars in the population cannot be decoupled from the problems of a small–number sample because of the assumption of the observed sample being representative of the population, implicit in the method.

## Upper Limits on the NS–NS Coalescence Rate

Observations of NS–NS systems and isolated pulsars related to NS–NS formation allow us to obtain upper limits on their Galactic coalescence rate in a number of different ways. Depending on how their value compares to the Galactic rate required for a LIGO II detection rate of 1 event per year, such limits can in principle provide us with valuable information about the prospects of gravitational–wave detection.

The absence of any young pulsars detected in NS–NS systems was used to obtain a rough upper limit to the rate of $`10^5`$ yr-1 B96 . Recently the same basic argument was reexamined in more detail and a more robust upper limit of $`10^4`$ yr-1 was derived A98 .

An upper bound to the NS–NS coalescence rate can also be obtained by combining our theoretical understanding of orbital dynamics (for supernovae with NS kicks in binaries) with empirical estimates of the birth rates of other types of pulsars related to NS–NS formation KL00 . Progenitors of NS–NS systems experience two supernova explosions. The second supernova explosion (forming the NS that is not observed as a pulsar) provides a unique tool for the study of NS–NS formation, since the post–supernova evolution of the system is simple, driven only by gravitational–wave radiation. There are three possible outcomes after the second supernova: (i) a coalescing NS–NS is formed (CB), (ii) a wide NS–NS (with a coalescence time longer than the Hubble time) is formed (WB), or (iii) the binary is disrupted (D) and a single pulsar similar to the ones seen in NS–NS systems is ejected. Based on supernova orbital dynamics we can accurately calculate the probability branching ratios for these three outcomes, $`P_{\mathrm{CB}}`$, $`P_{\mathrm{WB}}`$, and $`P_\mathrm{D}`$. For a given kick magnitude, we can calculate the maximum ratio $`(P_{\mathrm{CB}}/P_\mathrm{D})^{\mathrm{max}}`$ for the complete range of pre-supernova parameters defined by the necessary constraint $`P_{\mathrm{CB}}0`$ (Figure 2). Given that the two types of systems have a common parent progenitor population, the ratio of probabilities is equal to the ratio of the birth rates $`(BR_{\mathrm{CB}}/BR_\mathrm{D})`$.

We can then use (i) the absolute maximum of the probability ratio ($`0.27`$ from Figure 2) and (ii) an empirical estimate of the birth rate of single pulsars similar to those in NS–NS based on the current observed sample to obtain an upper limit to the coalescence rate. The selection of this sample involves some subtleties KL00 , and the analysis results in $`BR_{\mathrm{CB}}<1.5\times 10^5`$ yr-1. Note that this number could be increased because of the small–number sample and luminosity bias, which this time affects the empirical estimate of $`BR_\mathrm{D}`$ by a factor of $`26`$. Such an upward correction can bring the upper limit in the range $`39\times 10^5`$ yr-1.

This is an example of how we can use observed systems other than NS–NS to improve our understanding of their coalescence rate. A similar calculation can be done using the wide NS–NS systems instead of the single pulsars KL00 .

## Non–recycled double neutron stars

We have already pointed out that the empirical methods employed to obtain rate estimates for NS–NS coalescence include the implicit assumption that the properties of the observed sample are representative of the Galactic NS–NS population. This assumption extends to the pulsar properties and their evolutionary history of recycling (spin–up by accretion). Consistent with the pulsars observed in the detected NS–NS systems, it turns out that so far theoretical studies of NS–NS formation have considered systems where one of the neutron stars had the opportunity to be recycled, at least in principle (through stellar winds, Roche–lobe overflow accretion, or even possibly in a common–envelope phase).

Here, we report on a new evolutionary path leading to the formation of close NS–NS binaries, with the unique characteristic that none of the two NS ever had the chance to be recycled by accretion. As we will discuss in more detail, such NS–NS systems have a negligible probability of being detected as binary pulsars, and could represent a “dormant” NS–NS population in galaxies with important implications for gravitational–wave detection of NS–NS inspiral events. The existence of this recently identified BK00  evolutionary channel stems from the evolution of helium–rich stars (cores of massive NS progenitors), which has been neglected in most previous studies of double compact object formation. We find that these non–recycled NS–NS binaries are formed from bare carbon–oxygen cores in tight orbits, with formation rates comparable to or maybe even higher than those of recycled NS–NS binaries.

### The Method

We study NS–NS binaries formed through a multitude of evolutionary sequences that are not predefined, but instead are realized in Monte Carlo population synthesis calculations.

To describe the evolution of single stars (hydrogen– and helium–rich) from the zero age main sequence (ZAMS) to carbon–oxygen (CO) core formation, we employ analytical formulae from stellar evolution fits H00  However, we have adopted a prescription for the masses of compact objects formed at core–collapse events, based on the relation between CO core masses and final FeNi core masses W86 .

Concerning the evolution of interacting binaries, we model the changes of mass and orbital parameters taking into account mass and angular momentum transfer between the stars or loss from the system during Roche–lobe overflow, tidal circularization, rejuvenation of stars due to mass accretion, wind mass loss from massive and/or evolved stars, dynamically unstable mass transfer episodes leading to common–envelope (CE) evolution and spiral–in of the stars. We also account for the possibility of hyper-critical accretion onto compact objects during CE phases B95  and effects of asymmetric supernovae (SN) on a binary orbit (mass loss and a kick velocity a newly born compact object receives in SN). More details about the treatment of various evolutionary processes are presented elsewhere BKB00 .

In the synthesis calculations, we typically evolve a few million of primordial binaries to satisfy the requirement that the statistical (Poisson) fractional errors ($`1/\sqrt{N}`$) of the final NS–NS population are lower than 10%. The formation rates are calibrated using the latest Type II SN empirical rates and normalized to our Galaxy CET99 .

In our standard model, primordial binaries follow given distributions: for primary masses ($`5100`$ M), $`M_1^{2.7}dM_1`$; for mass ratios ($`0<q<1`$), $`dq`$; for orbital separations (from a minimum, so both ZAMS stars fit within their Roche lobes, up to $`10^5`$ R), $`dA/A`$; for eccentricities, $`2e`$. Each of the models is also characterized by a set of assumptions, which, for our standard model, are: (1) Kick velocities. We use a weighted sum of two Maxwellian distributions with $`\sigma =175`$ km s-1 (80%) and $`\sigma =700`$ km s-1 (20%) CC97 ; (2) Maximum NS mass. We adopt a conservative value of $`M_{\mathrm{max}}=3`$ M KB96 ; (3) Common envelope efficiency. We assume $`\alpha _{\mathrm{CE}}\times \lambda =1.0`$, where $`\alpha `$ is the efficiency with which orbital energy is used to unbind the stellar envelope, and $`\lambda `$ is a measure of the central concentration of the giant; (4) Non–conservative mass transfer. In cases of dynamically stable mass transfer between non–degenerate stars, we allow for mass and angular momentum loss from the binary P92 , assuming that half of the mass lost from the donor is also lost from the system ($`1f_\mathrm{a}=0.5`$) with specific angular momentum equal to $`\beta 2\pi `$A2/P ($`\beta =1`$); (5) Star formation history. We assume that star formation has been continuous in the disk of our Galaxy for the last 10 Gyr G01 .

### Results

We use our population synthesis models to investigate all possible formation channels of NS–NS binaries realized in the simulations. We find that a significant fraction of coalescing NS–NS systems are formed through a new, previously not identified evolutionary path. The evolution along this new channel begins with two phases of Roche–lobe overflow. The first, from the primary to the secondary, involves non–conservative but dynamically stable mass transfer and ends when the hydrogen envelope is consumed. The second, from the initial secondary to the helium core of the initial primary, involves dynamically unstable mass transfer, i.e., CE evolution. The post–CE binary consists of two bare helium stars of relatively low masses. As they evolve through core and shell helium burning, the two stars develop ‘giant–like” structures, with clear CO cores and convective envelopes. Their radial expansion eventually brings them into contact and the system evolves through a double CE phase (similar to Brown , for hydrogen–rich stars). During this double CE phase, the combined helium envelopes are ejected at the expense of orbital energy. The tight, post–CE system consists of two CO cores, which eventually end their lives as Type Ic supernovae leaving double neutron star system.

The unique qualitative characteristic of this NS–NS formation path is that both NS have avoided recycling. Based on comparison of non–recycled NS–NS relative to that of recycled pulsars, for each of our models, we derive a correction factor for empirical estimates of the Galactic NS–NS coalescence rate. Since these estimates account only for NS–NS systems with recycled pulsars, they must be increased to include any non–recycled systems formed. We have performed an extensive parameter study to assure robustness of our results. In Table 1 we present the formation rates of non–recycled NS–NS binaries and the total NS–NS population with merger times shorter than 10 Gyr, along with the upwards correction factor for the Galactic empirical rate estimates. Results are shown only for models where the derived factor differs from our standard model by more than 25%. We find that these factors are typically $`1.53`$ but can be higher for some models.

We note that the identification of the formation path for non–recycled NS–NS binaries stems entirely from accounting for the evolution of helium stars and for the possibility of double CE phases, both of which have typically been ignored in previous calculations.

## Conclusions

The current theoretical estimates of NS–NS coalescence rates appear to have a rather limited predictive power. They cover a range of values in excess of 3 orders of magnitude. Most importantly, this range includes the value of $`5\times 10^7`$ yr-1 required for a LIGO II detection rate of 1 event per year. This means that at the two edges of the range the conclusion swings from no detection to many per month, and therefore the detection prospects of NS–NS coalescence cannot be assessed firmly. On the other hand, empirical estimates based on the observed sample of coalescing NS–NS systems appear to be more robust. Taking into account recent empirical estimates and the associated uncertainties K00 , we find the Galactic NS–NS inspiral rate in the range $`10^65\times 10^4`$ yr-1. If we also include the independently derived upper limit of $`10^4`$ yr-1, we expect a detection rate of $`2300`$ events per year for LIGO II.

It is important to note here that another implicit assumption in derivation of the empirical estimates is that all NS–NS binaries have at some point in their lifetime contained a recycled pulsar with rather long lifetimes ($`10^9`$ yr). However, recent models of NS–NS formation BK00  show that there may exist a significant NS–NS population with neutron that never had the chance to be recycled and therefore have very short lifetimes (by 2-3 orders of magnitude, thus preventing their detection). For a variety of population synthesis models, the birth rate of this separate population of coalescing NS–NS binaries is typically comparable or higher than that of the systems with one recycled NS. The total number of coalescing NS–NS systems could be higher by factors of at least 50%, and up to 10 or even higher. Such an increase has important implications for prospects of gravitational wave detection by ground–based interferometers. Using the recent results on the empirical NS–NS coalescence rate K00 , we find that the most optimistic prediction for the LIGO I detection rate could be raised to at least 1 event per 2–3 years, and the most pessimistic LIGO II detection rate could be raised to 3–6 events per year or even higher.

Estimates of the coalescence rate of BH–NS and BH–BH systems rely solely on our theoretical understanding of their formation. As in the case of NS–NS binaries, the model uncertainties are significant and the ranges extend to more than 2 orders of magnitude. However, the requirement on the Galactic rate is less stringent for 10 M BH–BH binaries, only $`5\times 10^9`$ yr-1. Therefore, even with the pessimistic estimates for BH–BH coalescence rates ($`10^7`$ yr-1), we would expect at least a few to several detections per year (with LIGO II), which is quite encouraging. We also point out that that a recent examination of formation of close BH-BH through dynamical processes (stellar interactions) in globular clusters leads to detection rates as high as a few per day for LIGO II and 1 event per 2 years for LIGO I PZ00 .


text:
# 0 Introduction

## 0 Introduction

The 3-dimensional incompressible Euler equation (“3D Euler”) is as follows:

$$\begin{array}{ccc}\left(\frac{}{t}+u_x\right)u=_xp\hfill & & \left(xR^3,t0\right)\hfill \\ & & \\ _xu=\mathrm{\hspace{0.17em}0}\hfill & & (xR^3,t0)\hfill \\ & & \\ u(x,0)=u^0(x)\hfill & & (xR^3)\hfill \end{array}$$

with $`u^0`$ a given, smooth, divergence-free, rapidly decreasing vector field on $`R^3`$. Here, $`u(x,t)`$ and $`p(x,t)`$ are the unknown velocity and pressure for an ideal, incompressible fluid flow at zero viscosity. An outstanding open problem is to determine whether a 3D Euler solution can develop a singularity at a finite time $`T`$. A classic result of Beale-Kato-Majda  asserts that, if a singularity forms at time $`T`$, then the vorticity $`\omega (x,t)=_x\times u(x,t)`$ grows so rapidly that

$$\underset{0}{\overset{T}{}}\underset{x}{sup}|\omega (x,t)|dt=\mathrm{}.$$

In , Constantin-Fefferman-Majda showed that, if the velocity remain bounded up to the time $`T`$ of singularity formation, then the vorticity direction $`\omega (x,t)/|\omega (x,t)|`$ cannot remain uniformly Lipschitz continuous up to time $`T`$.

One scenario for possible formation of a singularity in a 3D Euler solution is a constricting vortex tube. Recall that a vortex line in a fluid is an arc on an integral curve of the vorticity $`\omega (x,t)`$ for fixed $`t`$, and a vortex tube is a tubular neighborhood in $`R^3`$ arising as a union of vortex lines. In numerical simulations of 3D Euler solutions, one routinely sees that vortex tubes grow longer and thinner, while bending and twisting. If the thickness of a piece of a vortex tube becomes zero in finite time, then one has a singular solution of 3D Euler. It is not known whether this can happen.

Our purpose here is to adapt our work ,  on two-dimensional flows to three dimensions, for application to 3D Euler. We introduce below the notion of a “regular tube”. Under the mild assumption that

$$\underset{0}{\overset{T}{}}\underset{x}{sup}|u(x,t)|dt<\mathrm{},$$

we show that a regular tube cannot reach zero thickness at time $`T`$. In particular, for 3D Euler solutions, a vortex tube cannot reach zero thickness in finite time, unless it bends and twists so violently that no part of it forms a regular tube. This significantly sharpens the conclusion of  for possible singularities of 3D Euler solutions arising from vortex tubes. On the other hand,  applies to arbitrary singularities of 3D Euler solutions, while our results apply to “regular tubes”.

Although we are mainly interested in 3D Euler solutions, our result is stated for arbitrary incompressible flows in 3 dimensions. The proof is simple and elementary. The main novelty for readers familiar with ,  is that we can adapt the ideas of  to three dimensions, even though there is no scalar that plays the rôle of the stream function on $`R^2`$.

## 1  Regular Tubes

Let $`Q=I_1\times I_2\times I_3R^3`$ be a closed rectangular box (with $`I_j`$ a bounded interval), and let $`T>0`$ be given.

A regular tube is an open set $`\mathrm{\Omega }_tQ`$ parametrized by time $`t[0,T)`$, having the form

$$\mathrm{\Omega }_t=\{(x_1,x_2,x_3)Q:\theta (x_1,x_2,x_3,t)<\mathrm{\hspace{0.17em}0}\}$$
(1)
with

$$\theta C^1(Q\times [0,T)),$$
(2)
and satisfying the following properties:

$$|_{x_1,x_2}\theta |0\mathrm{for}(x_1,x_2,x_3,t)Q\times [0,T),\theta (x_1,x_2,x_3,t)=\mathrm{\hspace{0.17em}0};$$
(3)
$$\mathrm{\Omega }_t(x_3):=\{(x_1,x_2)I_1\times I_2:(x_1,x_2,x_3)\mathrm{\Omega }_t\}\mathrm{is}\mathrm{non}\mathrm{empty},$$
(4)
for all $`x_3I_3,t[0,T)`$;

$$\mathrm{closure}\left(\mathrm{\Omega }_t(x_3)\right)\mathrm{interior}(I_1\times I_2)$$
(5)
for all $`x_3I_3`$, $`t[0,T)`$.

For example, a thin tubular neighborhood of a curve $`\mathrm{\Gamma }`$ forms a regular tube, provided the tangent vector $`\mathrm{\Gamma }^{}`$ stays transverse to the $`(x_1,x_2)`$ plane.

Let $`u(x,t)=(u_k(x,t))_{1k3}`$ be a $`C^1`$ velocity field defined on $`Q\times [0,T)`$. We say that the regular tube $`\mathrm{\Omega }_t`$ moves with the velocity field $`u`$, if we have

$$\left(\frac{}{t}+u_x\right)\theta =0\mathrm{whenever}(x,t)Q\times [0,T),\theta (x,t)=0.$$
(6)
It is well-known that a vortex tube arising from a 3D Euler solution moves with the fluid velocity.

## 2 Statement of the Main Result

Theorem: Let $`\mathrm{\Omega }_tQ(t[0,T))`$ be a regular tube that moves with a $`C^1`$, divergence free velocity field $`u(x,t)`$.

If

$$\underset{0}{\overset{T}{}}\underset{xQ}{sup}|u(x,t)|dt<\mathrm{}$$
(7)
then

$$\underset{tT}{lim\; inf}\mathrm{Vol}(\mathrm{\Omega }_t)>\mathrm{\hspace{0.17em}0}.$$
(8)
## 3 Calculus Formulas for Regular Tubes

Let $`\mathrm{\Omega }_t`$ be a regular tube, as in (1)$`\mathrm{}`$(5). Recall that

$$\mathrm{\Omega }_t(x_3)=\{(x_1,x_2)I_1\times I_2:\theta (x_1,x_2,x_3,t)<\mathrm{\hspace{0.17em}0}\}.$$
(9)
Define also

$$S_t(x_3)=\{(x_1,x_2)\mathrm{interior}(I_1\times I_2):\theta (x_1,x_2,x_3,t)=0\}\mathrm{for}x_3I_3,t[0,T).$$
(10)
Also, for intervals $`II_3`$, and for $`t[0,T)`$, define

$$\mathrm{\Omega }_t(I)=\{(x_1,x_2,x_3)Q:x_3I\mathrm{and}\theta (x_1,x_2,x_3,t)<\mathrm{\hspace{0.17em}0}\},\mathrm{and}$$
(11)
$$S_t(I)=\{(x_1,x_2,x_3)Q:x_3I\mathrm{and}(x_1,x_2)S_t(x_3)\}.$$
(12)
Let $`\nu `$ denote the outward-pointing unit normal to $`S_t(I_3)`$, and let $`\stackrel{~}{\nu }=(\stackrel{~}{\nu }_1,\stackrel{~}{\nu }_2,0)`$, where $`(\stackrel{~}{\nu }_1,\stackrel{~}{\nu }_2)`$ is the outward-pointing unit normal to $`S_t(x_3)`$. Thus, $`\nu `$ and $`\stackrel{~}{\nu }`$ are continuous vector-valued functions, defined on $`𝒮=\{(x_1,x_2,x_3,t)Q\times [0,T):(x_1,x_2,x_3)S_t(I_3)\}`$. Define also scalar-valued functions $`\sigma ,\stackrel{~}{\sigma }`$ on $`𝒮`$ by requiring that

$$\left(\frac{}{t}+\sigma \nu _x\right)\theta =\left(\frac{}{t}+\stackrel{~}{\sigma }\stackrel{~}{\nu }_x\right)\theta =\mathrm{\hspace{0.17em}0}\mathrm{for}xS_t(I_3).$$
(13)
Again, $`\sigma `$ and $`\stackrel{~}{\sigma }`$ are well-defined and continuous on $`𝒮`$, thanks to (3). Let $`F`$ be any continuous function on $`Q`$. We will establish the following elementary formulas.

$$\frac{d}{dt}\left[\underset{\mathrm{\Omega }_t(x_3)}{}F(x_1,x_2,x_3)d(\mathrm{Area})\right]=\underset{(x_1,x_2)S_t(x_3)}{}F(x_1,x_2,x_3)\stackrel{~}{\sigma }(x_1,x_2,x_3,t)d(\mathrm{length})$$
(14)
$$\underset{S_t(I)}{}Fd(\mathrm{Area})=\underset{x_3I}{}\left\{\underset{S_t(x_3)}{}\frac{F}{\nu \stackrel{~}{\nu }}d(\mathrm{length})\right\}𝑑x_3.$$
(15)
To check (14) and (15), we may use a partition of unity to reduce to the case in which $`F`$ is supported in a small neighborhood $`U`$. Also, we may restrict attention to a small time interval $`J`$. In a small enough $`U\times J`$, we may assume that $`S_t(I_3)`$ is given by the graph of a $`C^1`$ function $`\psi `$, thanks to (3). Thus, without loss of generality, we may suppose that, in $`U\times J`$, we have

$$(x_1,x_2,x_3)\mathrm{\Omega }_t(I_3)\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}x_1<\psi (x_2,x_3,t),\mathrm{and}$$
(16)
$$(x_1,x_2,x_3)S_t(I_3)\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}x_1=\psi (x_2,x_3,t).$$
(17)
From (16), (17) we have also

$$(x_1,x_2)\mathrm{\Omega }_t(x_3)\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}x_1<\psi (x_2,x_3,t),\mathrm{and}$$
(18)
$$(x_1,x_2)S_t(x_3)\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}x_1=\psi (x_2,x_3,t).$$
(19)
In view of (16) $`\mathrm{}`$ (19), we have

$$\nu =(1,\frac{}{x_2}\psi ,\frac{\psi }{x_3})/\sqrt{1+\left(\frac{}{x_2}\psi \right)^2+\left(\frac{}{x_3}\psi \right)^2}\mathrm{and}$$
(20)
$$\stackrel{~}{\nu }=(1,\frac{}{x_2}\psi ,\mathrm{\hspace{0.17em}0})/\sqrt{1+\left(\frac{}{x_2}\psi \right)^2}\mathrm{on}𝒮(U\times J).$$
(21)
From (10), (17), we have $`\theta (\psi (x_2,x_3,t),x_2,x_3,t)=0`$ on $`𝒮(U\times J)`$. Differentiating in $`t`$, we obtain

$$\left(\frac{}{t}+(\frac{\psi }{t}(x_2,x_3,t),0,0)_x\right)\theta =0\mathrm{on}𝒮(U\times J).$$

Subtracting this from (13), we find that $`\left[(\frac{\psi }{t}(x_2,x_3,t),0,0)\sigma \nu \right]`$ is orthogonal to $`_x\theta `$, hence also to $`\nu `$. This yields the formula

$$\sigma =(\frac{\psi }{t}(x_2,x_3,t),0,0)\nu =\left(\frac{\psi }{t}\right)/\sqrt{1+\left(\frac{\psi }{x_2}\right)^2+\left(\frac{\psi }{x_3}\right)^2}.$$
(22)
Similarly, subtracting the two equations (13), we learn that $`(\sigma \nu \stackrel{~}{\sigma }\stackrel{~}{\nu })`$ is orthogonal to $`_x\theta `$, hence also to $`\nu `$. Therefore,

$$\sigma =\stackrel{~}{\sigma }(\nu \stackrel{~}{\nu })\mathrm{on}𝒮.$$
(23)
From (20) $`\mathrm{}`$ (23), we obtain

$$\stackrel{~}{\sigma }=\left(\frac{\psi }{t}\right)/\sqrt{1+\left(\frac{\psi }{x_2}\right)^2}\mathrm{on}𝒮(U\times J).$$
(24)
Now we can read off (14) and (15). In fact, (18) gives

$$\frac{d}{dt}\left[\underset{\mathrm{\Omega }_t(x_3)}{}Fd(\mathrm{Area})\right]=\frac{d}{dt}\left[_{_{_{x_1<\psi (x_2,x_3,t)}}}F𝑑x_1𝑑x_2\right]=F(\psi ,x_2,x_3)\frac{\psi }{t}𝑑x_2,$$

and (19), (24) yield

$$\underset{S_t(x_3)}{}F\stackrel{~}{\sigma }d(\mathrm{length})=F\stackrel{~}{\sigma }\sqrt{1+\left(\frac{\psi }{x_2}\right)^2}𝑑x_2=F(\psi ,x_2,x_3)\frac{\psi }{t}𝑑x_2,$$

proving (14). Similarly, (17) gives

$$\underset{S_t(I)}{}Fd(\mathrm{Area})=F(\psi ,x_2,x_3)\sqrt{1+\left(\frac{}{x_2}\psi \right)^2+\left(\frac{}{x_3}\psi \right)^2}𝑑x_2𝑑x_3$$

while (20) and (21) imply

$$\nu \stackrel{~}{\nu }=\sqrt{1+\left(\frac{}{x_2}\psi \right)^2}/\sqrt{1+\left(\frac{}{x_2}\psi \right)^2+\left(\frac{}{x_3}\psi \right)^2},$$

so that (19) yields

$$\underset{x_3I}{}\left\{\underset{S_t(x_3)}{}\frac{F}{\nu \stackrel{~}{\nu }}d(\mathrm{length})\right\}𝑑x_3=\underset{x_3I}{}\left\{\left(\frac{F}{\nu \stackrel{~}{\nu }}\right)\sqrt{1+\left(\frac{}{x_2}\psi \right)^2}𝑑x_2\right\}𝑑x_3$$

$$=F(\psi ,x_2,x_3)\sqrt{1+\left(\frac{}{x_2}\psi \right)^2+\left(\frac{}{x_3}\psi \right)^2}𝑑x_2𝑑x_3.$$

This completes the proof of (15). Note that (23) allows us to rewrite (15) in the form

$$\underset{S_t(I)}{}Fd(\mathrm{Area})=\underset{x_3I}{}\left\{\underset{S_t(x_3)}{}F\frac{\stackrel{~}{\sigma }}{\sigma }d(\mathrm{length})\right\}𝑑x_3.$$
(25)
## 4 Proof of the Theorem

We retain the notation of the previous sections. We will define a time-dependent interval

$$J_t=[A(t),B(t)]I_3$$
(26)
and establish an obvious formula for the time derivative of Vol $`\mathrm{\Omega }_t(J_t)`$. We assume that the endpoints $`A(t)`$, $`B(t)`$ are $`C^1`$ functions of $`t`$. We have

$$\mathrm{Vol}\mathrm{\Omega }_t(J_t)=\underset{x_3J_t}{}\mathrm{Area}\mathrm{\Omega }_t(x_3)𝑑x_3,\mathrm{so}\mathrm{that}$$

$$\frac{d}{dt}\mathrm{Vol}\mathrm{\Omega }_t(J_t)=B^{}(t)\mathrm{Area}\mathrm{\Omega }_t(B(t))A^{}(t)\mathrm{Area}\mathrm{\Omega }_t(A(t))+\underset{x_3J_t}{}\frac{}{t}\mathrm{Area}\mathrm{\Omega }_t(x_3)𝑑x_3.$$

Applying (14) with $`F1`$, we find that

$$\frac{d}{dt}\mathrm{Vol}\mathrm{\Omega }_t(J_t)=B^{}(t)\mathrm{Area}\mathrm{\Omega }_t(B(t))A^{}(t)\mathrm{Area}\mathrm{\Omega }_t(A(t))+\underset{x_3J_t}{}\left\{\underset{S_t(x_3)}{}\stackrel{~}{\sigma }d(\mathrm{length})\right\}𝑑x_3.$$

In view of (25) (with $`F\sigma `$ on $`𝒮`$), this is equivalent to

$$\frac{d}{dt}\mathrm{Vol}\mathrm{\Omega }_t(J_t)=B^{}(t)\mathrm{Area}\mathrm{\Omega }_t(B(t))A^{}(t)\mathrm{Area}\mathrm{\Omega }_t(A(t))+\underset{S_t(J_t)}{}\sigma d(\mathrm{Area}).$$
(27)
Now we bring in the hypothesis that $`\mathrm{\Omega }_t`$ moves with a divergence-free $`C^1`$ velocity field $`u`$. From (6) and (13), we see that $`(\sigma \nu u)_x\theta =0`$ on $`S_t(J_t)`$. Thus $`(\sigma \nu u)`$ is orthogonal to $`\nu `$, so that $`\sigma =u\nu `$ on $`S_t(J_t)`$, and (27) may be rewritten as

$$\frac{d}{dt}\mathrm{Vol}\mathrm{\Omega }_t(J_t)=B^{}(t)\mathrm{Area}\mathrm{\Omega }_t(B(t))A^{}(t)\mathrm{Area}\mathrm{\Omega }_t(A(t))+\underset{S_t(J_t)}{}u\nu d(\mathrm{Area}).$$
(28)
On the other hand, since $`u=(u_1,u_2,u_3)`$ is divergence-free, the divergence theorem yields

$$0=\underset{\mathrm{\Omega }_t(J_t)}{}(_xu)d(\mathrm{Vol})=\underset{S_t(J_t)}{}u\nu d(\mathrm{Area})+\underset{\mathrm{\Omega }_t(B(t))}{}u_3d(\mathrm{Area})\underset{\mathrm{\Omega }_t(A(t))}{}u_3d(\mathrm{Area}).$$

Hence, (28) may be rewritten in the form

$$\frac{d}{dt}\mathrm{Vol}\mathrm{\Omega }_t(J_t)=\underset{\mathrm{\Omega }_t(B(t))}{}[B^{}(t)u_3(x,t)]d(\mathrm{Area})\underset{\mathrm{\Omega }_t(A(t))}{}[A^{}(t)u_3(x,t)]d(\mathrm{Area}).$$
(29)
This is our final formula for the time derivative of Vol $`\mathrm{\Omega }_t(J_t)`$. It is intuitively clear.

We now pick the time-dependent interval $`J_t=[A(t),B(t)]I_3`$. Let $`I_3=[a,b]`$, and let $`t_0(0,T)`$ be a time to be picked below. We define

$$B(t)=b\underset{t}{\overset{T}{}}\underset{xQ}{\mathrm{max}}|u(x,\tau )|d\tau ,$$
(30)
and

$$A(t)=a+\underset{t}{\overset{T}{}}\underset{xQ}{\mathrm{max}}|u(x,\tau )|d\tau .$$
(31)
We are assuming that $`u(x,\tau )`$ is continuous on $`Q\times [0,T)`$, and that $`\underset{0}{\overset{T}{}}\underset{xQ}{\mathrm{max}}|u(x,\tau )|d\tau <\mathrm{}`$. It follows that $`A(t),B(t)`$ are $`C^1`$ functions on $`[0,T)`$, and that

$$aA(t)<B(t)b\mathrm{for}t[t_0,T),$$
(32)
provided we pick $`t_0`$ close enough to $`T`$. We pick $`t_0`$ so that (32) holds. Thus, $`\mathrm{\Omega }_t(J_t)Q`$ for $`t[t_0,T)`$. Immediately from (30), (31), we obtain

$$B^{}(t)=A^{}(t)=\underset{xQ}{\mathrm{max}}|u(x,t)|\underset{x\mathrm{\Omega }_t(A(t))\mathrm{\Omega }_t(B(t))}{\mathrm{max}}|u_3(x,t)|$$
(33)
(recall $`u=(u_1,u_2,u_3))`$.

From (29) and (33) we see at once that

$$\frac{d}{dt}\mathrm{Vol}\mathrm{\Omega }_t(J_t)\mathrm{\hspace{0.17em}0}\mathrm{for}t[t_0,T).$$
(34)
On the other hand, (4) and (32) show that $`\mathrm{Vol}\mathrm{\Omega }_{t_0}(J_{t_0})>\mathrm{\hspace{0.17em}0}`$.

Consequently,

$$\underset{tT}{lim\; inf}\mathrm{Vol}\mathrm{\Omega }_t\underset{tT}{lim\; inf}\mathrm{Vol}\mathrm{\Omega }_t(J_t)\mathrm{Vol}\mathrm{\Omega }_{t_0}(J_{t_0})>\mathrm{\hspace{0.17em}0}.$$

The proof of our theorem is complete. $`\mathrm{}`$

## 5  References

1. J. Beale, T. Kato, and A. Majda, “Remarks on the breakdown of smooth solutions for the 3D Euler equations,” Comm. Math. Phys., 94, pp.61-64, (1984).
2. P. Constantin, C. Fefferman, and A. Majda, “Geometric constraints on potentially singular solutions for the 3D Euler equations,” Commun. Part. Diff. Eq., 21, pp.559-571, (1996).
3. D. Cordoba and C. Fefferman, “Scalars convected by a 2D incompressible flow,” (to appear).
4. D. Cordoba and C. Fefferman, “Behavior of several 2D fluid equations in singular scenarios,” submitted to Proc. Nat. Acad. Sci., U.S.A.