text:
# High frequency sources of gravitational waves

## 1 Introduction

The new generation of gravitational wave (GW) detectors is already collecting data by improving the sensitivity by at least one order of magnitude compared to the operating resonant detectors. Broadband GW detectors are sensitive to frequencies between 50 and a few hundred Hz. The next generation will broaden the bandwidth but it will still not be sensitive enough to frequencies over 500-600Hz, unless they are operated as narrow-band detectors . There are also suggestions for wide band resonant detectors in the kHz band .

In this short review we will discuss some of the sources that are in the frequency band ($`500600`$Hz), where the interferometers are sensitive enough only if they are narrowbanded. Since there is a variety of GW sources with very interesting physics associated to them, this high-frequency window deserves special attention. If either resonant or narrow-band interferometers achieve the needed sensitivity, there is a plethora of unique information that can be collected.

## 2 Gravitational collapse

One of the most spectacular events in the Universe is the supernova (SN) collapse to create a neutron star (NS) or a black hole (BH). Core collapse is a very complicated event and a proper study of the event demands a deeper understanding of neutrino emission, amplification of the magnetic fields, angular momentum distribution, pulsar kicks, etc. There are many viable explanations for each of the above issues but it is still not possible to combine all of them together into a consistent theory. Gravitational waves emanating from the very first moments of the core collapse might shed light on all the above problems and help us understand the details of this dramatic event. Gravitational collapse compresses matter to nuclear densities, and is responsible for the core bounce and the shock generation. The event proceeds extremely fast, lasting only less than a second, and the dense fluid undergoes motions with relativistic speeds ($`v/c0.20.4`$). Even small deviations from spherical symmetry during this phase can generate copious amounts of GWs. However, the size of these asymmetries is not known. From observations in the electromagnetic spectrum we know that stars more massive than $`8M_{}`$ end in core collapse and that $`90\%`$ of them are stars with masses $`820M_{}`$. During the collapse most of the material is ejected and if the progenitor star has a mass $`M20M_{}`$ it leaves behind a neutron star. If $`M20M_{}`$ more than 10% falls back and pushes the proto-neutron-star (PNS) above the maximum NS mass leading to the formation of black holes (type II collapsars). Finally, if the progenitor star has a mass $`M40M_{}`$ no supernova is launched and the star collapses to form a BH (type I collapsars).

A significant amount of the ejected material can fall back, subsequently spinning up and reheating the nascent NS. Instabilities can be excited again during such a process. If a BH formed its quasi-normal modes (QNM) can be excited for as long as the process lasts. “Collapsars” accrete material during the very first few seconds, at rates $`12M_{}`$/sec. Later the accretion rate is reduced by an order of magnitude but still material is accreted for a few tenths of seconds. Typical frequencies of the emitted GWs are at the range 1-3kHz for $`310M_{}`$ BHs. If the disk around the central object has a mass $`1M_{}`$ self-gravity becomes important and gravitational instabilities (spiral arms, bars) might develop and radiate GWs. There is also the possibility that the collapsed material might fragment into clumps, which orbit for some cycles like a binary system (fragmentation instability).

The supernova event rate is 1-2 per century per galaxy  and about 5-40% of them produce BHs through the fall back material . Conservation of angular momentum suggests that the final objects should rotate close to the mass shedding limit, but this is still an open question, since there is limited knowledge of the initial rotation rate of the final compact object. Pulsar statistics suggest that the initial periods are probably considerably shorter than $`20`$ms. This strong increase of rotation during the collapse has been observed in many numerical simulations (see e.g. ).

Core collapse as a potential source of GWs has been studied for more than three decades (some of the most recent calculations can be found in ). All these numerical calculations show that signals from Galactic supernova ($`d10`$kpc) are detectable even with the initial LIGO/Virgo sensitivity at frequencies $``$1kHz. Advanced LIGO can detect signals from distances of 1Mpc but it will be difficult with the designed broadband sensitivity to resolve signals from the Virgo cluster ($``$15Mpc). The typical GW amplitude from the 2D numerical simulations  for an observer located on the equatorial plane of the source is

$$h9\times 10^{21}\epsilon \left(\frac{10\mathrm{k}\mathrm{p}\mathrm{c}}{d}\right)$$
(1)
where $`\epsilon 1`$ is the normalized GW amplitude. The total energy radiated in GWs during the collapse is $`10^610^8M_{}c^2`$. These numerical estimates are not yet conclusive, important aspects such as 3D hydrodynamics combined with proper spacetime evolution have been neglected. The influence of the magnetic fields have been ignored in most calculations. The proper treatment of these issues might not change the above estimations by orders of magnitude but it will provide a conclusive answer. There are also issues that need to be understood such as the pulsar kicks (velocities even higher than 1000 km/s) which suggest that in a fraction of newly-born NSs (and BHs) the process may be strongly asymmetric, or the polarization of the light spectra in supernovae which also indicate significant asymmetries . Better treatment of the microphysics and construction of accurate progenitor models for the angular momentum distributions are needed. All these issues are under investigation by many groups.

## 3 Rotational instabilities

Newly born neutron stars are expected to rotate rapidly, being subject to rotation induced instabilities. These arise from non-axisymmetric perturbations having angular dependence $`e^{im\varphi }`$. Early Newtonian estimates have shown that a dynamical bar-mode ($`m=2`$) instability is excited if the ratio $`\beta =T/W`$ of the rotational kinetic energy $`T`$ to the gravitational binding energy $`W`$ is larger than $`\beta _{\mathrm{dyn}}=0.27`$. The instability develops on a dynamical time scale (the time that a sound wave needs to travel across the star) which is about one rotation period and may last from 1 to 100 rotations depending on the degree of differential rotation in the PNS. Another class of instabilities are those driven by dissipative effects such as fluid viscosity or gravitational radiation. Their growth time is much longer (many rotational periods) but they can be excited for significantly lower rotational rates, $`\beta 0.14`$.

### 3.1 Bar-mode instability

The bar-mode instability can be excited in a hot PNS, a few milliseconds after the core-bounce, given a sufficiently large $`\beta `$. It might also be excited a few tenths of seconds later, when the NS cools enough due to neutrino emission and contracts still further ($`\beta 1/R`$). The amplitude of the emitted gravitational waves can be estimated as $`hMR^2\mathrm{\Omega }^2/d`$, where $`M`$ is the mass of the body, $`R`$ its size, $`\mathrm{\Omega }`$ the rotational rate and $`d`$ the distance from Earth. This leads to an estimation of the GW amplitude

$$h9\times 10^{23}\left(\frac{ϵ}{0.2}\right)\left(\frac{f}{3\mathrm{k}\mathrm{H}\mathrm{z}}\right)^2\left(\frac{15\mathrm{M}\mathrm{p}\mathrm{c}}{d}\right)M_{1.4}R_{10}^2.$$
(2)
where $`ϵ`$ measures the ellipticity of the bar. Note that the GW frequency $`f`$ is twice the rotational frequency $`\mathrm{\Omega }`$. Such a signal is detectable only from sources in our galaxy or the nearby ones (our Local Group). If the sensitivity of the detectors is improved in the kHz region, then signals from the Virgo cluster will be detectable. If the bar persists for many ($``$ 10-100) rotation periods, then even signals from distances considerably larger than the Virgo cluster will be detectable. The event rate is of the same order as the SN rate (a few events per century per galaxy): this means that given the appropriate sensitivity at frequencies between 1-3kHz we might be able to observe a few events per year. Bars can be also created during the merging of NS-NS, BH-NS, BH-WD and even type II collapsars (see discussion in ).

The above estimates in general, rely on Newtonian hydrodynamics calculations; GR enhances the onset of the instability, $`\beta _{\mathrm{dyn}}0.24`$  and $`\beta _{\mathrm{dyn}}`$ may be even lower for large values of the compactness (larger $`M/R`$). The bar-mode instability may be excited for significantly smaller $`\beta `$ if centrifugal forces produce a peak in the density off the source s rotational center. Rotating stars with a high degree of differential rotation are also dynamically unstable for significantly lower $`\beta _{\mathrm{dyn}}0.01`$ . According to this scenario the unstable neutron star settles down to a non-axisymmetric quasi-stationary state which is a strong emitter of quasi-periodic gravitational waves

$$h_{\mathrm{eff}}3\times 10^{22}\left(\frac{R_{\mathrm{eq}}}{30\mathrm{k}\mathrm{m}}\right)\left(\frac{f}{800\mathrm{H}\mathrm{z}}\right)^{1/2}\left(\frac{100\mathrm{M}\mathrm{p}\mathrm{c}}{d}\right)M_{1.4}^{1/2}.$$
(3)
The bar-mode instability of differentially rotating neutron stars is an excellent source of gravitational waves but it is based on the assumption that the dissipation of non-axisymmetric perturbations by viscosity and magnetic fields is negligible. Magnetic fields might actually enforce the uniform rotation of the star on a dynamical timescale and the persistent non-axisymmetric structure might not have time to develop at all.

Numerical simulations have shown that the $`m=1`$ one-armed spiral mode might become dynamically unstable for considerably lower rotational rates . The $`m=1`$ instability depends critically on the softness of the equation of state (EoS) and the degree of differential rotation.

### 3.2 CFS instability, f and r-modes

After the initial bounce, neutron stars may maintain a considerable amount of deformation. They settle down to an axisymmetic configuration mainly due to emission of GWs, viscosity and magnetic fields. During this phase QNMs are excited. Technically speaking, an oscillating non-rotating star has equal values $`\pm |\sigma |`$ (the frequency of a mode) for the forward and backward propagating modes (corresponding to $`m=\pm |m|`$). Rotation changes the mode frequency by an amount $`\delta \sigma m\mathrm{\Omega }`$ and both the prograde and retrograde modes will be dragged forward by the stellar rotation. If the star spins sufficiently fast, this mode will appear moving forwards in the inertial frame (an observer at infinity), but still backwards in the rotating frame (an observer rotating with the star). Thus, an inertial observer sees GWs with positive angular momentum emitted by the retrograde mode, but since the perturbed fluid rotates slower than it would in absence of the perturbation, the angular momentum of the mode itself is negative. The emission of GWs consequently makes the angular momentum of the mode increasingly negative thus leading to an instability. From the above, one can easily conclude that a mode will be unstable if it is retrograde in the rotating frame and prograde for a distant observer measuring a mode frequency $`\sigma m\mathrm{\Omega }`$ i.e. the criterion will be $`\sigma (\sigma m\mathrm{\Omega })<0`$.

This class of frame-dragging instabilities is usually referred to as Chandrasekhar-Friedman-Schutz (CFS) instabilities. For the high frequency ($`f`$ and $`p`$) modes this is possible only for large values of $`\mathrm{\Omega }`$ or for quite large $`m`$. In general, for every mode there will always be a specific value of $`m`$ for which the mode will become unstable, although only modes with $`|m|<5`$ have an astrophysically significant growth time. The CFS mechanism is not only active for fluid modes but also for the spacetime or the so called w-modes. It is easy to see that the CFS mechanism is not unique to gravitational radiation: any radiative mechanism will have the same effect.

