text:
# Multiband superconductivity in NbSe2 from heat transport

## Abstract

The thermal conductivity of the layered $`s`$-wave superconductor NbSe2 was measured down to $`T_c`$/100 throughout the vortex state. With increasing field, we identify two regimes: one with localized states at fields very near $`H_{c1}`$ and one with highly delocalized quasiparticle excitations at higher fields. The two associated length scales are most naturally explained as multi-band superconductivity, with distinct small and large superconducting gaps on different sheets of the Fermi surface.

, , , , , , , , , ,

thanks: Permanent address: Inst. Surf. Chem., N.A.S. Ukraine.thanks: Current address: Department of Physics, University of Sherbrooke, Sherbrooke, Canada J1K 2R1

In the classical theory of superconductivity all electrons on the Fermi surface contribute equally to the superconducting pairing, giving a constant superconducting gap $`\mathrm{\Delta }`$. The difference of $`\mathrm{\Delta }`$ on different sheets of the Fermi surface, or multiband superconductivity (MBSC), considered theoretically back in the 50s , has emerged recently as a possible explanation for the unusual properties of MgB2 .

Based on the observation of a sizable difference of $`\mathrm{\Delta }`$ on two Fermi surface sheets by angle resolved photoemission spectroscopy, it has been suggested that the layered superconductor NbSe2 is also host to MBSC . However, these surface-sensitive measurements were performed only down to 5.3 K, close to $`T_c`$=7.0 K. We present evidence of bulk MBSC at low temperatures in this compound.

The thermal conductivity $`\kappa `$ of pure samples of NbSe2 (residual resistivity $`\rho _0`$=3 $`\mu \mathrm{\Omega }`$ cm) was measured upon warming in a magnetic field ($`Hc`$ axis) . In Fig. 1, $`\kappa `$ is presented as $`\kappa /T`$ vs $`T^2`$, enabling a separation of the electronic contribution, $`\kappa _0`$ $`T`$, and the phononic contribution, $`\kappa ^gT^3`$. In zero field, $`\kappa _0/T`$ goes to zero in the $`T`$=0 limit as expected for $`s`$-wave superconductors, thus ruling out the possibility of nodes in the gap at any point on the Fermi surface. With the application of a magnetic field, $`\kappa _0/T`$ rapidly increases (Fig. 2, main panel) up to $`H_{c2}`$, following closely the increase of the electronic specific heat $`\gamma `$ . This behavior is very different from that expected for conventional superconductors. There, electronic quasiparticles remain localized inside the vortex cores, such that the increase in quasiparticle density (seen as an increase of $`\gamma `$ with $`H`$) is not followed by an increase of $`\kappa _0`$. This is shown for V3Si in the upper panel of Fig. 2 .

A closer examination of the field dependence of $`\kappa _0`$ and of $`\gamma `$ at fields close to $`H_{c1}`$ is consistent with the presence of localized states, another indication of a gap without nodes. Indeed, in a limited field range, $`\kappa _0`$ shows an activated increase in field whereas the specific heat increases rapidly above $`H_{c1}`$ (lower inset in Fig. 2).

In Fig. 3, we show a more detailed comparison of $`\kappa `$ and $`\gamma `$ by plotting the ratio $`\kappa _0/\gamma `$. In a rough sense, it represents the ”mobility” of quasiparticle excitations in the vortex state. It is clear that the field dependence of $`\kappa _0/\gamma `$ shows a strong change of slope at $`H^{}H_{c2}/9`$, supporting the existence of a field scale which plays the same role for quasiparticle localization as $`H_{c2}`$ does in standard $`s`$-wave superconductors. In the inset of Fig. 3, we show that the low field behavior of NbSe2 can be matched to that of V3Si by using $`H^{}H_{c2}`$/9.

While $`\xi `$ is associated with the upper critical field ($`\xi ^2=\mathrm{\Phi }_0/2\pi H_{c2}`$), a second length scale $`\xi ^{}`$ must be associated with $`H^{}`$. The crossover between $`\xi `$ and $`\xi ^{}`$ with $`H`$ can naturally explain the shrinking of the vortex cores observed by muon spin relaxation in NbSe2 . Moreover, the superconducting coherence length is is related to $`\mathrm{\Delta }`$ via $`\xi v_F/\mathrm{\Delta }`$, where $`v_F`$ is the Fermi velocity. In NbSe2, $`v_F`$ is approximately the same for different sheets of the Fermi surface, such that the ratio $`\xi ^{}/\xi `$ 3 gives a ratio of associated superconducting gaps $`\mathrm{\Delta }/\mathrm{\Delta }^{}`$ 3. This value is consistent with previous heat capacity and tunneling results .

In conclusion, we have identified the anomalous evolution of thermal conductivity in NbSe2 versus magnetic field with the existence of two length scales in the superconducting state. This finds a natural explanation in a model of multi-band superconductivity, assuming a ratio of larger gap to smaller gap of about 3.

This work was supported by the Canadian Institute for Advanced Research and funded by NSERC.


format: text

text:
# Rotational states of Bose gases with attractive interactions in anharmonic traps

Present address: \]Department of Physics, KTH, SE-10691 Stockholm

## Abstract

A rotated and harmonically trapped Bose gas with attractive interactions is expected to either remain stationary or escape from the trap. Here we report that, on the contrary, in an anharmonic trapping potential the Bose gas with attractive interactions responds to external rotation very differently, namely through center-of-mass motion or by formation of vortices, depending on the strengths of the interaction and the anharmonic term in the trapping potential.

Trapped ultracold gases of atoms with attractive s-wave interactions offer quite different physics compared to the more common case of repulsive interactions pethicksmith . The behavior under rotation is an important example. While repulsive Bose-Einstein condensed gases respond to rotation by forming quantized vortices, such states are not necessarily energetically favorable in attracting gases. In fact, Wilkin et al. wilkin  have studied a harmonically trapped gas of attracting bosons under rotation in the limit where the interactions can be treated perturbatively. They find that in the lowest-energy many-body state for a given angular momentum (the yrast state) the center of mass (CM) carries all the angular momentum. According to Wilkin et al., for this state the single-particle density matrix has more than one eigenvalue of macroscopic size, which means that more than one single-particle state is macroscopically occupied. Such a state is called a fragmented condensate nozieres . The form of the lowest-energy state has since been proved in greater generality hussein . Pethick and Pitaevskii pethickpitaevskii  made the important observation that the rotational state in question can be seen as a Bose-Einstein condensate if it is viewed in the CM reference frame. Nevertheless, it seems that this interesting state of motion is never thermodynamically stable, because the critical frequency is equal to the frequency for destabilization of the cloud. On the contrary, in this Letter we show that an anharmonic trapping potential makes the state of CM rotation, as well as vortex states, attainable.

We consider a trapped, attractively interacting Bose gas at zero temperature, subject to a rotational force with a given angular velocity, and study the states that minimize the energy under such conditions. To this end, we first derive an equation of motion governing such a gas, taking into account the CM motion. We then study numerically and variationally the behavior in an anharmonic trapping potential, and show that the parameter space is divided into regions of CM motion and vortex formation.

The system under study consists of $`N`$ bosons subject to the many-body Hamiltonian

$$H=\frac{\mathrm{}^2}{2m}\underset{i=1}{\overset{N}{}}\frac{d^2}{d𝐫_i^2}+\underset{i=1}{\overset{N}{}}V(𝐫_i)+\frac{U_0}{2}\underset{ij}{}\delta (𝐫_i𝐫_j).$$
(1)
The external potential $`V(𝐫)`$ can in experiment be of either magnetic or optical origin, and is usually taken to be a harmonic-oscillator well. But the harmonic potential has special properties that simplify the problem but in fact prevent rotation, and we therefore study the more general case of a cylindrically symmetric trap that is anharmonic in the radial direction:

$$V(r,\theta ,z)=\frac{1}{2}m\left[\omega ^2(r^2+\lambda \frac{r^4}{a_{\mathrm{osc}}^2})+\omega _z^2z^2\right].$$
(2)
The oscillator length $`a_{\mathrm{osc}}`$ is defined as $`a_{\mathrm{osc}}=(\mathrm{}/m\omega )^{1/2}`$. The third term in the Hamiltonian (1) is the inter-atomic interaction potential which is here replaced by its usual s-wave approximation; the strength $`U_0=4\pi \mathrm{}^2a/m`$ depends on the scattering length $`a`$ which in this study will be assumed negative, implying an attractive interaction.

For the case of purely harmonic trapping ($`\lambda =0`$), the state of CM rotation was discussed in Refs. wilkin ; hussein ; mottelson . The energy of $`L`$ quanta of angular momentum carried by the CM is exactly $`L\mathrm{}\omega `$, independent of interactions (fetterrokhsar , Appendix A; cf. kohn ; elliottskyrme ). This has the important consequence that the gas cannot be rotated, because the critical frequency for the excitation of any rotation is equal to the frequency where the cloud is destabilized. The energy of a system in CM rotation as seen in a frame rotating with angular velocity $`\mathrm{\Omega }`$ is, namely,

