text:
# 𝑟-modes in Relativistic Superfluid Stars

## I Introduction

The $`r`$-mode instability A98 ; FM98  belongs to the secular instability attributable to the so-called CFS mechanism C70 ; JF78 ; FS78  that drives unstable various modes of oscillation in a rotating star. The $`r`$-modes excited in neutron stars are expected to work as a mechanism of decelerating the spin velocity of neutron stars by emitting gravitational waves, which will be detectable by LIGO and other detectors LOM98 ; AKS99 ; OEL98 ; LB98 ; AKSt99 . For a recent review, see, e.g., AK01 .

In old and cold neutron stars having the interior temperature below the superfluid transition temperatures $`T10^9\mathrm{K}`$, neutrons in the inner crust and neutrons and protons in the core are believed to be in superfluid states ST83 . It is therefore likely that most of observed neutron stars, which have cooled down very quickly via emission of neutrinos BP79 , have a core with superfluids. One of the important effects of rotating superfluids in cold neutron stars is mutual friction, which provides a strong dissipation mechanism for the oscillation modes mutual . The mutual friction is produced by scattering of normal fluid particles off the vortices in the rotating superfluids. Lindblom and Mendell LM00  calculated the mutual friction damping associated with the $`r`$-modes in neutron stars with a superfluid core and concluded that the mutual friction in the core is ineffective to suppress the $`r`$-mode instability. Recently, Lee and Yoshida LY03  have studied the $`r`$-modes in neutron stars with the superfluid core using a different numerical method, confirming most of the results obtained by LM00 . Lee and Yoshida LY03  have also shown that the $`r`$-modes are split into ordinary fluid-like $`r`$-modes and superfluid-like $`r`$-modes, as first suggested by Andersson and Comer AC01 , and that the instability caused by the superfluid-like $`r`$-modes is extremely weak and easily damped by dissipation processes in the interior. Very recently, Yoshida and Lee YL03  have calculated inertial modes of superfluid neutron stars, and shown that the inertial modes are also split into two families, namely, ordinary fluid-like inertial modes and superfluid-like inertial modes. It is quite reasonable since the $`r`$-modes belong to a subclass of the inertial modes. Yoshida and Lee YL03  found that all the inertial modes, except for the ordinary fluid-like $`r`$-modes, are strongly damped by the dissipation due to the mutual friction unless the entrainment effects between superfluids are extremely weak.

Most of the numerical studies on the $`r`$-modes and the inertial modes in superfluid neutron stars have been done within the framework of Newtonian dynamics. Since neutron stars are general relativistic objects having the typical relativistic factor $`M(R)/R`$ of order of $`10^1`$ with $`M(R)`$ and $`R`$ being the gravitational mass and radius of the star, respectively, it is necessary to extend the Newtonian analyses to general relativistic ones. In this paper, employing the two-constituent model developed by Cater and his co-workers C89 ; CL98 ; LSC98 ; CLL99 , and the formalism devised by Lockitch, Friedman, and Andersson LAF01 ; LFA02  for inertial modes of relativistic normal fluid stars, we derive a general relativistic formulation for inertial mode oscillations in superfluid neutron stars. In Sec. II we present the basic equations employed in this paper for the dynamics of relativistic superfluids, and in Sec. III numerical results are given, and Sec. IV is for discussions and conclusions. In this paper, we employ geometric units, given by $`c=G=1`$, where $`c`$ and $`G`$ are the speed of light and the gravitational constant, respectively, and sign conventions used in MTW .

## II Formulation

### II.1 Two-constituent formalism for superfluid dynamics

We assume neutron stars are composed of superfluid neutrons and a mixture of superfluid protons and normal fluid electrons. We assume perfect charge neutrality between the protons and the electrons because the plasma frequency of the mixture is much higher than the oscillation frequency considered in this study GM91 . The electrons therefore co-move with the protons. In other words, we do not need equations for describing the dynamics of the electrons. In the following, we loosely call the mixture of the protons and electrons “proton”. To describe the dynamics of the superfluids in a neutron star, we employ the relativistic two-constituent formalism, which has been developed by Carter and his co-workers C89 ; CL98 ; LSC98 ; CLL99 . The fundamental quantity of Carter’s superfluid formalism is the so-called master function. Although we have several choices for the master function, we take, following Comer, Langlois, and Lin CLL99 , the total thermodynamical energy density $`\mathrm{\Lambda }`$ as the master function. The master function $`\mathrm{\Lambda }`$ is assumed to depend on the three scalars $`n^2=n_\alpha n^\alpha `$, $`p^2=p_\alpha p^\alpha `$, and $`x^2=p_\alpha n^\alpha `$, where $`n^\alpha `$ and $`p^\alpha `$ are the conserved number density current of the neutrons and the protons, respectively.

Once an explicit functional form of the master function is given, the energy-momentum tensor $`T_\beta ^\alpha `$ is given by

$$T_\beta ^\alpha =\mathrm{\Psi }\delta _\beta ^\alpha +p^\alpha \chi _\beta +n^\alpha \mu _\beta ,$$
(1)
where the scalar $`\mathrm{\Psi }`$ denotes the generalized pressure, defined by

$$\mathrm{\Psi }=\mathrm{\Lambda }n^\alpha \mu _\alpha p^\alpha \chi _\alpha ,$$
(2)
and the one-forms $`\mu _\alpha `$ and $`\chi _\alpha `$ mean the chemical potential covectors, defined by

$`\mu _\alpha =n_\alpha +𝒜p_\alpha ,\chi _\alpha =𝒜n_\alpha +𝒞p_\alpha ,`$ (3)
where

$`𝒜={\displaystyle \frac{\mathrm{\Lambda }}{x^2}},=2{\displaystyle \frac{\mathrm{\Lambda }}{n^2}},𝒞=2{\displaystyle \frac{\mathrm{\Lambda }}{p^2}}.`$ (4)
The one-forms $`\mu _\alpha `$ and $`\chi _\alpha `$ are conjugate momenta to $`n^\alpha `$ and $`p^\alpha `$, respectively. The quantities $`\mu =(\mu _\alpha \mu ^\alpha )^{1/2}`$ and $`\chi =(\chi _\alpha \chi ^\alpha )^{1/2}`$ are interpreted as the chemical potentials of the neutrons and the protons, respectively. From equation (3), we can see that the thermodynamical quantity $`𝒜`$ determines the strength of the entrainment effects between the two superfluids. In other words, if we have $`𝒜=0`$, there is no entrainment effect. The system of the dynamical equations for the two superfluids is composed of two continuity equations, given by

$$_\alpha n^\alpha =0,_\alpha p^\alpha =0.$$
(5)
and two Euler equations, given by

$$n^\alpha _{[\alpha }\mu _{\beta ]}=0,p^\alpha _{[\alpha }\chi _{\beta ]}=0,$$
(6)
where $`_\alpha `$ means the covariant derivative associated with the metric tensor, and the square brackets denote antisymmetrised averaging. Since $`n^\alpha `$ and $`p^\alpha `$ are the conserved currents, it is convenient to introduce two time-like unit vectors $`u^\alpha `$ and $`v^\alpha `$ defined by

$$n^\alpha =nu^\alpha ,p^\alpha =pv^\alpha ,$$
(7)
where $`n`$ and $`p`$ denote the number density of the neutrons and the protons, respectively.

### II.2 Equilibrium configurations

Equilibrium states of a slowly rotating star are described by the stationary axisymmetric spacetime, given by the line element

$`ds^2=g_{\alpha \beta }dx^\alpha dx^\beta =e^{2\nu (r)}dt^2+e^{2\lambda (r)}dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2)2\omega r^2\mathrm{sin}^2\theta dtd\phi .`$ (8)
Here, we have taken account of the rotational effects up to the first order of the stellar rotation angular frequency. In this study, we assume that the neutron and proton superfluids in equilibrium are in the same rotational motion with the uniform rotation angular frequency $`\mathrm{\Omega }`$. The four-velocity of the two superfluids is therefore given by

$$u^\alpha =v^\alpha =\gamma (t^\alpha +\mathrm{\Omega }\phi ^\alpha )=e^{\nu (r)}(t^\alpha +\mathrm{\Omega }\phi ^\alpha ),$$
(9)
where $`t^\alpha `$ and $`\phi ^\alpha `$ are the time-like and rotational Killing vectors of the spacetime. In virtue of equation (9), the conjugate momentum one-forms in equilibrium reduce to

$`\mu _\alpha `$ $`=`$ $`(n+𝒜p)u_\alpha =\mu (n,p)u_\alpha ,`$
$`\chi _\alpha `$ $`=`$ $`(𝒜n+𝒞p)u_\alpha =\chi (n,p)u_\alpha ,`$ (10)
and the energy-momentum tensor of an equilibrium star is given by

$`T_\beta ^\alpha =\mathrm{\Lambda }u^\alpha u_\beta +\mathrm{\Psi }q_\beta ^\alpha =\mathrm{\Lambda }u^\alpha u_\beta +(\mathrm{\Lambda }+\mu n+\chi p)q_\beta ^\alpha ,`$ (11)
where $`q_\beta ^\alpha `$ denotes the projection tensor associated with the four-velocity $`u^\alpha `$, defined by $`q_\beta ^\alpha =\delta _\beta ^\alpha +u^\alpha u_\beta `$.

Substituting equations (10) into Euler equations (6), we obtain the hydrostatic equations given by

$`{\displaystyle \frac{1}{\mu }}{\displaystyle \frac{d\mu }{dr}}={\displaystyle \frac{d\nu }{dr}},{\displaystyle \frac{1}{\chi }}{\displaystyle \frac{d\chi }{dr}}={\displaystyle \frac{d\nu }{dr}},`$ (12)
which leads to

$`{\displaystyle \frac{\mu }{\chi }}=C,`$ (13)
where $`C`$ is an integral constant. In this study, we take $`C=1`$, assuming that the two superfluids are in chemical equilibrium in the equilibrium states. The metric coefficients are determined from the Einstein equations, which are written as

$`{\displaystyle \frac{d\nu }{dr}}=e^{2\lambda }\left({\displaystyle \frac{M}{r^2}}+4\pi r\mathrm{\Psi }\right),{\displaystyle \frac{dM}{dr}}=4\pi r^2\mathrm{\Lambda },`$ (14)
$`{\displaystyle \frac{d^2\varpi }{dr^2}}+\left({\displaystyle \frac{4}{r}}{\displaystyle \frac{d\nu }{dr}}{\displaystyle \frac{d\lambda }{dr}}\right){\displaystyle \frac{d\varpi }{dr}}{\displaystyle \frac{4}{r}}\left({\displaystyle \frac{d\nu }{dr}}+{\displaystyle \frac{d\lambda }{dr}}\right)\varpi =0,`$ (15)
where

$`M(r)={\displaystyle \frac{r}{2}}(1e^{2\lambda }),\varpi =\mathrm{\Omega }\omega .`$ (16)
In this paper, we neglect a normal fluid envelope (e.g., CLL99 ), and we require that the two superfluids have the common outer surface whose radius $`R`$ is defined by equations $`\mu (R)=m_n`$ and $`\chi (R)=m_p`$, where $`m_n`$ and $`m_p`$ stand for the rest mass of the neutron and the proton, respectively. The surface boundary conditions for the equilibrium structure are given as follows:

$`\nu ={\displaystyle \frac{1}{2}}\mathrm{ln}\left(1{\displaystyle \frac{2M(R)}{R}}\right),\varpi +{\displaystyle \frac{R}{3}}{\displaystyle \frac{d\varpi }{dr}}=\mathrm{\Omega },\mathrm{at}r=R,`$ (17)
where $`M(R)`$ stands for the gravitational mass of the star. The boundary conditions at the stellar center are the regularity condition for all the physical quantities.

### II.3 Perturbation equations for the two-fluid model

To describe fluid perturbations in a star, it is convenient to introduce two kinds of changes in physical quantities, called the Eulerian and the Lagrangian changes JF78 . An Eulerian change $`\delta Q`$ is the difference between the quantities $`Q`$ in the perturbed and the unperturbed states at a spacetime point. The relation between the Lagrangian changes and the Eulerian changes are given by

$`\mathrm{\Delta }_nQ=\delta Q+L_{\xi _n}Q,\mathrm{\Delta }_pQ=\delta Q+L_{\xi _p}Q,`$ (18)
where $`L_k`$ stands for the Lie derivative along the vector $`k^\alpha `$, and we have introduced the two distinct Lagrangian displacement vectors $`\xi _n^\alpha `$ and $`\xi _p^\alpha `$ to describe the perturbed motion of the neutron and proton superfluids because the two superfluids can move independently. The Lagrangian displacement is regarded as a vector that connects fluid elements in the unperturbed state to the corresponding elements in the perturbed state.

We suppose that perturbations of any physical quantity can be expressed in terms of the Lagrangian displacements $`\xi _n^\alpha `$ and $`\xi _p^\alpha `$, and the Eulerian changes in the metric $`h_{\alpha \beta }=\delta g_{\alpha \beta }`$. According to similar consideration to that in JF78 , we can write the Lagrangian changes of the fluid velocities of the neutron and proton in terms of the Lagrangian change in the metric:

$`\mathrm{\Delta }_nu^\alpha ={\displaystyle \frac{1}{2}}u^\alpha u^\mu u^\nu \mathrm{\Delta }_ng_{\mu \nu },\mathrm{\Delta }_pv^\alpha ={\displaystyle \frac{1}{2}}u^\alpha u^\mu u^\nu \mathrm{\Delta }_pg_{\mu \nu },`$ (19)
where

$`\mathrm{\Delta }_ng_{\alpha \beta }=h_{\alpha \beta }+_\alpha \xi _{n\beta }+_\beta \xi _{n\alpha },\mathrm{\Delta }_pg_{\alpha \beta }=h_{\alpha \beta }+_\alpha \xi _{p\beta }+_\beta \xi _{p\alpha }.`$ (20)
Here, we have used the relation $`u^\alpha =v^\alpha `$ in the equilibrium. In this study, we consider no transfusion between the neutrons and the protons. The particle numbers of the neutrons and the protons are therefore conserved separately. Then, the conservation equations are given by