In GR, the f-mode ($`l=m=2`$) becomes unstable for $`\beta 0.060.08`$ . If the star has significant differential rotation the instability is excited for somewhat higher values of $`\beta `$ (see e.g. ). The $`f`$-mode instability is an excellent source of GWs. After the brief dynamical phase, the PNS becomes unstable and the instability deforms the star into a non-axisymmetric configuration via the $`l=2`$ bar mode. Since the star loses angular momentum, it spins-down, and the GW frequency sweeps from 1kHz down to about 100Hz. Such a signal if properly modelled can be detected from a distance of 100Mpc (if the mode grows to a large nonlinear amplitude).

Rotation does not only shifts the spectra of these modes; it also gives rise to a new type of restoring force, the Coriolis force, with an associated new family of rotational  or inertial modes. Inertial modes are primarily velocity perturbations and of special interest is the quadrupole inertial mode ($`r`$-mode) with $`l=m=2`$. The frequency of the $`r`$-mode in the rotating frame of reference is $`\sigma =2\mathrm{\Omega }/3`$. Using the CFS criterion for stability we can easily show that the $`r`$-mode is unstable for any rotation rate of the star. For temperatures between $`10^710^9`$K and rotation rates larger than 5-10% of the Kepler limit, the growth time of the unstable mode is smaller than the damping times of the bulk and shear viscosity. The mode grows until it saturates due to non-linear effects. The amplitude of the emitted GWs depends on $`\alpha `$. Mode coupling might not allow the growth of the instability to high amplitudes $`\alpha 10^210^3`$ . The existence of a crust or of hyperons in the core  and strong magnetic fields, affect the efficiency of the instability (for an extended review see ). For newly-born neutron stars the amplitude of GWs might not be large enough and the signals will be detectable only from the local group of galaxies ( $`d<1`$Mpc)

$$h(t)10^{21}\alpha \left(\frac{\mathrm{\Omega }}{1\mathrm{k}\mathrm{H}\mathrm{z}}\right)\left(\frac{100\mathrm{k}\mathrm{p}\mathrm{c}}{d}\right)$$
(4)
If the compact object is a strange star, then the instability will not reach high amplitudes ($`\alpha 10^310^4`$) but it will persist for a few hundred years and in this case there might be up to ten unstable stars per galaxy at any time . Integrating data for a few weeks can lead to an effective amplitude $`h_{\mathrm{eff}}10^{21}`$ for galactic signals at frequencies $`7001000`$Hz. The frequency of the signal changes only slightly on a timescale of a few months, so the radiation is practically monochromatic.

Old accreting neutron stars, radiating GWs due to the $`r`$-mode instability, at frequencies 400-700Hz, are probably a better source . Still, the efficiency and the actual duration of the process depends on the saturation amplitude $`\alpha `$. If the accreting compact object is a strange star then it might be a persistent source which radiates GWs for as long as accretion lasts.

## 4 Oscillations of black holes and neutron stars

Black-hole ringing. If the collapse produces a black hole (collapsar type I or II) the black hole will ring until it settles down to the stationary Kerr state. Although the ringing phase does not last too long (a few tenths of msecs), the ringing due to the excitation by the fallback material might last for secs. The frequency and the damping time of the oscillations for the $`l=m=2`$ mode can be estimated via the relations

$`\sigma `$ $``$ $`3.2\mathrm{kHz}M_{10}^1\left[10.63(1a/M)^{3/10}\right]`$ (5)
$`Q`$ $`=`$ $`\pi \sigma \tau 2\left(1a\right)^{9/20}`$ (6)
These relations together with similar ones either for the 2nd QNM or the $`l=2`$, $`m=0`$ can uniquely determine the mass $`M`$ and angular momentum $`a`$ of the BH if the frequency and the damping time of the signal have been accurately extracted . The amplitude of the ring-down waves depends on the BH’s initial distortion. If the excitation of the BH is due to falling material then the energy is roughly $`\mathrm{\Delta }Eϵ\mu c^2(\mu /M)`$ where $`ϵ0.01`$ . This leads to an effective GW amplitude

$$h_{\mathrm{eff}}2\times 10^{21}\left(\frac{ϵ}{0.01}\right)\left(\frac{10\mathrm{M}\mathrm{p}\mathrm{c}}{d}\right)\left(\frac{\mu }{M_{}}\right)$$
(7)
Neutron star ringing. If the collapse leaves behind a compact star, various types of oscillation modes might be excited which can help us estimate parameters of the star such as radius, mass, rotation rate and EoS. This gravitational wave asteroseismology is a unique way to find the radius and the EoS of compact stars. One can derive approximate formulas in order to connect the observable frequencies and damping times of the various stellar modes to the stellar parameters. For example, for the fundamental oscillation ($`l=2`$) mode ($`f`$-mode) of non-rotating stars we get

$`\sigma (\mathrm{kHz})`$ $``$ $`0.8+1.6M_{1.4}^{1/2}R_{10}^{3/2}+\delta _1m\overline{\mathrm{\Omega }}`$ (8)
$`\tau ^1(\mathrm{secs}^1)`$ $``$ $`M_{1.4}^3R_{10}^4\left(22.914.7M_{1.4}R_{10}^1\right)+\delta _2m\overline{\mathrm{\Omega }}`$ (9)
where $`\overline{\mathrm{\Omega }}`$ is the normalized rotation frequency of the star, and $`\delta _1`$ and $`\delta _2`$ are constants estimated by sampling data from various EoS. The typical frequencies of the NS modes are larger than 1kHz. On the other hand, 2D simulations of rotating core-collapse have shown that if a rapidly rotating NS is created, then the dominant mode is the quasi-radial mode (“$`l=0`$”), radiating through its $`l=2`$ piece at frequencies $``$ 800Hz-1kHz . Since each type of mode is sensitive to the physical conditions where the amplitude of the mode is greatest, the more information we get from the various types of modes the better we understand the details of the star.

Concluding, we should mention that the tidal disruption of a NS by a BH  or the merging of two NSs  may give valuable information for the radius and the EoS if we can recover the signal at frequencies higher than 1 kHz.

This work has been supported by the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137).

## References


format: text

text:
# Unmasking the nuclear matter equation of state

## Abstract

Accurately calibrated (or “best fit”) relativistic mean-field models are used to compute the distribution of isoscalar monopole strength in 90Zr and 208Pb, and the isovector dipole strength in 208Pb using a continuum random-phase-approximation approach. It is shown that the distribution of isoscalar monopole strength in 208Pb—but not in 90Zr—is sensitive to the density dependence of the symmetry energy. This sensitivity hinders the extraction of the compression modulus of symmetric nuclear matter from the isoscalar giant monopole resonance (ISGMR) in 208Pb. Thus, one relies on 90Zr, a nucleus with both a small neutron-proton asymmetry and a well developed ISGMR peak, to constrain the compression modulus of symmetric nuclear matter to the range $`K=(248\pm 6)\mathrm{MeV}`$. In turn, the sensitivity of the ISGMR in 208Pb to the density dependence of the symmetry energy is used to constrain its neutron skin to the range $`R_nR_p0.22`$ fm. The impact of this result on the enhanced cooling of neutron stars is briefly addressed.

Constraining the equation of state (EoS) of neutron-rich matter remains a fundamental problem in nuclear physics and astrophysics. The stability of neutron-rich nuclei Todd and Piekarewicz (2003), the dynamics of heavy-ion collisions Danielewicz et al. (2002); Li (2002), the structure of neutron stars Lattimer and Prakash (2001), and the simulation of core-collapse supernova Sumiyoshi et al. (1995); Buras et al. (2003), all depend sensitively on the EoS. Unfortunately, our window to the EoS is limited by terrestrial experiments that have, until now, probed stable nucleonic matter at (or close to) nuclear-matter saturation density. Fortunately, dramatic improvements are unfolding on several fronts. First, the commissioning of new radioactive-beam facilities all over the world will probe the EoS at large neutron-proton asymmetries. By defining the limits of nuclear existence, these exotic nuclei will constrain the EoS of neutron-rich matter at (and below) normal nuclear densities. Second, space-based telescopes have started to place important constrains on the high-density component of the EoS Pons et al. (2002); Walter and Lattimer (2002). New telescopes operating at a variety of wavelengths are turning neutron stars from theoretical curiosities into powerful diagnostic tools.

The nuclear matter equation of state is conveniently parametrized in terms of the energy of symmetric nuclear matter ($`/𝒜`$) and the symmetry energy ($`𝒮/𝒜`$) in the following form:

$`/𝒜(k_\mathrm{F},b)M=/𝒜(k_\mathrm{F})+b^2𝒮/𝒜(k_\mathrm{F})+𝒪(b^4)=`$
$`\left(ϵ_0+{\displaystyle \frac{1}{2}}K\xi ^2+\mathrm{}\right)+b^2\left(J+L\xi +{\displaystyle \frac{1}{2}}K_{\mathrm{sym}}\xi ^2+\mathrm{}\right).`$ (1)
Here the deviation from the equilibrium Fermi momentum is denoted by $`\xi (k_\mathrm{F}k_\mathrm{F}^0)/k_\mathrm{F}^0`$, the neutron-proton asymmetry by $`b(NZ)/A`$, and the various coefficients ($`K`$, $`J`$, $`L`$, $`K_{\mathrm{sym}}`$) parametrize the density dependence of the EoS around saturation density.