$$E\mathrm{\Omega }L_z=E_0+\mathrm{}(\omega \mathrm{\Omega })L,$$
(3)
where $`E_0`$ is the energy of the cloud in the absence of rotation. The energy, Eq. (3), is minimized for $`L=0`$ when $`\mathrm{\Omega }<\omega `$ and for $`L\mathrm{}`$ when $`\mathrm{\Omega }>\omega `$. Consequently, a harmonically trapped attracting Bose gas cannot be set into steady rotational motion: it either stays still, for slower stirring, or escapes from the trap for faster stirring.

In order to stabilize a finite-angular momentum state one may add an anharmonic term to the trapping potential. In such a potential, the single-particle energy as a function of angular momentum increases faster than linear, and the energy in a rotating frame may therefore have a minimum for a finite value of $`L`$. However, there is a priori no guarantee that the lowest-energy state is that of CM rotation when the trap is anharmonic. In fact, for a noninteracting gas in an anharmonic trap it is known that the lowest-energy state for given angular momentum $`L_z=Nm_l\mathrm{}`$ is the state where all particles are in the same angular momentum eigenstate, in other words, a Bose-Einstein condensate containing a vortex with quantum number $`m_l`$ lundh2002 . This is so because the nonlinear dependence of the energy on angular momentum implies that an angular-momentum eigenstate has lower energy than any superposition of different states with the same expectation value of angular momentum. An intuitive picture can be constructed as follows: in a frame rotating with a fixed angular velocity $`\mathrm{\Omega }`$, the sum of the centrifugal and potential energies create a Mexican-hat shaped potential if the angular velocity is high enough. In the absence of interactions, the lowest-energy state is one of toroidal shape lying along the minimum of the effective potential. When the interactions are attractive, however, the tendency for the density to be maximized may favor the state of CM rotation. It remains to investigate the crossover between these two types of rotational state.

We now put up an ansatz that allows for the derivation of a simple set of equations of motion of an attractive gas. In a harmonically confined gas with isotropic interactions, the CM and relative motions decouple and therefore the wave function of the $`N`$-body system can be written in the form

$$\mathrm{\Psi }(𝐫_1,\mathrm{},𝐫_N)=\psi _C(𝐑)\psi _R(𝐫_1^{},\mathrm{},𝐫_N^{}),$$
(4)
where $`𝐑=(𝐫_1+\mathrm{}+𝐫_N)/N`$ is the CM coordinate and $`𝐫_i^{}=𝐫_i𝐑`$ are the coordinates of each particle relative to the CM, $`i=1,\mathrm{},N`$ pethickpitaevskii ; elliottskyrme . Furthermore, the relative part $`\psi _R`$ exhibits Bose-Einstein condensation, so that the associated density matrix has one macroscopic eigenvalue pethickpitaevskii . Consequently, the many-body wave function can be written on the form

$$\mathrm{\Psi }(𝐫_1,\mathrm{},𝐫_N)=\psi _C(𝐑)\underset{i=1}{\overset{N}{}}\chi (𝐫_i^{}),$$
(5)
as long as the depletion of the condensate is neglected. Although Eqs. (4-5) hold strictly only in the case of a harmonic trap, we shall assume that the decoupling between CM and internal motion holds approximately also in the anharmonic case and use Eq. (5) as an ansatz for the derivation of an equation of motion that can be used to study the rotational properties of the attractive gas. We expect the resulting equation to be approximately valid when $`\lambda `$ is small so that the trap does not deviate much from harmonic.

The interaction energy is independent of the center-of-mass motion, simply because

$$𝐫_i𝐫_j=𝐫_i^{}𝐫_j^{}.$$
(6)
The harmonic part of the trap potential energy can be expressed in the relative coordinates according to

$$\underset{i}{}x_i^2=NX^2+\underset{i}{}x_{i}^{}{}_{}{}^{2},$$
(7)
and equivalently for the other Cartesian components. The results (6-7) together imply that in the harmonic case, the CM motion is decoupled from the internal motion, as discussed above. The quartic part of the potential expressed in CM and relative coordinates is

$`{\displaystyle \underset{i=1}{\overset{N}{}}}r_i^4=NR^4+{\displaystyle \underset{i=1}{\overset{N}{}}}\left[2R^2r_{i}^{}{}_{}{}^{2}+4(Xx_i^{}+Yy_i^{})^2+4r_{i}^{}{}_{}{}^{2}(Xx_i^{}+Yy_i^{})+r_{i}^{}{}_{}{}^{4}\right].`$ (8)
The equations of motion are to be obtained in a variational procedure, slightly generalizing the derivation given in Refs. pethicksmith ; leggett . We express the energy of the $`N`$-body system in terms of the trial wave function, Eq. (5),

$$E=𝑑𝐫_1\mathrm{}𝑑𝐫_N\mathrm{\Psi }^{}(𝐫_1,\mathrm{},𝐫_N)H\mathrm{\Psi }(𝐫_1,\mathrm{},𝐫_N),$$
(9)
where $`H`$ is the Hamiltonian, Eq. (1). The transformation to the new system of $`N+1`$ coordinates, $`𝐑,𝐫_1^{},\mathrm{},𝐫_N^{}`$, introduces a Jacobian $`J=N^d`$ in $`d`$ dimensions and a delta function $`\delta (𝐫_1^{}+𝐫_2^{}+\mathrm{}+𝐫_N^{})`$ in the integral, because the relative coordinates are not independent. Neglecting this interdependence of the coordinates introduces only errors of order $`N^1`$, however, so that the $`N+1`$ variables can in effect be regarded as independent and the functions $`\psi _C`$ and $`\chi `$ can be assumed to be individually normalized to unity. The trap potential term is rewritten according to Eqs. (7-8), and the action of the kinetic-energy operator on the trial wave function is readily calculated. With a number of superfluous factors $`\chi (𝐫_i^{})`$ integrated out, the energy functional is

$`E={\displaystyle }Jd𝐑d𝐫^{}\psi _C^{}(𝐑)\chi ^{}(𝐫^{})[{\displaystyle \frac{\mathrm{}^2}{2Nm}}{\displaystyle \frac{^2}{𝐑^2}}+{\displaystyle \frac{Nm}{2}}(\omega ^2R^2+\omega _z^2Z^2)+`$
$`N\lambda {\displaystyle \frac{m\omega ^2}{2a_{\mathrm{osc}}^2}}\left(R^4+2R^2r_{}^{}{}_{}{}^{2}+4(Xx^{}+Yy^{})^2+4r_{}^{}{}_{}{}^{2}(Xx^{}+Yy^{})+r_{}^{}{}_{}{}^{4}\right)`$
$`(N1){\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2}{𝐫^2}}+N{\displaystyle \frac{m}{2}}(\omega ^2r_{}^{}{}_{}{}^{2}+\omega _z^2z_{}^{}{}_{}{}^{2})+N(N1){\displaystyle \frac{U_0}{2}}|\chi (𝐫^{})|^2]\psi _C(𝐑)\chi (𝐫^{}),`$ (10)
where $`J=N^3`$ as noted above. The energy is now functionally minimized with respect to the two wave functions $`\psi _C`$ and $`\chi `$. The result is the two coupled equations

$`\left\{{\displaystyle \frac{\mathrm{}^2}{2M}}{\displaystyle \frac{^2}{𝐑^2}}+V_C(𝐑)\right\}\psi _C(𝐑)=E_C\psi _C(𝐑),`$ (11)
$`\left\{{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2}{𝐫^2}}+V_R(𝐫^{})+U_0|\mathrm{\Phi }(𝐫)|^2\right\}\mathrm{\Phi }(𝐫)=\mu \mathrm{\Phi }(𝐫).`$ (12)
We have introduced here the condensate wave function $`\mathrm{\Phi }=\sqrt{N}\chi `$ and dropped the prime on the relative coordinate $`𝐫`$. The effective potentials for the center-of-mass and relative motion are defined as

$`V_C(R)={\displaystyle \frac{1}{2}}M\omega ^2(1+4\lambda {\displaystyle \frac{r^2}{a_{\mathrm{osc}}^2}})R^2+{\displaystyle \frac{\lambda M\omega ^2}{2a_{\mathrm{osc}}^2}}R^4,`$ (13)
$`V_R(r)={\displaystyle \frac{1}{2}}m\omega ^2(1+4\lambda {\displaystyle \frac{R^2}{a_{\mathrm{osc}}^2}})r^2+{\displaystyle \frac{\lambda m\omega ^2}{2a_{\mathrm{osc}}^2}}r^4,`$ (14)
where $`R^2=J𝑑𝐑R^2|\psi _C(𝐑)|^2`$ and $`r^2=N^1𝑑𝐫r^2|\mathrm{\Phi }(𝐫)|^2`$. Averages of odd powers in the relative coordinate vanish because the relative coordinates are centered around the CM. The quantities $`E_C`$ and $`\mu `$ are introduced as Lagrange multipliers to ensure separate normalization; $`E_C`$ is actually the center-of-mass energy and $`\mu `$ is the chemical potential as usual, and $`M=Nm`$ is the total mass of the gas. For harmonic trapping ($`\lambda =0`$), the equations decouple: the CM behaves like a single particle of mass $`M`$, and the condensate wave function for the internal motion obeys the Gross-Pitaevskii equation. The same effective potential for the internal motion was derived in Ref. zimmermann ; however, the Gross-Pitaevskii equation was there taken as a starting point.