$`{\displaystyle \frac{\mathrm{\Delta }_nn}{n}}={\displaystyle \frac{1}{2}}q^{\alpha \beta }\mathrm{\Delta }_ng_{\alpha \beta },{\displaystyle \frac{\mathrm{\Delta }_pp}{p}}={\displaystyle \frac{1}{2}}q^{\alpha \beta }\mathrm{\Delta }_pg_{\alpha \beta }.`$ (21)
The Eulerian perturbations of the velocities and the number densities can be expressed in terms of $`\xi _n^\alpha `$, $`\xi _p^\alpha `$, and $`h_{\alpha \beta }`$ as

$`\delta u^\alpha `$ $`=`$ $`q_\beta ^\alpha L_u\xi _n^\beta +{\displaystyle \frac{1}{2}}u^\alpha u^\mu u^\nu h_{\mu \nu }=\delta \widehat{u}^\alpha +{\displaystyle \frac{1}{2}}u^\alpha u^\mu u^\nu h_{\mu \nu },`$
$`\delta v^\alpha `$ $`=`$ $`q_\beta ^\alpha L_u\xi _p^\beta +{\displaystyle \frac{1}{2}}u^\alpha u^\mu u^\nu h_{\mu \nu }=\delta \widehat{v}^\alpha +{\displaystyle \frac{1}{2}}u^\alpha u^\mu u^\nu h_{\mu \nu },`$ (22)
$`\delta n={\displaystyle \frac{n}{2}}q^{\alpha \beta }h_{\alpha \beta }q_\beta ^\alpha _\alpha (n\xi _n^\beta ),\delta p={\displaystyle \frac{p}{2}}q^{\alpha \beta }h_{\alpha \beta }q_\beta ^\alpha _\alpha (p\xi _p^\beta ).`$ (23)
According to C89  and assuming $`u^\alpha =v^\alpha `$ in the equilibrium, we can write the Eulerian perturbations of the conjugate momenta $`\mu _\alpha `$ and $`\chi _\alpha `$ as

$`\delta \mu _\alpha `$ $`=`$ $`\mu \left(\delta \widehat{u}_\alpha +h_{\alpha \beta }u^\beta +{\displaystyle \frac{1}{2}}u_\alpha u^\mu u^\nu h_{\mu \nu }\right)+p𝒜(\delta \widehat{v}_\alpha \delta \widehat{u}_\alpha )+\left\{(𝒜\overline{𝒜})\delta p+(\overline{})\delta n\right\}u_\alpha ,`$
$`\delta \chi _\alpha `$ $`=`$ $`\chi \left(\delta \widehat{v}_\alpha +h_{\alpha \beta }u^\beta +{\displaystyle \frac{1}{2}}u_\alpha u^\mu u^\nu h_{\mu \nu }\right)+n𝒜(\delta \widehat{u}_\alpha \delta \widehat{v}_\alpha )+\left\{(𝒞\overline{𝒞})\delta p+(𝒜\overline{𝒜})\delta n\right\}u_\alpha ,`$ (24)
where

$`\overline{𝒜}`$ $`=`$ $`2np{\displaystyle \frac{}{p^2}}2n^2{\displaystyle \frac{𝒜}{n^2}}2p^2{\displaystyle \frac{𝒜}{p^2}}np{\displaystyle \frac{𝒜}{x^2}},`$
$`\overline{}`$ $`=`$ $`2n^2{\displaystyle \frac{}{n^2}}4np{\displaystyle \frac{𝒜}{n^2}}p^2{\displaystyle \frac{𝒜}{x^2}},`$
$`\overline{𝒞}`$ $`=`$ $`2p^2{\displaystyle \frac{𝒞}{p^2}}4np{\displaystyle \frac{𝒜}{p^2}}n^2{\displaystyle \frac{𝒜}{x^2}}.`$ (25)
Note that the second terms on the right-hand side of equations (24) represent the non-dissipative drag force between the two superfluids, which is proportional to both the function $`𝒜`$ and the velocity difference $`\delta \widehat{v}_\alpha \delta \widehat{u}_\alpha `$. From equation (6), the perturbed Euler equations are given by

$`q_\alpha ^\beta L_u\delta \mu _\beta `$ $``$ $`q_\alpha ^\beta _\beta (u^\rho \delta \mu _\rho )+\mu \delta \widehat{u}^\beta (_\beta u_\alpha _\alpha u_\beta )=0,`$
$`q_\alpha ^\beta L_u\delta \chi _\beta `$ $``$ $`q_\alpha ^\beta _\beta (u^\rho \delta \chi _\rho )+\chi \delta \widehat{v}^\beta (_\beta u_\alpha _\alpha u_\beta )=0,`$ (26)
where $`_\alpha `$ means the partial derivative.

It is convenient to introduce the vorticity equation when pulsations in a rotating star are considered. Making use of the conjugate momenta as dynamical variables, the vorticity equations can be obtained as

$`L_ud\mu _{\alpha \beta }=0,L_vd\chi _{\alpha \beta }=0,`$ (27)
where $`d\mu _{\alpha \beta }`$ and $`d\chi _{\alpha \beta }`$ are the exterior differentiation of the one-forms $`\mu _\beta `$ and $`\chi _\beta `$, defined by

$`d\mu _{\alpha \beta }=_\alpha \mu _\beta _\beta \mu _\alpha ,d\chi _{\alpha \beta }=_\alpha \chi _\beta _\beta \chi _\alpha .`$ (28)
Then, the Lagrangian variation of these equations can be derived straightforwardly:

$`\mathrm{\Delta }_nL_ud\mu _{\alpha \beta }=L_u\mathrm{\Delta }_nd\mu _{\alpha \beta }=0,\mathrm{\Delta }_pL_vd\chi _{\alpha \beta }=L_v\mathrm{\Delta }_pd\chi _{\alpha \beta }=0,`$ (29)
where $`\mathrm{\Delta }_nd\mu _{\alpha \beta }`$ and $`\mathrm{\Delta }_pd\chi _{\alpha \beta }`$ can be written in terms of the Lagrangian changes in the conjugate momenta as

$`\mathrm{\Delta }_nd\mu _{\alpha \beta }=d\mathrm{\Delta }_n\mu _{\alpha \beta }=_\alpha \mathrm{\Delta }_n\mu _\beta _\beta \mathrm{\Delta }_n\mu _\alpha ,\mathrm{\Delta }_pd\chi _{\alpha \beta }=d\mathrm{\Delta }_p\chi _{\alpha \beta }=_\alpha \mathrm{\Delta }_p\chi _\beta _\beta \mathrm{\Delta }_p\chi _\alpha .`$ (30)
Note that since the equations $`u^\beta d\mathrm{\Delta }_n\mu _{\alpha \beta }=0`$ and $`_{[\alpha }d\mathrm{\Delta }_n\mu _{\beta \gamma ]}=0`$ ($`v^\beta d\mathrm{\Delta }_p\chi _{\alpha \beta }=0`$ and $`_{[\alpha }d\mathrm{\Delta }_p\chi _{\beta \gamma ]}=0`$) are satisfied, $`d\mathrm{\Delta }_n\mu _{\alpha \beta }`$ ($`d\mathrm{\Delta }_p\chi _{\alpha \beta }`$) has only two independent components.

The metric perturbations $`h_{\alpha \beta }`$ are determined by the linearized Einstein equations, given by

$`\delta G_\beta ^\alpha =8\pi \delta T_\beta ^\alpha ,`$ (31)
where $`\delta G_\beta ^\alpha `$ denotes the linearized Einstein tensor. Here, $`\delta T_\beta ^\alpha `$ means the linearized energy-momentum tensor, and is given by

$`\delta T_\beta ^\alpha =\delta \mathrm{\Psi }\delta _\beta ^\alpha +\delta n^\alpha \mu _\beta +\delta p^\alpha \chi _\beta +n^\alpha \delta \mu _\beta +p^\alpha \delta \chi _\beta ,`$ (32)
where

$`\delta \mathrm{\Psi }={\displaystyle \frac{1}{2}}(n^\alpha \mu ^\beta +p^\alpha \chi ^\beta )h_{\alpha \beta }n^\alpha \delta \mu _\alpha p^\alpha \delta \chi _\alpha .`$ (33)
### II.4 Equations of state: Expanded master function

In this study, we make use of the same master function as that introduced in ACL02 , which is generally given by

$$\mathrm{\Lambda }(n^2,p^2,x^2)=\underset{i=0}{\overset{\mathrm{}}{}}\lambda _i(n^2,p^2)(x^2np)^i.$$
(34)
This expansion of the master function may be justified because $`x^2np=0`$ in equilibrium and the deviation from $`x^2np=0`$ for perturbed states is the same order of the perturbations. Although Comer and Joynt CJ02  have recently discussed in greater detail the entrainment effects in general relativistic superfluids, we use the master function given above for simplicity. In terms of the master function (34), we can obtain the thermodynamical functions in pulsation equations as follows:

$`𝒜`$ $`=`$ $`\lambda _1(n^2,p^2),={\displaystyle \frac{p}{n}}𝒜{\displaystyle \frac{1}{n}}{\displaystyle \frac{\lambda _0}{n}},𝒞={\displaystyle \frac{n}{p}}𝒜{\displaystyle \frac{1}{p}}{\displaystyle \frac{\lambda _0}{p}},`$ (35)
$`\overline{𝒜}`$ $`=`$ $`𝒜+{\displaystyle \frac{^2\lambda _0}{np}},\overline{}=+{\displaystyle \frac{^2\lambda _0}{n^2}},\overline{𝒞}=𝒞+{\displaystyle \frac{^2\lambda _0}{p^2}},`$ (36)
where we have used the relation $`x^2=np`$ in the equilibrium. Note that $`\overline{𝒜}`$ becomes equal to $`𝒜`$ if all the expansion coefficients $`\lambda _i`$ are separable in $`n`$ and $`p`$ in the sense that $`\lambda _i=f_i(n)+g_i(p)`$ for appropriate functions $`f_i(n)`$ and $`g_i(p)`$. Using the chemical potentials in the equilibrium given by

$`\mu ={\displaystyle \frac{\lambda _0}{n}},\chi ={\displaystyle \frac{\lambda _0}{p}},`$ (37)
we have

$`𝒜\overline{𝒜}={\displaystyle \frac{\mu }{p}}={\displaystyle \frac{\chi }{n}},\overline{}={\displaystyle \frac{\mu }{n}},𝒞\overline{𝒞}={\displaystyle \frac{\chi }{p}},`$ (38)
and, assuming all the expansion coefficients $`\lambda _i`$ are separable in $`n`$ and $`p`$, we obtain

$`𝒜=\overline{𝒜},\overline{}={\displaystyle \frac{d\mu }{dn}},𝒞\overline{𝒞}={\displaystyle \frac{d\chi }{dp}}.`$ (39)
### II.5 Oscillation equations for inertial modes

If we assume the star in the equilibrium is stationary and axisymmetric, the time and azimuthal dependence of the perturbations can be given by $`\mathrm{exp}(i\sigma t+im\phi )`$, where $`\sigma `$ is the oscillation frequency observed by an inertial observer at spatial infinity, and $`m`$ is the azimuthal wave number. Because of the rotation effects, it is generally impossible to achieve the separation of variables for the perturbed quantities. We expand the perturbations in terms of the tensor spherical harmonics with different $`l`$’s for a given $`m`$. In this paper, we consider the oscillation modes associated with $`m2`$ because we are interested in the modes that are unstable against gravitational radiation reactions. Thus, we can select the so-called Regge-Wheeler gauge in order to fix the gauge freedom for the metric perturbations RW . Then, the metric perturbations are expanded as

$`h_{\mu \nu }`$ $`=`$ $`{\displaystyle \underset{l|m|}{}}\left(\begin{array}{cccc}e^{2\nu }H_{0l}(r)& H_{1l}(r)& 0& 0\\ \mathrm{sym}& e^{2\lambda }H_{2l}(r)& 0& 0\\ \mathrm{sym}& \mathrm{sym}& K_l(r)r^2& 0\\ \mathrm{sym}& \mathrm{sym}& \mathrm{sym}& r^2\mathrm{sin}^2\theta K_l(r)\end{array}\right)Y_l^m(\theta ,\phi )e^{i\sigma t},`$ (40)
$`+`$ $`{\displaystyle \underset{l^{}|m|}{}}\left(\begin{array}{cccc}0& 0& ih_{0l^{}}(r)\frac{1}{\mathrm{sin}\theta }_\phi Y_l^{}^m(\theta ,\phi )& ih_{0l^{}}(r)\mathrm{sin}\theta _\theta Y_l^{}^m(\theta ,\phi )\\ \mathrm{sym}& 0& ih_{1l^{}}(r)\frac{1}{\mathrm{sin}\theta }_\phi Y_l^{}^m(\theta ,\phi )& ih_{1l^{}}(r)\mathrm{sin}\theta _\theta Y_l^{}^m(\theta ,\phi )\\ \mathrm{sym}& \mathrm{sym}& 0& 0\\ \mathrm{sym}& \mathrm{sym}& \mathrm{sym}& 0\end{array}\right)e^{i\sigma t},`$
where $`Y_l^m(\theta ,\phi )`$ is the spherical harmonic function, and $`l=|m|+2(k1)`$ and $`l^{}=l+1`$ for even modes, and $`l=|m|+2k1`$ and $`l^{}=l1`$ for odd modes, where $`k=1,2,3,\mathrm{}`$. The even and odd modes are characterized by symmetry and antisymmetry of the eigenfunctions with respect to the equatorial plane. The Lagrangian displacements $`\xi _{n\alpha }`$ and $`\xi _{p\alpha }`$ can be expanded in terms of the vector spherical harmonics as