Seven decades of nuclear physics have placed important constraints on the nuclear matter equation of state. Indeed, the energy systematics of medium to heavy nuclei, when combined with accurately calibrated models, place the saturation point of symmetric nuclear matter at a density of $`\rho _00.15\mathrm{fm}^3`$ ($`k_F^01.3\mathrm{fm}^1`$) and a binding-energy per nucleon of $`ϵ_016`$ MeV. It should be noted that one of the main virtues of the above Taylor-series expansion around saturation density \[Eq. (1)\] is that the linear term in $`\xi `$ for symmetric nuclear matter (i.e., the pressure) automatically vanishes. Yet no such special saturation point exists in the case of the symmetry energy. Indeed, the symmetry energy at saturation density is not well known. Rather, it is the symmetry energy at the lower density of $`\stackrel{~}{\rho }_00.10\mathrm{fm}^3`$ ($`\stackrel{~}{k}_F^01.15\mathrm{fm}^1`$) that seems to be accurately constrained (to within 1 MeV) by available ground-state observables Brown (2000); Furnstahl (2002). It should be emphasized that present-day experiments can fix only one isovector quantity. If one insists—and one should not—on constraining the parameters of the symmetry energy at saturation density, then one would find a strong correlation among its parameters ($`J`$, $`L`$, $`K_{\mathrm{sym}}`$, $`\mathrm{}`$Furnstahl (2002). For example, relativistic models consistently predict larger values for both the symmetry-energy coefficient $`J`$ and the slope $`L`$ at saturation density relative to nonrelativistic Skyrme models. This must be so if all models are to reproduce the value of the symmetry energy at the lower Fermi momentum of $`\stackrel{~}{k}_F^0`$. Thus, in the present contribution we adopt the following convention: the symmetry energy is expanded around $`\stackrel{~}{k}_F^01.15\mathrm{fm}^1`$ and the value of the symmetry energy at $`\stackrel{~}{k}_F^0`$ is fixed at $`\stackrel{~}{J}=25.67`$ MeV. That is,

$$𝒮/𝒜(k_\mathrm{F})=\stackrel{~}{J}+\stackrel{~}{L}\stackrel{~}{\xi }+\frac{1}{2}\stackrel{~}{K}_{\mathrm{sym}}\stackrel{~}{\xi }^2+\mathrm{}$$
(2)
where $`\stackrel{~}{\xi }=(k_\mathrm{F}\stackrel{~}{k}_\mathrm{F}^0)/\stackrel{~}{k}_\mathrm{F}^0`$. Note that henceforth “tilde” quantities refer to parameters of the symmetry energy at $`\stackrel{~}{k}_F^01.15\mathrm{fm}^1`$.

Having established that existing ground-state observables accurately determine the binding energy per nucleon $`ϵ_0`$ at $`k_F^0`$ and the symmetry-energy coefficient $`\stackrel{~}{J}`$ at $`\stackrel{~}{k}_F^0`$, how can one constrain any further the density dependence of the equation of state? In the case of symmetric nuclear matter, the dynamics of small density fluctuations around the saturation point is controlled by the compression modulus $`K`$. The isoscalar giant monopole resonance (ISGMR) in heavy nuclei has long been regarded as the optimal observable from which to determine the compression modulus Blaizot (1980). This is especially true now that the breathing mode has been measured on a variety of nuclei with unprecedented accuracy Youngblood et al. (1999). In contrast, the density dependence of the symmetry energy is poorly constrained. Indeed, one may fit a variety of ground-state observables (such as charge densities, binding energies, and single-particle spectra) using accurately-calibrated models that, nevertheless, predict a wide range of values for the neutron skin of 208Pb Horowitz and Piekarewicz (2001a). As the neutron skin of a heavy nucleus is strongly correlated to the slope of the symmetry energy Brown (2000); Furnstahl (2002), measuring the skin thickness of a single heavy nucleus will constrain the density dependence of the symmetry energy. The Parity Radius Experiment (PREX) at the Jefferson Laboratory aims to measure the neutron radius of 208Pb accurately (to within $`0.05`$ fm) and model independently via parity-violating electron scattering Horowitz et al. (2001); Michaels et al. (2002). This experiment should provide a unique observational constraint on the density dependence of the symmetry energy.

While the above arguments suggest a clear path toward constraining the density dependence of the EoS, theoretical uncertainties have clouded these issues. First and foremost is the apparent discrepancy between nonrelativistic and relativistic predictions for the value of the compressional modulus of symmetric nuclear matter required to reproduce the ISGMR in 208Pb. While nonrelativistic models predict $`K220235`$ MeV Colò et al. (1992); Blaizot et al. (1995); Hamamoto et al. (1997), relativistic models argue for a significantly larger value $`K250270`$ MeV Lalazissis et al. (1997); Vretenar et al. (2000, 2003). Further, relativistic models systematically predict larger values for the neutron skin of 208Pb relative to nonrelativistic Skyrme models. One of the goals of this contribution is to show that these two points are related. Indeed, the aim of this contribution is twofold. First, to vindicate—through the exclusive use of accurately-calibrated models—our previous assertion that the distribution of ISGMR in heavy nuclei, and therefore the inferred value of $`K`$, is sensitive to the density dependence of the symmetry energy Piekarewicz (2002). Second, to rely on existing data on the isoscalar giant-monopole resonance in 90Zr and 208Pb Youngblood et al. (1999), and on the isovector giant-dipole resonance in 208Pb Ritman and et al. (1993), to set limits—simultaneously—on the compression modulus of symmetric nuclear matter and on the neutron skin of 208Pb. Note that since first proposed Piekarewicz (2002), other groups have addressed the possible impact of the density dependence of the symmetry energy on the ISGMR in 208Pb Vretenar et al. (2003); Agrawal et al. (2003); Colò and Van Gai (2003).

The starting point for the calculations is an interacting Lagrangian density of the following form Horowitz and Piekarewicz (2001a, b):

$`_{\mathrm{int}}=\overline{\psi }\left[g_\mathrm{s}\varphi (g_\mathrm{v}V_\mu +{\displaystyle \frac{g_\rho }{2}}𝝉𝐛_\mu +{\displaystyle \frac{e}{2}}(1+\tau _3)A_\mu )\gamma ^\mu \right]\psi `$
$`{\displaystyle \frac{\kappa }{3!}}\mathrm{\Phi }^3{\displaystyle \frac{\lambda }{4!}}\mathrm{\Phi }^4+\mathrm{\Lambda }_\mathrm{v}\left(W_\mu W^\mu \right)\left(𝐁_\mu 𝐁^\mu \right),`$ (3)
where $`\mathrm{\Phi }=g_\mathrm{s}\varphi `$, $`W_\mu =g_\mathrm{v}V_\mu `$, and $`𝐁_\mu =g_\rho 𝐛_\mu `$. The Lagrangian density includes Yukawa couplings of the nucleon field to a scalar ($`\varphi `$) and to three vector ($`V^\mu `$, $`𝐛^\mu `$, and $`A^\mu `$) fields. In addition to the Yukawa couplings, the Lagrangian is supplemented by three nonlinear meson interactions. The inclusion of scalar-meson interactions (via $`\kappa `$ and $`\lambda `$) is used to soften the equation of state of symmetric nuclear matter, while the mixed isoscalar-isovector coupling ($`\mathrm{\Lambda }_\mathrm{v}`$) modifies the density dependence of the symmetry energy—without affecting well known ground-state properties. Note that this last term was absent from Ref. Piekarewicz (2002), so the softening of the symmetry energy had to be done artificially. This drawback has now been corrected.

The relativistic mean-field models employed here are motivated by the enormous success of the NL3 parametrization Lalazissis et al. (1997). Indeed, the NL3000 set used here (having $`\mathrm{\Lambda }_\mathrm{v}=0`$) is practically identical to the original NL3 model. The other sets are obtained by adding an isoscalar-isovector coupling $`\mathrm{\Lambda }_\mathrm{v}0`$, while at the same time re-adjusting the strength of the $`NN\rho `$ coupling constant ($`g_\rho `$) to maintain the symmetry-energy coefficient fixed at $`\stackrel{~}{J}=25.67`$ MeV (see discussion above). The aim of this added coupling is to change the neutron density of heavy nuclei, while leaving intact ground-state observables that are well constrained experimentally. One should stress that the addition of $`\mathrm{\Lambda }_\mathrm{v}`$ has no impact on the properties of symmetric nuclear matter, so the saturation properties remain unchanged. In summary, all the models used in this contribution share the following properties with the original NL3 model of Ref. Lalazissis et al. (1997): for symmetric nuclear matter, a Fermi momentum at saturation of $`k_F^0=1.30\mathrm{fm}^1`$ with a binding-energy per nucleon of $`ϵ_0=16.24`$ MeV, and a compression modulus of $`K=271`$ MeV. For the symmetry energy, a symmetry-energy coefficient of $`\stackrel{~}{J}=25.67`$ MeV at a Fermi momentum of $`\stackrel{~}{k}_F^0=1.15\mathrm{fm}^1`$.

While the success of the NL3 interaction in reproducing ground-state properties (such as binding energies, charge radii, energy separations, etc.) for a variety of nuclei all throughout the periodic table is well documented, we display in Table 1 ground-state properties for only the two nuclei of relevance to this contribution, namely, 90Zr and 208Pb. However, even these accurately calibrated models predict a wide range of values for the neutron skin of 208Pb, confirming that the neutron-skin of a heavy nucleus is not tightly constrained by known nuclear observables. The fifth column in the table displays the compression modulus of asymmetric nuclear matter with a neutron-proton asymmetry corresponding to 90Zr ($`b=0.111`$) and 208Pb ($`b=0.212`$). It is this quantity—not the compression modulus of symmetric nuclear matter—that is constrained by the breathing mode of nuclei. This simple fact makes the connection between the measured ISGMR and the compression modulus of symmetric nuclear matter sensitive to the density dependence of the symmetry energy. Recall that the compression modulus of symmetric nuclear matter was fixed in all models at $`K=271`$ MeV, yet for ($`b=0.212`$) asymmetric nuclear matter the compression modulus ranges from $`243`$ MeV (for the stiffest symmetry energy) all the way up to $`260`$ MeV (for the softest symmetry energy). Finally, the last column in the table shows peak and centroid energies for the ISGMR in 90Zr and 208Pb computed in a relativistic random-phase-approximation (RPA); Figs. 1 and 2 display the corresponding distribution of strength. As expected, there is a strong correlation between the centroid energies and the compression modulus of asymmetric nuclear matter. One should note in passing that the continuum RPA formalism employed here, but reported elsewhere Piekarewicz (2000, 2001), respects important symmetries of nature, such as translational invariance (in the form of Thouless’ theorem Thouless (1961); Dawson and Furnstahl (1990)) and the conservation of the vector current.

The great advantage of a nucleus such as 90Zr is that it has both a well developed isoscalar-monopole peak and a small neutron-proton asymmetry ($`b=0.111`$). The latter manifests itself into the near collapse of all curves in Fig. 1 into a single one, so that the former may directly constrain the compression modulus of symmetric nuclear matter. In contrast to 90Zr, the distribution of ISGMR strength in 208Pb is sensitive to the density dependence of the symmetry energy. While this sensitivity should be sufficient to constrain the density dependence of the symmetry energy, one could do even better. Indeed, one may constrain the density dependence of the symmetry energy by demanding that both the ISGMR and the isovector giant-dipole resonance (IVGDR) in 208Pb be simultaneously reproduced. The distribution of isovector dipole strength in 208Pb is displayed in Fig. 3. We note that the isovector-dipole response gets hardened as the symmetry energy is softened. As all models share the same value of the symmetry-energy coefficient at $`\stackrel{~}{k}_F^0=1.15\mathrm{fm}^1`$, the hardening of the response follows as a result of the symmetry energy being higher at the (low) densities relevant to the isovector-dipole mode Vretenar et al. (2003).

To constrain simultaneously the compression modulus of symmetric nuclear matter and the neutron radius of 208Pb, one starts by noticing that the theoretical centroid energy of the ISGMR in 90Zr overestimates the experimental value by about $`1`$ MeV. Although a proper adjustment of $`K`$ should be done through a re-calibration of parameters, a simple, yet accurate estimate may be obtained via the following scaling relation: $`E_{\mathrm{ISGMR}}\sqrt{K}`$ Blaizot et al. (1995). Using this relation and accounting for experimental uncertainties, an adjustment of about $`20`$ MeV in $`K`$ is required to reproduce the ISGMR in 90Zr. That is, $`K=271\mathrm{MeV}K=(248\pm 6)\mathrm{MeV}`$. This adjustment in $`K`$ induces a corresponding correction in the calculated values of the ISGMR in 208Pb. The centroid energies after correction, together with the peak energies of the IVGDR in 208Pb, are displayed in Fig. 4, alongside the experimental values Youngblood et al. (1999); Ritman and et al. (1993). The numbers in parenthesis indicate the predicted values for the neutron skin in 208Pb. The figure suggests that models with neutron skins in 208Pb larger than $`R_nR_p0.22`$ fm may have an unrealistically stiff symmetry energy.

In summary, relativistic mean-field models have been used to compute the distribution of isoscalar-monopole strength in 90Zr and 208Pb, and of isovector-dipole strength in 208Pb using a continuum RPA approach. It was demonstrated—using exclusively accurately calibrated models—that the distribution of isoscalar monopole strength in 208Pb is sensitive to the density dependence of the symmetry energy. Further, existing experimental data were used to set limits on both the compression modulus of symmetric nuclear matter and on the neutron skin of 208Pb. It appears that medium-mass nuclei, having a well-developed ISGMR peak and a small neutron-proton asymmetry (such as 90Zr but not 208Pb) allow for the best determination of the compression modulus of symmetric nuclear matter. In turn, the sensitivity of the ISGMR and the IVGDR in 208Pb to the density dependence of the symmetry energy may be used to impose constraints on the neutron skin of 208Pb. From the present analysis, a compression modulus of $`K=(248\pm 6)\mathrm{MeV}`$ and a neutron skin in 208Pb of $`R_nR_p0.22`$ fm were obtained. These values appear closer to those predicted in nonrelativistic studies.

We conclude with a comment on the impact of these results on the cooling of neutron stars. In earlier publications we have demonstrated how improved values for neutron radii could have a widespread impact on the structure and dynamics of neutron stars Horowitz and Piekarewicz (2001a, b); Horowitz and Piekarewicz (2002); Carriere et al. (2003). In particular, we suggested that the enhanced cooling of the neutron star in 3C58 Slane et al. (2002) may be due to the conventional URCA process—provided the symmetry energy is stiff enough to generate a neutron skin in 208Pb larger than $`R_nR_p0.24`$ fm. In view of our present findings, this now seems unlikely. Thus, the possibility that 3C58 harbors an exotic star, such as a quark star, looms large.

###### Acknowledgements.

The author is grateful to the organizers of the International Workshop on Nuclear Response under Extreme Conditions in Trento and to the ECT\* for their support and hospitality. This work was supported in part by the U.S. Department of Energy under Contract No.DE-FG05-92ER40750.


text:
# Specific Heat of Disordered Superconductors Induced by Negative Centers

## Abstract

We show that superconductors with inhomogeneous order parameters can show similar features as anisotropic ones. In this paper we study the low temperature specific heat dependence in such a system and we show that the disorder associated with randomly distributed attractive centers convert the BCS temperature exponential behaviour into a power law formula. To describe superconductivity we used a random version of the negative U Hubbard model, while the disorder was treated by means of the Coherent Potential Approximation (CPA).

thanks: Tel.: +48- 81- 5381573; Fax: +48- 81- 5241004; E-mail: litak@archimedes.pol.lublin.pl

According to the Anderson theorem the effect of non-magnetic disorder on s-wave superconductors can be neglected unless spatial fluctuations of the order parameter are negligible . This theorem for conventional superconductors can also be applied to other superconductors with a long coherence length where the fluctuations of the order parameter induced by disorder  are small and the quasiparticle density of states shows a clear gap. As a consequence, superconducting properties, such as the specific heat and the magnetic penetration depth, feature the typical exponential BCS dependences. The situation is quite different for exotic superconductors. There the order parameter is modulated by a wave vector k and even disappears (for d-wave symmetry of the order parameter) at Fermi surface line nodes.

In the case of the recently discovered MgB2  with a high critical temperature $`T_C40`$ K and a relatively short coherence length $`\xi `$ (from 25 to 120Å , depending on the sample ) the Anderson assumption of a uniform order parameter is also not satisfied. In fact, there is a lot of evidence  that in this compound the multi-gap scenario with a spatially fluctuating s-wave order parameter is realized. Moreover Sharma et al.  found a percolative transition of superconductivity in Mg1-xB2.

Prompted by these findings we examine the effects of a low temperature specific heat dependence using the random negative U Hubbard model proposed to describe percolation in superconductors . We study the case where $`U_i`$ is $`|U|`$ and 0 with probability $`c`$ and $`1c`$, respectively on a square lattice whose sites are labelled $`i`$ using the Gorkov decoupling  and the Coherent Potential Approximation (CPA) .

For the above assumptions we have done calculations of the order parameter $`\mathrm{\Delta }_A`$ for various concentrations of negative centers $`c`$ in a binary alloy AcB1-c . Here, for simplicity, we neglected the Hartree term $`Un/2`$. This term may be important in some other contexts, as it is known to lead to a band splitting phenomenon and consequently to the appearance of the critical concentration of a superconducting phase $`c_0`$ .

The temperature dependences of $`\mathrm{\Delta }_A`$ for $`U_A=8t`$ are presented in Fig. 1a, while examples of quasiparticle densities of states (QDOS) for different temperatures are depicted in Fig. 1b. Note that the pairing potential $`\mathrm{\Delta }_A`$ for $`c1`$ differs slightly from the uniform case (c=1), however the QDOS with states inside the gap is fundamentally different.

Now the low temperature specific heat can be easily evaluated using the following formula :

$$C=2T_{\mathrm{}}^{\mathrm{}}D(\omega )\frac{\omega }{T}\frac{\mathrm{d}}{\mathrm{d}T}\left(\frac{1}{\mathrm{e}^{\omega /T}+1}\right),$$
(1)
where $`D(\omega )`$ denotes QDOS.

Fig. 1c shows the results for the specific heat for $`c=1.0`$, 0.7, 0.5 and 0.3. Note that, in contrast to a clean case (’1’) where the corresponding plot can be associated with exponential formula $`C/T\mathrm{exp}(\alpha /T)`$ for constant $`\alpha `$, the next curves (’2-4’) mimic a distinct power law behaviour (excluding some lowest temperatures).

Thus, the behaviour of superconductors with randomly distributed negative centers showing fluctuations of the pairing potential in real space resembles, in some sense, the behaviour of anisotropic but clean systems with a pairing potential varying in k space. Namely, a wave vector modulation of the superconducting order parameter in the case of s-wave and most of all line nodes in the case of d-wave superconductors lead to similar effects with a low temperature unusual dependence of the specific heat.

This work has been partially supported by the KBN Grant No. 3P03B06225.


text:
# Zero Temperature Glass Transition in the Two-Dimensional Gauge Glass Model

## I Introduction

Despite intense studies of the superconducting glass transition in vortex systems, the precise nature of the transition remains controversial. The presence of quenched disorder in a strongly type II superconductor in an applied magnetic field may create a superconducting vortex glass phase with zero linear resistance, where the vortices are localized in random positions given by the disorder and interactions Fisher (1989). Considerable theoretical and simulation efforts have been focused on whether a superconducting glass state exists at finite temperature, or alternatively, only at zero temperature. In particular the situation in two dimensions still remains debated. In this paper we present results from Monte Carlo simulation and finite size scaling analysis for the glass transition in two dimensions. We focus mainly on the dynamic scaling properties, which have not been studied in detail before, using a dynamic scaling ansatz for a zero temperature transition. From our static and dynamic scaling analysis, we find evidence for a zero temperature glass transition, and present various predictions for the physical properties at the transitions.

The gauge glass (GG) model is an XY model with a random quenched vector potential on the links of the lattice; see Eq. (1) below. The GG model contains important aspects of order parameter symmetry, interactions, and quenched disorder, and has been frequently studied as a model for a vortex glass transition. In a real glass, the system may fall out of equilibrium at the transition, but here we focus on the equilibrium critical behavior of the model. This enables the study of universal scaling laws at the transition, that depend only on the universality class. It has not been settled if the GG model represents the correct universality class of the vortex glass transition, and it is thus motivated to obtain the precise critical behavior of the GG model. The GG model is not realistic in a few important aspects: disorder is included only in terms of a random vector potential, instead of a random pinning energy for the vortex cores, and it has zero applied magnetic field. Alternatively, a suitable modification of the GG model as a model for the elastic energy of a dislocated vortex lattice is discussed in Ref. Lidmar, 2003. Another worry is the use of a classical model for a zero temperature transition, where quantum effects may be important. However, since vortices localize in the glass phase, their quantum fluctuations are small and can be neglected, and a classical treatment is a reasonable starting point, also at $`T=0`$. This neglects the possibility of quantum tunneling of vortices, which may actually be important.

Considerable simulation efforts have been invested in trying to locate the lower critical dimension $`d_c`$ of the vortex glass transition, i.e., the limiting dimensionality such that a glass phase exists at finite temperature only for $`d>d_c`$. In the 3D case, several simulations Reger et al. (1991); Olson and Young (2000); Katzgraber and Young (2002) show that a vortex glass phase does exist at finite temperatures in the absence of magnetic screening. In 2D the situation is less clear. It has been shown rigorously that there can be no vortex-glass order in 2D for $`T>0`$ Nishimori (1994). For a finite screening length it is expected that no glass phase can exist at $`T>0`$ Bokil and Young (1995). Akino and Kosterlitz computed the domain wall stiffness exponent $`\theta `$ at zero temperature, and obtained a negative value, $`\theta =0.36`$, indicating that $`T_c=0`$ Akino and Kosterlitz (2002). Holme and Olsson computed the energy barrier for vortex motion across the system and argue that $`T_c>0`$ Holme and Olsson (2002). Katzgraber and Young recently obtained $`T_c=0`$ and $`\theta =1/\nu =0.39\pm 0.03`$ from finite temperature simulations Katzgraber and Young (2002); Katzgraber (2003). In contrast, several other studies of static quantities obtained $`T_c0.2`$ Holme et al. (2003); Choi and Park (1999). Dynamic scaling in the 2D GG model has been studied in several papers Li (1992); Kim (2000); Choi and Park (1999); Hyman et al. (1995); Wengel and Young (1996); Katzgraber and Campbell (2004), with conflicting results. It is important to consider dynamic scaling in order to obtain information about physically measurable quantities like the resistance. Previous studies of zero temperature scaling of dynamic properties considered nonlinear $`IV`$ characteristics Hyman et al. (1995), and in addition the linear resistance was found to be roughly consistent with a simple Arrhenius form Hyman et al. (1995); Wengel and Young (1996), $`Re^{K/T}`$, but no further support for zero temperature dynamic scaling has been obtained. In contrast, Refs. Li, 1992; Kim, 2000; Choi and Park, 1999 obtained $`T_c>0`$ from vortex dynamics simulations. Hence the behavior of the GG model remain controversial, in particular for the dynamic properties.

In this paper we study critical scaling properties at the glass transition for both dynamic and static quantities, assuming a zero temperature transition. Our simulations are motivated in part by new methods that have recently been developed. At low temperature the correlation times in the simulation grow very long and efficient convergence acceleration methods become highly desirable. The temperature exchange method Hukushima and Nemoto (1996); Marinari (1998) has proven very useful for such problems, where configurations are exchanged between different temperatures, in order to avoid getting stuck in metastable states. This considerably speeds up the convergence of the simulation, and thereby enables simulation at larger system sizes and lower temperatures. Our scaling results for static quantities are in good agreement with some of the previous results in the literature for the correlation length exponent and the domain wall energy exponent $`\theta `$ Hyman et al. (1995); Wengel and Young (1996); Katzgraber and Young (2002); Katzgraber (2003); Akino and Kosterlitz (2002). For our dynamics simulations we also get valuable input from the exchange method, even though exchange is not directly applicable for dynamics, by using typical equilibrium configurations from the exchange runs as input to the dynamics runs. This considerably extends the regime of temperatures and system sizes accessible for the dynamics simulation. The main new results of our work are for the dynamic finite size scaling analysis of the linear resistance and of the dynamic correlation function for the order parameter at the zero temperature glass transition. We obtain correlation times from finite size scaling of our dynamic MC data that show novel behavior. We find evidence for two distinct diverging correlation times, one associated with phase slip fluctuations, and the other associated with local phase fluctuations that are independent of the phase slip processes.

The organization of the paper is as follows: In Sec. II we introduce the GG model and describe our simulation methods. In Sec. III we discuss the quantities calculated in the simulation and the finite size scaling methods. In Sec. IV the static results are presented, and Sec. V contains the results for the dynamic quantities. Section VI contains discussion and conclusions.

## II Model and Monte Carlo simulation details

The gauge glass (GG) model Shih et al. (1984); Huse and Seung (1990); Fisher et al. (1991) is defined by the Hamiltonian

$$H=\underset{i,j}{}\mathrm{cos}(\theta _i\theta _jA_{ij})$$
(1)
where $`\theta _i`$ is the phase of the superconducting order parameter on site $`i`$ of a square lattice, and $`_{i,j}`$ denotes summation over all nearest-neighbor pairs in the lattice. The vector potential enters through $`A_{ij}=\frac{2\pi }{\varphi _0}_i^j𝑑𝐫𝐀`$, which is a quenched random variable on each link in the lattice, with a uniform probability distribution in $`[0,2\pi )`$, corresponding to fixed random magnetic induction $`𝐁=\times 𝐀`$ through the plaquettes of the lattice. The partition function is $`Z=\left(_j_0^{2\pi }𝑑\theta _j\right)\mathrm{exp}(H/T)`$, where $`T`$ is the temperature. The simulation uses square lattices of size $`L\times L`$ with periodic boundary conditions (PBC) in both the $`x`$ and $`y`$-directions.

We also consider the GG model with different boundary conditions, where the system repeats periodically, but with fluctuating phase twist variables $`\mathrm{\Delta }_x,\mathrm{\Delta }_y`$ added to the respective phase differences across the boundary of the system. This allows for phase twist fluctuations across the system, and enables calculation of the resistance of the model, as is discussed below. In this case the Hamiltonian is

$$H=\underset{i,j}{}\mathrm{cos}(\theta _i\theta _jA_{ij}\delta _{x_i,L}\mathrm{\Delta }_x\delta _{y_i,L}\mathrm{\Delta }_y)$$
(2)
where $`\mathrm{\Delta }_x,\mathrm{\Delta }_y`$ are integrated over $`(\mathrm{},\mathrm{})`$ in the partition function. In a static simulation, the partition function is invariant if the the phase twist at the boundary is spread out in the bulk as a constant extra phase difference $`\mathrm{\Delta }/L`$ on each link. However, dynamically these both prescriptions are different. We define the phase twist fluctuations located at the boundary in order to maintain local dynamics.

In the Monte Carlo simulations the trial moves are attempts to assign the phase, at a randomly chosen site, a random value from a uniform distribution in $`[0,2\pi )`$. The moves are accepted with probability $`1/(1+\mathrm{exp}\mathrm{\Delta }E/T)`$, where $`\mathrm{\Delta }E`$ is the energy difference for the trial move. One sweep is defined as on average one attempt to update each phase. In the simulations of static quantities we use a temperature exchange method Hukushima and Nemoto (1996); Marinari (1998), that considerably enhances the convergence of the simulations at low temperature and big system sizes. We use ten sweeps of local MC moves, followed by one exchange sweep. In the exchange sweep the trial moves are attempts to exchange configurations between neighboring temperatures, where up to 39 closely spaced temperatures, with $`T_n=T_{n1}(11/L)`$, are simulated in parallel, sharing the same disorder realization. The exchange moves are accepted with probability $`\mathrm{exp}[(\beta _n\beta _{n+1})(E_{n+1}E_n)]`$, where $`\beta =1/T`$. To reach equilibrium we discard up to $`610^4`$ MC sweeps, followed by equally many sweeps where data is collected. The results are averaged over up to $`10^4`$ realizations of the random vector potential $`A_{ij}`$. To avoid systematic errors in the calculation of squares of expectation values, two replicas of the system with the same disorder are used. We use various tests to verify that sufficient equilibration is achieved, which will be presented in the results section below. In this paper we denote thermal averages by $`\mathrm{}`$, and disorder averages by $`[\mathrm{}]`$.

In the simulation of dynamic quantities the exchange MC method does not apply directly, since we must use a method that respects the local nature of the dynamics of the system. Instead we use the usual Monte Carlo scheme with local updates of the phase. However, also for the dynamics runs we take advantage of the exchange algorithm in two ways. First, we use equilibrium configurations from an exchange run as initial configurations for the dynamics runs. Second, during the dynamics simulation we also recalculate the static quantities, and compare to the same quantities from the exchange runs. This gives an additional test of convergence for the dynamics simulation. We define one time step to be one MC sweep through the system with local updates. In the simulation with fluctuating twists, one time step corresponds to one sweep with local phase updates, followed by one update of $`\mathrm{\Delta }_x`$ and of $`\mathrm{\Delta }_y`$. The phase twist updates are attempts to add a uniformly distributed random number in $`(\pi ,\pi )`$. As initial configurations in the dynamics runs, we use configurations generated in an exchange run, followed by up to $`410^5`$ sweeps of local MC moves where data is taken. Occasional runs with even more sweeps were also used to test convergence. Disorder averages are formed over up to $`10^4`$ realizations of the disorder. We stress that there are many dynamic universality classes, and here we only consider that of local relaxation of the phase, which has proven useful in previous studies of the superconducting phase transition Hyman et al. (1995); Lidmar et al. (1998).

## III Calculated quantities and finite size scaling relations

The basic quantity in our analysis of thermodynamic properties of the glass transition is the root mean square (RMS) current Reger et al. (1991); Olson and Young (2000). The current in the $`x`$-direction is defined as the derivative of the free energy $`F`$ with respect to a phase twist $`\mathrm{\Delta }_x`$ on the boundary in the $`x`$-direction, for $`\mathrm{\Delta }_x=0`$:

$$I=\frac{F}{\mathrm{\Delta }_x}|_{\mathrm{\Delta }_x=0}=\frac{1}{L}\underset{j}{}\mathrm{sin}(\theta _{j+\widehat{x}}\theta _jA_{j,j+\widehat{x}})$$
(3)
where $`j+\widehat{x}`$ denotes the nearest neighbor (NN) site next to $`j`$ in the $`x`$-direction. Here the phase twist $`\mathrm{\Delta }_x`$ at the boundary has been absorbed into the bulk by adding a uniform phase twist $`A=\mathrm{\Delta }_x/L`$ on each link in the $`x`$-direction, which gives the $`1/L`$ factor in the last equation. The RMS current is defined as

$$I_{\mathrm{rms}}=[I^aI^b]^{1/2}$$
(4)
where $`a,b`$ denote two replicas of the system sharing identical disorder.

In order to use the scaling properties of MC data to discriminate between the possibilities that $`T_c=0`$ and $`T_c>0`$, we must formulate suitable scaling assumptions for each case Hyman et al. (1995); Katzgraber and Young (2002). At the glass transition the correlation length is assumed to diverge as $`\xi |TT_c|^\nu `$ if $`T_c>0`$, and as $`\xi T^\nu `$ if $`T_c=0`$. In the case $`T_c>0`$, we assume that close to $`T_c`$ the RMS current obeys the finite size scaling relation Olson and Young (2000)

$$I_{\mathrm{rms}}(T,L)=\stackrel{~}{I}(L^{1/\nu }(TT_c))$$
(5)
where $`\stackrel{~}{I}`$ is a universal scaling function. This scaling form implies that the RMS current is independent of system size $`L`$ at $`T=T_c`$, which gives a convenient criterion for locating $`T_c`$. In contrast, if $`T_c=0`$, the scaling relation becomes modified, since for $`T_c=0`$ we have $`\xi T^\nu =\beta ^\nu `$, and hence $`\beta \xi ^{1/\nu }`$, which must be included in the scaling relation. This gives $`\beta I=\beta F/\mathrm{\Delta }\xi ^0`$. The scaling ansatz for the RMS current valid for the case $`T_c=0`$ becomes

$$I_{\mathrm{rms}}(T,L)=L^{1/\nu }\stackrel{~}{I}(L^{1/\nu }T)$$
(6)
Below we will compare both these scaling forms to our data, and find that the scaling form with $`T_c>0`$ is not fulfilled by our MC data, while the scaling form with $`T_c=0`$ works well.

Next we discuss the quantities that we calculate in the dynamic simulation, and their respective scaling relations. The aim of our analysis is to obtain scaling relations from a dynamic scaling ansatz that assumes $`T_c=0`$, and here we will only formulate the scaling relations for this case. The scaling relations will first be formulated in terms of the correlation time, $`\tau `$, whose scaling will in turn be discussed below.

The resistance can be calculated by using fluctuating twist boundary conditions instead of PBC. It is obtained from the Kubo formula for the voltage-voltage correlation function. We define a time dependent resistance function as

$$R(t,T,L)=\frac{1}{2T}\underset{t^{}=t}{\overset{t}{}}\mathrm{\Delta }t^{}\left[V(t^{})V(0)\right]$$
(7)
where the voltage at MC timestep $`t`$ is given by the Josephson relation $`V(t)=\frac{\mathrm{}}{2e}\frac{d\mathrm{\Delta }_x}{dt}`$ (we henceforth assume units such that $`\mathrm{}/2e=1`$). Here $`d\mathrm{\Delta }_x/dt`$ is the rate of change of the phase shift across the system in the $`x`$-direction, which gives the resistance in the $`x`$-direction. The usual Kubo formula Young (1994) for the resistance is obtained from this function in the limit when the summation over MC timesteps $`t\mathrm{}`$. In practice this limit means that the summation time has to be long enough that the resistance is independent of $`t`$. However, since this time is going to be very long in the simulation, it is quite useful to consider $`R`$ as a function of $`t`$ and study the scaling properties as function of $`t`$. Since the voltage scales as $`V1/\tau `$, we make the finite size scaling ansatz

$$R(t,T,L)=\frac{1}{\tau T}\stackrel{~}{R}(L^{1/\nu }T,t/\tau )$$
(8)
where $`\stackrel{~}{R}`$ is a scaling function.

Another useful quantity is the dynamic correlation function for the real part of the superconducting order parameter, $`\mathrm{\Psi }^{}=_j\mathrm{cos}\theta _j`$, which is gauge dependent and hence not directly observable in a superconductor, but gives a convenient measure of the dynamic correlations of the local order in the system. We define the time dependent summed autocorrelation function

$$G(t,T,L)=\underset{t^{}=t}{\overset{t}{}}\mathrm{\Delta }t^{}\left[\frac{\mathrm{\Psi }^{}(t^{})\mathrm{\Psi }^{}(0)}{\mathrm{\Psi }^2}\right]$$
(9)
For $`t\mathrm{}`$, $`G`$ measures the correlation time $`\tau `$ of the fluctuations in $`\mathrm{\Psi }^{}`$, and we have

$$G(t,T,L)=\tau \stackrel{~}{G}(L^{1/\nu }T,t/\tau )$$
(10)
Note that all the (unknown) dynamic scaling functions contain two arguments, which complicates their analysis. We will make use of two possible ways to obtain a function of one variable only, and this function can then be conveniently analyzed. One way is to obtain data for $`t\tau `$, and the second way is to fix the first argument to a constant value by using different temperatures for each system size so that $`L^{1/\nu }T=const`$ Hyman et al. (1995), which then permits scaling analysis of the entire $`t/\tau `$ dependence.

In order to quantitatively fit the dynamic scaling formulas to MC data, we need an ansatz for the form of $`\tau `$. This is related to the divergence of the free energy barriers at the glass transition. The thermodynamic free energy barrier exponent $`\theta `$ is defined from the interfacial free energy RMS difference, $`U_L`$, for changing the boundary condition in one direction from PBC to anti-PBC. For a system of size $`L`$ we have $`U_L=KL^\theta `$, where $`K`$ is a constant. Now, a change in free energy from a change in the boundary conditions must scale in the same way as the current in Eqs. (3) and (6), which gives $`\theta =1/\nu `$. Thus $`\theta `$ can be obtained from analysis of $`I`$ at finite temperatures.

A dynamic barrier exponent $`\psi `$ can be defined by $`V_L=KL^\psi `$, which is the free energy cost to move a vortex a distance $`L`$ using only local moves. In general $`V_L`$ and $`U_L`$ do not need to coincide, and the dynamic barrier for vortex motion across the system can be very big. We assume that the correlation time scales in an activated form given by

$$\tau \mathrm{exp}(V_L/T)\mathrm{exp}(C/T^{1+\nu \psi })$$
(11)
where $`C`$ is a constant. If $`\psi =0`$ this reduces to a simple Arrhenius form, $`\tau \mathrm{exp}C/T`$. We also consider a few other possibilities. The usual finite temperature form of the critical slowing down is $`\tau L^z`$. Another possible form is a logarithmic barrier, $`\tau \mathrm{exp}K(\mathrm{ln}L)/T`$.

## IV Results for static quantities

We first study the approach to equilibration of the exchange MC simulation. This is necessary in order to avoid systematic errors due to insufficient equilibration. We consider the RMS current, $`I_{\mathrm{rms}}`$, defined in Eq. (4). Figure 1 shows some of our data for the RMS current as a function of simulation warmup time. The data points in the figure are calculated by first discarding $`t`$ MC sweeps, and then forming averages over equally many sweeps. In Fig. 1 we see that the equilibration time, where the curves become independent of $`t`$, becomes very long already for moderate system sizes $`L`$, despite the considerable convergence speed-up due to the exchange method.

Next we consider finite size scaling properties of the RMS current. Figure 2 shows MC data for the RMS current vs. temperature. According to Eq. (5) data curves for $`I_{\mathrm{rms}}`$ should become system size independent and thus intersect at $`T=T_c`$, if $`T_c>0`$. The absence of any such intersection in the figure indicates that $`T_c`$ is much smaller than the temperatures where we have data, i.e., $`T_c0.1`$. Henceforth we will therefore assume that $`T_c=0`$.

To scale MC data for $`I_{\mathrm{rms}}`$ assuming that $`T_c=0`$, Eq. (6) shows that a plot of $`L^{1/\nu }I_{\mathrm{rms}}`$ vs. $`L^{1/\nu }T`$ for different $`T,L`$ should collapse onto a system size independent curve, for the correct value of $`1/\nu `$. Scaling plots are shown in Fig. 3 for the choices (a) $`1/\nu =0.39`$ and (b) $`1/\nu =0.36`$. In (a) we see that $`1/\nu =0.39`$ gives the best overall data collapse, but the fit is somewhat poor for the smallest values of $`L^{1/\nu }T`$, which for $`T_c=0`$ is expected to be the most significant part of the data collapse. We therefore drop data points for large $`L^{1/\nu }T`$ and only fit the value of $`\nu `$ to the data points with smallest $`L^{1/\nu }T`$. Fitting the data points at the two lowest temperatures for sizes $`L=10,12,16`$ to a straight line, and minimizing the RMS fit error, gives $`\nu =0.36\pm 0.03`$. The error bar is obtained by the bootstrap method Newman and Barkema (1999). The resulting data collapse is shown in Fig. 3 (b), and we see that, for small $`L^{1/\nu }T`$, a better fit is obtained than in (a). The estimate $`\nu =0.36\pm 0.03`$ is insensitive within the error bar to varying the precise details of the fit, such as including the three lowest temperatures for each $`L`$.

## V Results for dynamic quantities

The dynamics simulations contain only local MC moves, and the averages computed in the dynamics runs correspond to time integration of the observables. The statics simulation contains both local MC moves and nonlocal temperature exchange moves that take big steps in configuration space, which corresponds to ensemble averaging. If ergodicity is broken these different averages may disagree Holme and Olsson (2002). To investigate this issue we plot the RMS current calculated in a dynamics run in Fig. 2, together with the corresponding result from the static simulation. The figure shows that the results from static and dynamic runs coincide within error bars.

We will make use of the following input from the statics simulations. The static scaling results gives the value of the correlation length exponent, which is needed also in the dynamic scaling analysis. We will assume the value $`1/\nu =0.39`$, but make occasional comparison with the corresponding results for $`1/\nu =0.36`$. Furthermore, the static scaling results also give an indication of the range of system sizes and temperatures where we can expect dynamic scaling of the data to hold. From Fig. 3 we estimate that dynamic scaling can be expected when $`L^{1/\nu }T`$ is smaller than about 0.8.

We expect that the longest emergent time scale comes from phase slip events corresponding to vortex motion across the sample. An illustration of these long correlation times is seen in Fig. 4, which shows phase slip events in a dynamic simulation where fluctuations in the boundary conditions are included, for a few different realizations of the quenched disorder. Phase slip events are often as far as about $`10^5`$ sweeps apart.

Figure 5 shows the time dependent resistance function given by Eq. (7) for $`L^{0.39}T=0.8`$. The data curves show unexpectedly complicated behavior, with a power-law regime at short times, and a time-independent limit at long times, with various crossover times in between. We will now analyze the scaling behavior of these different time scales.

We first consider the long time limit, $`t\tau `$, of the data curves in Fig. 5, which gives the resistance $`R`$. According to Eq. (8), the correct value of $`\tau `$ for each system size and temperature, such that $`L^{1/\nu }T=A`$ is a constant, should make data curves for $`RT\tau `$ plotted vs. $`t/\tau `$ collapse onto a common curve for $`t\tau `$. Thus we adjust $`\tau `$ for each data curve to the value that gives the best collapse. The result is shown in Fig. 6. We obtain a data collapse to a common curve at large $`t`$, for each set of data curves belonging to the same value for $`A=L^{1/\nu }T`$. At shorter times the collapse breaks down, indicating a different short time behavior. By this procedure we determine an estimate of the correlation time $`\tau (T,L)`$, up to an unknown multiplicative constant that depends on $`A=L^{1/\nu }T`$. The scaling result shown in this figure demonstrates that, up to the level of statistical accuracy of the simulation data, the resistance obeys a dynamic scaling relation which assumes $`T_c=0`$. This result is independent of any assumption about the functional form of $`\tau (L,T)`$.

To proceed further, we need to fit the values for $`\tau `$ and we then need an assumption of a functional form for $`\tau `$. All the functional forms $`\tau \mathrm{exp}K(\mathrm{ln}L)/T`$, $`\tau \mathrm{exp}KL^\psi /T`$, and $`\tau L^z`$ approximately fit the data, and our present data is therefore not sufficient to clearly distinguish between these functional forms. Larger system sizes and lower temperatures are necessary in order to discriminate between them. The form $`\tau \mathrm{exp}K(\mathrm{ln}L)/T`$ gives a good fit to the data, with $`K0.56`$ independent of $`T`$ and $`L`$, shown in Fig. 7. However, this should only be considered as one possible parameterization of the correlation time.

Next we consider the scaling behavior of $`R(t,T,L)`$ for intermediate times. Surprisingly the scaling found for long times in Fig. 6 is not obeyed for intermediate times where the data collapse breaks up. We therefore attempt another scaling relation to apply in this regime. According to Eq. (8), we again set $`L^{1/\nu }T=A`$ to a constant value, and plot $`RT\tau `$ vs. $`t/\tau `$. We can then obtain an approximate data collapse in the intermediate power-law regime in the figure for $`\tau L^z`$ with $`z=2.7\pm 0.2`$, where the error bar is estimated from the interval outside of which a significantly worse fit is obtained. We can also fit the same data with other functional forms for $`\tau `$, but we get the best fit with $`\tau L^z`$. These results show that the resistance has two separate dynamic scaling regimes, characterized by two separate diverging correlation times. It is natural to associate the longest of these correlation times with the phase slip processes that make the resistance finite at any finite temperature, since the correlation time roughly matches the time between phase slips in Fig. 4. We will now study these different dynamic scaling regimes in more detail.

To reach further insight into the scaling properties at intermediate times we consider the order parameter autocorrelation function $`G(t,T,L)`$, defined in Eq. (9). Figure 9 shows MC data for $`G(t,T,L)`$ scaled according to Eq. (10) with $`\tau L^z`$. The data in the figure has $`L^{1/\nu }T=0.8`$, and the data collapse is obtained for $`z=2.74`$. The data in the figure is for the case of fluctuating phase twists at the boundaries. We also calculated the same quantity with periodic boundary conditions, and obtained an equally good fit with the same exponent. Approximately the same exponent also gives a collapse for other values of $`A`$, but a slight drift in $`z`$ is noticeable, just as the data collapse in Fig. 3 for the RMS current is somewhat dependent on $`A`$. Including this uncertainty we estimate $`z=2.7\pm 0.2`$. We conclude that the dynamics measured by $`G`$ is independent of the phase twist fluctuations. The dynamic scaling relation $`\tau L^{2.7}`$, instead comes from local dynamic fluctuations other than the global phase slips. This means that we can interpret the intermediate resistance scaling regime in Fig. 6 as coming from local fluctuations rather than global phase slips.

## VI Discussion and conclusions

This paper considers finite size scaling properties of Monte Carlo data for the glass transition in the two-dimensional gauge glass model. Equilibration of our simulations at low temperatures is achieved by a temperature exchange simulation method. Up to the level of statistical uncertainty, system sizes and temperatures reached in the simulations, our finite size scaling results for both static and dynamic properties show scaling properties that are clearly consistent with a glass transition occurring only at zero temperature, and thus indicates the absence of a superconducting glass state at any finite temperature in 2D.

Our results for static scaling properties reproduces, and adds some details to, some of the results of Refs. Katzgraber and Young, 2002; Katzgraber, 2003, in addition to providing a useful starting point for our dynamics simulations. From the scaling of the RMS current we find the correlation length exponent given by $`1/\nu 0.39`$, obtained by optimizing the quality of a finite size scaling data collapse over a finite temperature interval close to $`T=0`$. This value agrees with that of Refs. Katzgraber and Young, 2002; Katzgraber, 2003. If we instead optimize the fit at the lowest temperatures where we have data, where scaling is also expected to work best, we obtain $`1/\nu =0.36\pm 0.03`$. It is interesting to note that this is in very good agreement with the result for the ground state domain wall exponent $`\theta =1/\nu `$ in Ref. Akino and Kosterlitz, 2002. To obtain an even more accurate estimate of $`\nu `$, data for lower temperatures and bigger system sizes is necessary. To converge such data would require a considerable increase in simulation efforts.

For the dynamic properties, we consider the linear resistance computed from the diffusion of the phase slip across the boundary of the system. Since the correlation times grow very long at low temperature, it is not possible to converge the simulation at as low temperature as for the static data where the temperature exchange algorithm helps considerably, which makes it more difficult to study dynamic scaling than static scaling. A further complication is that our data for the linear resistance shows more than one characteristic time scale. We identify the longest correlation time as the characteristic time scale for phase slip fluctuations across the sample. Unfortunately it is not possible to determine the correlation time uniquely from the data, since the interval of sizes and temperatures where we have data is too limited. All the forms $`\tau \mathrm{exp}KL^\psi /T`$, $`\tau \mathrm{exp}K\mathrm{ln}L/T`$, and $`\tau L^z`$ can be made to approximately fit the data. To discriminate between the possible forms, data at lower $`T`$ and larger $`L`$ is needed, which is beyond the present calculation. However, a logarithmically growing phase slip barrier, with $`\tau \mathrm{exp}K\mathrm{ln}L/T`$, gives a good fit to the data with $`K0.56`$. A qualitatively similar result was obtained in a zero temperature calculation of the phase slip barrier energy in Ref. Holme and Olsson, 2002. Since $`R\mathrm{exp}KL^\psi /T`$, of which the logarithmic barrier is a special case with $`\psi =0`$, a superconducting glass with zero resistance is obtained only for $`T=0`$. In other words the dynamic glass transition occurs at zero temperature.

Surprisingly, data for the time dependent resistance function $`R(t,T,L)`$ does not scale with the longest correlation time for intermediate times. Instead, in this regime scaling is obtained by assuming a new correlation time of the form $`\tau L^z`$ with $`z=2.7\pm 0.2`$. Hence the scaling behavior of our data implies two different diverging correlation times with quite different properties at the zero temperature glass transition. A further indication of the smaller correlation time is given by the scaling of the order parameter autocorrelation function, which scales over the entire time interval where we have data, including the longest times, quite in contrast to the resistance. The interpretation of this shorter correlation time is that it gives the characteristic time scale for fluctuations other than the global phase slip fluctuations that lead to a finite resistance. Similar results showing different dynamics for local fluctuations and phase twist fluctuations have been obtained previously at the finite temperature transition in XY models without disorder Jensen et al. (2000). We also note that the scaling relation $`\tau L^z`$ can be written $`\tau \mathrm{exp}z\mathrm{ln}L\mathrm{exp}zL^0\mathrm{exp}C/T^{1+\psi \nu }`$, which means that $`\psi =\theta =1/\nu `$. Hence the relation $`\tau L^z`$ suggests that the static and dynamic free energy barrier exponents are the same for the local fluctuations of the order parameter.

The two different diverging correlation times suggested here may explain why signatures of a glass phase with zero resistance at finite temperatures can easily be obtained in simulations. For the shorter correlation time we have $`\psi =\theta =1/\nu `$, which gives $`R(T)\mathrm{exp}C/T^{1+\psi \nu }=const`$, for the resistance at small temperatures, i.e. the resistance stays at the $`T=0`$ value $`R=0`$ for small $`T`$. This shows that in a regime of intermediate time scales, zero resistance can be obtained at finite temperature. However, at much longer time scales that may be hard to reach in a simulation, phase slip fluctuation set in, and qualitatively different behavior follows with a superconducting glass transition only at zero temperature.

In summary, from finite size scaling of a time dependent resistance function and an order parameter autocorrelation function, we constructed two different diverging correlation times from the $`T_c=0`$ scaling behavior of Monte Carlo data. For further work, it would be of interest to further clarify the precise details of dynamic scaling at the zero temperature glass transition, and to include quantum fluctuations in the model. It would furthermore be interesting to search for the proposed novel dynamic scaling behavior in experiments.

We acknowledge very helpful discussions with Petter Minnhagen, Helmut Katzgraber, Peter Olsson, Stephen Teitel, and Peter Young. This work was supported by the Swedish Research Council, PDC, NSC, and the Göran Gustafsson foundation.


text:
# Search for Lepton-Flavor Violation in the Decay 𝜏⁻→ℓ⁻ℓ⁺ℓ⁻

BABAR-PUB-03/044

SLAC-PUB-10256

thanks: Deceased

The BABAR Collaboration

## Abstract

A search for the lepton-flavor-violating decay of the tau into three charged leptons has been performed using 91.5$`\text{ fb}^1`$ of data collected at an $`e^+e^{}`$ center-of-mass energy of 10.58$`\mathrm{GeV}`$ with the BABAR detector at the PEP-II storage ring. In all six decay modes considered, the numbers of events found in data are compatible with the background expectations. Upper limits on the branching fractions are set in the range $`(13)\times 10^7`$ at 90% confidence level.

preprint: BABAR-PUB-03/044preprint: SLAC-PUB-10256

Lepton-flavor violation (LFV) involving charged leptons has never been observed, and stringent experimental limits exist from muon branching fractions: $`(\mu e\gamma )<1.2\times 10^{11}`$ brooks99  and $`(\mu eee)<1.0\times 10^{12}`$ sindrum88  at 90% confidence level (CL). Recent results from neutrino oscillation experiments neut  show that LFV does indeed occur, although the branching fractions expected in charged lepton decay due to neutrino mixing alone are probably no more than $`10^{14}`$ pham98 .

In tau decays, the most stringent limit on LFV is $`(\tau \mu \gamma )<3.1\times 10^7`$ at 90% CL belle03 . Many extensions to the Standard Model (SM), particularly models seeking to describe neutrino mixing, predict enhanced LFV in tau decays over muon decays with branching fractions from $`10^{10}`$ up to the current experimental limits ma02 . Observation of LFV in tau decays would be a clear signature of non-SM physics, while improved limits will provide further constraints on theoretical models.

This analysis is based on data recorded by the BABAR detector at the PEP-II asymmetric-energy $`e^+e^{}`$ storage ring operated at the Stanford Linear Accelerator Center. The data sample consists of 81.9$`\text{ fb}^1`$ recorded at $`\sqrt{s}=10.58\mathrm{GeV}`$ and 9.6$`\text{ fb}^1`$ recorded at $`\sqrt{s}=10.54\mathrm{GeV}`$. With an expected cross section for tau pairs at the luminosity-weighted $`\sqrt{s}`$ of $`\sigma _{\tau \tau }=(0.89\pm 0.02)`$ nb kk , this data sample contains over 160 million tau decays.

The BABAR detector is described in detail in Ref. detector . Charged-particle (track) momenta are measured with a 5-layer double-sided silicon vertex tracker and a 40-layer drift chamber inside a 1.5-T superconducting solenoidal magnet. The transverse momentum resolution parameterized as $`\sigma _{p_T}/p_T=(0.13p_T/[\mathrm{GeV}/c]+0.45)\%`$ is achieved. An electromagnetic calorimeter consisting of 6580 CsI(Tl) crystals is used to identify electrons and photons, a ring-imaging Cherenkov detector is used to identify charged hadrons, and the instrumented magnetic flux return (IFR) is used to identify muons. Particle attributes are reconstructed in the laboratory frame and then boosted to the $`e^+e^{}`$ center-of-mass (CM) frame using the measured asymmetric beam energies.

This paper presents a search for LFV in the neutrinoless decay $`\tau ^{}\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$. All possible lepton combinations consistent with charge conservation are considered, leading to six distinct decay modes ($`e^{}e^+e^{}`$, $`\mu ^+e^{}e^{}`$, $`\mu ^{}e^+e^{}`$, $`e^+\mu ^{}\mu ^{}`$, $`e^{}\mu ^+\mu ^{}`$, $`\mu ^{}\mu ^+\mu ^{}`$) cc . The signature of this process is three charged particles, each identified as either an electron or muon, with an invariant mass and energy equal to that of the parent tau lepton. Candidate signal events in this analysis are required to have a “1-3 topology,” where one tau decay yields three charged particles (3-prong), while the second tau decay yields one charged particle (1-prong). Four well reconstructed tracks are required with zero net charge, pointing towards a common region consistent with $`\tau ^+\tau ^{}`$ production and decay. One of these tracks must be separated from the other three by at least 90 in the CM frame. The plane perpendicular to this isolated track divides the event into two hemispheres and defines the 1-3 topology. Pairs of oppositely charged tracks identified as photon conversions in the detector material with an $`e^+e^{}`$ invariant mass below 30$`\mathrm{MeV}/c^2`$ are ignored.

Each of the charged particles found in the 3-prong hemisphere must be identified as either an electron or muon candidate. Electrons are identified using the ratio of calorimeter energy to track momentum $`(E/p)`$, the ionization loss in the tracking system $`(\mathrm{d}E/\mathrm{d}x)`$, and the shape of the shower in the calorimeter. Muons are identified by hits in the IFR and small energy deposits in the calorimeter. Muons with momentum less than $`0.5\mathrm{GeV}/c`$ cannot be identified because they do not penetrate far enough into the IFR.

The particle identification (PID) requirements are not sufficient to suppress certain backgrounds, particularly those from higher-order radiative Bhabha and $`\mu ^+\mu ^{}`$ events that can have four leptons in the final state. To reduce these backgrounds, additional selection criteria are applied to the six different decay modes. For all decay modes, the momentum of the 1-prong track is required to be less than 4.8 GeV/c in the CM frame. For the $`e^{}e^+e^{}`$ and $`e^{}\mu ^+\mu ^{}`$ decay modes, the charged particle in the 1-prong hemisphere must not be identified as an electron, while for the $`\mu ^{}e^+e^{}`$ and $`\mu ^{}\mu ^+\mu ^{}`$ decay modes it must not be a muon. For all four of these decay modes, the angle $`\theta _{13}`$ between the 1-prong momentum and the vector sum of the 3-prong momenta in the CM frame must satisfy $`\mathrm{cos}\theta _{13}>0.9999`$, while the net transverse momentum of the four tracks must be greater than 0.1$`\mathrm{GeV}/c`$. Additional requirements are imposed to reduce the $`q\overline{q}`$ and SM $`\tau ^+\tau ^{}`$ backgrounds. Events in the four decay modes specified above are required to have no unassociated calorimeter clusters (photons) in the 3-prong hemisphere with energy greater than $`100\mathrm{MeV}`$ in the laboratory frame, while events in all six decay modes are required to have no track in the 3-prong hemisphere that is also consistent with being a kaon.

To reduce backgrounds further, signal events are required to have an invariant mass and total energy in the 3-prong hemisphere consistent with a parent tau lepton. These quantities are calculated from the observed track momenta assuming the corresponding lepton masses for each decay mode. The energy difference is defined as $`\mathrm{\Delta }EE_{\mathrm{rec}}^{}E_{\mathrm{beam}}^{}`$, where $`E_{\mathrm{rec}}^{}`$ is the total energy of the tracks observed in the 3-prong hemisphere and $`E_{\mathrm{beam}}^{}`$ is the beam energy, both in the CM frame. The mass difference is defined as $`\mathrm{\Delta }MM_{\mathrm{rec}}m_\tau `$ where $`M_{\mathrm{rec}}`$ is the reconstructed invariant mass of the three tracks and $`m_\tau =1.777\mathrm{GeV}/c^2`$ is the tau mass bes .

The signal distributions in the $`(\mathrm{\Delta }M,\mathrm{\Delta }E)`$ plane are broadened by detector resolution and radiative effects. The radiation of photons from the incoming $`e^+e^{}`$ particles before annihilation affects all decay modes, leading to a tail at low values of $`\mathrm{\Delta }E`$. Radiation from the final-state leptons is more likely for electrons than muons, and produces a tail at low values of $`\mathrm{\Delta }M`$ as well. Rectangular signal regions are defined separately for each decay mode as follows. For all six decay modes, the upper right corner of the signal region is fixed at $`(30\mathrm{MeV}/c^2,50\mathrm{MeV})`$, while the lower left corner is at $`(70,120)`$ for the $`e^{}e^+e^{}`$ and $`\mu ^{}e^+e^{}`$ decay modes, $`(100,200)`$ for $`\mu ^+e^{}e^{}`$, $`(50,200)`$ for $`e^+\mu ^{}\mu ^{}`$, $`(50,150)`$ for $`e^{}\mu ^+\mu ^{}`$, and $`(30,150)`$ for $`\mu ^{}\mu ^+\mu ^{}`$. All values are given in units of $`(\mathrm{MeV}/c^2,\mathrm{MeV})`$. These signal region boundaries are chosen to provide the smallest expected upper limits on the branching fractions in the background-only hypothesis. These expected upper limits are estimated using only Monte Carlo (MC) simulations and data control samples, not candidate signal events. Figure 1 shows the observed data in the $`(\mathrm{\Delta }M,\mathrm{\Delta }E)`$ plane, along with the signal region boundaries and the expected signal distributions. To avoid bias, a blinded analysis procedure was adopted with the number of data events in the signal region remaining unknown until the selection criteria were finalized and all cross checks were performed.

The efficiency of the selection for signal events is estimated with a MC simulation of LFV tau decays. Simulated tau-pair events including higher-order radiative corrections are generated using KK2f kk  with one tau decaying to three leptons with a 3-body phase space distribution, while the other tau decays according to measured rates PDG  simulated with Tauola tauola . Final state radiative effects are simulated for all decays using Photos photos . The detector response is simulated with GEANT4 geant , and the simulated events are then reconstructed in the same manner as data.

About 50% of the MC signal events pass the 1-3 topology requirement. The lepton identification efficiencies and misidentification probabilities are measured using tracks in kinematically-selected data samples (radiative Bhabha, radiative $`\mu ^+\mu ^{}`$, two-photon $`e^+e^{}\mathrm{}^+\mathrm{}^{}`$, and $`J/\psi \mathrm{}^+\mathrm{}^{}`$) and parameterized as a function of particle momentum, polar angle, and azimuthal angle in the laboratory frame. These data-derived efficiencies are then used to give the probability that a simulated MC particle will be identified (or misidentified) as an electron or a muon. For the lepton momentum spectrum predicted by the signal MC, the electron and muon identification requirements are found to have an average efficiency per lepton of 91% and 63%, respectively. The probability for a hadron to be misidentified as an electron in SM 3-prong tau decays is 2.2%, while the probability to be misidentified as a muon is 4.8% pid . The final efficiency for signal events to be found in the signal region is shown in Table 1 for each decay mode and ranges from 7% to 12%. This efficiency includes the 85% branching fraction for 1-prong tau decays.

There are three main classes of background remaining after the selection criteria are applied: low multiplicity $`q\overline{q}`$ events (mainly continuum light-quark production), QED events (Bhabha and $`\mu ^+\mu ^{}`$), and SM $`\tau ^+\tau ^{}`$ events. These three background classes have distinctive distributions in the $`(\mathrm{\Delta }M,\mathrm{\Delta }E)`$ plane: $`q\overline{q}`$ events tend to populate the plane uniformly, while QED backgrounds are restricted to a narrow band at positive values of $`\mathrm{\Delta }E`$, and $`\tau ^+\tau ^{}`$ backgrounds are restricted to negative values of both $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }M`$. A negligible two-photon background remains.

The expected background rates for each decay mode are determined by fitting a set of probability density functions (PDFs) to the observed data in the $`(\mathrm{\Delta }M,\mathrm{\Delta }E)`$ plane in a grand sideband (GS) region. The GS region, shown in Fig. 1, is defined as the rectangle bounded by the points $`(600\mathrm{MeV}/c^2,700\mathrm{MeV})`$ and $`(400\mathrm{MeV}/c^2,400\mathrm{MeV})`$, excluding the signal region. For both the $`q\overline{q}`$ and $`\tau ^+\tau ^{}`$ backgrounds, an analytic PDF is constructed from the product of two PDFs $`P_M`$ and $`P_E`$, where $`P_M(\mathrm{\Delta }M)`$ is the sum of two Gaussians with a common mean and $`P_E(\mathrm{\Delta }E)=(1x/\sqrt{1+x^2})(1+ax+bx^2+cx^3)`$ with $`x=(\mathrm{\Delta }Ed)/e`$ opal00 . The shapes of these PDFs are described by a total of nine free parameters, which are determined by fits to MC $`q\overline{q}`$ and $`\tau ^+\tau ^{}`$ background samples for each decay mode.

For the QED backgrounds, an analytic PDF is constructed from the product of a Crystal Ball function CBF  in $`\mathrm{\Delta }E^{}`$ and a linear function in $`\mathrm{\Delta }M^{}`$, where the $`(\mathrm{\Delta }M^{},\mathrm{\Delta }E^{})`$ axes have been rotated slightly from $`(\mathrm{\Delta }M,\mathrm{\Delta }E)`$ to fit the observed distribution. The six parameters of this PDF, including the rotation angle, are obtained by fitting control samples with a 1-3 topology that are enhanced in Bhabha or $`\mu ^+\mu ^{}`$ events by requiring that the particle in the 1-prong hemisphere is identified as an electron or muon. Any value for $`\mathrm{cos}\theta _{13}`$ is allowed, but the control sample events otherwise pass the selection criteria.

With the shapes of the three background PDFs determined, an unbinned maximum likelihood fit to the data in the GS region is used to find the expected rate of each background type in the signal region, as shown in Table 1. The PDF shape determinations and background fits are performed separately for each of the six decay modes. Figure 2 shows the data and the background PDFs for values of $`\mathrm{\Delta }E`$ in the signal range.

The largest systematic uncertainty in the signal efficiency is due to the uncertainty in measuring the PID efficiencies. This uncertainty is determined from the statistical precision of the PID control samples, and ranges from $`0.7\%`$ for $`e^{}e^+e^{}`$ to $`6.2\%`$ for $`\mu ^{}\mu ^+\mu ^{}`$relative to the efficiency uncertain . The modeling of the tracking efficiency contributes an additional 2% uncertainty, as does the statistical limitation of the MC signal sample. All other sources of uncertainty are found to be small, including the modeling in the generator of radiative effects, track momentum resolution, trigger performance, observables used in the selection criteria, and knowledge of the tau 1-prong branching fractions. The efficiency has been estimated using a 3-body phase space model and no uncertainty is assigned for possible model dependence. The selection efficiency is found to be uniform within 10% across the Dalitz plane, provided the invariant mass for any pair of leptons is less than $`1.4\mathrm{GeV}/c^2`$.

Since the background levels are extracted directly from the data, systematic uncertainties on the background estimation are directly related to the background parameterization and the fit technique used. The finite data available in the GS region to determine the background rates is the largest uncertainty and varies from 10% to 25% depending upon the decay mode. Additional uncertainties are estimated by varying the fit procedure and changing the functional form of the background PDFs. Cross checks of the background estimation were performed by considering the number of events expected and observed in sideband regions immediately neighboring the signal region for each decay mode.

The numbers of events observed ($`N_{\mathrm{obs}}`$) and the background expectations ($`N_{\mathrm{bgd}}`$) are shown in Table 1, with no significant excess found in any decay mode. Upper limits on the branching fractions are calculated according to $`_{\mathrm{UL}}^{90}=N_{\mathrm{UL}}^{90}/(2\epsilon \sigma _{\tau \tau })`$, where $`N_{\mathrm{UL}}^{90}`$ is the 90% CL upper limit for the number of signal events when $`N_{\mathrm{obs}}`$ events are observed with $`N_{\mathrm{bgd}}`$ background events expected. The values $`\epsilon `$, $``$, and $`\sigma _{\tau \tau }`$ are the selection efficiency, luminosity, and $`\tau ^+\tau ^{}`$ cross section, respectively. The estimates of $`=91.5\text{ fb}^1`$ and $`\sigma _{\tau \tau }=0.89`$ nb are correlated lum , and the uncertainty on the product $`\sigma _{\tau \tau }`$ is 2.3%. The branching fraction upper limits have been calculated including all uncertainties using the technique of Cousins and Highland cousins92  following the implementation of Barlow barlow02 . The 90% CL upper limits on the $`\tau ^{}\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$ branching fractions, shown in Table 1, are in the range $`(13)\times 10^7`$.

We are grateful for the excellent luminosity and machine conditions provided by our PEP-II colleagues, and for the substantial dedicated effort from the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and kind hospitality. This work is supported by DOE and NSF (USA), NSERC (Canada), IHEP (China), CEA and CNRS-IN2P3 (France), BMBF and DFG (Germany), INFN (Italy), FOM (The Netherlands), NFR (Norway), MIST (Russia), and PPARC (United Kingdom). Individuals have received support from the A. P. Sloan Foundation, Research Corporation, and Alexander von Humboldt Foundation.