In the absence of rotation, the spatial extent of the CM wave function is on the order of the effective CM oscillator length $`a_C=(\mathrm{}/M\omega )^{1/2}`$ while the size of the condensate is of order $`a_{\mathrm{osc}}`$, which is a factor $`N^{1/2}`$ larger. As a result, the CM motion has a negligible effect on the motion of the condensate, and the usual Gross-Pitaevskii equation is retained. However, we are interested in the case of highly excited CM motion, where the total angular momentum $`L`$ carried by the CM is comparable to the number of particles $`N`$, and in that limit, the effects of the CM wave function cannot be neglected.

Numerical solutions to the coupled equations of motion (11-14) are obtained with the steepest descent method. The $`z`$ dependence of the CM wave function $`\psi _C`$ decouples and is exactly Gaussian, while for the condensate wave function $`\mathrm{\Phi }`$ both the radial and axial dependence is subject to numerical minimization. We have chosen the harmonic part of the potential to be isotropic, $`\omega _z=\omega `$. For large values of $`L`$, the CM wave function becomes very peaked and difficult to handle numerically. It has therefore been necessary to restrict the number of particles $`N`$ to $`N15`$, and the possible angular momenta $`L`$ to a few times $`N`$.

The result for the total angular momentum is shown in Fig. 1. We have chosen $`Na/a_{\mathrm{osc}}=0.1`$, and $`\lambda =0.05`$ and $`\lambda =0.1`$ as two illustrative examples. Both these parameter values lie in the region where for slow rotation the CM carries the motion, hence the total angular momentum per particle $`L_z/N`$ increases only in steps of $`\mathrm{}N^1`$, but for faster rotation the angular momentum is carried by a vortex in the condensate wave function, and the angular momentum per particle is integer. For sufficiently large $`\lambda `$ and weak attraction, vortex states are energetically favorable for all $`\mathrm{\Omega }`$, because the anharmonic potential favors vortex states where the density distribution is toroidal. Conversely, for large enough $`Na/a_{\mathrm{osc}}`$ and small enough $`\lambda `$, the attraction between bosons favors the state of CM rotation where the density distribution is maximally concentrated.

The numerical results give at hand that a mixed state, where both $`m_l`$ and $`L`$ are nonzero, is never the energetically favorable situation. This is clear, since the density distribution of such a state is a toroid that is not centered at the middle of the trap, but occupies regions of high potential energy in the Mexican-hat shaped effective potential formed by the anharmonic and centrifugal energies.

In order to be able to cover larger parameter spaces, we must consider a variational solution of the equations of motion (11-14). We choose trial functions with Gaussian envelope for the CM and condensate wave functions:

$`\psi _C(𝐑)`$ $`=`$ $`𝒩_C(X+iY)^Le^{R^2/2ba_C^2(\omega _z/\omega )Z^2/2a_C^2},`$
$`\mathrm{\Phi }(𝐫)`$ $`=`$ $`𝒩_R(x+iy)^{m_l}e^{r^2/2ca_{\mathrm{osc}}^2(\omega _z/\omega )z^2/2(da_{\mathrm{osc}})^2}.`$ (15)
That this is a good ansatz has been confirmed by inspecting the shape of the numerically calculated wave functions. $`𝒩_C`$ and $`𝒩_R`$ are normalization constants. As already noted, the $`z`$ dependence of the CM wave function is not subject to variational minimization, but is a Gaussian with the fixed width $`a_C\sqrt{\omega /\omega _z}`$. On the other hand, the widths $`b`$, $`c`$ and $`d`$ are variational parameters. The ansatz allows for the angular momentum to be shared between the CM and internal motion and thus incorporates both the known limits of a BEC in a vortex state for $`U_0=0`$, and CM rotation for $`\lambda =0`$. The total angular momentum is $`L_z=(L+Nm_l)\mathrm{}`$; we have already seen that the lowest-energy states in a rotating frame have either $`L=0`$ or $`m_l=0`$. Defining $`q=L/N`$, the energy is

$`E={\displaystyle \frac{N}{2}}q\mathrm{}\omega \left({\displaystyle \frac{1}{b}}+b+4(m_l+1)\lambda bc+\lambda qb^2\right)`$
$`+{\displaystyle \frac{N}{2}}(m_l+1)\mathrm{}\omega \left[{\displaystyle \frac{1}{c}}+c+\lambda (m_l+2)c^2\right]+{\displaystyle \frac{N}{4}}\mathrm{}\omega _z\left({\displaystyle \frac{1}{d^2}}+d^2\right)`$
$`+N\mathrm{}\omega \left({\displaystyle \frac{\omega }{2\pi \omega _z}}\right)^{1/2}{\displaystyle \frac{Na}{a_{\mathrm{osc}}}}{\displaystyle \frac{(2m_l)!}{(m_l!)^22^{2m_l}}}{\displaystyle \frac{1}{cd}}.`$ (16)
Terms of order unity are discarded. The energy is now minimized with respect to the parameters $`b`$, $`c`$ and $`d`$, and the energetically optimal way of distributing the total angular momentum between the CM and internal motion is determined. The result is plotted as solid lines in Fig. 1. It should be noted that the variational solution is carried out for the limit of large $`N`$ while the full numerical solution is restricted to finite values of $`N`$.

Fig. 2 shows the boundary between the CM-rotation and vortex regimes for a few fixed angular momenta as calculated in the variational scheme. As expected, a large $`\lambda `$, large angular momentum and weak coupling favors the vortex state. The coupling at which the cloud becomes unstable against collapse is also shown; the apparent shrinking of the stable regime with increasing $`\lambda `$ is due to the decrease of the effective size of the cloud, and therefore an increase in density, when $`\lambda `$ is increased at fixed $`\omega `$. Included is also the result of a perturbative calculation, treating the exact many-body wavefunctions to first order in $`\lambda `$ and $`Na/a_{\mathrm{osc}}`$ (cf. wilkin ; mottelson ); it is seen that perturbation theory only works when these quantities are very small.

In conclusion, a Bose gas with negative s-wave scattering length confined in an anharmonic trap will responds to a rotational force by either exciting center-of-mass rotational motion or forming a multiply quantized vortex, depending on the strength of the interaction and the anharmonic term. The wave function of the CM position and the condensate wave function associated with the internal motion are described by two equations that in the general anharmonic case are coupled. These results show that the behavior of condensates with attractive interactions is markedly dependent on the harmonicity of the potential, and that an anharmonic trapping potential is expected to allow for vortex creation. Reference dalibard  reported experiments using an anharmonic potential with a small quartic term ($`\lambda 0.001`$), created by imposing a Gaussian optical potential over a magnetic trap. A similar setup, possibly with a stronger anharmonicity through increased laser power or decreased beam width, could be employed in order to explore the rotational states discussed in this paper.

While preparing this manuscript we became aware of the related work by Ghosh ghosh .

Discussions with C. J. Pethick are gratefully acknowledged. The authors acknowledge support from the Academy of Finland (Project 50314), the European Network “Cold Atoms and Ultra-Precise Atomic Clocks” (CAUAC), and the Magnus Ehrnrooth foundation.


text:
# Breached pairing superfluidity: Possible realization in QCD

## BCS superfluidity.

For $`\delta p_F\frac{2\mathrm{\Delta }}{\sqrt{V^uV^s}}`$, $`E_𝐩^+`$ never intersects $`E_𝐩^+=0`$. Then $`E_𝐩^+>0`$ and $`E_𝐩^{}<0`$ for all $`𝐩`$. Pairing is possible in the whole range $`(p_0\lambda )|𝐩|(p_0+\lambda )`$. We arrive at the BCS solution for the gap

$$\mathrm{\Delta }=\lambda (V^u+V^s)\mathrm{exp}\left(\frac{1}{N_+(0)g}\right),$$
(11)
where $`N_+(0)=\frac{1}{(2\pi )^3}d^3𝐩\delta (ϵ_𝐩^+)`$ is the density of states, and $`\lambda `$ is the ultraviolet cutoff. Since $`N_+(0)p_0^2`$, the gap parameter decreases as a function of increasing $`\delta p_F`$ with $`p_F^u`$ fixed.

## Breached pairing (BP) superfluidity.

This occurs when $`\delta p_F>\frac{2\mathrm{\Delta }}{\sqrt{V^uV^s}}`$. $`E_𝐩^+=0`$ has two roots. We integrate below $`p_\mathrm{\Delta }^{}`$ and above $`p_\mathrm{\Delta }^+`$, but exclude the region $`p_\mathrm{\Delta }^{}|𝐩|<p_\mathrm{\Delta }^+`$. A solution of the gap equation exists only if $`g>g_c`$, with the critical coupling constant