$`i\kappa \mathrm{\Omega }\xi _{nt}=0,i\kappa \mathrm{\Omega }\xi _{pt}=0,`$ (41)
$$i\kappa \mathrm{\Omega }\xi _{nr}=\underset{l|m|}{}\frac{e^{2\lambda }}{r}W_{nl}(r)Y_l^m(\theta ,\phi )e^{i\sigma t},i\kappa \mathrm{\Omega }\xi _{pr}=\underset{l|m|}{}\frac{e^{2\lambda }}{r}W_{pl}(r)Y_l^m(\theta ,\phi )e^{i\sigma t},$$
(42)
$`i\kappa \mathrm{\Omega }\xi _{n\theta }`$ $`=`$ $`{\displaystyle \underset{l,l^{}|m|}{}}\left(V_{nl}(r)_\theta Y_l^m(\theta ,\phi )iU_{nl^{}}{\displaystyle \frac{1}{\mathrm{sin}\theta }}_\phi Y_l^{}^m(\theta ,\phi )\right)e^{i\sigma t},`$
$`i\kappa \mathrm{\Omega }\xi _{p\theta }`$ $`=`$ $`{\displaystyle \underset{l,l^{}|m|}{}}\left(V_{pl}(r)_\theta Y_l^m(\theta ,\phi )iU_{pl^{}}{\displaystyle \frac{1}{\mathrm{sin}\theta }}_\phi Y_l^{}^m(\theta ,\phi )\right)e^{i\sigma t},`$ (43)
$`i\kappa \mathrm{\Omega }\xi _{n\phi }`$ $`=`$ $`{\displaystyle \underset{l,l^{}|m|}{}}\left(V_{nl}(r)_\phi Y_l^m(\theta ,\phi )+iU_{nl^{}}\mathrm{sin}\theta _\theta Y_l^{}^m(\theta ,\phi )\right)e^{i\sigma t},`$
$`i\kappa \mathrm{\Omega }\xi _{p\phi }`$ $`=`$ $`{\displaystyle \underset{l,l^{}|m|}{}}\left(V_{pl}(r)_\phi Y_l^m(\theta ,\phi )+iU_{pl^{}}\mathrm{sin}\theta _\theta Y_l^{}^m(\theta ,\phi )\right)e^{i\sigma t},`$ (44)
where $`\kappa \mathrm{\Omega }\sigma +m\mathrm{\Omega }`$ is the oscillation frequency in the corotating frame of the star. Note that the gauge freedom for the Lagrangian displacement is fixed so as to satisfy $`\xi _{nt}=\xi _{pt}=0`$. The Eulerian perturbations of the number densities for the neutrons and protons are given by

$`\delta n={\displaystyle \underset{l|m|}{}}\delta n_l(r)Y_l^m(\theta ,\phi )e^{i\sigma t},\delta p={\displaystyle \underset{l|m|}{}}\delta p_l(r)Y_l^m(\theta ,\phi )e^{i\sigma t}.`$ (45)
When we are interested in inertial modes in a slowly rotating superfluid star in the lowest order in the rotation frequency $`\mathrm{\Omega }`$, it is rather easy to generalize the formulation devised by Lockitch et al. LAF01  for relativistic inertial modes in normal fluid neutron stars. According to LAF01 , we assume that the perturbations satisfy the following ordering laws in the slow rotation limit of $`\mathrm{\Omega }0`$:

$`H_1`$ $`=`$ $`O(1),W_n=O(1),W_p=O(1),V_n=O(1),V_p=O(1),`$ (46)
$`H_0`$ $`=`$ $`O(\mathrm{\Omega }),H_2=O(\mathrm{\Omega }),K=O(\mathrm{\Omega }),\delta n=O(\mathrm{\Omega }),\delta p=O(\mathrm{\Omega }),`$ (47)
$`h_0=O(1),U_p=O(1),U_n=O(1),`$ (48)
and $`\kappa =O(1)`$. Comer GC02  showed that solutions obeying these ordering laws are allowed in the perturbation equations for rotating stars, exploring the so-called zero-frequency subspace of eigensolutions to perturbation equations in a non-rotating relativistic superfluid star. If we consider no rotational effects, the solutions subject to the ordering laws (47) are interpreted as infinitely degenerate $`g`$-modes of zero-frequency, while the solutions subject to the ordering laws (48) are consider to be a relativistic and superfluid counterpart of the trivial toroidal mode in a non-rotating Newtonian normal fluid star Uetal .

Substituting the perturbed quantities into linearized equations (23), (29), and (31) and assuming the ordering laws for the eigenfunctions (47) and (48), we can obtain a system of infinitely coupled ordinary differential equations for inertial modes in a slowly rotating superfluid star. To write down the oscillation equations for the inertial modes, it is convenient to use vector notation for the eigenfunctions $`𝐰_n`$, $`𝐰_p`$, $`𝐯_n`$, $`𝐯_p`$, $`𝐮_n`$, $`𝐮_p`$, $`𝐡`$, and $`𝐇`$, whose components are given by

$`w_{n,k}=W_{nl},w_{p,k}=W_{pl},v_{n,k}=V_{nl},v_{p,k}=V_{pl},u_{n,k}=U_{nl^{}},u_{p,k}=U_{pl^{}},h_k=h_{0l^{}},H_k=H_{1l},`$ (49)
where $`l=|m|+2(k1)`$, $`l^{}=l+1`$ for even modes and $`l=|m|+2k1`$, $`l^{}=l1`$ for odd modes, where $`k=1,2,3,\mathrm{}`$. The particle number conservation equations (23) for the neutron and proton are then written as

$`e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐰_n)}{dr}}`$ $`=`$ $`\left(1+{\displaystyle \frac{r}{n}}{\displaystyle \frac{dn}{dr}}+r{\displaystyle \frac{d\lambda }{dr}}+2r{\displaystyle \frac{d\nu }{dr}}\right)𝐰_n+\mathrm{\Lambda }_0𝐯_n,`$
$`e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐰_p)}{dr}}`$ $`=`$ $`\left(1+{\displaystyle \frac{r}{p}}{\displaystyle \frac{dp}{dr}}+r{\displaystyle \frac{d\lambda }{dr}}+2r{\displaystyle \frac{d\nu }{dr}}\right)𝐰_p+\mathrm{\Lambda }_0𝐯_p.`$ (50)
Independent components of the vorticity equations (29) lead to

$`L_1𝐮_nqM_0𝐯_n+\overline{q}K_0𝐰_n+𝒜_1(𝐮_p𝐮_n)+𝐡=0,`$
$`L_1𝐮_pqM_0𝐯_p+\overline{q}K_0𝐰_p+𝒜_2(𝐮_n𝐮_p)+𝐡=0,`$ (51)
$`L_0e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐯_n)}{dr}}`$ $``$ $`qM_1e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐮_n)}{dr}}m\overline{q}\mathrm{\Lambda }_0^1e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐰_n)}{dr}}\left(e^{2\lambda }+mr{\displaystyle \frac{d\overline{q}}{dr}}\mathrm{\Lambda }_0^1\right)𝐰_n`$
$`+`$ $`\left(\overline{q}K_1\mathrm{\Lambda }_1r{\displaystyle \frac{dq}{dr}}M_1\right)𝐮_n+m\left(\overline{q}r{\displaystyle \frac{dq}{dr}}\mathrm{\Lambda }_0^1\right)𝐯_nr𝐇`$
$`+`$ $`r{\displaystyle \frac{d𝒜_1}{dr}}(𝐯_p𝐯_n)+𝒜_1\left(e^{2\lambda }(𝐰_n𝐰_p)+e^{2\nu }r{\displaystyle \frac{d}{dr}}(e^{2\nu }𝐯_pe^{2\nu }𝐯_n)\right)=0,`$
$`L_0e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐯_p)}{dr}}`$ $``$ $`qM_1e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐮_p)}{dr}}m\overline{q}\mathrm{\Lambda }_0^1e^{2\nu }r{\displaystyle \frac{d(e^{2\nu }𝐰_p)}{dr}}\left(e^{2\lambda }+mr{\displaystyle \frac{d\overline{q}}{dr}}\mathrm{\Lambda }_0^1\right)𝐰_p`$ (52)
$`+`$ $`\left(\overline{q}K_1\mathrm{\Lambda }_1r{\displaystyle \frac{dq}{dr}}M_1\right)𝐮_p+m\left(\overline{q}r{\displaystyle \frac{dq}{dr}}\mathrm{\Lambda }_0^1\right)𝐯_pr𝐇`$
$`+`$ $`r{\displaystyle \frac{d𝒜_2}{dr}}(𝐯_n𝐯_p)+𝒜_2\left(e^{2\lambda }(𝐰_p𝐰_n)+e^{2\nu }r{\displaystyle \frac{d}{dr}}(e^{2\nu }𝐯_ne^{2\nu }𝐯_p)\right)=0,`$
where

$`q=2\kappa ^1{\displaystyle \frac{\varpi }{\mathrm{\Omega }}},\overline{q}=\kappa ^1{\displaystyle \frac{e^{2\nu }}{r}}{\displaystyle \frac{d}{dr}}\left(e^{2\nu }r^2{\displaystyle \frac{\varpi }{\mathrm{\Omega }}}\right),𝒜_1={\displaystyle \frac{p}{\mu }}𝒜,𝒜_2={\displaystyle \frac{n}{\chi }}𝒜.`$ (53)
The linearized Einstein equations are reduced to

$`r𝐇=\mathrm{\Lambda }_0^116\pi e^{2\lambda }r^2(n\mu 𝐰_n+p\chi 𝐰_p),`$ (54)
$`r{\displaystyle \frac{d}{dr}}\left(r{\displaystyle \frac{d}{dr}}𝐡\right)\left(1+r{\displaystyle \frac{d\nu }{dr}}+r{\displaystyle \frac{d\lambda }{dr}}\right)r{\displaystyle \frac{d}{dr}}𝐡\left(e^{2\lambda }\mathrm{\Lambda }_1+2\left(1e^{2\lambda }+r{\displaystyle \frac{d\nu }{dr}}+r{\displaystyle \frac{d\lambda }{dr}}\right)\mathrm{𝟏}\right)𝐡=16\pi e^{2\lambda }r^2(n\mu 𝐮_n+p\chi 𝐮_p),`$ (55)
where 1 stands for the unit matrix. Here, the matrices $`L_0`$, $`L_1`$, $`K_0`$, $`K_1`$, $`\mathrm{\Lambda }_0`$, $`\mathrm{\Lambda }_1`$, $`M_0`$, and $`M_1`$ are defined in Appendix. The oscillation equations we use are given as a set of linear ordinary differential equations for the variables $`𝐰_n`$, $`𝐰_p`$, $`𝐮_n`$, $`𝐮_p`$, $`𝐡`$, and $`r\frac{d𝐡}{dr}`$, which are obtained by eliminating the vectors $`𝐯_n`$, $`𝐯_p`$, and $`𝐇`$ in equations (50), (52), and (55) using the algebraic equations (51) and (54).

The boundary conditions imposed at the center of the star are regularity condition that all the perturbation functions are regular and do not diverge at $`r=0`$. This implies that the eigenfunctions we solve must vanish at $`r=0`$. In order to determine the boundary conditions at the stellar surface, on the other hand, we follow the arguments similar to those given in CLL99 . In our specific equilibrium models, for which we have assumed an equation of state that is separable in $`n`$ and $`p`$, the generalized pressure $`\mathrm{\Psi }`$ of the fluids can be written as $`\mathrm{\Psi }(n,p)=\mathrm{\Psi }_n(n)+\mathrm{\Psi }_p(p)`$, where $`\mathrm{\Psi }_n(n)`$ and $`\mathrm{\Psi }_p(p)`$ represent the pressures due to the neutron and proton superfluids, respectively. Since the two superfluids can flow independently, as appropriate surface boundary conditions we require that the Lagrangian changes in each of the pressures vanish at the free surface of the star, which leads to

$`\mathrm{\Delta }_n\mathrm{\Psi }_n={\displaystyle \frac{d\mathrm{\Psi }_n}{dn}}\mathrm{\Delta }_nn=0,\mathrm{\Delta }_p\mathrm{\Psi }_p={\displaystyle \frac{d\mathrm{\Psi }_p}{dp}}\mathrm{\Delta }_np=0.`$ (56)
The conditions $`\mathrm{\Delta }_nn=0`$ and $`\mathrm{\Delta }_pp=0`$ at the surface therefore result in the surface boundary conditions given by $`W_{nl}=0`$ and $`W_{pl}=0`$ for inertial mode solutions in the lowest order of $`\mathrm{\Omega }`$. In the exterior of the star, the only non-trivial metric perturbations $`h_{0l^{}}`$ are determined by equation (55), which has two independent solutions, one is regular at spatial infinity, the other singular. The regular solution can be given as a series expansion LAF01 ; LFA02 :

$`h_{0l^{}}={\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}\widehat{h}_{l^{},s}\left({\displaystyle \frac{R}{r}}\right)^{l^{}+s},`$ (57)
where

$`\widehat{h}_{l^{},s}={\displaystyle \frac{(l^{}+s2)!(l^{}+s+1)!(2l^{}+1)!}{s!(l^{}2)!(l^{}+1)!(2l^{}+s+1)!}}\left({\displaystyle \frac{2M(R)}{R}}\right)^s\widehat{h}_{l^{},0},`$ (58)
and $`\widehat{h}_{l^{},0}`$ is an arbitrary constant. If we require that the interior solution $`h_{0l^{}}`$ must be continuous with the exterior solution (57) at the stellar surface because the spacetime must be regular everywhere, we obtain the boundary conditions for the metric perturbations at the stellar surface given by

$`\underset{ϵ0}{lim}\left[h_{0l^{}}(R+ϵ)r{\displaystyle \frac{d}{dr}}h_{0l^{}}(Rϵ)r{\displaystyle \frac{d}{dr}}h_{0l^{}}(R+ϵ)h_{0l^{}}(Rϵ)\right]=0.`$ (59)
For numerical computation, oscillation equations of a finite dimension are obtained by disregarding the terms with $`l`$ larger than $`l_{max}`$ in the expansions of the perturbations, where $`l_{max}`$ is determined so that the eigenfrequency and the eigenfunctions are well converged as $`l_{max}`$ increases. We solve the oscillation equations of a finite dimension as an eigenvalue problem with the scaled oscillation frequency $`\kappa `$ using a Henyey type relaxation method (see, e.g., Uetal ).

## III Numerical results

### III.1 Equilibrium models

For the equilibrium structure with $`x^2np=0`$, we use a generalized polytropic equation of state, introduced by Comer et al. CLL99 , which is given by

$`\lambda _0(n^2,p^2)=m_nn\sigma _nn^{\beta _n}m_pp\sigma _pp^{\beta _p},`$ (60)
where $`\sigma _n`$, $`\sigma _p`$, $`\beta _n`$, and $`\beta _p`$ are constants. Note that this master function is separable in $`n`$ and $`p`$. The generalized pressure and chemical potentials in the equilibrium can then be written as

$`\mathrm{\Psi }`$ $`=`$ $`\sigma _n(\beta _n1)n^{\beta _n}+\sigma _p(\beta _p1)p^{\beta _p},`$ (61)
$`\mu `$ $`=`$ $`m_n+\sigma _n\beta _nn^{\beta _n1},`$ (62)
$`\chi `$ $`=`$ $`m_p+\sigma _p\beta _pp^{\beta _p1}.`$ (63)
Because chemical equilibrium is assumed in the equilibrium state, we can analytically write $`p`$ in terms of $`n`$ as

$`p=\left({\displaystyle \frac{\sigma _n\beta _n}{\sigma _p\beta _p}}\right)^{1/(\beta _p1)}n^{(\beta _n1)/(\beta _p1)},`$ (64)
where we have assumed $`m_p=m_n`$. We confirm from equations (63) and (64) that two superfluids in an equilibrium state given by equation (60) have the common outer surface if $`\beta _n`$ and $`\beta _p`$ satisfy the conditions $`\beta _n1`$ and $`\beta _p1`$. In order to introduce the entrainment effects in the perturbed states with $`x^2np0`$, we include, according to ACL02 , the expansion coefficient $`\lambda _1`$, given by

$`\lambda _1=\eta {\displaystyle \frac{m_nm_p}{m_pp+\eta (m_nn+m_pp)}},`$ (65)
where $`\eta `$ is a parameter whose appropriate range is considered to be $`0.04\eta 0.2`$ S76 ; BJK96 ; ACL02 . In this paper, we call $`\eta `$ the entrainment parameter. Making use of the formulas obtained in the preceding sections, we can explicitly write the relevant thermodynamical coefficients as follows:

$`𝒜=\overline{𝒜}=\eta {\displaystyle \frac{m_nm_p}{m_pp+\eta (m_nn+m_pp)}}.`$ (66)
In order to avoid unnecessary singular behaviors of $`p`$ and $`n`$ at the stellar surface, we consider only the case of $`\beta _n=\beta _p`$. Thus, all the models we use are non-stratified models in the sense that $`n/p=\mathrm{constant}`$ inside the star. The constant parameters $`\sigma _n`$, $`\sigma _p`$, $`\beta _n`$, and $`\beta _p`$ in $`\lambda _0`$ for the equilibrium equation of state (see Table 1) are the same as those used for model 1 calculated in CLL99 . In this paper, we consider two equilibrium models, whose physical parameters are tabulated in Table I. The model I is almost Newtonian whose relativistic factor $`M(R)/R`$ is $`M(R)/R=10^3`$, while the model II is a relativistic one with the relativistic factor $`M(R)/R=0.15`$, the value of which is similar to those for the models used in LFA02 . Although the model I is not appropriate as a physical model of neutron stars, we use it to compare with completely Newtonian calculations LY03 ; YL03 .

### III.2 $`r`$-mode oscillations

We calculated the $`r`$-modes of the two equilibrium models for a wide range of the entrainment parameter. It is important to note that the $`l^{}=m`$ fundamental $`r`$-modes, whose eigenfunction $`U_{nm}`$ and/or $`U_{pm}`$ is dominating and has no nodes in radial direction, are the only $`r`$-modes we find in this study. This situation is quite similar to that for the $`r`$-modes in a barotropic ordinary fluid star, in which neither the overtone $`r`$-modes with $`l^{}=m`$ nor the $`r`$-modes with $`l^{}m`$ are found LF99 ; YL00a ; YL00b ; LFA02 . It is generally found that the $`r`$-modes in a relativistic superfluid star are split into two families, which we call ordinary fluid $`r`$-modes ($`r^o`$-modes) and superfluid $`r`$-modes ($`r^s`$-modes), according to the notation of Lee and Yoshida LY03 . Because of doubling of the dynamical degrees of freedom for the system of two superfluids, the mode splitting of this kind also appears in other oscillation modes essentially associated with fluid motions E88 ; LM94 ; L95 ; AC01 ; ACL02 ; PR02 ; LY03 ; YL03 . In Figures 1 and 2, we plot the scaled frequencies $`\kappa \sigma /\mathrm{\Omega }+m`$ of the $`r`$-modes associated with $`m=2`$ and $`3`$ as a function of the entrainment parameter $`\eta `$. As shown by the figures, although $`\kappa `$ of the $`r^s`$-modes increases almost linearly with increasing $`\eta `$, $`\kappa `$ of the $`r^o`$-modes is almost independent of $`\eta `$. Comparing the two frequency curves for the models I and II, it is found that both the $`r^o`$-mode and $`r^s`$-mode frequencies are strongly dependent on the relativistic factor $`M(R)/R`$ of the equilibrium models. This is because the frequency of the $`r`$-modes is roughly proportional to an $`average`$ of the function $`\varpi `$ within the star, which reflects the general relativistic effects such as the rotational frame dragging and the gravitational redshift. Similar dependence of $`\kappa `$ on the parameter $`M(R)/R`$ for the $`r`$-modes has been found for relativistic ordinary fluid stars YK98 ; LAF01 ; SY01 ; RK02 ; LFA02 ; YL02 .

From Table II, where $`\kappa `$’s for the $`r^o`$\- and $`r^s`$-modes associated with $`m=2`$ and $`3`$ are tabulated for several values of $`\eta `$, we can confirm that $`\kappa `$’s for the $`r^o`$-modes for model I are nearly equal to $`2m/l^{}(l^{}+1)`$, the value found for the $`r`$-modes in Newtonian ordinary fluid stars. It is important to note that $`\kappa `$’s for the $`r^o`$\- and $`r^s`$-modes at $`\eta =0`$ are not equal to each other, and the difference between the $`\kappa `$’s increases as the relativistic factor $`M(R)/R`$ is increased. In a Newtonian star, however, $`\kappa `$’s for the $`r^o`$\- and $`r^s`$-modes at $`\eta =0`$ are equal to each other in the lowest order of $`\mathrm{\Omega }`$, as shown in LY03 . The reason for the difference between the Newtonian and relativistic $`r`$-modes at $`\eta =0`$ is that even at $`\eta =0`$ there still exist couplings between two superfluid motions through the Einstein equations in the relativistic case, but no such couplings exist in the Newtonian case because all the Newtonian gravitational perturbations vanish in the lowest order of $`\mathrm{\Omega }`$.

In order to illustrate how the eigenfunctions of the $`r^o`$\- and $`r^s`$-modes in a relativistic superfluid star behave, it is convenient to introduce a new set of variables defined by

$`W_{+l}`$ $`=`$ $`{\displaystyle \frac{nW_{nl}+pW_{pl}}{n+p}},W_l=W_{nl}W_{pl},`$
$`V_{+l}`$ $`=`$ $`{\displaystyle \frac{nV_{nl}+pV_{pl}}{n+p}},V_l=V_{nl}V_{pl},`$ (67)
$`U_{+l^{}}`$ $`=`$ $`{\displaystyle \frac{nU_{nl^{}}+pU_{pl^{}}}{n+p}},U_l^{}=U_{nl^{}}U_{pl^{}}.`$
In Figures 3 and 4, we display four dominant coefficients $`U_{+2}`$, $`W_{+3}`$, $`U_{+4}`$, and $`h_{02}`$ for the $`r^o`$-mode, and $`U_2`$, $`W_3`$, $`U_4`$, and $`h_{02}`$ for the $`r^s`$-mode for the case of $`m=2`$ and $`\eta =0.1`$ for model II. Here, the normalization condition has been given by $`U_{nm}=1`$ at $`r=R`$. The coefficients that are not displayed in these figures have very small and negligible amplitudes. Figures 3 and 4 represent that the basic properties of the eigenfunctions are quite similar to those of the $`r^o`$\- and $`r^s`$-modes in a Newtonian star with superfluidity LY03 ; YL03 . The neutrons and protons co-move for the $`r^o`$-modes, while they counter-move for the $`r^s`$-modes. It is also found in Figure 4 that the metric perturbation $`h_{02}`$ almost completely vanishes for the $`r^s`$-modes, which is consistent with the fact that $`U_{+m}0`$. This means that gravitational radiations due to $`r^s`$-mode oscillations are negligible at least in the lowest order in $`\mathrm{\Omega }`$. It is important to note that the other coefficients $`W_{\pm 3}`$ and $`U_{\pm 4}`$ are not necessarily negligible compared with $`U_{\pm 2}`$ because of the general relativistic effects LAF01 ; YL02 ; LFA02 .

## IV Discussions and Conclusions

In this paper we have studied the modal properties of the $`r`$-modes of relativistic superfluid neutron stars, taking account of the entrainment effects between the neutron and proton superfluids. To describe the general relativistic dynamics of the superfluids, we employed the two-constituent formalism developed by Carter and his co-workers C89 ; CL98 ; LSC98 ; CLL99 . We derived the perturbation equations for the relativistic inertial modes in neutron stars filled with the superfluids by generalizing the formalism employed by Lockitch, Andersson, and Friedman LAF01 ; LFA02  for the relativistic inertial modes in normal fluid neutron stars. We found that the basic properties of the $`r`$-modes in a relativistic star with the two superfluids are very similar to those in a Newtonian superfluid star LY03 ; YL03 . We confirmed that the $`r`$-modes of relativistic superfluid stars are split into two families, ordinary fluid-like $`r`$-modes ($`r^o`$-mode) and superfluid-like $`r`$-modes ($`r^s`$-mode). The two superfluids counter-move for the $`r^s`$-modes, while they co-move for the $`r^o`$-modes. The dimensionless frequency $`\kappa `$ for the $`r^o`$-modes is almost independent of the entrainment parameter $`\eta `$. For the $`r^s`$-modes, on the other hand, $`\kappa `$ almost linearly increases with $`\eta `$. The gravitational radiation driven instability due to the $`r^s`$-modes is much weaker than that of the $`r^o`$-modes because the matter current associated with the axial parity perturbations vanish almost completely for the former.

A solid crust near the surface of the neutron star has a significant influence on the modal properties of the $`r`$-modes in the superfluid core, since the solid crust supports its own oscillation modes MVH88 ; LS96 ; YL02b ; MPS01 , and resonance phenomena are expected between the $`r`$-modes in the superfluid core and the torsional sound waves in the solid crust even in the general relativistic context YL01 ; LU01 . Dissipation in the viscous boundary layer at the interface between the fluid core and the solid crust is another important issue for the $`r`$-mode instability, as shown by ekman . How the dissipation in the core-crust interface can be important for the $`r`$-mode instability will be a quite interesting problem when superfluid neutrons and protons in the core and superfluid neutrons in the inner crust are taken into account simultaneously (see, e.g., LCHE93 ).

Newtonian superfluid neutron stars can support infinite number of inertial modes, which are also split into two families, ordinary fluid-like inertial modes ($`i^o`$-mode) and superfluid-like inertial modes ($`i^s`$-mode) LY03 ; YL03 . Non-linear couplings between the $`r`$-modes and inertial modes will be important to limit the amplitude growth of the $`r`$-mode instability, as recently shown by Arras et al. AEL  for Newtonian normal fluid stars. To investigate the relativistic inertial modes in superfluid neutron stars will be one of our future studies.

###### Acknowledgements.

SY acknowledges financial support from Fundação para a Ciência e Tecnologia (FCT) through project SAPIENS 36280/99.

## Appendix: Matrices used in pulsation equations

The components of the matrices, $`L_0`$, $`L_1`$, $`K_0`$, $`K_1`$, $`\mathrm{\Lambda }_0`$, $`\mathrm{\Lambda }_1`$, $`M_0`$, and $`M_1`$ are given as follows:

For even modes,

$$(L_0)_{i,i}=1\frac{mq}{l(l+1)},(L_1)_{i,i}=1\frac{mq}{(l+1)(l+2)},$$

$$(K_0)_{i,i}=\frac{J_{l+1}^m}{l+1},(K_0)_{i,i+1}=\frac{J_{l+2}^m}{l+2},$$

$$(K_1)_{i,i}=\frac{J_{l+1}^m}{l+1},(K_1)_{i+1,i}=\frac{J_{l+2}^m}{l+2},$$

$$(\mathrm{\Lambda }_0)_{i,i}=l(l+1),(\mathrm{\Lambda }_1)_{i,i}=(l+1)(l+2),$$

$$(M_0)_{i,i}=\frac{l}{l+1}J_{l+1}^m,(M_0)_{i,i+1}=\frac{l+3}{l+2}J_{l+2}^m,$$

$$(M_1)_{i,i}=\frac{l+2}{l+1}J_{l+1}^m,(M_1)_{i+1,i}=\frac{l+1}{l+2}J_{l+2}^m,$$

where $`l=|m|+2i2`$ for $`i=1,2,3,\mathrm{}`$.

For odd modes,

$$(L_0)_{i,i}=1\frac{mq}{l(l+1)},(L_1)_{i,i}=1\frac{mq}{l(l1)},$$