$$g_c=\left[N_+(0)\mathrm{ln}\left(\frac{\lambda }{\delta p_F}\frac{(V^u+V^s)}{\sqrt{V^uV^s}}\right)\right]^1.$$
(12)
For $`V^s/V^u1`$ with fixed Fermi momenta, $`g_c0`$ becomes arbitrarily weak.

Expanding the gap equation (10) for small $`\mathrm{\Delta }<\delta p_F`$, we find the critical behavior

$$\mathrm{\Delta }(gg_c)^{\frac{1}{2}}$$
(13)
with $`g_c`$ defined in Eq. (12).

The quasiparticle excitations of the breached pairing phase differ qualitatively from those of BCS. In the BCS phase, $`E_𝐩^\pm `$ branches are both gapped with their minimum energies given by $`|ϵ_{p_0}^{}\pm 2\mathrm{\Delta }\sqrt{V_uV_s}/(V_u+V_s)|0`$, respectively. By contrast, in the breached pairing phase $`E_𝐩^+`$ intersects zero energy. Only the $`E_p^{}`$ branch is fully gapped.

For finite $`\mathrm{\Delta }`$, the gap equation (10) needs to be solved self-consistently, implementing the restriction in integration region, $`|𝐩|p_\mathrm{\Delta }^{}`$ and $`|𝐩|p_\mathrm{\Delta }^+`$, which depends on $`\mathrm{\Delta }`$. To this end, we transform the integral gap equation into an algebraic one by introducing a (BCS) gap paramter, $`\mathrm{\Delta }_0\mathrm{\Delta }(\delta p_F=0)`$. The BCS gap $`\mathrm{\Delta }_0`$ implicitly defines the interaction strength, $`1/g=N_+(0)\mathrm{ln}\left[\lambda (V^u+V^s)/\mathrm{\Delta }_0\right]`$. We shall use $`\mathrm{\Delta }_0`$ to parameterize the interaction strength below. Performing the integrals

$`{\displaystyle _{p_0\lambda }^{p_0+\lambda }}{\displaystyle \frac{dp}{\sqrt{\epsilon _𝐩^{+\mathrm{\hspace{0.17em}2}}+\mathrm{\Delta }_0^2}}}`$ (14)
$`=`$ $`{\displaystyle _{p_0\lambda }^{p_0+\lambda }}{\displaystyle \frac{dp}{\sqrt{\epsilon _𝐩^{+\mathrm{\hspace{0.17em}2}}+\mathrm{\Delta }^2}}}{\displaystyle _{p_\mathrm{\Delta }^{}}^{p_\mathrm{\Delta }^+}}{\displaystyle \frac{dp}{\sqrt{\epsilon _𝐩^{+\mathrm{\hspace{0.17em}2}}+\mathrm{\Delta }^2}}}`$
where $`\epsilon _𝐩^+=\frac{1}{2}(V^u+V^s)(pp_0)`$, we obtain

$$\frac{\mathrm{\Delta }_0^2}{\mathrm{\Delta }^2}=\frac{\epsilon _{p_\mathrm{\Delta }^+}^s}{\epsilon _{p_\mathrm{\Delta }^{}}^s}$$
(15)
for $`\lambda \mathrm{}`$. After some elementary algebra, the above equation reduces to

$$\frac{\mathrm{\Delta }_0^2}{\mathrm{\Delta }^2}=\frac{\delta p_F+\sqrt{\delta p_F^2\frac{4\mathrm{\Delta }^2}{V^uV^s}}}{\delta p_F\sqrt{\delta p_F^2\frac{4\mathrm{\Delta }^2}{V^uV^s}}}.$$
(16)
It has a (stable) solution only if

$$\delta p_F<\frac{\mathrm{\Delta }_0}{\sqrt{V^uV^s}},$$
(17)
which is equivalent to $`g>g_c`$ (see Eq. (12)). In solving Eq. (16), we distinguish three cases:

## I: Chemical potentials fixed.

Then $`\delta p_F`$ remains unchanged before and after pairing, and the gap equation permits two solutions:

$$\mathrm{\Delta }=\{\begin{array}{cc}\mathrm{\Delta }_0,\hfill & \text{(BCS)}\hfill \\ \left[\mathrm{\Delta }_0\sqrt{V^uV^s}\left(\delta p_F\delta p_{\mathrm{th}}\right)\right]^{\frac{1}{2}},\hfill & \text{(unstable BP)}\hfill \end{array}$$
(18)
where $`\delta p_{\mathrm{th}}=\frac{\mathrm{\Delta }_0}{\sqrt{V^uV^s}}`$ (threshold value). In this case, the BP state always has higher energy than the normal state, so it is unstable (see Eq. (24)). Fig. 1a makes it further clear that even if it is a metastable solution at $`\delta p_F>\delta p_{\mathrm{th}}`$, the BP state always has a smaller $`\mathrm{\Delta }`$ than BCS. It is therefore energetically disfavored. This case was considered by Sarma and others Sar  a long time ago.

## II: Overall quark density fixed.

In this case one has to adjust the overall chemical potential, $`\mu _{oa}\frac{1}{2}(V^up_F^u+V^sp_F^s)=\frac{1}{2}(V^u+V^s)p_0`$, in order to accommodate a fixed overall density when pairing occurs. We hold the relative chemical potential, $`\delta \mu =V^up_F^uV^sp_F^s`$ fixed. The mismatched Fermi momentum is then altered, $`\delta p_F\delta \stackrel{~}{p}_F`$, with ‘$``$’ denoting the new Fermi momenta $`\stackrel{~}{p}_F^{u,s}`$ after pairing.

We find $`\delta \stackrel{~}{p}_F`$ is determined by

$$p_0^3p_0^2\frac{\epsilon _{p_F^u}^{}\epsilon _{p_F^s}^{}}{V^u+V^s}=\stackrel{~}{p}_0^3\stackrel{~}{p}_0^2\frac{\stackrel{~}{\epsilon }_{\stackrel{~}{p}_\mathrm{\Delta }^+}^{}\stackrel{~}{\epsilon }_{\stackrel{~}{p}_\mathrm{\Delta }^{}}^{}}{V^u+V^s},$$
(19)
where $`\stackrel{~}{p}_0=p_0+\frac{2V^uV^s}{V^{u\mathrm{\hspace{0.17em}2}}V^{s\mathrm{\hspace{0.17em}2}}}(\delta p_F\delta \stackrel{~}{p}_F).`$ For large Fermi momenta compared to the mismatch in Fermi surfaces, $`\delta p_Fp_F^u,p_F^s`$, we find

$$\delta \stackrel{~}{p}_F=\frac{\delta p_F}{12\alpha ^2\frac{\mathrm{\Delta }^2}{\mathrm{\Delta }_0^2+\mathrm{\Delta }^2}},$$
(20)
where $`\alpha (V^uV^s)/(V^u+V^s)`$. The gap equation (16) (with $`\delta p_F\delta \stackrel{~}{p}_F`$) and Eq. (20) are solved by

$$\mathrm{\Delta }=\{\begin{array}{cc}\mathrm{\Delta }_0,\hfill & \text{(BCS)}\hfill \\ \left[\mathrm{\Delta }_0\sqrt{V^uV^s}\left(\frac{\delta p_{\mathrm{th}}\delta p_F}{2\alpha ^21}\right)\right]^{\frac{1}{2}}.\hfill & \text{(BP)}\hfill \end{array}$$
(21)
The breached pairing phase is realized for $`2\alpha ^21>0`$, i.e., $`V^s/V^u<32\sqrt{2}`$. (See Fig. 1b.) The dashed line shows a continuous phase transition between BCS ($`\delta p_F\delta p_F^c`$) and breached pairing ($`\delta p_F>\delta p_F^c`$), where $`\delta p_F^c=\frac{2\mathrm{\Delta }_0}{\sqrt{V^uV^s}}(1\alpha ^2)`$. There is no overlap region where both states are solutions. The breached pairing state is stable and has lower energy than that of the unpaired matter (see (24) below).

For higher values of $`V^s/V^u`$ (but smaller than one), the BP solution becomes unstable, being of the same class as (18). Very close to $`2\alpha ^21=0`$ one finds from Eq. (21) that $`\delta p_F^c`$ approached $`\delta p_{\mathrm{th}}`$. The value $`32\sqrt{2}`$ is model dependent, and will be different when a more realistic model and a more accurate approximation scheme are used. In the limit $`V^s/V^u0`$ (so $`\alpha 1`$), the gap equation (16) exhibits the breached pairing phase for all $`\delta p_F>0`$; BCS occurs only at the point $`\delta p_F=0`$.

In mean field theory there is a single second-order phase transition between the BCS and BP phases (Fig.1b). At fixed mismatch in Fermi momenta $`\delta p_F`$ and as $`\mathrm{\Delta }_0`$ increases, one first encounters BP and then the BCS phase. The transition between the BP and BCS phases occurs at $`\delta p_F=\delta p_F^c`$. Since (21) is symmetric with respect to interchange of $`V^u`$ and $`V^s`$, the expression for critical points holds equally for both interior and “exterior” gap limits.

## III: Relative quark density fixed.

In this case, we also hold the overall chemical potential fixed while allowing the overall density to vary. In the limit $`p_0\lambda \mathrm{\Delta },\delta p_F`$, we find

$$\delta \stackrel{~}{p}_F=\sqrt{\delta p_F^2+\frac{4\mathrm{\Delta }^2}{V^uV^s}}.$$
(22)
The gap equation (16) (with $`\delta p_F\delta \stackrel{~}{p}_F`$) has no BCS but BP solution for any $`\delta p_F>0`$,

$$\mathrm{\Delta }=\left[\mathrm{\Delta }_0\sqrt{V^uV^s}(\delta p_{\mathrm{th}}\delta p_F)\right]^{\frac{1}{2}}.\text{(BP)}$$
(23)
The BCS solution only occurs at a single point, $`\delta p_F=0`$.

## Condensation energy.

To leading order in $`\mathrm{\Delta }/\delta p_F`$, the energy difference between the breached pairing and normal ($`\mathrm{\Delta }=0`$) phases is

$$E_{BP}E_N=\{\begin{array}{cc}+\frac{N_+(0)\mathrm{\Delta }^4}{2\delta p_F^2}\frac{1}{V^uV^s},\hfill & \text{(I)}\hfill \\ \frac{N_+(0)\mathrm{\Delta }^4}{2\delta p_F^2}\left(\frac{2\alpha ^21}{V^uV^s}\right),\hfill & \text{(II)}\hfill \\ \frac{N_+(0)\mathrm{\Delta }^4}{2\delta p_F^2}\frac{1}{V^uV^s}.\hfill & \text{(III)}\hfill \end{array}$$
(24)
The BP state is energetically favorable in case II if $`V^s/V^u<32\sqrt{2}`$ (with overall density and $`\delta \mu `$ fixed), and it is always favorable in case III (with the relative density and overall $`\mu `$ fixed).

In the breached pairing scenario pairing occurs at zero total mometum and peaks near the light species Fermi surface. It resembles standard BCS in that the pairing occurs at all angles, there is no necessity for translation or rotation symmetry breaking, and there is a branch of gapped excitations. There are, nevertheless, also gapless modes at two new “effective” Fermi surfaces, between which pairing is suppressed. For small $`\mathrm{\Delta }`$ the energy gain of the breached pairing phase with respect to unpaired normal matter is $`\mathrm{\Delta }^4`$, parametrically less than $`\mathrm{\Delta }^2`$ in the BCS case.

## Discussion.

In the context of QCD, enforcing charge neutrality will make case III relevant for the $`u`$ and $`d`$ quarks. For pairing between strange and light quarks, case II is relevant. Thus, breached pairing allows all flavors to form condensates $`ud_{BP},us_{BP},ds_{BP}`$ of breached pairing type. We suspect that there is a density range where it is energetically favorable compared to unpaired matter or the $`2`$SC phase with only $`ud_{BCS}0`$ condensation. Alford and Rajagopal Alford and Rajagopal (2002) showed the absence of $`2`$SC phase in compact stars, and suggested that there should be some non-BCS pattern of pairing including possible unpaired matter in the range of densities between hadronic and CFL matter. The breached pairing version of CFL appears to be a serious candidate for these intermediate densities. Since there is a simple second order phase transition from BP to BCS superfluidity, as a function of increasing density, it is natural to expect a single phase transition from breached pairing CFL to standard CFL quark matter. Detailed calculations, including shifts of Fermi surfaces due to electric and color neutrality, will be required to determine the range, if any, over which breached pairing CFL is favored (in preparation). At strong coupling one could even speculate about breached pairing CFL states at zero strangeness, where the strange chemical potential falls between $`\pm m_s`$.

The possibility that one might rearrange the free-particle Fermi surfaces to take advantage of opportunities for pairing was considered long ago in condensed matter literature by Sarma and others Sar . Superfluidity with quasiparticle dispersion relation crossing zero, thus leading to gapless modes, was encountered in a model resembling ours in Ref. Alford et al. (2000). In neither case, however, was the possibility of a stable breached pairing state recognized. (See also a recent study in cold atoms Wu- .) As this work was being completed reference Ref. Sho  appeared, in which a related problem, imposing a charge neutrality condition directly, was analyzed.

We have phrased our discussion in terms of QCD and quarks, but one can easily imagine closely related universality classes for condensed matter or cold atom systems.

The authors thank Jeff Bowers for useful discussions and valuable suggestions, Wolfram Schroers for help in numerical analyses; and M. Forbes, J. Kundu, Y. Schroeder and I. Sigalov. This work is supported in part by funds provided by the U.S.Department of Energy (D.O.E.) under cooperative research agreement DF-FC02-94ER40818.


text:
# All-sky upper limit for gravitational radiation from spinning neutron stars

## 1 Introduction

A unique property of the gravitational-wave detectors is that with a single observation resulting in one time series one can search for signals coming from all sky directions. In the case of other instruments like optical and radio telescopes, to cover the whole or even part of the sky requires a large amount of expensive telescope time as each sky position needs to be observed independently. The difficulty in the search for gravitational-wave signals is that they are very weak and they are deeply buried in the noise of the detector. Consequently the detection of these signals and interpretation of data analysis results is a delicate task. In this paper we present results of an all-sky search for continuous sources of gravitational radiation. A prime example of such a source is a spinning neutron star. A signal form such a source has definite characteristics that make it suitable for application of the optimal detection techniques based on matched filtering. Moreover such signals are stable as a result of the stability of the rotation of the neutron star and they will be present in the data for time periods much longer than the observational interval. This enables a reliable verification of the potential candidates by repeating the observations both by the same detector and by different detectors. To perform our all-sky search we have used the data of the EXPLORER  resonant bar detector. The directed search of the galactic center with the EXPLORER detector has already been carried out and an upper limit for the amplitude of the gravitational waves has been established .

Our paper is divided into two parts. In the first part we summarize the theoretical tools that we use in our analysis and in the second part we present results of our all-sky search. The main result of our analysis is an upper limit of $`2\times 10^{23}`$ for dimensionless amplitude of gravitational waves originating from continuous sources located in any position in the sky and in the frequency band from $`921.00`$ Hz to $`921.76`$ Hz.

The data analysis was performed by a team consisting of Pia Astone, Kaz Borkowski, Piotr Jaranowski, Andrzej Królak and Maciej Pietka and was carried out on the basis of Memorandum of Understanding between the ROG collaboration and Institute of Mathematics of Polish Academy of Sciences. More details about the search can be found on the website: http://www.astro.uni.torun.pl/`~`kb/AllSky/AllSky.html.

## 2 Data analysis methods

In this section we give a summary of data analysis techniques that we used in the search. The full exposition of our method is given in Ref. .

### 2.1 Response of a bar detector to a continuous gravitational-wave signal

Dimensionless noise-free response function $`h`$ of a resonant bar gravitational-wave detector to a weak plane gravitational wave in the long wavelength approximation \[i.e., when the size of the detector is much smaller than the reduced wavelength $`\lambda /(2\pi )`$ of the wave\] can be written as a linear combination of the wave polarization functions $`h_+`$ and $`h_\times `$:

$$h(t)=F_+(t)h_+(t)+F_\times (t)h_\times (t),$$
(1)
where $`F_+`$ and $`F_\times `$ are called the beam-pattern functions and can be written as

$`F_+(t)`$ $`=`$ $`a(t)\mathrm{cos}2\psi +b(t)\mathrm{sin}2\psi ,`$
$`F_\times (t)`$ $`=`$ $`b(t)\mathrm{cos}2\psi a(t)\mathrm{sin}2\psi ,`$ (2)
where $`\psi `$ is the polarization angle of the wave and

$`a(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\mathrm{cos}^2\gamma \mathrm{sin}^2\gamma \mathrm{sin}^2\varphi )(1+\mathrm{sin}^2\delta )\mathrm{cos}[2(\alpha \varphi _r\mathrm{\Omega }_rt)]`$
$`+{\displaystyle \frac{1}{2}}\mathrm{sin}2\gamma \mathrm{sin}\varphi (1+\mathrm{sin}^2\delta )\mathrm{sin}[2(\alpha \varphi _r\mathrm{\Omega }_rt)]`$
$`{\displaystyle \frac{1}{2}}\mathrm{sin}^2\gamma \mathrm{sin}2\varphi \mathrm{sin}2\delta \mathrm{cos}(\alpha \varphi _r\mathrm{\Omega }_rt)`$
$`+{\displaystyle \frac{1}{2}}\mathrm{sin}2\gamma \mathrm{cos}\varphi \mathrm{sin}2\delta \mathrm{sin}(\alpha \varphi _r\mathrm{\Omega }_rt)`$
$`+{\displaystyle \frac{1}{2}}(13\mathrm{sin}^2\gamma \mathrm{cos}^2\varphi )\mathrm{cos}^2\delta ,`$
$`b(t)`$ $`=`$ $`\mathrm{sin}2\gamma \mathrm{sin}\varphi \mathrm{sin}\delta \mathrm{cos}[2(\alpha \varphi _r\mathrm{\Omega }_rt)]`$ (3)
$`+(\mathrm{cos}^2\gamma \mathrm{sin}^2\gamma \mathrm{sin}^2\varphi )\mathrm{sin}\delta \mathrm{sin}[2(\alpha \varphi _r\mathrm{\Omega }_rt)]`$
$`\mathrm{sin}2\gamma \mathrm{cos}\varphi \mathrm{cos}\delta \mathrm{cos}(\alpha \varphi _r\mathrm{\Omega }_rt)`$
$`\mathrm{sin}^2\gamma \mathrm{sin}2\varphi \mathrm{cos}\delta \mathrm{sin}(\alpha \varphi _r\mathrm{\Omega }_rt).`$
In Eqs. (2.1) and (2.1) the angles $`\alpha `$ and $`\delta `$ are respectively right ascension and declination of the gravitational-wave source. The geodetic latitude of the detector’s site is denoted by $`\varphi `$, the angle $`\gamma `$ determines the orientation of the bar detector with respect to local geographical directions: $`\gamma `$ is measured counter-clockwise from East to the bar’s axis of symmetry. The rotational angular velocity of the Earth is denoted by $`\mathrm{\Omega }_r`$, and $`\varphi _r`$ is a deterministic phase which defines the position of the Earth in its diurnal motion at $`t=0`$ (the sum $`\varphi _r+\mathrm{\Omega }_rt`$ essentially coincides with the local sidereal time of the detector’s site, i.e., with the angle between the local meridian and the vernal point).

We are interested in a continuous gravitational wave, which is described by the wave polarization functions of the form

$$h_+(t)=h_{0+}\mathrm{cos}\mathrm{\Psi }(t),h_\times (t)=h_{0\times }\mathrm{sin}\mathrm{\Psi }(t),$$
(4)
where $`h_{0+}`$ and $`h_{0\times }`$ are independent constant amplitudes of the two wave polarizations. These amplitudes depend on the physical mechanisms generating the gravitational wave. In the case of a wave originating from a spinning neutron star these amplitudes can be estimated by

$$h_o=4.23\times 10^{23}\left(\frac{ϵ}{10^5}\right)\left(\frac{I}{10^{45}\text{g cm}^2}\right)\left(\frac{\text{1 kpc}}{r}\right)\left(\frac{f}{1\mathrm{kHz}}\right)^2,$$
(5)
where $`I`$ is the neutron star moment of inertia with respect to the rotation axis, $`ϵ`$ is the poloidal ellipticity of the star, and $`r`$ is the distance to the star. The value of $`10^5`$ of the parameter $`ϵ`$ in the above estimate should be treated as an upper bound. In reality it may be several orders of magnitude less.

The phase $`\mathrm{\Psi }`$ of the wave is given by

$`\mathrm{\Psi }(t)`$ $`=`$ $`\mathrm{\Phi }_0+\mathrm{\Phi }(t),`$
$`\mathrm{\Phi }(t)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{s_1}{}}}\omega _k{\displaystyle \frac{t^{k+1}}{(k+1)!}}+{\displaystyle \frac{𝐧_0𝐫_{\mathrm{SSB}}(t)}{c}}{\displaystyle \underset{k=0}{\overset{s_2}{}}}\omega _k{\displaystyle \frac{t^k}{k!}}.`$ (6)
In Eq. (2.1) the parameter $`\mathrm{\Phi }_0`$ is the initial phase of the wave form, $`𝐫_{\mathrm{SSB}}`$ is the vector joining the solar system barycenter (SSB) with the detector, $`𝐧_0`$ is the constant unit vector in the direction from the SSB to the source of the gravitational-wave signal. We assume that the gravitational-wave signal is almost monochromatic around some angular frequency $`\omega _0`$ which we define as instantaneous angular frequency evaluated at the SSB at $`t=0`$, $`\omega _k`$ ($`k=1,2,\mathrm{}`$) is the $`k`$th time derivative of the instantaneous angular frequency at the SSB evaluated at $`t=0`$. To obtain formulas (2.1) we model the frequency of the signal in the rest frame of the neutron star by a Taylor series. For the detailed derivation of the phase model (2.1) see Sec. II B and Appendix A of Ref. . The integers $`s_1`$ and $`s_2`$ are the number of spin downs to be included in the two components of the phase. We need to include enough spin downs so that we have a sufficiently accurate model of the signal to extract it from the noise.

The vector $`𝐫_{\mathrm{SSB}}`$ is the sum of two vectors: the vector $`𝐫_{\mathrm{ES}}`$ joining the solar system barycenter and the Earth barycenter and vector $`𝐫_\mathrm{E}`$ joining the Earth barycenter and the detector. The vector $`𝐫_{\mathrm{ES}}`$ can be obtained form the JPL Planetary and Lunar Ephemerides: “DE405/LE405”. They are available via the Internet at anonymous ftp: `navigator.jpl.nasa.gov`, the directory: `ephem/export`. Whereas the vector $`𝐫_\mathrm{E}`$ can be accurately computed using the International Earth Rotation Service (IERS) tables available at `http://hpiers.obspm.fr/iers/eop/eopc04`. The codes to read the above files and to calculate the vector $`𝐫_{\mathrm{SSB}}`$ are described in Ref. .

It is convenient to write the response of the gravitational-wave detector given above in the following form

$$h(t)=A_1a(t)\mathrm{cos}\mathrm{\Phi }(t)+A_2b(t)\mathrm{cos}\mathrm{\Phi }(t)+A_3a(t)\mathrm{sin}\mathrm{\Phi }(t)+A_4b(t)\mathrm{sin}\mathrm{\Phi }(t),$$
(7)
where $`A_i`$ are four constant amplitudes, the functions $`a`$ and $`b`$ are given by Eqs. (2.1), and $`\mathrm{\Phi }`$ is the phase given by Eq. (2.1). The modulation amplitudes $`a`$ and $`b`$ depend on the right ascension $`\alpha `$ and the declination $`\delta `$ of the source (they also depend on the angles $`\varphi `$ and $`\gamma `$). The phase $`\mathrm{\Phi }`$ depends on the frequency $`\omega _0`$, $`s`$ spin-down parameters $`\omega _k`$ $`(k=1,\mathrm{},s)`$, and on the angles $`\alpha `$, $`\delta `$. We call the parameters $`\omega _0`$, $`\omega _k`$, $`\alpha `$, $`\delta `$ the intrinsic parameters and the remaining ones the extrinsic parameters. Moreover the phase $`\mathrm{\Phi }`$ depends on the latitude $`\varphi `$ of the detector. The whole signal $`h`$ depends on $`s+5`$ unknown parameters: $`h_{0+}`$, $`h_{0\times }`$, $`\alpha `$, $`\delta `$, $`\omega _0`$, $`\omega _k`$.

### 2.2 Optimal data analysis method

We assume that the noise $`n`$ in the detector is an additive, stationary, Gaussian, and zero-mean random process. Then the data $`x`$ (when the signal $`h`$ is present) can be written as

$$x(t)=n(t)+h(t).$$
(8)
The log likelihood function has the form

$$\mathrm{log}\mathrm{\Lambda }=(x|h)\frac{1}{2}(h|h),$$
(9)
where the scalar product $`(|)`$ is defined by

$$(x|y):=4\mathrm{}_0^{\mathrm{}}\frac{\stackrel{~}{x}(f)\stackrel{~}{y}^{}(f)}{S_h(f)}𝑑f.$$
(10)
In Eq. (10) $`\mathrm{}`$ denotes the real part, tilde denotes Fourier transform, asterisk is complex conjugation, and $`S_h`$ is the one-sided spectral density of the detector’s noise. Assuming that over the bandwidth of the signal $`h`$ the spectral density $`S_h(f)`$ is nearly constant and equal to $`S_h(f_0)`$, where $`f_0`$ is the frequency of the signal measured at the SSB at $`t=0`$, the log likelihood ratio from Eq. (9) can be approximated by

$$\mathrm{ln}\mathrm{\Lambda }\frac{2T_o}{S_h(f_0)}\left(xh\frac{1}{2}h^2\right).$$
(11)
where $``$ denotes time averaging over the observational interval $`[0,T_o]`$:

$$x:=\frac{1}{T_o}_0^{T_o}x(t)𝑑t.$$
(12)
The signal $`h`$ depends linearly on four amplitudes $`A_i`$. The equations for the maximum likelihood (ML) estimators $`\widehat{A}_i`$ of the amplitudes $`A_i`$ are given by

$$\frac{\mathrm{ln}\mathrm{\Lambda }}{A_i}=0,i=1,\mathrm{},4.$$
(13)
One can easily find the explicit analytic solution to Eqs. (13). To simplify formulas we assume that the observation time $`T_o`$ is an integer multiple of one sidereal day, i.e., $`T_o=n(2\pi /\mathrm{\Omega }_r)`$ for some positive integer $`n`$. Then the time average of the product of the functions $`a`$ and $`b`$ vanishes, $`ab=0`$, and the analytic formulas for the ML estimators of the amplitudes are given by

$$\widehat{A}_12\frac{xh_1}{a^2},\widehat{A}_22\frac{xh_2}{b^2},\widehat{A}_32\frac{xh_3}{a^2},\widehat{A}_42\frac{xh_4}{b^2}.$$
(14)
The reduced log likelihood function $``$ is the log likelihood function where amplitude parameters $`A_i`$ were replaced by their estimators $`\widehat{A}_i`$. By virtue of Eqs. (14) from Eq. (11) one gets

$$\frac{2}{S_h(f_0)T_o}\left(\frac{|F_a|^2}{a^2}+\frac{|F_b|^2}{b^2}\right),$$
(15)
where

$`F_a`$ $`:=`$ $`{\displaystyle _0^{T_o}}x(t)a(t)\mathrm{exp}[i\mathrm{\Phi }(t)]𝑑t,`$
$`F_b`$ $`:=`$ $`{\displaystyle _0^{T_o}}x(t)b(t)\mathrm{exp}[i\mathrm{\Phi }(t)]𝑑t.`$ (16)
The ML estimators of the signal parameters are obtained in two steps. Firstly, the estimators of the frequency, the spin-down parameters, and the angles $`\alpha `$ and $`\delta `$ are obtained by maximizing the functional $``$ with respect to these parameters. Secondly, the estimators of the amplitudes $`A_i`$ are calculated from the analytic formulas (14) with the correlations $`xh_i`$ evaluated for the values of the parameters obtained in the first step.

In order to calculate the optimum statistics $``$ efficiently we introduce an approximation to the phase of the signal that is valid for observation times short comparing to the period of 1 year. The phase of the gravitational-wave signal contains terms arising from the motion of the detector with respect to the SSB. These terms consist of two contributions, one which comes from the motion of the Earth barycenter with respect to the SSB, and the other which is due to the diurnal motion of the detector with respect to the Earth barycenter. The first contribution has a period of one year and thus for shorter observation times can be well approximated by a few terms of the Taylor expansion. The second term has a period of 1 sidereal day and to a very good accuracy can be approximated by a circular motion. We find that for two days of observation time that we used in the analysis of the EXPLORER data presented in Section 3 we can truncate the expansion at terms that are quadratic in time. Moreover in the Taylor expansion of the frequency parameter it is sufficient to include terms up to the first spin down. We thus propose the following approximate simple model of the phase of the gravitational-wave signal:

$$\mathrm{\Psi }_s(t)=p+p_0t+p_1t^2+A\mathrm{cos}(\mathrm{\Omega }_rt)+B\mathrm{sin}(\mathrm{\Omega }_rt),$$
(17)
where $`\mathrm{\Omega }_r`$ is the rotational angular velocity of the Earth. The characteristic feature of the above approximation is that the phase is a linear function of the parameters of the signal. The parameters $`A`$ and $`B`$ can be related to the right ascension $`\alpha `$ and the declination $`\delta `$ of the gravitational-wave source through the equations

$`A`$ $`=`$ $`{\displaystyle \frac{\omega _0r}{c}}\mathrm{cos}\delta \mathrm{cos}(\alpha \varphi _r),`$ (18)
$`B`$ $`=`$ $`{\displaystyle \frac{\omega _0r}{c}}\mathrm{cos}\delta \mathrm{sin}(\alpha \varphi _r),`$ (19)
where $`\omega _0`$ is the angular frequency of the gravitational-wave signal and $`r`$ is the equatorial component of the detector’s radius vector. The parameters $`p`$, $`p_0`$, and $`p_1`$ contain contributions both form the intrinsic evolution of the gravitational-wave source and the modulation of the signal due to the motion of the Earth around the Sun. With this approximation the integrals given by Eqs. (2.2) and needed to compute $``$ become Fourier transforms and they can be efficiently calculated using the FFT algorithm. Thus the evaluation of $``$ consists of correlation of the data with two linear filters depending on parameters $`p_1,A,B`$ followed by FFTs. In Ref.  we have verified that for the case of our search the linear approximation to the phase does not lead to the loss of signal-to-noise ratio of more than 5%.

In order to identify potential gravitational-wave signals we apply a two step procedure consisting of a coarse search followed by a fine search. The coarse search consists of evaluation of $``$ on a discrete grid constructed in such a way that the loss of the signal-to-noise is minimized and comparison of the obtained values of $``$ with a predefined threshold $`_o`$. The parameters of the nodes of the grid for which the threshold is crossed are registered as potential gravitational-wave signals. The threshold is calculated from a chosen false alarm probability which is defined as the probability that $``$ crosses the threshold when no signal is present and the data is only noise. The fine search consists of finding a local maximum of $``$ using a numerical implementation of the Nelder-Mead algorithm, where coordinates of the starting point of the maximization procedure are the parameter values of the coarse search.

To calculate the false alarm probability as a function of a threshold and to construct a grid in the parameter space we introduce yet another approximation of our signal. Namely we use a signal with a constant amplitude and the phase given by Eq. (17). In paper  we have shown that the Fisher matrix for the exact model with amplitude modulations given by Eq. (7) can be accurately reproduced by the Fisher matrix of a constant amplitude model. As the calculations of the false alarm probability and construction of a grid in the parameter space are based on the Fisher matrix we expect that the constant amplitude model is a good approximation for the purpose of the above calculations. We stress that in the search of real data we used the full model with amplitude modulations.

We shall now summarize the statistical properties of the functional $``$. We first assume that the phase parameters are known and consequently that $``$ is only a function of the random data $`x`$. When the signal is absent $`2`$ has a $`\chi ^2`$ distribution with four degrees of freedom and when the signal is present it has a noncentral $`\chi ^2`$ distribution with four degrees of freedom and noncentrality parameter equal to the *optimal signal-to-noise ratio* $`d`$ defined as

$$d:=\sqrt{(h|h)}.$$
(20)
Consequently the probability density functions $`p_0`$ and $`p_1`$ when respectively the signal is absent and present are given by

$`p_0()`$ $`=`$ $`\mathrm{exp}(),`$ (21)
$`p_1(d,)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{d}}I_1\left(d\sqrt{2}\right)\mathrm{exp}\left({\displaystyle \frac{1}{2}}d^2\right),`$ (22)
where $`I_1`$ is the modified Bessel function of the first kind and order 1. The false alarm probability $`P_F`$ is the probability that $``$ exceeds a certain threshold $`_o`$ when there is no signal. In our case we have

$$P_F(_o):=__o^{\mathrm{}}p_0()𝑑=(1+_o)\mathrm{exp}(_o).$$
(23)
The probability of detection $`P_D`$ is the probability that $``$ exceeds the threshold $`_o`$ when the signal-to-noise ratio is equal to $`d`$:

$$P_D(d,_o):=__o^{\mathrm{}}p_1(d,)𝑑.$$
(24)
The integral in the above formula cannot be evaluated in terms of known special functions.

We see that when the noise in the detector is Gaussian and the phase parameters are known the probability of detection of the signal depends on a single quantity: the optimal signal-to-noise ratio $`d`$. When the phase parameters are unknown the optimal statistics $``$ depends not only on the random data $`x`$ but also on the phase parameters that we shall denote by $`\xi `$. Such an object is called in the theory of stochastic processes a random field. Let us consider the correlation function of the random field

$$C(\xi ,\xi ^{}):=E\left\{[(\xi )m(\xi )][(\xi ^{})m(\xi ^{})]\right\},$$
(25)
where

$$m(\xi ):=E\left\{(\xi )\right\}.$$
(26)
For the case of the constant amplitude, linear phase model we have

$$m(\xi )=1,$$
(27)
and

$`C(\xi ,\xi ^{})`$ $`=C(\xi \xi ^{})`$ (28)
$`=\mathrm{cos}[\mathrm{\Delta }p_0t+\mathrm{\Delta }p_1t^2+\mathrm{\Delta }A\mathrm{cos}(\mathrm{\Omega }_rt)+\mathrm{\Delta }B\mathrm{sin}(\mathrm{\Omega }_rt)]^2`$
$`+\mathrm{sin}[\mathrm{\Delta }p_0t+\mathrm{\Delta }p_1t^2+\mathrm{\Delta }A\mathrm{cos}(\mathrm{\Omega }_rt)+\mathrm{\Delta }B\mathrm{sin}(\mathrm{\Omega }_rt)]^2,`$
where $`\mathrm{\Delta }`$ denotes the difference between the parameter values. Thus the correlation function $`C`$ depends only on the difference $`\xi \xi ^{}`$ and not on the values of the parameters themselves. In this case the random field $``$ is called homogeneous. For such fields we have developed  a method to calculate the false alarm probability that we only summarize here. The main idea is to divide the space of the intrinsic parameters $`\xi `$ into elementary cells. The size of the cell is determined by the *characteristic correlation hypersurface* of the random field $``$. The correlation hypersurface is defined by the condition that at this hypersurface the correlation $`C`$ equals to half its maximum value. Assuming that $`C`$ attains its maximum value when $`\xi \xi ^{}=`$ 0 the equation of the characteristic correlation hypersurface reads

$$C(\tau )=\frac{1}{2}C(0),$$
(29)
where we have introduced $`\tau :=\xi \xi ^{}`$. Let us expand the left-hand side of Eq. (29) around $`\tau =`$ 0 up to terms of the second order in $`\tau `$. We arrive at the equation

$$\underset{i,j=1}{\overset{M}{}}G_{ij}\tau _i\tau _j=1,$$
(30)
where $`M`$ is the dimension of the intrinsic parameter space and the matrix $`G`$ is defined as follows

$$G_{ij}:=\frac{1}{C(0)}\frac{^2C(\tau )}{\tau _i\tau _j}|_{\tau =\mathrm{𝟎}}.$$
(31)
Equation (30) defines the boundary of an $`M`$-dimensional hyperellipsoid which we call the *correlation hyperellipsoid*. The $`M`$-dimensional hypervolume $`V_{cell}`$ of the hyperellipsoid defined by Eq. (30) equals

$$V_{cell}=\frac{\pi ^{M/2}}{\mathrm{\Gamma }(M/2+1)\sqrt{detG}},$$
(32)
where $`\mathrm{\Gamma }`$ denotes the Gamma function. We estimate the number $`N_c`$ of elementary cells by dividing the total Euclidean volume $`V_{total}`$ of the parameter space by the volume $`V_{cell}`$ of the correlation hyperellipsoid, i.e. we have

$$N_c=\frac{V_{total}}{V_{cell}}.$$
(33)
We approximate the probability distribution of $`(\xi )`$ in each cell by the probability distribution $`p_0()`$ when the parameters are known \[in our case it is the one given by Eq. (21)\]. The values of the statistics $``$ in each cell can be considered as independent random variables. The probability that $``$ does not exceed the threshold $`_o`$ in a given cell is $`1P_F(_o)`$, where $`P_F(_o)`$ is given by Eq. (23). Consequently the probability that $``$ does not exceed the threshold $`_o`$ in *all* the $`N_c`$ cells is $`[1P_F(_o)]^{N_c}`$. The probability $`P_F^T`$ that $``$ exceeds $`_o`$ in *one or more* cells is thus given by

$$P_F^T(_o)=1[1P_F(_o)]^{N_c}.$$
(34)
This is the false alarm probability when the phase parameters are unknown.

The basic quantity to consider in the construction of the grid of templates in the parameter space is the expectation value $`E_1\left\{\right\}`$ of the statistics when the signal is present. For the constant amplitude linear phase model we have that

$$E_1\left\{(\xi )\right\}=1+\frac{1}{2}d^2C(\xi ),$$
(35)
where $`C`$ is the correlation function of the random field introduced above and $`d`$ is the signal-to-noise ratio of the signal present in the data. We choose the grid of templates in such a way that the correlation between any potential signal present in the data and the nearest point of the grid never falls below a certain value. In the case of the approximate model of the signal that we use the grid is unform and consists of regular polygons in the space parametrized by $`p_0,p_1,A,B`$. The construction of the grid is described in detail in Section VIIA of . For the grid that we used in the search of the EXPLORER data the correlation function $`C`$ for any signal present in the data was greater than 0.77.

## 3 An all-sky search of the EXPLORER data

We have implemented the theoretical tools presented in Section 2 and we have performed an all-sky search for continuous sources of gravitational waves in the data of the resonant bar detector EXPLORER<sup>1</sup><sup>1</sup>1The EXPLORER detector is operated by the ROG collaboration located in Italian Istituto Nazionale di Fisica Nucleare (INFN); see http://www.roma1.infn.it/rog/explorer/explorer.html. . The EXPLORER detector has collected many years of data with a high duty cycle (e.g. in 1991 the duty cycle was 75%). The EXPLORER detector was most sensitive over certain two narrow bandwidths (called minus and plus modes) of about 1 Hz wide at two frequencies around 1 kHz. To make the search computationally manageable we analyzed two days of data in the narrow band were the detector had the best sensitivity. By narrowing the bandwidth of the search we can shorten the length of the data to be analyzed as we need to sample the data at only twice the bandwidth. For the sake of the FFT algorithm it is best to keep the length of the data to be a power of $`2`$. Consequently we have chosen the number of data points to analyze to be $`N=2^{18}`$. Thus for $`T_o=2`$ days of observation time the bandwidth $`\mathrm{\Delta }\nu `$ was $`\mathrm{\Delta }\nu =N/(2T_o)0.76`$ Hz. We have also chosen to analyze the data for the plus mode. As a result we searched the bandwidth from 921.00 Hz to 921.76 Hz. We have used the filters with the phase linear in the parameters given by Eq. (17). In the filters we have included the amplitude modulation. The number of cells $`N_c`$ calculated from Eq. (33) was around $`1.6\times 10^{12}`$. Consequently from Eq. (34) the threshold signal-to-noise ratio for 1% false alarm probability was equal to $`8.3`$. In the search that we have performed we have used a lower threshold signal-to-noise of $`6.7`$. The aim of lowering the threshold was to make up for the loss of the signal-to-noise ratio due the discreteness of the grid of templates and due to the use of filters that only approximately matched the true signal. The number of points in the grid over which we had to calculate the statistics $``$ turned out to be $`\mathrm{183\hspace{0.17em}064\hspace{0.17em}440}`$. This number involved $`\mathrm{63\hspace{0.17em}830}`$ positions in the sky and $`\mathrm{2\hspace{0.17em}868}`$ spin down values for each sky position. We have carried out the search on a network of PCs and workstations. We had around two dozens of processors at our disposal. The data analysis started in September 2001 and was completed in November 2002.

The two-day stretch of data that we analyzed was taken from a larger set of $`13`$ days of data taken by the EXPLORER detector in November 1991. We have chosen the 2-day stretch of data on the basis of conformity of the data to the Gaussian random process. We have divided the data into $`2^{16}`$ points sections corresponding to around $`11`$ hours of data. For each stretch we have performed the Kolmogorov-Smirnov (KS) test. The results of the KS test are presented in Figure 1. For the case of KS test the higher the probability the more Gaussian the data are. From the above test we conclude that large parts of data are approximately Gaussian. On the basis of the above analysis we have chosen the two day data stretch to begin at the first sample of the 7th data stretch corresponding to Modified Julian Date of $`48580.7909`$. In Figure 2 we have presented the spectral density of the two day stretch of data that we analyzed. We see that the minimum spectral density was close to $`10^{21}/\sqrt{\text{Hz}}`$.

We have obtained $`\mathrm{22\hspace{0.17em}295}`$ threshold crossings for the Northern hemisphere and $`\mathrm{44\hspace{0.17em}701}`$ for the Southern hemisphere. We consider all candidates contained within a single cell of the parameter space defined in Section 2.2 as dependent and we choose as an independent candidate the candidate that has the highest signal-to-noise ratio within one cell. Consequently we have obtained $`\mathrm{11\hspace{0.17em}703}`$ independent candidates for the Northern hemisphere and $`\mathrm{18\hspace{0.17em}702}`$ for the Southern hemisphere. In Figures 3 and 4 we have plotted the histograms of the values of the statistics $``$ for the independent candidates and we have compared it with the theoretical distribution for $``$ when no signal is present in the data. A good agreement with the theoretical distribution is another indication of Gaussianity of the data. It also reveals that there are no obvious populations of the continuous gravitational-wave sources at the level of sensitivity of our search. In the search no event has crossed our $`99\%`$ confidence threshold of $`8.3`$. The strongest signal obtained by the coarse search had the signal-to-noise ratio of $`7.9`$ and the fine search increased that value to $`8.2`$.

As we do not have a detection of a gravitational-wave signal we can make a statement about the upper bound for the gravitational-wave amplitude. To do this we take our strongest candidate of signal-to-noise ratio $`d_o`$ and we suppose that it resulted from a gravitational-wave signal. Then, using formula (24), we calculate the signal-to-noise $`d_{ul}`$ of the gravitational-wave signal so that there is $`1\%`$ probability that it crosses the threshold $`_o`$ corresponding to $`d_o`$, where $`_o=2+\frac{1}{2}d_o^2`$. The $`d_{ul}`$ is the desired $`99\%`$ confidence upper bound. For $`d_o=8.2`$, which corresponds to the signal-to-noise ratio of our strongest candidate, we find that $`d_{ul}=5.9`$. For the EXPLORER detector this corresponds to the dimensionless amplitude of the gravitational-wave signal equal to $`2\times 10^{23}`$. Thus we have the following conclusion:

In the frequency band from 921.00 Hz to 921.76 Hz and for signals coming from any sky direction the dimensionless amplitude of the gravitational-wave signal from a continuous source is less than $`\mathrm{𝟐}\times \mathrm{𝟏𝟎}^{\mathrm{𝟐𝟑}}`$ with $`99\%`$ confidence.

## Acknowledgments

The work of P. Jaranowski, A. Królak, K. M. Borkowski, and M. Pietka was supported in part by the KBN (Polish State Committee for Scientific Research) Grant No. 2P03B 094 17. The part of A. Królak research described in this paper was performed while he held an NRC-NASA Resident Research Associateship at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We would also like to thank Interdisciplinary Center for Mathematical and Computational Modeling of Warsaw University for computing time.

## References