$$(K_0)_{i,i}=\frac{J_l^m}{l},(K_0)_{i+1,i}=\frac{J_{l+1}^m}{l+1},$$

$$(K_1)_{i,i}=\frac{J_l^m}{l},(K_1)_{i,i+1}=\frac{J_{l+1}^m}{l+1},$$

$$(\mathrm{\Lambda }_0)_{i,i}=l(l+1),(\mathrm{\Lambda }_1)_{i,i}=l(l1),$$

$$(M_0)_{i,i}=\frac{l+1}{l}J_l^m,(M_0)_{i+1,i}=\frac{l}{l+1}J_{l+1}^m,$$

$$(M_1)_{i,i}=\frac{l1}{l}J_l^m,(M_1)_{i,i+1}=\frac{l+2}{l+1}J_{l+1}^m,$$

where $`l=|m|+2i1`$ for $`i=1,2,3,\mathrm{}`$. Here, $`J_l^m`$ is the function of $`m`$ and $`l`$, defined by $`J_l^m=[(l^2m^2)/(4l^21)]^{1/2}`$, and $`q=2\varpi /\kappa \mathrm{\Omega }`$.


format: text

text:
# Properties of the Ideal Ginzburg-Landau Vortex Lattice

## I Introduction

Since Abrikosov’s 1  prediction of the flux-line lattice in Type-II superconductors from Ginzburg-Landau (GL) theory 2 , several approximate formulas for the magnetization $`M=B/\mu _0H`$versus the applied magnetic field $`H`$ or average induction $`B`$ have been published 1 ; 3 ; 4 ; 5 ; 6 ; 7 . In these papers and below, the basic situation is considered where a macroscopically large, homogeneous and isotropic, long superconductor is exposed to a uniform parallel field $`H`$. In this ideal case demagnetization effects, flux-line pinning, and surface effects may be disregarded, and thus the flux-lines are straight lines forming an ideal periodic lattice. These results are easily extended to anisotropic superconductors (where an anisotropic effective-mass tensor is introduced into the GL theory) by defining an effective GL parameter $`\stackrel{~}{\kappa }`$ that depends on the orientation of the flux lines; this transformation works when $`H`$ is along a principle symmetry axis 8 ; 9 ; 10 . Generalization to geometries where demagnetization effects occur, are possible by introduction of a demagnetizing factor; but this concept works only for homogeneous specimens with the shape of an ellipsoid. In this case the flux lines in the bulk are still straight and form an ideal flux-line lattice (FLL). For other specimen shapes the FLL is distorted, i.e., the orientation and density of the FLL varies spatially and can be calculated only numerically 11 ; 12 .

The aim of the present paper is to compare the widely used approximate expressions for $`M(H,\kappa )`$ with the exact value obtained numerically and to give useful general analytic interpolation formulas valid in the entire ranges of $`H`$ and $`\kappa `$ where the FLL exists, namely, $`H_{c1}HH_{c2}`$ for $`H`$, or $`0BB_{c2}=\mu _0H_{c2}`$, and $`1/\sqrt{2}\kappa <\mathrm{}`$ for $`\kappa `$, where $`H_{c1}(T)`$ and $`H_{c2}(T)`$ are the lower and upper critical fields and $`\kappa `$ is the GL parameter. Interestingly, such general formulas have not been published yet, and thus the accuracy of the commonly used expressions is not known, probably due to the difficulty of the numerical solution of the complex-valued GL equations. Early numerics 13  used the circular cell method (CCM), which approximates the hexagonal unit cell of the triangular FLL (or the quadratic unit cell of the square FLL) by a circle and the two-dimensional (2D) solution by the 1D rotationally symmetric solution inside this circular cell; both the GL function and the magnetic field are forced to have vanishing slope on this circular boundary, as the exact solution has on the boundary of the Wigner-Seitz cell. This method yields the exact $`H_{c1}`$ and is expected to be best at low inductions $`BB_{c2}`$ where the flux lines are well separated. But surprisingly, the circular cell approximation gives very good magnetization curves at all $`B`$ (see Fig. 1) and even yields the exact value of the upper critical field $`H_{c2}`$. Some more exact results of the CCM are listed below. Another method 7  uses a similar circular symmetric GL order parameter and a linear superposition of circular symmetric magnetic fields to obtain excellent approximate $`M(H,\kappa )`$, see also Ref. 14, .

An in principle exact numerical method 15  uses periodic real trial functions for the squared GL function $`|\psi (x,y)|^2`$ and magnetic field $`B(x,y)`$ and minimizes the resulting free energy functional with respect to a finite number of Fourier coefficients. The same method was later applied 16  to solve the microscopic BCS-Gor’kov theory for the properties of the FLL in the entire temperature interval $`0TT_c`$ where $`T_c`$ is the superconducting transition temperature (GL theory, strictly spoken, applies only close to $`T_c`$). Recently this variational method was improved 17  by keeping the same periodic trial functions but now solving the GL equations iteratively; this iteration works much faster and allows to use many more Fourier coefficients (many thousands instead of only five in Ref. 15, ). I shall use this 2D iterative precision method of Ref. 17,  for the calculation of the FLL at $`B>0`$. At low inductions $`BB_{c2}`$ this 2D method is supplemented by an iterative circular cell method presented in Appendix A. This 1D method yields accurate values of $`h_{c1}(\kappa )=H_{c1}/H_{c2}`$, which then can be used in interpolation formulas. For convenience, I introduce the reduced fields $`b=B/B_{c2}`$, $`h=H/H_{c2}`$, $`m=M/H_{c2}`$, such that one has $`m=bh`$, $`h_{c1}h1`$, $`0b1`$, and $`h_{c1}m0`$.

For completeness it should be mentioned that the isolated vortex 18  and the FLL 19  have also been computed from BCS theory (valid at all temperatures) using the quasiclassical Eilenberger theory based on energy-integrated Green functions. This method was recently extended to compute the FLL structure and local density of states for s-wave 19 ; 20 ; 21 , d-wave 21 ; 22 , and chiral p-wave 23  superconductors. Very recently the GL method 17  was generalized phenomenologically to lower temperatures and to charged vortices 24 .

## II Triangular and square flux-line lattices and the circular cell method

The properties of the FLL within GL theory are calculated by minimizing the GL free energy of the superconductor with respect to the complex GL function $`\psi (𝐫)`$ and to the vector potential $`𝐀(𝐫)`$ of the local magnetic field $`𝐁(𝐫)=\times 𝐀`$. In the usual reduced units 1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7  (length $`\lambda `$, magnetic field $`\sqrt{2}H_c`$, energy density $`\mu _0H_c^2`$, where $`H_c`$ is the thermodynamic critical field) the spatially averaged free energy $`F`$ of the GL theory referred to the Meissner state ($`\psi =1`$, $`𝐁=0`$) reads

$$F=\frac{(1|\psi |^2)^2}{2}+\left|\left(\frac{}{i\kappa }𝐀\right)\psi \right|^2+𝐁^2.$$
(1)
Here $`\mathrm{}=(1/V)d^3r\mathrm{}`$ means spatial averaging over the superconductor of volume $`V`$. Introducing the supervelocity $`𝐐(𝐫)=𝐀\phi /\kappa `$ and the magnitude $`f(𝐫)=|\psi |`$ of $`\psi (𝐫)=f(𝐫)\mathrm{exp}[i\phi (𝐫)]`$ one may express $`F`$ as a functional of the real and gauge invariant functions $`f`$ and $`𝐐`$,

$$F=\frac{(1f^2)^2}{2}+\frac{(f)^2}{\kappa ^2}+f^2Q^2+(\times 𝐐)^2.$$
(2)
In the presence of vortices $`𝐐(𝐫)`$ has to be chosen such that $`\times 𝐐`$ has the appropriate singularities along the vortex cores, see e.g. Eq.(B4) in App. B.

In this paper I consider the ideal periodic FLL in a homogeneous (pin-free) large superconductor in a uniform magnetic field $`H`$ along $`z`$. In this 2D situation one has $`f=f(x,y)`$, $`𝐐=𝐐(x,y)`$, and $`𝐁=\widehat{𝐳}B(x,y)`$. Within GL theory in reduced units the properties of this ideal FLL depend only on two parameters: the GL parameter $`\kappa `$ and the average induction $`B=B(x,y)`$. The equilibrium magnetic field $`H`$, and the magnetization $`M=B/\mu _0H`$, are obtained either from the definition $`H=F/B`$ or, more elegantly, from the virial theorem discovered by Doria, Gubernatis, and Rainer, 25  which in reduced units reads

$$H=\frac{f^2f^4+2B(x,y)^2}{2B}.$$
(3)
Some of the properties of the FLL, and all properties of the isolated flux line, may be calculated in an elegant way by the circular cell approximation 7 ; 13 ; 14  as described in App. A. In the circular cell method the hexagonal Wigner-Seitz cell around each flux line is replaced by a circle with radius $`R`$ and same area $`\pi R^2=\mathrm{\Phi }_0/B`$ if each flux line carries one quantum of flux $`\mathrm{\Phi }_0=h/2e=2.0710^^5`$ Tm2. In reduced units one has $`\mathrm{\Phi }_0=2\pi /\kappa `$ and $`R/\lambda =R=(2/b\kappa ^2)^{1/2}`$ with $`b=B/B_{c2}`$. The boundary conditions on the CCM circle $`r=R`$ are $`df/dr=dB/dr=0`$. I find that the free energy of the triangular FLL, $`F_{tr}`$, and its magnetization, $`M_{tr}`$, are reproduced by the CCM with high accuracy in the entire ranges of $`\kappa `$ and $`B`$, $`1/\sqrt{2}\kappa <\mathrm{}`$ and $`0b<1`$. In particular, the CCM not only yields $`H_{c1}`$ (in the limit $`R\mathrm{}`$) but it also reproduces the exact upper critical field $`H_{c2}(\kappa )`$, and in the special case $`\kappa =1/\sqrt{2}`$ even the exact result $`H(B)=\mathrm{const}=H_c=H_{c1}=H_{c2}`$. These somewhat surprising features of this approximation are related to the facts that $`H_{c2}`$, and in the case $`\kappa =1/\sqrt{2}`$ even the entire curve $`H(B)`$, are independent of the detailed arrangement of the flux lines, i.e., they are the same for triangular and square or honey-comb FLLs and for any other arrangement of single or multiple quanta flux lines. Another surprising finding is that the virial theorem, Eq. (3), works perfectly in the CCM. Figure 1 shows the magnetization curves $`M(H)`$ and the equilibrium field $`H(B)`$ of the superconductor obtained by the CCM for $`\kappa =0.85`$, 1, 1.2, 1.5, 2, 3, 5, 7, 10, and 20.

In the limit $`b0`$ the CCM yields the lower critical field $`H_{c1}`$, which with high accuracy is fitted by the formula

$`\mu _0H_{c1}={\displaystyle \frac{\mathrm{\Phi }_0}{4\pi \lambda ^2}}[\mathrm{ln}\kappa +\alpha (\kappa )],`$
$`h_{c1}={\displaystyle \frac{H_{c1}}{H_{c2}}}={\displaystyle \frac{\mathrm{ln}\kappa +\alpha (\kappa )}{2\kappa ^2}},`$
$`\alpha (\kappa )=\alpha _{\mathrm{}}+\mathrm{exp}[c_0c_1\mathrm{ln}\kappa c_2(\mathrm{ln}\kappa )^2]\pm ϵ`$ (4)
with $`\alpha _{\mathrm{}}=0.49693`$, $`c_0=0.41477`$, $`c_1=0.775`$, $`c_2=0.1303`$, and $`ϵ0.00076`$. This expression yields at $`\kappa =1/\sqrt{2}`$ the correct value $`h_{c1}=1`$ and for $`\kappa 1`$ it has the limit $`\alpha =0.49693`$. A simpler expression for $`\alpha (\kappa )`$, yielding an $`h_{c1}`$ with error still less than 1% and with the correct limits at $`\kappa =1/\sqrt{2}`$ and $`\kappa 1`$, is

$`\alpha (\kappa )=0.5+{\displaystyle \frac{1+\mathrm{ln}2}{2\kappa \sqrt{2}+2}}.(4\mathrm{a})`$
The CCM in principle cannot yield properties related to the different symmetries of the FLL, or to its shear modulus, and it cannot give the form factors (Fourier coefficients) of the magnetic field $`B(x,y)`$ that may be measured by neutron scattering. These subtle properties can be computed by the 2D method presented in Ref. 17,  and in App. B. This effective numerical method expresses the smooth functions $`f(x,y)^2`$ and $`B(x,y)`$ as 2D Fourier series and determines the Fourier coefficients by iteration.

Figure 2 (top) shows the difference of the free energy densities of the triangular ($`F_{tr}`$) and square ($`F_{sq}`$) FLLs. This difference is proportional to the shear modulus $`c_{66}`$ of the triangular FLL (the shear modulus of the unstable square FLL is negative within GL theory) by the relation 17

$$c_{66}=(3\pi ^2/2)(F_{sq}F_{tr}).$$
(5)
Note that this difference is very small, $`0<(F_{sq}F_{tr})/(\mu _0H_c^2)<0.0018`$. Even smaller (by ten times) is the difference between the free energy densities of the CCM ($`F_{cc}`$) and of the triangular FLL plotted in Fig. 2 (bottom). One has $`0<(F_{cc}F_{tr})/(\mu _0H_c^2)<0.00020`$. This result shows that the CCM is an excellent approximation for global properties of the FLL. Both differences are largest for large $`\kappa `$ and have a maximum near $`b0.3`$. The finding $`F_{sq}>F_{tr}`$ means that the triangular FLL is stable for all $`\kappa >1/\sqrt{2}`$. Note that for $`\kappa =1/\sqrt{2}`$ one has exactly $`F_{sq}=F_{cc}=F_{tr}=0`$ for all $`b`$.

Figure 3 (top) shows the difference between the magnetizations $`M_{sq}`$ of the square FLL and $`M_{tr}`$ of the triangular FLL. Again, this difference is small, $`0.0008<(M_{sq}M_{tr})/H_{c2})0.00014`$ and the relative difference has the limits $`0.018<(M_{sq}M_{tr})/M_{tr}0.0095`$. Figure 3 (bottom) shows the difference between the magnetization $`M_{cc}`$ obtained by the CCM (see Fig. 1) and the exact value $`M_{tr}`$ of the triangular lattice. Like with the free energy, this difference is again smaller by a factor of ten than the difference between two lattice symmetries, $`0.00016<(M_{cc}M_{tr})/H_{c2}0.00008`$ and $`0.0011<(M_{cc}M_{tr})/M_{tr}0.0017`$. The differences vanish exactly at $`\kappa =1/\sqrt{2}`$, and also at $`\kappa \mathrm{}`$, since there $`m=M/H_{c2}0`$. The relative differences (insets in Fig. 3) are maximum at $`\kappa 1`$.

The smallness of these differences explains why in Fig. 1 the magnetization curves for all three cases $`M_{tr}`$, $`M_{sq}`$, and $`M_{cc}`$ coincide within line thickness.

Figure 4 shows an example ($`b=0.3`$, $`\kappa =1.5`$) comparing the spatial functions $`f`$ and $`B`$ of the triangular FLL and from the CCM. Shown are the cross sections $`f(x,0)`$ along the nearest neighbor direction $`x`$ and $`f(0,y)`$ perpendicular to this, and $`f(r)`$ from the CCM \[ $`a`$ is the vortex spacing, $`a^2/\lambda ^2=4\pi /(\sqrt{3}b\kappa ^2)`$ \]. It is seen that $`f(x,0)`$ and $`f(r)`$, and also $`B(x,0)`$ and $`B(r)`$, coincide closely; at lower $`b<0.3`$ the difference is smaller than the line thickness. The solutions for the square FLL deviate more from the circular cell solutions.

The maximum, minimum, and saddle-point fields of the triangular FLL, $`B_{max}=B(0,0)`$, $`B_{min}=B(0,a/\sqrt{3})`$, and $`B_{sad}=B(a/2,0)`$, depend on $`b`$ and $`\kappa `$. $`B_{max}`$ is only slightly above the equilibrium field $`H`$, and $`B_{sad}`$ and $`B_{min}`$ are close to each other and lie somewhat below the average field $`B`$. $`B_{max}`$ and $`B_{min}`$ are shown in Fig. 3 of Ref. 17,  as functions of $`b`$ for several $`\kappa =0.707\mathrm{}5`$. In Fig. 5 the small differences $`B_{max}H`$, $`B_{sad}B`$, and $`B_{sad}B_{min}`$ are plotted versus $`b`$, in units $`B_{c2}`$ and multiplied by a function of $`\kappa `$ such that the curves for all $`\kappa 1/\sqrt{2}`$ collapse at $`b1`$. One finds for all $`\kappa `$ near $`b=1`$:

$`{\displaystyle \frac{B_{max}H}{B_{c2}}}0.0351{\displaystyle \frac{\kappa ^20.5}{(\kappa ^20.069)^2}}(1b)^2,`$ (6)
$`{\displaystyle \frac{B_{sad}B}{B_{c2}}}0.146{\displaystyle \frac{1b}{\kappa ^20.069}},`$ (7)
$`{\displaystyle \frac{B_{sad}B_{min}}{B_{c2}}}0.0526{\displaystyle \frac{(1b)}{\kappa ^20.069}}.`$ (8)
The factor 0.069 in Eqs. (6)–(8) is $`0.50.5/\beta _A=0.0688`$ where $`\beta _A=1.1596`$ is the Abrikosov parameter of the triangular FLL. Plots of $`B_{cc}(R)B_{min}`$ where $`B_{cc}(R)`$ is the field value at the boundary of the circular cell in the CCM, look similar to the plots of $`B_{sad}B_{min}`$ in Fig. 5 (lower panel), since the value $`B_{cc}(R)`$ lies approximately in the middle between $`B_{min}`$ and $`B_{sad}`$, see Fig. 4. Since for $`\kappa 1`$ and $`b1/\kappa ^2`$ the field in the vortex center equals $`B_{max}=2H_{c1}`$, one has $`B_{max}HH_{c1}`$, and thus the function plotted in Fig. 5 (upper panel) for $`b0`$ tends to the limit $`(b_{max}h)\times \kappa ^2h_{c1}\kappa ^2\frac{1}{2}(\mathrm{ln}\kappa +0.50)`$, cf. Eq. (4).

The variance of the magnetic field is

$`\sigma =[B(x,y)B]^2=B(x,y)^2B^2={\displaystyle \underset{𝐊0}{}}B_𝐊^2,`$ (9)
where $`B_𝐊`$ are the Fourier coefficients of $`B(x,y)=_𝐊B_𝐊\mathrm{cos}\mathrm{𝐊𝐫}`$ and $`𝐊`$ the vectors of the reciprocal lattice of the FLL (App. B). Near $`b=1`$ the Abrikosov solution of the linearized GL theory 11 ; 26  yields for all $`\kappa `$ values 27

$`\sigma =7.5210^4{\displaystyle \frac{\mathrm{\Phi }_0^2}{\lambda ^4}}{\displaystyle \frac{\kappa ^4(1b)^2}{(\kappa ^20.069)^2}},`$
$`S{\displaystyle \frac{\sqrt{\sigma }}{B_{c2}}}=0.172{\displaystyle \frac{1b}{\kappa ^20.069}}.`$ (10)
The functions $`S`$ and $`S/(1b)`$ are plotted in Fig. 6 versus $`\sqrt{b}`$ for various $`\kappa `$. It can be seen that Eq. (10) is a rather good approximation for the large range $`0.25<b<1`$. At smaller $`b`$ the variance $`\sigma (b)`$ has a maximum and then goes to zero again at $`b=0`$.

For small inductions $`b1`$ and large $`\kappa `$ one can use the London approximation $`B_𝐊=B/(1+K^2\lambda ^2)`$. For the appropriate cutoff at large magnitudes $`K\xi ^1=\kappa /\lambda `$ see Ref. 28, ; 29,  and below. In the range $`0.13/\kappa ^2b1`$ the unity in the denominator of $`B_𝐊`$ may be disregarded since $`K^2\lambda ^2(4\pi /\sqrt{3})b\kappa ^2=7.255b\kappa ^2`$. Thus $`B`$ drops out and $`\sigma `$ becomes independent of $`b`$: 27

$`\sigma =0.00371{\displaystyle \frac{\mathrm{\Phi }_0^2}{\lambda ^4}},S={\displaystyle \frac{0.383}{\kappa ^2}}.`$ (11)
This often used formula corresponds to the upper axis in Fig. 6. One can see that this approximation is good only for very large $`\kappa 70`$ and only in the range of $`b`$ near the maximum of $`\sigma `$. At very small $`b0.13/\kappa ^2`$ both $`\sigma (b)`$ and $`S(\sqrt{b})`$ drop linearly to zero when $`b0`$. In this range the sum in Eq. (9) can be evaluated as an integral, yielding

$`\sigma ={\displaystyle \frac{b\kappa ^2}{8\pi ^2}}{\displaystyle \frac{\mathrm{\Phi }_0^2}{\lambda ^4}},S={\displaystyle \frac{\sqrt{b/2}}{\kappa }}.`$ (12)
This approximation is good for $`\kappa 5`$ and very small $`b`$ ($`b<0.01/\kappa ^2`$ for $`\kappa =5`$; $`b<0.04/\kappa ^2`$ for $`\kappa 10`$), see the two straight lines in Fig. 6. For $`\kappa 5`$ a good approximation, with less than 5% error from $`b=1`$ down to the value $`b0.25/\kappa ^{1.3}`$ where the maximum of $`\sigma `$ occurs, is

$$S\frac{\sqrt{\sigma }}{B_{c2}}0.172\frac{1b}{\kappa ^2}[\mathrm{\hspace{0.17em}1}+1.21(1\sqrt{b})^3].$$
(13)
This approximation is much better that the interpolation, Eq. (3) of Ref. 27, .

## III Magnetization curves

This section presents analytic expressions which approximate the computed magnetization $`m=M/H_{c2}=bh`$ (Fig. 1) as a function of the induction $`b=B/B_{c2}`$ or of the thermodynamic field $`h=H/H_{c2}`$. We distinguish approximations working at high or low inductions.

### III.1 Approximation for high inductions

The linearized GL theory yields for $`1b1`$ Abrikosov’s $`B_{c2}`$ solution 1 ; 11

$$mm_A=\frac{1b}{(2\kappa ^21)\beta _A+1},$$
(14)
where $`\beta _A=\omega _A^2/\omega _A^2=1+_{m,n}\mathrm{exp}[K_{mn}^2S/(4\pi )]`$ (Ref. 11, ; 30,  and App. B) is the Abrikosov parameter, $`\beta _A=1.1595953`$ ($`\beta _A=1.1803406`$) for the triangular (square) FLL. The linear magnetization $`m_A(b,\kappa )`$ is a good approximation in the range $`0.5b1`$, see Fig. 1. This suggests the following fit to the exact $`m`$:

$`m(b,\kappa )=m_A(1b)^2\mathrm{exp}[f_1(b)]g_1(\kappa )+ϵ_1,`$
$`f_1(b)=2.50u^28.08u+0.39,u=(1b)^{0.41},`$
$`g_1(\kappa )=(1.133+1.926/\kappa ^{2.25})(2\kappa ^21)/(2\kappa ^4).`$ (15)
with relative error $`|ϵ_1/m|<0.0013`$ for $`b>0.5`$ for the triangular FLL. Formula (15) is a good approximation with relative error $`<1\%`$ for all $`\kappa `$ in the large range of fields $`1/(4\kappa ^2)+510^4b1`$, see Fig. 7.

The same expression (14) fits also the $`m(b,\kappa )`$ of the square FLL, with somewhat larger error if the same functions $`f_1(b)`$, $`g_1(\kappa )`$ are used rather than the optimally fitted ones. For the difference $`m_{\mathrm{tr}}m_{\mathrm{sq}}`$ see Fig. 3.

### III.2 Approximation for “intermediate fields”

For completeness I mention here also the London approximation 3  which was supposed to be good in the “intermediate field range” $`H_{c1}HH_{c2}`$ that exists only in superconductors with extremely large $`\kappa `$. Within London theory the induction is (see App. B)

$`B(x,y)=B{\displaystyle \underset{𝐊}{}}{\displaystyle \frac{\mathrm{cos}\mathrm{𝐊𝐫}}{1+K^2\lambda ^2}},`$ (16)
where the sum goes over all $`K`$-vectors with length from $`K=0`$ to some cutoff $`K\xi ^1`$. Inserting this into the London free energy density $`[B(𝐫)^2+\lambda ^2(\times 𝐁)^2]/(2\mu _0)`$ and averaging over the superconductor one gets

$`F={\displaystyle \underset{𝐊}{}}{\displaystyle \frac{(B^2/2\mu _0)}{1+K^2\lambda ^2}}{\displaystyle \frac{B^2}{2\mu _0}}+{\displaystyle \frac{B\mathrm{\Phi }_0}{2\mu _0}}{\displaystyle \frac{d^2k}{4\pi ^2}\frac{1}{k^2\lambda ^2}}.`$ (17)
The integral from $`k_{\mathrm{min}}^2(K_{10}/2)^2\pi ^2B/\mathrm{\Phi }_0`$ to $`k_{\mathrm{max}}^2\xi ^2=2\pi B_{c2}/\mathrm{\Phi }_0`$ equals $`(4\pi \lambda ^2)^1\mathrm{ln}(\gamma ^{}/b)`$ where $`\gamma ^{}`$ is some constant and $`b=B/B_{c2}`$ as above. This yields

$`M`$ $`=`$ $`H{\displaystyle \frac{B}{\mu _0}}={\displaystyle \frac{F}{B}}{\displaystyle \frac{B}{\mu _0}}={\displaystyle \frac{\mathrm{\Phi }_0}{8\pi \lambda ^2\mu _0}}\mathrm{ln}{\displaystyle \frac{\gamma }{b}},`$
$`m`$ $`=`$ $`{\displaystyle \frac{M}{B_{c2}}}={\displaystyle \frac{1}{4\kappa ^2}}\mathrm{ln}{\displaystyle \frac{0.358}{b}}`$ (18)
with constant $`\gamma =\gamma ^{}/e=0.3575\mathrm{}`$ obtained by our fit to the numerical $`m(b)`$ at $`\kappa =200`$. This old London approximation is shown in Fig. 7 as dotted lines. One sees that this fit works only at large $`\kappa 20`$ in the relatively small interval $`1/(2\kappa ^2)b0.01`$, i.e. at very low $`b`$ (but not too low $`b`$). It gives $`m=0`$ at $`b=\gamma `$ for all $`\kappa `$. This fit is slightly improved by replacing $`\mathrm{ln}(\gamma /b)`$ by $`\mathrm{ln}(1\gamma +\gamma /b)`$, which gives the correct $`m=0`$ at $`b=1`$.

A much better fit in the spirit of this logarithmic approximation is (see Fig. 8)

$`m`$ $`=`$ $`{\displaystyle \frac{1}{4\kappa ^2}}\mathrm{ln}\left[1+{\displaystyle \frac{1b}{b}}f_2(b)\right],`$
$`f_2(b)`$ $`=`$ $`0.357+2.890b1.581b^2.`$ (19)
This fit is good for $`\kappa 3`$ (error $`<3\%`$) and $`\kappa 5`$ (error $`<1\%`$ ) in the large ranges $`(\mathrm{ln}\kappa +1)/(10\kappa ^2)b1`$ for $`\kappa =3\mathrm{}200`$. These intervals of validity may also be expressed as $`M/H_{c1}=m/h_{c1}0.8`$ (0.85) for $`\kappa 20`$ ($`\kappa 50`$).

### III.3 Approximation for low inductions

All the above approximations do not describe the correct vertical slope of $`M(H)`$ at $`H=H_{c1}`$, or zero slope of $`H(B)`$ and unity slope of $`M(B)=BH`$ at $`B=0`$. This is achieved by the London approximation of pairwise interacting vortices described in App. C. For very small $`b1`$ one may account only for the nearest neighbor shell of six vortices in the triangular FLL of spacing $`a=c\lambda `$. With $`h(b)`$, Eq. (C8), this yields for $`m(b)=h(b)b`$:

$`mh_{c1}b+{\displaystyle \frac{3\sqrt{\pi c}}{2\kappa ^2}}e^c\left[1+{\displaystyle \frac{19}{8c}}{\displaystyle \frac{47}{128c^2}}\right],`$
$`c={\displaystyle \frac{a}{\lambda }}=\left({\displaystyle \frac{4\pi /\sqrt{3}}{b\kappa ^2}}\right)^{1/2}.`$ (20)
Formula (20) correctly describes the steep diverging slope of $`m(h)\mathrm{}`$ or slopes $`m(b)^{}1`$ and $`h(b)^{}0`$ as $`b0`$ and is valid for $`0b2.5/\kappa ^2`$ for $`\kappa 7`$, Accidentally it also fits well $`m(b)`$ for $`\kappa 2`$ and $`b0.2`$, see the dotted lines in Fig. 8. A smoother fit is obtained by the exact London expression (C6) if one or three neighbor shells are included in the sum. But taking more terms in the sum improves the fit only at large $`\kappa `$. Accounting for neighbors up to $`\nu =100`$ lattice spacings apart (about 5000 terms) one gets good approximation to $`m`$ and $`h`$ for $`0b0.01`$ (0.02, 0.05) if $`\kappa 20`$ ($`\kappa =7`$, $`\kappa =2`$), see Fig.8. In the limit $`\nu \mathrm{}`$ the infinite sum (C6) reproduced Eq. (18), i.e, the dashed curves in Fig. 8 for $`\kappa =50,\mathrm{\hspace{0.17em}200}`$ then will straighten and cut the axis $`M=0`$ at $`b=\gamma =0.358`$ ($`b^{1/2}=0.60`$ in Fig. 7, $`b^{1/3}=0.71`$ in Fig. 8).

### III.4 General interpolation

All the approximations for $`m(b)`$ and $`h(b)`$ known so far, including the above formulas, fit either the low or high field regions. The formulas (15) and (20) \[or better, Eq. (C6) with the sum taken over three shells\] have a small overlap for all $`\kappa `$ and thus, together, they fit the entire range $`0b1`$ \[though the good fit of the low-$`\kappa `$ data by the London expression (20) or (C6) is accidental\].

For practical purposes one may construct interpolation formulas that approximate the numerically obtained magnetization in the entire range $`0<b<1`$. They should satisfy the five conditions

$`h(0)=h_{c1},h^{}(0)=0,h(1)=1,h^{}(1)=1p,`$
$`h^{\prime \prime }(1)=0,p=m^{}(1)=[(2\kappa ^21)\beta _A+1]^1,`$ (21)
with $`h_{c1}(\kappa )`$ from Eq. (4). A simple expression that satisfies all these conditions is

$`m(b,\kappa )=hb=p(1b)+(h_{c1}p)(1b)^\eta `$ (22)
with $`\eta (\kappa )=(1p)/(h_{c1}p)`$. Formula (22) approximates the exact $`m(b)`$ well for $`\kappa 2`$ with relative deviation $`|ϵ|<3\%`$, for $`\kappa =3`$ with $`2\%<ϵ<6\%`$, and for $`\kappa =5`$ with $`1\%<ϵ<16\%`$, see the dotted lines in Fig. 1, top.

For large $`\kappa `$, general interpolation formula are more difficult to construct because of the nonanalytic limiting expression Eq. (20). One may, however, combine the $`m_{\mathrm{low}}`$ from Eq. (20) with the $`m_{\mathrm{high}}`$ from Eq. (19) using a smooth transition at $`b(2\kappa ^2)^1`$, e.g., with weights $`1w`$ and $`w=\frac{1}{2}+\frac{1}{2}\mathrm{tanh}[\mathrm{\hspace{0.17em}2.5}(2b\kappa ^21)]`$, or slightly better, $`w=\frac{1}{2}+\frac{1}{2}\mathrm{erf}[\mathrm{\hspace{0.17em}2}(2b\kappa ^21)]`$, yielding

$`m(b,\kappa )=(1w)m_{\mathrm{low}}+wm_{\mathrm{high}}.`$ (23)
This interpolation between expressions (19) and (20) works well for $`0<b<1`$ with relative error $`|ϵ|<2\%`$ for $`\kappa 5`$ and $`3.5\%<ϵ<2\%`$ for $`\kappa 1`$, see the dotted lines in Fig. 1 (bottom). Thus, $`m`$ in the entire ranges of $`b`$ and $`\kappa `$ may be approximated by Eq. (22) or Eq. (23).

## IV Conclusion

The properties of the ideal periodic flux-line lattice in superconductors are calculated from Ginzburg-Landau theory for the entire ranges of GL parameters $`1/\sqrt{2}\kappa <\mathrm{}`$ and inductions $`0b=B/B_{c2}<1`$. The differences between the free energies and magnetizations of the triangular and square vortex lattices and the values obtained by the circular cell approximation are investigated in detail. Approximate analytical expressions are given for the variance $`\sigma (b,\kappa )`$ of the periodic induction and for the magnetization $`m(b,\kappa )`$. These limiting and interpolation formulas should replace previous approximate expressions that have rather limited validity.

The numerical methods presented in the appendices, in principle may be applied also to theories going beyond the isotropic GL theory considered here.

## Appendix A Isolated vortex and circular cell method

The calculation of the isolated flux line and of the FLL within the circular cell method, is a cylindrically symmetric problem. The free energy depends on the magnitude of the GL function $`f(r)`$ and on the magnetic induction $`B(r)`$ (along $`z`$) related to the vector potential $`A(r)`$ and supervelocity $`Q(r)`$ (along $`\phi `$) by

$$B=\frac{(Ar)^{}}{r}=\frac{(Qr)^{}}{r},Q=A\frac{1}{\kappa r}.$$
(24)
In reduced units $`\sqrt{2}H_c=\mu _0H_c^2=\lambda =1`$, the free energy of a flux line or of the circular cell with radius $`R`$ ($`\pi R^2=\mathrm{\Phi }_0/B`$) averaged over this cell and referred to the Meissner state ($`f=1`$, $`B=0`$) reads

$$F_{cc}=_0^R\left[\frac{(1f^2)^2}{2}+\frac{(f^{})^2}{\kappa ^2}+f^2Q^2+B^2\right]\frac{2\pi rdr}{\pi R^2}$$
(25)
with $`f^{}=df/dr`$. Minimizing the functional (A2) with respect to $`f(r)`$ and $`Q(r)`$ we obtain the two GL equations, which may be written in the form

$`f^{\prime \prime }+\kappa ^2f=\kappa ^2(2ff^3Q^2f)+f^{}/r,`$ (26)
$`B^{}=f^2Q=j,`$ (27)
where $`j=B^{}`$ is the current density. In Eq. (A3) a term $`\kappa ^2f`$ was added on both sides to improve the convergence of the iteration below. The boundary conditions are

$$f(0)=f^{}(R)=j(0)=j^{}(R)=0.$$
(28)
An appropriate ansatz in terms of Fourier series is

$`f(r)={\displaystyle \underset{m=1}{\overset{M}{}}}f_G\mathrm{sin}Gr,G={\displaystyle \frac{\pi (2m1)}{2R}},`$ (29)
$`A(r)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K\mathrm{sin}Kr+{\displaystyle \frac{r}{2}}B,K={\displaystyle \frac{\pi n}{R}},`$ (30)
$`B(r)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K{\displaystyle \frac{\mathrm{sin}Kr+Kr\mathrm{cos}Kr}{r}}+B,`$ (31)
$`Q(r)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K\mathrm{sin}Kr{\displaystyle \frac{1r^2/R^2}{\kappa r}},`$ (32)
$`j(r)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K{\displaystyle \frac{Kr\mathrm{cos}Kr(1+K^2r^2)\mathrm{sin}Kr}{r^2}}.`$ (33)
For the equidistant grid $`r_i=(i\frac{1}{2})R/N_r`$, $`i=1,2,\mathrm{}N_r`$, one has the orthogonality relation

$$\underset{i=1}{\overset{N_r}{}}\mathrm{sin}Gr_i\mathrm{sin}G^{}r_i=\frac{1}{2}N_r\delta _{GG^{}}$$
(34)
and similar equations for $`\mathrm{sin}Kr_i`$ and $`\mathrm{cos}Kr_i`$. The GL equations (A3) and (A4) thus may be written in the form of equations for the Fourier coefficients $`f_G`$ and $`a_K`$:

$`f_G={\displaystyle \frac{1}{G^2+\kappa ^2}}{\displaystyle \frac{2}{N_r}}{\displaystyle \underset{i=1}{\overset{N_r}{}}}\mathrm{sin}Gr_i\times `$
$`[\kappa ^2(2ff^3Q^2f)+f^{}/r_i],`$ (35)
$`a_K={\displaystyle \frac{1}{K^2+1}}[a_K+{\displaystyle \frac{2}{N_r}}{\displaystyle \underset{i=1}{\overset{N_r}{}}}\mathrm{sin}Kr_i\times `$
$`({\displaystyle \underset{n^{}=1}{\overset{N}{}}}a_K^{}{\displaystyle \frac{Kr\mathrm{cos}Kr\mathrm{sin}Kr}{r^2}}f^2Q)].`$ (36)
These two equations may be used to obtain the $`f_G`$ and $`a_K`$ by iteration, starting with appropriate initial values. The iteration becomes more stable and faster if the value of the previous iteration step is added with a certain weight $`(1c)<1`$, e.g., c=0.6, according to the algorithm:

$`f_G(1c)f_G+cF_G\{f,Q\},`$ (37)
$`a_K(1c)a_K+cA_K\{f,Q\},`$ (38)
with the symbols $`F_G\{f,Q\}`$ and $`A_K\{f,Q\}`$ denoting the right hand sides of Eqs. (A12) and (A13), respectively. Rapid convergence is achieved by iterating equations (A14), (A15) alternately. The equilibrium magnetic field $`H`$ is then obtained from Eq. (3), and the magnetization from

$`M={\displaystyle \frac{2}{BR}}{\displaystyle _0^R}\left[{\displaystyle \frac{f^4f^2}{2}}+B^2B(r)^2\right]r𝑑r.`$ (39)
At very large $`\kappa `$ and very small $`b`$ a large number $`N_r`$ of grid points $`r_i`$ is needed to achieve high accuracy, $`N_rR/\xi =R\kappa =\sqrt{2/b}`$. In this case the accuracy with a limited number of grid points may be improved by choosing a nonequidistant grid, e.g., $`r_i=u_i^2`$ with equidistant $`u_i=(i\frac{1}{2})\sqrt{R}/N_r`$. To use the orthogonality relations one then has to express $`f`$, $`B`$, and $`Q`$ as Fourier series in the new variable $`u=r^2`$ and also write the two GL equations in terms of the variable $`u`$, using, e.g., $`f^{}(r)=f^{}(u)/2u`$ and $`f^{\prime \prime }(r)=f^{\prime \prime }(u)/4u^2f^{}(u)/4u^3`$. This yields

$`f^{\prime \prime }(u)=4u^2\kappa ^2(f+f^3+Q^2f)+f^{}/u,`$ (40)
$`B^{}(u)=2uf^2Q,`$ (41)
and the Fourier series

$`f(u)={\displaystyle \underset{m=1}{\overset{M}{}}}f_G\mathrm{sin}Gu,G={\displaystyle \frac{\pi (2m1)}{2R}},`$ (42)
$`A(u)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K\mathrm{sin}Ku+{\displaystyle \frac{u^2}{2}}B,K={\displaystyle \frac{\pi n}{R}},`$ (43)
$`B(u)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K{\displaystyle \frac{2\mathrm{sin}Ku+Ku\mathrm{cos}Ku}{2u^2}}+B,`$ (44)
$`Q(u)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K\mathrm{sin}Ku{\displaystyle \frac{1u^4/R^2}{\kappa u^2}},`$ (45)
$`j(u)={\displaystyle \underset{n=1}{\overset{N}{}}}a_K{\displaystyle \frac{Ku\mathrm{cos}Ku(4+K^2u^2)\mathrm{sin}Ku}{4u^4}}.`$ (46)
The equations for the new Fourier coefficients are

$`f_G={\displaystyle \frac{1}{G^2+4\kappa ^2}}[4\kappa ^2f_G+{\displaystyle \frac{2}{N_r}}{\displaystyle \underset{i=1}{\overset{N_r}{}}}\mathrm{sin}Gu_i\times `$
$`(4u_i^2\kappa ^2f(1f^2Q^2)+f^{}/u_i)],`$ (47)
$`a_K={\displaystyle \frac{1}{K^2+1}}[a_K+{\displaystyle \frac{2}{N_r}}{\displaystyle \underset{i=1}{\overset{N_r}{}}}\mathrm{sin}Ku_i\times `$
$`({\displaystyle \underset{n^{}=1}{\overset{N}{}}}a_K^{}{\displaystyle \frac{Ku_i\mathrm{cos}Ku_i4\mathrm{sin}Ku_i}{u_i^2}}4u^2Qf^2)].`$ (48)
For better convergence a term $`4\kappa f_G`$ was added on both sides of Eq. (A17) to yield (A24). The corresponding iteration scheme using (A14), (A15) needs a smaller weight $`c`$ and more iteration steps, but for large $`\kappa ^2/b`$ it is faster than the first scheme since it needs less grid points $`N_r`$ to reach the same accuracy.

## Appendix B Periodic vortex lattice

The properties of the ideally periodic FLL within GL theory may be calculated by minimizing the GL free energy of the superconductor, Eq. (2), with respect to appropriate periodic trial functions, e.g., Fourier series with a large number of terms. For the smooth function $`\omega =f^2(𝐫)`$ we write the ansatz

$`\omega (𝐫)=f^2`$ $`=`$ $`{\displaystyle \underset{𝐊}{}}a_𝐊(1\mathrm{cos}\mathrm{𝐊𝐫})`$ (49)
with $`𝐫=(x,y)`$, $`𝐊=(K_x,K_y)`$. In all sums here and below the term $`𝐊=0`$ is excluded. For vortex positions $`𝐑=𝐑_{mn}=(mx_1+nx_2,ny_2)`$ the reciprocal lattice vectors are $`𝐊=𝐊_{mn}=(2\pi /S)(my_2,nx_1+mx_2)`$ with $`S=x_1y_2=\mathrm{\Phi }_0/B`$ the unit cell area and $`m,n=0,\pm 1,\pm 2,\mathrm{}`$ . For the triangular lattice one has $`x_2=x_1/2`$, $`y_2=x_1\sqrt{3}/2`$, and for the square lattice $`x_2=0`$, $`y_2=x_1`$. For supervelocity $`𝐐`$ and induction $`𝐁=\times 𝐐=B(𝐫)\widehat{𝐳}`$ we choose

$`B(𝐫)`$ $`=`$ $`B+{\displaystyle \underset{𝐊}{}}b_𝐊\mathrm{cos}\mathrm{𝐊𝐫},`$ (50)
$`𝐐(𝐫)`$ $`=`$ $`𝐐_A(𝐫)+{\displaystyle \underset{𝐊}{}}b_𝐊{\displaystyle \frac{\widehat{𝐳}\times 𝐊}{K^2}}\mathrm{sin}\mathrm{𝐊𝐫}.`$ (51)
Here $`𝐐_A(x,y)`$ is the supervelocity of the Abrikosov $`B_{c2}`$ solution, which satisfies

$`\times 𝐐_A=\left[B\mathrm{\Phi }_0{\displaystyle \underset{𝐑}{}}\delta _2(𝐫𝐑)\right]\widehat{𝐳},`$ (52)
where $`\delta _2(𝐫)=\delta (x)\delta (y)`$ is the 2D delta function. This relation shows that $`𝐐_A`$ is the velocity field of a lattice of ideal vortex lines but with zero average rotation. Close to each vortex center one has $`𝐐_A(𝐫)\widehat{𝐳}\times 𝐫^{}/(2\kappa r^2)`$ and $`\omega (𝐫)r^2`$ with $`𝐫^{}=𝐫𝐑`$. In principle $`𝐐_A(𝐫)`$ may be expressed as a slowly converging Fourier series by integrating (B4) using $`\mathrm{div}𝐐=\mathrm{div}𝐐_A=0`$ as in Ref. 15, . But it is more convenient to take $`𝐐_A`$ from the exact relation

$`𝐐_A(𝐫)={\displaystyle \frac{\omega _A\times \widehat{𝐳}}{2\kappa \omega _A}},`$ (53)
where $`\omega _A(x,y)`$ is the Abrikosov $`B_{c2}`$ solution given by the rapidly converging series (B1) with coefficients 30 ; 31

$`a_𝐊^A=(1)^{m+mn+n}\mathrm{exp}[K_{mn}^2S/(8\pi )]`$ (54)
for general lattice symmetry, and $`a_𝐊^A=(1)^{\nu ^2}\mathrm{exp}(\pi \nu ^2/\sqrt{3})`$ ($`\nu ^2=m^2+mn+n^2`$) for the triangular lattice. This $`\omega _A`$ is normalized to $`\omega _A(x,y)=1`$, which means that $`_𝐊^{}a_𝐊^A=1`$ for any lattice symmetry. Another strange property of the Abrikosov solution (B6) is that $`(\omega _A/\omega _A)^2^2\omega _A/\omega _A=4\pi /S=\mathrm{const}`$, although both terms diverge at the vortex positions; this relation follows from (B4) and (B5) using $`B=\mathrm{\Phi }_0/S=2\pi /(\kappa S)`$. The useful formula (B5) may be proven via the complex $`B_{c2}`$ solution $`\psi _A(x,y)`$; it means that near $`B_{c2}`$ the third and fourth term in $`F`$, Eq. (2), are identical.

Approximate solutions $`\omega (𝐫)`$ and $`B(𝐫)`$ may be computed by using a finite number of Fourier coefficients $`a_𝐊`$ and $`b_𝐊`$ and minimizing the free energy $`F(B,\kappa ,a_𝐊,b_𝐊)`$ with respect to these coefficients 15 . However, a much faster and more accurate solution method 17  is to iterate the two GL equations $`\delta F/\delta \omega =0`$ and $`\delta F/\delta 𝐐=0`$ written in appropriate form. Namely, the iteration is stable and converges rapidly if one isolates a term $`(^2+\mathrm{const})(\omega ,𝐐)`$ on the l.h.s. and puts the remaining terms to the r.h.s. as a kind of “inhomogeneity” of such London-like equations, e.g.,

$`(^2+2\kappa ^2)\omega `$ $`=`$ $`2\kappa ^2(2\omega \omega ^2\omega Q^2g),`$ (55)
$`(^2+\overline{\omega })𝐐_b`$ $`=`$ $`\omega 𝐐_A(\omega \overline{\omega })𝐐_b,`$ (56)
with the abbreviations $`g(𝐫)=(\omega )^2/(4\kappa ^2\omega )`$, $`𝐐_b=𝐐𝐐_A`$, $`\times 𝐐_b=B(𝐫)B`$, and $`\overline{\omega }=\omega =_𝐊^{}a_𝐊`$. Equations (B7), (B8) introduced some “penetration depths” $`(2\kappa ^2)^{1/2}=\xi /\sqrt{2}`$ and $`\overline{\omega }^{1/2}=\lambda /\overline{\omega }^{1/2}`$ (in real units), which stabilize the convergence of the iteration. Acting on the Fourier series $`\omega `$ (B1) and $`𝐐_b`$ (B3) the Laplacian $`^2`$ yields a factor $`K^2`$; this facilitates the inversion of (B7) and (B8). Using the orthonormality

$`2\mathrm{cos}\mathrm{𝐊𝐫}\mathrm{cos}𝐊^{}𝐫=\delta _{\mathrm{𝐊𝐊}^{}}`$ (57)
(for $`𝐊0`$) one obtains from (B1), (B2) $`a_𝐊=2\omega (𝐫)\mathrm{cos}\mathrm{𝐊𝐫}`$ and $`b_𝐊=2B(𝐫)\mathrm{cos}\mathrm{𝐊𝐫}`$. The convergence of the iteration is considerably improved by adding a third equation which minimizes $`F`$, Eq. (2), with respect to the amplitude of $`\omega `$, i.e., $`F/\overline{\omega }=0`$. This step gives the largest decrease of $`F`$. The resulting three iteration equations for the parameters $`a_𝐊`$ and $`b_𝐊`$ then read 17

$`a_𝐊`$ $`:=`$ $`{\displaystyle \frac{4\kappa ^2(\omega ^2+\omega Q^22\omega +g)\mathrm{cos}\mathrm{𝐊𝐫}}{K^2+2\kappa ^2}},`$ (58)
$`a_𝐊`$ $`:=`$ $`a_𝐊\omega \omega Q^2g/\omega ^2,`$ (59)
$`b_𝐊`$ $`:=`$ $`{\displaystyle \frac{2[(\omega \overline{\omega })B(𝐫)+p]\mathrm{cos}\mathrm{𝐊𝐫}}{K^2+\overline{\omega }}},`$ (60)
with $`p=(\omega \times 𝐐)\widehat{𝐳}=Q_x\omega /yQ_y\omega /x`$ and $`g=(\omega )^2/(4\kappa ^2\omega )=(f)^2/\kappa ^2`$ as above.

The solutions $`\omega (𝐫)`$, $`𝐁(𝐫)`$, and $`𝐐(𝐫)`$ are obtained by starting, e.g., with $`a_𝐊=(1b)a_𝐊^A`$ and $`b_𝐊=0`$ and then iterating the three equations (B10), (B11), (B12) by turns until the coefficients do not change any more. After typically 25 such triple steps, the solution stays constant to all 15 digits and the GL equations are exactly satisfied. Since all terms in (B10) - (B12) are smooth periodic functions of $`𝐫`$, high accuracy is achieved by using a regular spatial 2D grid, e.g., $`x_i=(i1/2)x_1/N_x`$ ($`i=1\mathrm{}N_x`$) and $`y_j=(j1/2)y_2/(2N_y)`$ ($`j=1\mathrm{}N_y`$, $`2N_yN_xy_2/x_1`$) with constant weights $`x_1/N_x`$ and $`y_2/(2N_y)`$. These $`N=N_xN_y`$ = 100 to 5000 grid points fill the rectangular basic area $`0xx_1`$, $`0yy_2/2`$, which is valid for any unit cell with the shape of a parallelogram. Spatial averaging $`\mathrm{}`$ then just means summing $`N`$ terms and dividing by $`N`$.

Best accuracy is achieved by considering all $`𝐊_{mn}`$ vectors within a half circle $`|𝐊_{mn}|K_{max}`$, with $`K_{max}^220N/S`$ chosen such that the number of the $`𝐊_{mn}`$ is slightly less than the number $`N`$ of grid points. The high precision of this method may be checked with the identity $`B(x,y)/B_{c2}=1\omega (x,y)`$, which is valid at $`\kappa =1/\sqrt{2}`$ for all $`b`$. This relation is confirmed with an error $`<10^9`$. The equilibrium field $`H`$ or reversible magnetization $`M=BH`$ is computed from Doria’s virial theorem, Eq. (3).

## Appendix C London theory

The modified London equation for a lattice of straight vortex lines at regular positions $`𝐑=𝐑_{mn}`$ (App. B) is

$`(1\lambda ^2^2)B(x,y)=\mathrm{\Phi }_0{\displaystyle \underset{𝐑}{}}\delta (𝐫𝐑_{mn}),`$ (61)
where $`\delta (x,y)`$ is the 2D delta function. The solution for the magnetic field of one isolated vortex at $`𝐑=0`$ is

$`B_v(r)=(\mathrm{\Phi }_0/2\pi \lambda ^2)K_0(r/\lambda ).`$ (62)
The modified Bessel function

$`K_0(r/\lambda )={\displaystyle \frac{d^2k}{2\pi }\frac{\mathrm{cos}\mathrm{𝐤𝐫}}{\lambda ^2+k^2}}`$ (63)
has the derivative $`K_0(x)^{}=K_1(x)`$ with the limits $`K_0(x1)\mathrm{ln}x`$, $`K_1(x1)1/x`$, and for $`x1`$: 32

$`K_0(x)`$ $``$ $`\sqrt{{\displaystyle \frac{\pi }{2x}}}e^x\left(1{\displaystyle \frac{1}{8x}}+{\displaystyle \frac{9}{128x^2}}{\displaystyle \frac{225}{3972x^3}}\right),`$
$`K_1(x)`$ $``$ $`\sqrt{{\displaystyle \frac{\pi }{2x}}}e^x\left(1+{\displaystyle \frac{3}{8x}}{\displaystyle \frac{15}{128x^2}}+{\displaystyle \frac{315}{3972x^3}}\right).`$ (64)
For a periodic FLL one obtains the Fourier series $`B(x,y)`$, Eq. (16), which may also be written as a sum over isolated vortex fields, $`B(x,y)=_𝐑B_v(𝐫𝐑)`$. Similarly, the free energy of the FLL may be written as a sum of vortex self energies ($`\mathrm{\Phi }_0H_{c1}`$ per unit length) plus a double sum over all interactions between two vortices. The average energy density $`F`$, Eq.(17), then reads

$`F=BH_{c1}+{\displaystyle \frac{B\mathrm{\Phi }_0}{4\pi \lambda ^2\mu _0}}{\displaystyle \underset{𝐑}{}}K_0(R/\lambda ).`$ (65)
For the triangular vortex lattice we write $`R/\lambda =\nu c`$ with $`c=a/\lambda =(4\pi /\sqrt{3})^{1/2}(b\kappa ^2)^{1/2}`$ ($`a`$ = vortex spacing) and $`\nu ^2=m^2+mn+n^2=`$ 1, 3, 4, 7, 7, 9, $`\mathrm{}`$ . Taking the derivative $`H=F/B`$ one obtains for $`h=H/H_{c2}`$ with $`h_{c1}=H_{c1}/H_{c2}`$:

$`h=h_{c1}+{\displaystyle \frac{3}{\kappa ^2}}{\displaystyle \underset{\nu }{}}\left[K_0(\nu c)+{\displaystyle \frac{\nu c}{2}}K_1(\nu c)\right].`$ (66)
Here the sum is over $`\nu =1,\sqrt{3},2,\mathrm{}`$ , i.e. the number of six flux lines per shell is already accounted for. Equation (C6) is still exact. It works for $`b1`$ (i.e. for nonoverlapping vortex cores) and for $`\kappa >1.4`$ (i.e. when the long-range interaction of vortices is purely magnetic 11 ; 33 ). With the expansions (C4) one obtains for $`x=\nu c1`$:

$`hh_{c1}+{\displaystyle \frac{3\sqrt{\pi }}{2\kappa ^2}}{\displaystyle \underset{\nu }{}}e^x\sqrt{x}\left[1+{\displaystyle \frac{19}{8x}}{\displaystyle \frac{47}{128x^2}}\right].`$ (67)
At very small $`b`$, namely for $`c=a/\lambda 1`$, the sum may be restricted to the nearest neighbor shell, i.e. to the first term, $`\nu =1`$, yielding

$`hh_{c1}+{\displaystyle \frac{3\sqrt{\pi c}}{2\kappa ^2}}e^c\left[1+{\displaystyle \frac{19}{8c}}{\displaystyle \frac{47}{128c^2}}\right].`$ (68)