format: text

text:
# Detecting Microscopic Black Holes with Neutrino Telescopes

## I Introduction

The possibility that we live in $`D=4+n>4`$ spacetime dimensions has profound implications. In particular, if gravity propagates in these extra dimensions, the fundamental Planck scale $`M_D`$ at which gravity becomes comparable in strength to other forces may be far below $`M_410^{19}\text{GeV}`$, leading to a host of potential signatures for high energy physics Antoniadis:1998ig .

Among the most striking consequences of low-scale gravity is the possibility of black hole creation in high-energy particle collisions Amati:1987wq ; 'tHooft:rb ; Argyres:1998qn ; Banks:1999gd ; Emparan:2000rs ; Giddings:2000ay ; Giddings:2001bu ; Dimopoulos:2001hw . In most gravitational processes, such as those involving graviton emission and exchange, analyses rely on a perturbative description that breaks down for center-of-mass energies of $`M_D`$ and above. In contrast, black hole properties are best understood for energies above $`M_D`$, where semi-classical and thermodynamic descriptions become increasingly valid Hawking:1975sw . In principle, then, black holes provide a robust probe of extra dimensions and low-scale gravity, as long as particle collisions with center-of-mass energies above $`M_D1\text{TeV}`$ are available.

Nature provides interactions at the necessary energies in the form of cosmic rays with energies above $`10^{10}\text{GeV}`$. In collisions with nucleons, these cosmic rays probe center-of-mass energies exceeding 100 TeV, beyond both current man-made colliders and those of the foreseeable future. Cosmic neutrinos may create black holes deep in the Earth’s atmosphere, resulting in spectacular signals of giant air showers in ground arrays and fluorescence detectors Feng:2001ib ; Anchordoqui:2001ei ; Emparan:2001kf ; Ringwald:2001vk ; Anchordoqui:2001cg . With a handful of events, standard model (SM) and most alternative explanations may be excluded Feng:2001ib ; Anchordoqui:2001cg  by comparison with rates for Earth-skimming neutrinos Bertou:2001vm ; Feng:2001ue ; Kusenko:2001gj , and with more events, black holes may be identified through their shower characteristics Anchordoqui:2001ei . Bounds have been derived Ringwald:2001vk ; Anchordoqui:2001cg  from the absence of such showers in current data from Fly’s Eye Baltrusaitis:mt  and AGASA agasa . For conservative fluxes and the geometric black hole cross section, the AGASA data require $`M_D1.31.8\text{TeV}`$ for $`n4`$, the most stringent constraint to date Anchordoqui:2001cg . At the same time, the Auger Observatory, scheduled for completion by 2004, may observe tens of black hole events per year Feng:2001ib ; Ringwald:2001vk ; Anchordoqui:2001cg . Related phenomena related to $`p`$-brane production may also be observed in cosmic rays Ahn:2002mj ; Jain:2002kf .

Here we examine the possibility of detecting and studying black holes produced by cosmic neutrinos in neutrino telescopes. Several large-scale neutrino telescope projects are underway, including IceCube IceCube  in the Antarctic ice, and ANTARES ANTARES  and NESTOR NESTOR  in the Mediterranean. Among the many possible black hole signatures, such detectors are most sensitive to contained hadronic showers and through-going muons and taus from the evaporation of black holes produced in the Earth’s crust. The possibility of black hole detection in neutrino telescopes has recently been studied in Ref. Kowalski:2002gb ; where possible, we will compare our results to those of Ref. Kowalski:2002gb  below. For a preliminary study, see also Ref. Uehara:2001yk .

For TeV-scale gravity and conservative flux assumptions, we find that IceCube could detect several black holes per year. These rates may be enhanced by larger fluxes, and observable rates are possible even given a postulated exponential suppression factor in the black hole cross section Voloshin:2001vs ; Voloshin:2001fe . The relative event rates in the three channels may differ from the SM, and the energy and angle distributions of black hole events are also distinctive. These will not only help identify black holes, but may also constrain parameters such as $`n`$ and $`M_D`$, and determine if suppression factors in the cross section are present or absent. The search for black holes in neutrino telescopes therefore complements black hole searches in other cosmic ray detectors, as well as searches for the effects of perturbative gravity processes at center-of-mass energies below $`M_D`$ Nussinov:1999jt ; Jain:2000pu ; Tyler:2001gt ; Alvarez-Muniz:2001mk .

## II Black Hole Production and Evaporation

In $`D=4+n`$ dimensions, gravity is described by the Einstein action

$$S_E=\frac{1}{8\pi G_D}d^{4+n}x\sqrt{g}\frac{1}{2},$$
(1)
where $`G_D`$ is the $`D`$-dimensional Newton’s constant. To define the fundamental Planck scale $`M_D`$, we adopt the convention

$$\frac{1}{8\pi G_D}=\frac{M_D^{2+n}}{(2\pi )^n}.$$
(2)
In the most straightforward scenarios Antoniadis:1998ig  with flat extra dimensions of equal length, TeV-scale gravity is excluded for $`n=1`$ by solar system tests of Newtonian gravity. Astrophysical bounds Cullen:1999hc  on supernova cooling and neutron star heating provide the most stringent bounds for $`n=2`$ ($`M_D600\text{TeV}`$) and $`n=3`$ ($`M_D10\text{TeV}`$). For $`n4`$, the most stringent bounds are from collider searches for perturbative graviton effects Giudice:1998ck ; Pagliarone:2001ff ; Abbott:2000zb  and cosmic ray bounds on black hole production Ringwald:2001vk ; Anchordoqui:2001cg ; Baltrusaitis:mt ; agasa . These constraints are each subject to a variety of theoretical assumptions (for a comparison and discussion, see, e.g., Ref. Anchordoqui:2001cg ), but the most stringent of these require roughly $`M_D1\text{TeV}`$.

To determine event rates for neutrino telescopes, we must first model black hole production and evaporation. For production by cosmic neutrinos, we follow the analysis of Ref. Feng:2001ib . Black holes produced in parton collisions are typically far smaller than the length scales of the extra dimensions. These black holes are then well-approximated by (4+$`n`$)-dimensional solutions. The Schwarzschild radius for a (4+$`n`$)-dimensional black hole with mass $`M_{\text{BH}}`$ and vanishing charge and angular momentum is Myers:1986un

$$r_s(M_{\text{BH}}^2)=\frac{1}{M_D}\left[\frac{M_{\text{BH}}}{M_D}\right]^{\frac{1}{1+n}}\left[\frac{2^n\pi ^{\frac{n3}{2}}\mathrm{\Gamma }\left(\frac{3+n}{2}\right)}{2+n}\right]^{\frac{1}{1+n}}.$$
(3)
We assume that two partons $`i`$ and $`j`$ with center-of-mass energy $`\sqrt{\widehat{s}}`$ form a black hole of mass $`M_{\text{BH}}=\sqrt{\widehat{s}}`$ when they pass within a distance $`r_s(\widehat{s})`$, leading to a geometric cross section of Banks:1999gd ; Giddings:2001bu ; Dimopoulos:2001hw

$$\widehat{\sigma }(ij\text{BH})(\widehat{s})=\pi r_s^2(\widehat{s}).$$
(4)
Evidence from analyses of axisymmetric D'Eath:hb  and off-axis Eardley:2002re  classical collisions, an analysis in a simple model framework Solodukhin:2002ui , and a string calculation Dimopoulos:2001qe  suggests that this picture is valid semi-classically and is not subject to large corrections Giddings:2001ih . Modifications for non-vanishing angular momentum and spinning black holes have also been found to be small Anchordoqui:2001cg . However, Voloshin has argued Voloshin:2001vs ; Voloshin:2001fe  that the cross section could be suppressed by the factor $`e^I`$, where the action is

$$I=\frac{S}{n+1}=\frac{4\pi M_{\text{BH}}r_s}{(n+1)(n+2)},$$
(5)
with $`S`$ the black hole entropy. This implies vanishing cross sections in the classical limit, contrary to the evidence noted above. To explore the impact of modifications to the black hole production cross section, however, we consider cases both with and without this suppression.

The neutrino-nucleon scattering cross section is then Feng:2001ib

$$\sigma (\nu N\text{BH})=\underset{i}{}_{(x_{\text{min}}M_D)^2/s}^1𝑑x\widehat{\sigma }_i(xs)f_i(x,Q),$$
(6)
where

$$x_{\text{min}}M_{\text{BH}}^{\text{min}}/M_D1$$
(7)
parameterizes the minimal black hole mass for which we expect a semi-classical description to be valid, $`s=2m_NE_\nu `$, the sum is over all partons $`i`$, and $`f_i`$ are the parton distribution functions. The cross section is highly insensitive to choice of momentum transfer Feng:2001ib ; e.g., the choices $`Q=M_{\text{BH}}`$ Dimopoulos:2001hw  and $`Q=r_s^1`$ Emparan:2001kf , lead to cross section differences of only 10% to 20%. These cross sections are typically also insensitive to uncertainties at low $`x`$. For example, even extremely high energy neutrinos with $`E_\nu 10^{10}\text{GeV}`$ probe only $`x(1\text{TeV})^2/10^{10}\text{GeV}^210^4`$, within the range of validity of the CTEQ5 parton distribution functions we use Lai:2000wy .

Cross sections for black hole production by cosmic neutrinos are given in Fig. 1. The SM cross section for $`\nu N\mathrm{}X`$ is included for comparison (dotted curve). As a result of the sum over all partons and the lack of suppression from small perturbative couplings, the black hole cross section may exceed SM interaction rates by two or more orders of magnitude.

Note that in our conventions, the cross section rises for increasing $`n`$ and fixed $`M_D`$. In conventions where the fundamental Planck scale is taken to be $`M_{}`$ with

$$\frac{1}{G_D}=M_{}^{2+n},$$
(8)
this behavior is reversed: the cross section decreases for increasing $`n`$ and fixed $`M_{}`$ Feng:2001ib . The dependence of the cross section on $`n`$ for fixed Planck scale is convention-dependent and unphysical; we have adopted the $`M_D`$ convention to simplify comparison with existing collider bounds.

Assuming a constant density of $`1.0\text{g}/\text{cm}^3`$ for the Earth’s surface near the detector, as is valid for IceCube, the neutrino’s interaction length in Earth is

$$L=1.7\times 10^7\text{km}\left(\frac{\text{pb}}{\sigma }\right).$$
(9)
The center of IceCube is at a depth of roughly 1.8 km. A neutrino reaching this point horizontally passes through 150 km of Earth. The cross sections corresponding to neutrino interaction lengths equal to these two lengths, that is, the horizontal and vertical depths of IceCube, are also given in Fig. 1. We see that for the geometric cross section, $`M_D1\text{TeV}`$, and neutrino energies $`E_\nu 10^9\text{GeV}`$ where the cosmogenic flux peaks (see Sec. III), black hole production increases the probability of conversion in down-going neutrinos without increasing the cross section so much that vertical neutrinos are shadowed by the Earth. We therefore expect significantly enhanced rates in neutrino telescopes. (On the other hand, the upgoing event rates will be even more suppressed.) For the exponentially suppressed cross section, a similar enhancement is also possible for $`M_D=1`$ to $`2\text{TeV}`$.

Once produced, these tiny black holes evaporate with rest lifetime of order $`\text{TeV}^110^{27}\text{s}`$. Even though highly boosted, they decay before accreting matter. They evaporate in a thermal distribution with temperature $`T_H=(1+n)/(4\pi r_s)`$ Hawking:1975sw ; Myers:1986un  and average multiplicity Giddings:2001bu ; Dimopoulos:2001hw

$$N\frac{M_{\text{BH}}}{2T_H}=\frac{2\pi }{1+n}\left[\frac{M_{\text{BH}}}{M_D}\right]^{\frac{2+n}{1+n}}\left[\frac{2^n\pi ^{\frac{n3}{2}}\mathrm{\Gamma }\left(\frac{3+n}{2}\right)}{2+n}\right]^{\frac{1}{1+n}}.$$
(10)
Neglecting particle masses, the decay products are distributed according to the number of degrees of freedom Emparan:2000rs : quarks (72), gluons (16), charged leptons (12), neutrinos (6), $`W`$ and $`Z`$ bosons (9), photons (2), Higgs bosons (1), and gravitons (2). We neglect the possibility of other low mass degrees of freedom, such as right-handed neutrinos and supersymmetric particles. About 75% of the black hole’s energy is radiated in hadronic degrees of freedom, while the probability of any given decay particle being a muon (or a tau) is approximately 3%.

As through-going muons and taus will be a promising signal in neutrino telescopes, the typical multiplicity of black hole decays is of great importance. To quantify this, we define a weighted multiplicity

$$\overline{N}=\frac{1}{\sigma }\underset{i}{}_{(x_{\text{min}}M_D)^2/s}^1𝑑xN\widehat{\sigma }_i(xs)f_i(x,Q),$$
(11)
where $`\sigma `$ is given in Eq. (6), and $`N`$ is as in Eq. (10) with $`M_{\text{BH}}^2=xs`$. The weighted multiplicity $`\overline{N}`$ is given in Fig. 2 for various cases with $`n=6`$ and $`M_D=1\text{TeV}`$. For the geometric black hole cross section, these multiplicities may be substantially enhanced for lower $`n`$; for the exponentially suppressed cross section, the dependence on $`n`$ is slight. We find that $`\overline{N}𝒪(10)`$ is possible for ultra-high energy neutrinos. Note that while raising $`x_{\text{min}}`$ and including the exponential suppression factor both suppress the total cross section, they have opposite effects on $`\overline{N}`$: raising $`x_{\text{min}}`$ eliminates events with relatively low multiplicity, and so raises $`\overline{N}`$, while the exponential suppression is largest for events with large $`M_{\text{BH}}`$ and large multiplicity, and so suppresses $`\overline{N}`$.

## III Astrophysical Sources of Ultra-High Energy Neutrinos

Observations of black hole formation and decay are most easily made with neutrino fluxes that are large near or above the EeV scale. Several such sources have been proposed in the literature. While a detailed discussion of these sources is beyond the scope of this paper, we briefly describe here a few of them, including those used in our study.

First, neutrinos are almost certainly produced through pion decays in the scattering of protons off the cosmic microwave background, $`p\gamma _{\mathrm{CMB}}n\pi ^+`$ Stecker:1979ah . This cosmogenic flux is subject to a number of quantitative uncertainties, including cosmological source evolution. As representative fluxes, we consider the recent results presented in Fig. 4 of Ref. Engel:2001hd . These fluxes are shown in Fig. 3.

Second, gamma ray bursts have also been considered as a possible source of the highest energy cosmic rays. If this is the case, Fermi accelerated protons from shocks will generate extremely high energy neutrinos with energy spectrum $`d\mathrm{\Phi }_\nu /dE_\nu E_\nu ^2`$ Waxman:1997ti ; Waxman:2001tk ; Halzen:1999xc ; Alvarez-Muniz:2000st . This neutrino flux, as well as those from other compact sources, such as active galactic nuclei, is limited by the Waxman-Bahcall (WB) bound Waxman:1998yy . This constraint is valid for all astrophysical neutrino sources that are optically thin to $`p\gamma `$ and $`pp`$ interactions. We consider a conservative estimate of this bound, $`E_\nu ^2d\mathrm{\Phi }_\nu /dE_\nu =1\times 10^8\text{cm}^2\text{s}^1\text{sr}^1`$ Waxman:1998yy , where here $`\nu =\nu _e,\nu _\mu ,\overline{\nu }_\mu `$. This flux is also shown in Fig. 3. It is approximately equal to the cosmogenic flux for $`E_\nu 10^{10}\text{GeV}`$, but is much larger for lower (and higher) energies.

Third, if the highest energy cosmic rays observed are generated by the annihilation of superheavy dark matter particles or from the decay of topological defects, neutrinos will also be produced. Such fluxes have been described in Refs. Halzen:2001ec ; Sarkar:2001se ; Birkel:1998nx ; Berezinsky:1997hy . We will not discuss these sources further, but note that they may predict large neutrino fluxes at extremely high energies, enhancing the results given below.

In propagating to the Earth, the neutrino fluxes of Fig. 3 will mix Athar:2000yw . Given the solutions preferred by neutrino oscillation experiments and the enormous distances traveled, we take the neutrinos that reach the Earth to be in the ratio $`\nu _e:\nu _\mu :\nu _\tau =1:1:1`$ and similarly for anti-neutrinos.

## IV Event Simulation

With regard to signals at neutrino telescopes, black hole evaporation products may be divided into three categories: showers (hadronic and electromagnetic), muons, and taus. Black holes also decay to neutrinos, but we neglect this flux in this work.

The ranges of typical hadronic and electromagnetic showers are much less than the linear dimension of large-scale neutrino telescopes, and so, to first approximation, only contained showers from black holes produced inside the detector may be detected. The backgrounds for showers from black hole evaporation consist of hadronic showers from neutral and charged current neutrino events and electromagnetic showers from charged current electron neutrino events. At IceCube, hadronic shower energies should be measurable with an accuracy of about 30%.

In contrast to showers, at the typical energies of black hole events, muons travel several kilometers before losing a decade of energy. The dominant signal is therefore through-going muons. IceCube can measure the energy and direction of any observed muon. The angular resolution is about $`2^{}3^{}`$ while the energy resolution is approximately a factor of three. Signal and background muons may therefore be differentiated with an energy cutoff. As we will see, for down-going muons with energy above 500 TeV to 1 PeV, the black hole signal may be well above the SM background from atmospheric neutrinos. Note that for black holes produced sufficiently near the detector, muon events may be obscured by showers. This occurrence is rare, however, and we ignore this possibility below.

At high energies, when tau decay is sufficiently time dilated, taus have ranges as large or larger than muons, and so the dominant tau signal is from through-going taus. These events have a characteristic signature consisting of a “clean” track, i.e., a track without much energy lost through low energy cascades. We will consider all events in which a tau track passes through the detector. Taus can be differentiated from muons by the cleanliness of their tracks or through their decays in the detector leading to “lollipop” events IceCubetaus , described later in this paper. We assume that taus and muons are distinguishable at all energies, although it will be difficult to distinguish a slow muon of energy $`200\text{GeV}`$ from a very high energy tau, as they can both generate clean, through-going tracks in the detector.

Very massive black holes with large multiplicity decays may evaporate to several muons (or taus). However, these travel in coincidence through the detector. Their angular spread is only $`\mathrm{\Delta }\theta T_HN/E_\nu `$, much less than a degree, and so cannot be resolved with large Cerenkov detectors such as IceCube. Spectacular multi-lepton signatures are therefore not possible. Note, however, that more massive black holes are easier to detect, as they produce leptons with greater energy, which travel further before dropping below the cutoff energy.

To evaluate black hole detection prospects, it is essential to determine, in a unified framework, both the SM background and the black hole event rate. We now describe both of these calculations.

### IV.1 Standard Model Events

#### IV.1.1 Showers

In the SM, a general expression for the total number of shower events in an underground detector is

$`N_{\mathrm{sh}}`$ $`=`$ $`{\displaystyle \underset{i,j}{}}2\pi AT{\displaystyle d\mathrm{cos}\theta _z𝑑E_{\nu _i}\frac{d\mathrm{\Phi }_{\nu _i}}{dE_{\nu _i}}(E_{\nu _i})P_{\mathrm{surv}}}`$ (12)
$`\times {\displaystyle _{y_{\mathrm{min}}^{i,j}}^{y_{\mathrm{max}}^{i,j}}}dy{\displaystyle \frac{1}{\sigma _{\mathrm{SM}}^j(E_{\nu _i})}}{\displaystyle \frac{d\sigma _{\mathrm{SM}}^j}{dy}}(E_{\nu _i})P_{\mathrm{int}},`$
where $`\theta _z`$ is the zenith angle ($`\theta _z=0`$ is vertically downward), and the sums are over neutrino (and anti-neutrino) flavors $`i=e,\mu ,\tau `$ and interactions $`j=\mathrm{CC}`$ (charged current) and NC (neutral current). $`A`$ is the detector’s cross sectional area with respect to the $`\nu `$ flux, $`T`$ is its observation time, and $`d\mathrm{\Phi }_{\nu _i}/dE_{\nu _i}`$ is the differential neutrino flux that reaches the Earth. For $`i=\tau `$, Eq. (12) is modified to include the effects of regeneration, as discussed below.

$`P_{\mathrm{surv}}`$ is the probability that a neutrino survives to reach the detector. It is given by

$$P_{\mathrm{surv}}\mathrm{exp}[X(\theta _z)\sigma _{\mathrm{SM}}^{\mathrm{tot}}(E_{\nu _i})N_A],$$
(13)
where $`N_A6.022\times 10^{23}\text{g}^1`$, and the total neutrino interaction cross section is

$$\sigma _{\mathrm{SM}}^{\mathrm{tot}}=\sigma ^{\mathrm{CC}}+\sigma ^{\mathrm{NC}}.$$
(14)
Note that this is conservative, as it neglects the possibility of a neutrino interacting through a NC interaction and continuing on to create a contained shower. $`X(\theta _z)`$ is the column density of material the neutrino must traverse to reach the detector with zenith angle $`\theta _z`$. It depends on the depth of the detector and is given by

$$X(\theta _z)=_{\theta _z}\rho (r(\theta _z,l))𝑑l,$$
(15)
the path length along direction $`\theta _z`$ weighted by the Earth’s density $`\rho `$ at distance $`r`$ from the Earth’s center. For the Earth’s density profile, we adopt the piecewise continuous density function $`\rho (r)`$ of the Preliminary Earth Model preliminary .

$`P_{\mathrm{int}}`$ is the probability that the neutrino interacts in the detector. It is given by

$$P_{\mathrm{int}}=1\mathrm{exp}\left[\frac{L}{L_{\mathrm{SM}}^j(E_{\nu _i})}\right],$$
(16)
where, for showers, $`L`$ is the linear dimension of the detector, and $`L_{\mathrm{SM}}^j(E_{\nu _i})`$ is the mean free path for neutrino interaction of type $`j`$. For realistic detectors, $`LL_{\mathrm{SM}}^j(E_{\nu _i})`$, and so $`P_{\mathrm{int}}L/L_{\mathrm{SM}}^j(E_{\nu _i})`$. To an excellent approximation, then, the shower event rate scales linearly with detector volume $`V=AL`$, and we present results in units of events/volume/time.

Finally, the inelasticity parameter $`y`$ is the fraction of the neutrino energy carried away in hadrons. The limits of integration are determined by the interaction type and neutrino flavor. For NC $`\nu _e`$ interactions and all $`\nu _\mu `$ and $`\nu _\tau `$ interactions, $`y_{\mathrm{max}}=1`$ and $`y_{\mathrm{min}}=E_{\mathrm{sh}}^{\mathrm{thr}}/E_\nu `$, where $`E_{\mathrm{sh}}^{\mathrm{thr}}`$ is the threshold energy for shower detection. For CC $`\nu _e`$ interactions, the outgoing electron also showers, and so $`y_{\mathrm{max}}=1`$ and $`y_{\mathrm{min}}=0`$.

#### IV.1.2 Muons

Energetic through-going muons are produced only by $`\nu _\mu `$ CC interactions. For a muon to be detected, it must reach the detector with energy above some threshold $`E_\mu ^{\mathrm{thr}}`$. The expression of Eq. (12) then also describes the number of muon events with $`P_{\mathrm{surv}}`$ as before, but with

$$P_{\mathrm{int}}=1\mathrm{exp}\left[\frac{R_\mu }{L_{\mathrm{SM}}^{\mathrm{CC}}(E_{\nu _\mu })}\right],$$
(17)
where $`R_\mu `$ is the range of a muon with initial energy $`E_\mu =(1y)E_{\nu _\mu }`$ and final energy $`E_\mu ^{\mathrm{thr}}`$. We assume muons lose energy continuously according to

$$\frac{dE}{dX}=\alpha \beta E,$$
(18)
where $`\alpha =2.0\text{MeV}\text{cm}^2/\text{g}`$ and $`\beta =4.2\times 10^6\text{cm}^2/\text{g}`$ Dutta:2000hh . The muon range is then

$$R_\mu =\frac{1}{\beta }\mathrm{ln}\left[\frac{\alpha +\beta E_\mu }{\alpha +\beta E_\mu ^{\mathrm{thr}}}\right].$$
(19)
In this case, $`y_{\mathrm{max}}=1E_\mu ^{\mathrm{thr}}/E_\nu `$ and $`y_{\mathrm{min}}=0`$.

The event rate for muons is significantly enhanced by the possibility of muons propagating from several km into the detector. Note, however, that this enhancement is $`\theta _z`$-dependent: for nearly vertical down-going paths, the path length of the muon is limited by the amount of matter above the detector, not by the muon’s range. This is taken into account explicitly in the simulations, and its effect will be evident in the results presented in Sec. VI.

At extremely high neutrino energies, the approximate form $`P_{\mathrm{int}}R_\mu /L_{\mathrm{SM}}^{\mathrm{CC}}(E_{\nu _\mu })`$, often presented in the literature, is less accurate than the one used here, and the difference may be significant in scenarios in which the neutrino cross section is enhanced with respect to the SM value, such as the ones we explore in this paper. Note also that in this case, our expression for $`P_{\mathrm{surv}}`$ is conservative, as it demands that the neutrino survive all the way to the detector, neglecting the possibility that it may convert a significant distance from the detector and still produce a signal.

#### IV.1.3 Taus

Taus are produced only by CC $`\nu _\tau `$ interactions. This process differs significantly from the muon case, as tau neutrinos are regenerated by tau decay through $`\nu _\tau \tau \nu _\tau `$ Halzen:1998be . As a result, for tau neutrinos, CC and NC interactions do not deplete the $`\nu _\tau `$ flux, but serve only to soften the neutrino energy. We include this important effect by first performing a dedicated simulation that determines $`\overline{E}_{\nu _\tau }(E_{\nu _\tau },\theta _z)`$, the average energy a $`\nu _\tau `$ has when it reaches the detector, as a function of its initial energy $`E_{\nu _\tau }`$ and zenith angle $`\theta _z`$. The tau event rate is then given by

$`N_\tau `$ $`=`$ $`2\pi AT{\displaystyle d\mathrm{cos}\theta _z𝑑E_{\nu _\tau }\frac{d\mathrm{\Phi }_{\nu _\tau }}{dE_{\nu _\tau }}(E_{\nu _\tau })_{y_{\mathrm{min}}}^{y_{\mathrm{max}}}𝑑y\frac{1}{\sigma _{\mathrm{SM}}^{\mathrm{CC}}(\overline{E}_{\nu _\tau })}\frac{d\sigma _{\mathrm{SM}}^{\mathrm{CC}}}{dy}(\overline{E}_{\nu _\tau })}`$ (20)
$`\times \left[1\mathrm{exp}\left({\displaystyle \frac{R_\tau ((1y)\overline{E}_{\nu _\tau })}{L_{\mathrm{SM}}^{\mathrm{CC}}(\overline{E}_{\nu _\tau })}}\right)\right]\mathrm{\Theta }((1y)\overline{E}_{\nu _\tau }E_\tau ^{\mathrm{thr}}),`$
where $`R_\tau ((1y)\overline{E}_{\nu _\tau })`$ is the range of the produced tau, but evaluated at the energy of the tau neutrino after regeneration. $`R_\tau `$ is given by Eq. (19) but now with $`\beta =3.6\times 10^7\text{cm}^2/\text{g}`$ Dutta:2000hh . The last factor takes into account the requirement that the tau track be long enough to be identified in the detector. We require $`E_\tau ^{\mathrm{thr}}2.5\times 10^6\text{GeV}`$ so that the tau decay length is above 125 m, the string separation length in IceCube. It is not clear, at this time, whether through-going tau events will be separable from less energetic muon events. Those tau events that include one (lollipop events) or two (double bang events) showers in the detector volume will be identifiable, however. The rate of down-going lollipop events in a $`\text{km}^3`$ neutrino telescope is expected to be of the order of the rate of down-going shower events, probably slightly smaller. Double bang events will be mostly observed for neutrino energies in a limited range between roughly 10 and 100 PeV IceCubetaus . Our results assume that all tau events can be distinguished from muon events, although this may be difficult to realize.

As with muons, at very high energies taus can travel several kilometers before decaying or suffering significant energy loss. The enhancement to tau event rates from this effect is $`\theta _z`$-dependent as discussed above for muons.

### IV.2 Black Hole Events

Black hole event rates may be determined with only minor modifications. For showers, the corresponding expression is

$`N_{\mathrm{BH}}`$ $`=`$ $`{\displaystyle \underset{i}{}}2\pi AT{\displaystyle d\mathrm{cos}\theta _z𝑑E_{\nu _i}\frac{d\mathrm{\Phi }_{\nu _i}}{dE_{\nu _i}}(E_{\nu _i})P_{\mathrm{surv}}}`$ (21)
$`\times {\displaystyle _0^1}dy{\displaystyle \frac{1}{\sigma _{\mathrm{BH}}(E_{\nu _i})}}{\displaystyle \frac{d\sigma _{\mathrm{BH}}}{dy}}(E_{\nu _i})P_{\mathrm{int}}P_{\mathrm{BH}}(E_{\nu _i}).`$
The survival and interaction probabilities are now

$`P_{\mathrm{surv}}`$ $`=`$ $`\mathrm{exp}[X(\theta _z)\sigma ^{\mathrm{tot}}(E_{\nu _i})N_A]`$ (22)
$`P_{\mathrm{int}}`$ $`=`$ $`1\mathrm{exp}\left[{\displaystyle \frac{L}{L_{\mathrm{BH}}(E_{\nu _i})}}\right],`$ (23)
where

$$\sigma ^{\mathrm{tot}}=\sigma _{\mathrm{SM}}^{\mathrm{tot}}+\sigma _{\mathrm{BH}},$$
(24)
with $`\sigma _{\mathrm{BH}}`$ the cross section for black hole production, and $`L_{\mathrm{BH}}(E_{\nu _i})`$ is the neutrino mean free path for black hole production. As in the SM, $`L`$ is the linear dimension of the detector in the case of shower events. We assume $`d\sigma _{\mathrm{BH}}/dy\delta (y0.75)`$, i.e., that 75% of the black hole energy is carried away by showers, and impose that the generated shower has energy above $`E_{\mathrm{sh}}^{\mathrm{thr}}`$ to be detected. $`P_{\mathrm{BH}}(E_\nu )=1`$ for shower events.

In the case of muon and tau events, $`L`$ is the corresponding range for the average muon or tau energy. We assume $`E_{\mu ,\tau }=E_\nu /\overline{N}`$, where $`\overline{N}`$ is the weighted multiplicity discussed in Eq. (11). The true $`E_{\mu ,\tau }`$ distribution is essentially flat with endpoints 0 and $`2E_\nu /\overline{N}`$. Our simplification is valid except for lepton energies near threshold. In addition, spreading the lepton energy has two compensating effects, as some leptons below threshold become detectable, and some above threshold drop below threshold. We have checked that the error made is insignificant at the $`10\%`$ level.

To account for the branching fraction to leptons, we include

$$P_{\mathrm{BH}}(E_{\nu _i})=1\mathrm{exp}\left[\frac{\overline{N}(E_{\nu _i})}{30}\right],$$
(25)
the probability of obtaining at least one muon (or tau) in the decay of a black hole when the expected muon (or tau) multiplicity is $`\overline{N}/30`$.

## V Rates

The integrated event rates for showers, muons, and taus are shown in Tables 1, 2, and 3, respectively. We consider IceCube, with a representative depth of $`1.8\text{km}`$. We give results for various $`n`$, $`M_D`$, $`x_{\text{min}}`$, and with and without exponential suppression in the parton cross section. Only down-going rates are presented. Up-going rates are, as expected, extremely suppressed in the presence of black hole production, as will be seen in Sec. VI.

For showers, the geometric cross section $`\widehat{\sigma }=\pi r_s^2`$, and $`M_D=1\text{TeV}`$, we find a few events per year in each channel for cosmogenic fluxes, and as many as tens of events per year for the WB flux. These event rates are far above SM background. For the WB flux, bounds from AGASA and Fly’s Eye will imply limits above $`M_D=1\text{TeV}`$, but even for $`M_D=2\text{TeV}`$, we find reasonable rates. In this case, however, the SM background will be important. For the exponentially suppressed parton cross section, event rates are suppressed, but not drastically so. In fact, for $`x_{\text{min}}=1`$, $`𝒪(1)`$ event per year is possible given the cosmogenic flux. For $`x_{\text{min}}=3`$, the exponential suppression is large, and event rates are highly suppressed.

Generally, event rates are on more solid footing for larger $`x_{\text{min}}`$, where both the semi-classical production cross section and the assumption of thermal decay are more reliable. For showers, however, our event rate calculation relies essentially only on the requirement that black holes decay to visible showers and is insensitive to the details of evaporation. At the same time, while the production cross section for black holes (or their stringy Planck mass progenitors Dimopoulos:2001qe ) is subject to significant quantum corrections for $`M_{\text{BH}}M_D`$, there is no reason to expect it to vanish or be greatly suppressed. For showers, then, we find the requirement $`x_{\text{min}}=1`$ reasonable.

Qualitatively similar conclusions apply for muons and taus. While these event rates are suppressed relative to shower rates by the branching ratio for black hole decay to leptons, they are enhanced by the possibility of muons and taus propagating from many km away into the detector. We find that these effects effectively balance each other. In the case of muons and taus, we have assumed thermal decay distributions, which may not be accurate for $`x_{\text{min}}1`$. The event rates are hardly reduced for $`x_{\text{min}}=3`$, however; as noted in Sec. II, the more stringent $`x_{\text{min}}`$ requirement preferentially eliminates events with low multiplicities, and so has little impact on lepton rates.

The muon flux is significantly larger for the WB flux than for the cosmogenic flux. This is straightforward to understand — the WB flux is larger at low energies. The tau event rate is not as greatly enhanced, however. The tau decay length is $`4.9\text{km}(E_\tau /10^8\text{GeV})`$. While the WB flux enhances fluxes at low energies, the resulting taus decay before they lose energy, and the tau range is therefore diminished. Larger low energy fluxes therefore do not enhance the tau rate as significantly.

In addition to enhancing the SM event rates by more than an order of magnitude, Tables 1-3 show that black hole production may also change the relative event rates of the various channels. The ability of neutrino telescopes, unique among cosmic ray experiments, to differentiate showers, muons, and possibly taus, allows one to measure these relative event rates. Given sufficient statistics, this may provide an important signal of physics beyond the SM.

## VI Angle and Energy Distributions

Neutrino telescopes may also measure angle and energy distributions, providing additional opportunities to distinguish black hole events from SM or other possible physics and for constraining black hole properties.

The zenith angle distributions of SM and black hole events for showers, muons, and taus are given in Figs. 4, 5, and 6, respectively. Black hole production makes the Earth even more opaque to ultra-high energy neutrinos, reducing the up-going rate, while increasing the probability of interaction, and thereby the event rate, of down-going neutrinos. In all cases, up-going events rates are below SM contributions. No sensitivity to black hole production is then expected when looking for up-going events.

The characteristic features of the zenith angle distribution of showers compared to the distributions of muons and taus are also noteworthy. For showers, vertical fluxes are less attenuated than horizontal ones, and so the angular distribution is maximized for vertically down-going events. In contrast, as noted in Sec. IV.1, vertically down-going muon and tau events cannot benefit from muon and tau ranges beyond the depth of the detector ($`2`$ km). The optimal direction for muon and tau events is therefore closer to the horizon. In scenarios with enhanced cross sections such as the one we are exploring in this paper, the neutrino flux is attenuated even in the horizontal direction as can be seen in Figs. 5 and 6. As a consequence, the quasi-horizontal direction with small but positive $`\mathrm{cos}\theta _z`$, where the effects of attenuation are maximally offset by lepton range, leads to the largest rates. For muons, this optimal zenith angle is around $`\mathrm{cos}\theta _z0.2`$. For taus, with even longer ranges, it is $`\mathrm{cos}\theta _z0.4`$. These are general considerations that apply to the detection of ultra-high energy neutrinos in any scenario predicting enhanced cross sections with respect to the SM cross sections.

Energy distributions for the various channels are given in Figs. 7, 8, and 9. The energy distributions are clearly sensitive to $`x_{\text{min}}`$ and the presence or absence of exponential suppression in the parton cross section. Different types of events therefore probe different and complementary aspects of black hole production and decay. By measuring the down-going energy distributions with reasonable statistics, the IceCube detector may be able to discriminate between the different possibilities. Note that to facilitate comparison with conventional presentations of SM rates at neutrino telescopes, we have plotted distributions in initial neutrino energy. Shower energies are related to these neutrino energies in a fairly direct way, and the shower energy distributions will have roughly the same shape. The energy distributions of through-going leptons may be significantly distorted, however, as the lepton energy is a function of black hole multiplicity as well as the distance the lepton propagates before reaching the detector. Note also that, while shower energies should be well-measured, lepton energy measurements present significant challenges, especially for taus. Lollipop tau events will be a clear tau signature with well-measured energy. However, determining the energy of through-going tau events will be difficult in IceCube as they typically do not radiate and, therefore, appear similar to minimum ionizing muons with energy $`200\text{GeV}`$. More detailed detector simulations will be necessary to effectively use the information from through-going taus described in this paper.

## VII Prospects and Conclusions

In the presence of TeV-scale gravity and extra dimensions, ultra-high energy cosmic neutrinos may produce microscopic black holes in the surface of the Earth. We have explored the prospects for detecting such events in neutrino telescopes. We considered contained showers, through-going muons, and also through-going taus.

The rates for a $`\text{km}^3`$ detector, such as IceCube, are of the order of a few per year in each channel for $`M_D`$ near current bounds. These rates are well above SM backgrounds, and provide a significant opportunity to observe black hole production. In several years of data taking at a kilometer-scale neutrino telescope, it should be possible to probe most of the cases described in this paper: values of $`M_D`$ well above 1 TeV may be within reach, and observable event rates are even possible for exponentially suppressed parton cross sections. At the same time, these results imply that black hole observation at ANTARES and NESTOR appears unlikely, given their effective areas of $`0.1\text{km}^2`$. The RICE experiment Kravchenko:2001id , which aims to detect coherent Cerenkov radiation from electromagnetic and hadronic showers at frequencies of 100 MHz to 1 GHz, provides another interesting and complementary probe with effective volume for shower detection reaching $`1\text{km}^3`$ above a few PeV.

The expected rates for black hole events are, of course, highly sensitive to the choice of neutrino flux. As this work was being completed, a paper also discussing black holes at neutrino telescopes appeared Kowalski:2002gb . There, results for two fluxes are presented: a “limit from hidden sources” flux and a cosmogenic flux. The first flux is roughly 3 times higher than the limit placed by the AMANDA experiment on diffuse fluxes Andres:2001ty . The WB flux we use is much more conservative and is roughly 2 orders of magnitude smaller. A variety of cosmogenic fluxes Yoshida:pt ; Protheroe:1995ft  are also considered in Ref. Kowalski:2002gb . Our ratios of black hole to SM rates are in agreement, but here, too, it appears we have chosen a more conservative representative flux, leading to a factor of 3 to 5 fewer shower and muon events from cosmogenic neutrinos. The possibility of tau events, exponentially suppressed cross sections, and angle and energy distributions were not addressed in Ref. Kowalski:2002gb .

Relative to IceCube, the prospects for black hole detection appear to be slightly brighter at the Auger Observatory, expected to be completed by 2004, where tens of events per year may be discovered Feng:2001ib ; Ringwald:2001vk ; Anchordoqui:2001cg . IceCube is scheduled to reach its ultimate goal of $`1\text{km}^3`$ in 2009, but it will collect data continuously as it grows in stages, beginning with the existing AMANDA II detector, roughly of area $`0.1\text{km}^2`$, and reaching a volume of around $`0.5\text{km}^3`$ in 2006. It is difficult to model this time-varying detection volume, but by 2007, an integrated exposure of roughly $`1\text{km}^3\text{yr}`$ will have been achieved. Given the rates of Tables 1-3, we conclude that, for scenarios with $`M_D`$ near current bounds, IceCube may detect black hole events along with the Auger Observatory before the LHC begins operation.

Of course, enhanced event rates alone are not necessarily a signal of new physics. In particular, the uncertainty of astrophysical neutrino fluxes makes it difficult to determine whether observed events result from the SM or new physics based on counting experiments alone. Comparisons of quasi-horizontal showers and Earth-skimming neutrino rates provide a powerful discriminant in ground arrays Feng:2001ib ; Anchordoqui:2001cg . However, the possibility of measuring relative event rates in three separate channels at IceCube provides a direct, unique, and complementary tool for distinguishing black hole events from the SM and other new physics possibilities. As noted in Sec. VI, angular and energy distributions will also be helpful to separate the signal of a large astrophysical flux from the signal of new physics. The separation of these channels in ground arrays and fluorescence detectors is extremely difficult.

Given sufficient statistics, then, the information provided by IceCube may provide information on black hole branching ratios, and the various angle and energy distributions may also help distinguish various properties of black hole production and decay. Black hole production may be the best prospect for discovering extra dimensions in the near future, and neutrino telescopes may be an excellent tool for studying them.

###### Acknowledgements.

We thank Gonzalo Parente for invaluable help with implementing the CTEQ parton distribution functions in our simulation, Ralph Engel for providing the cosmogenic neutrino flux, and Alfred Shapere for helpful conversations. J.A.-M. is supported by NASA Grant NAG5-7009. The work of J.L.F. was supported in part by the Department of Energy under cooperative research agreement DF-FC02-94ER40818. This work was supported in part by DOE grant No. DE-FG02-95ER40896 and in part by the Wisconsin Alumni Research Foundation.


text:
# Classification of three-dimensional exceptional log canonical hypersurface singularities II

## Introduction

In this paper the classification of three-dimensional exceptional hypersurface singularities started in  is completed. The main classification theorem is the following one:

###### Main theorem.

Let $`(X,0)(^4,0)`$ be a three-dimensional exceptional canonical $`(`$respectively strictly log canonical$`)`$ hypersurface singularity defined by a polynomial $`f`$. Then there exists a biholomorphic coordinate change $`\psi :(^4,0)(_{t,z,x,y}^4,0)`$ and unique primitive vector $`𝐩N_{}`$ such that just one of the following two possibilities holds:

1. The quasihomogeneous polynomial $`\stackrel{~}{f}_𝐩=(f\psi )_𝐩`$ defines an exceptional canonical $`(`$respectively strictly log canonical and canonical outside 0$`)`$ singularity $`(X_𝐩,0)(_{t,z,x,y}^4,0)`$. In this case $`𝐩`$-blow-up of $`^4`$ induces purely log terminal blow-ups $`\phi :(Y,E)(X,0)`$ and $`\phi _𝐩:(Y_𝐩,E_𝐩)(X_𝐩,0)`$, where $`(E,\mathrm{Diff}_E(0))=(E_𝐩,\mathrm{Diff}_{E_𝐩}(0))`$. That is, these singularities have the same type and in particular the same complement index.

The canonical singularities satisfying the condition of well-formedness — $`\mathrm{Diff}_{E/(𝐩)}(0)=0`$ are classified in the theorems 3.24, 3.28 of  and in the tables of chapter \[4, §4\]. The polynomial $`\stackrel{~}{f}_𝐩`$; $`(E,\mathrm{Diff}_E(0))`$; minimal complement index are written in the tables.

The canonical singularities not satisfying the condition of well-formedness are classified in the tables of chapter §3. The polynomial $`\stackrel{~}{f}_𝐩`$; $`(E,\mathrm{Diff}_E(0))`$; minimal complement index are written in the tables.

The strictly log canonical and canonical outside 0 quasihomogeneous singularities are always exceptional $`(`$in any dimension$`)`$ by theorem \[4, 2.13\]. In the three-dimensional case $`(E,\mathrm{Diff}_E(0))=(\stackrel{~}{f}_𝐩(𝐩),0)`$ is $`K3`$ surface with Du Val singularities and $`(X,0)`$ is 1-complementary.
2. $`\stackrel{~}{f}_𝐩=t^3+g_2^2(z,x,y)`$, where $`g_2`$ is an irreducible homogeneous polynomial of degree two. In this case the purely log terminal blow-ups are constructed in the theorems \[4, 3.3, 3.5\]. Also it was obtained the similar classification depending on the type of jets $`\stackrel{~}{f}_5`$ and $`\stackrel{~}{f}_6`$.

###### Proof.

According to the main theorem proved in  it is enough to classify the three-dimensional exceptional canonical quasihomogeneous singularities which are not well-formed. The required classification is given in the tables of chapter §3. ∎

It follows from the classification that there are only finite number of types of three-dimensional exceptional log canonical hypersurface singularities. There is a conjecture that the similar result about the finite number of types is true in any dimension for the exceptional singularities. See also examples \[4, 1.19, 1.20\].

The present paper is the direct continuation of . Therefore see  in connection with the notations, definitions and etc.

The paper is organized in the following way. In chapter §1 the purely log terminal blow-ups of not well-formed singularities are studied. In chapter §2 the three-dimensional quasihomogeneous not well-formed singularities are considered. In chapter §3 the summarizing tables are written.

I am grateful to Professor Yu.G.Prokhorov for useful discussions and valuable remarks.

The research was partially supported by a grant 99-01-01132 from the Russian Foundation of Basic Research and a grant INTAS-OPEN 2000#269.

## 1. Structure of purely log terminal blow-up of not well-formed singularity

Now we explain in detail the condition of well-formedness of quasihomogeneous singularity.

Let $`(X,0)(^n,0)`$ be a quasihomogeneous log canonical, canonical outside 0 singularity with the integer weights $`𝐩=(p_1,\mathrm{},p_n)`$ and quasihomogeneous degree $`d`$ defined by a polynomial $`g(x_1,\mathrm{},x_n)`$. Without loss of generality it can be assumed that $`(p_1,\mathrm{},p_n)=1`$.

Consider $`𝐩`$-blow-up $`\psi :^n(𝐩)^n`$. The exceptional divisor of $`\psi `$ is the weighted projective space $`(𝐩)`$. This blow-up induces a plt blow-up $`\psi |_{X𝐩}:(X(𝐩),E)(X,0)`$ (see theorem 2.13 ). Put $`q_i=(p_1,\mathrm{},\widehat{p}_i,\mathrm{},p_n)`$. Then $`E=(\stackrel{~}{g}(x_1,\mathrm{},x_n)(\stackrel{~}{p}_1,\mathrm{},\stackrel{~}{p}_n))=(g(x_1^{1/q_1},\mathrm{},x_n^{1/q_n})(p_1q_1/(q_1\mathrm{}q_n),\mathrm{},p_nq_n/(q_1\mathrm{}q_n))`$. Put $`q_{ij}=(\stackrel{~}{p}_1,\mathrm{},\widehat{\stackrel{~}{p}}_i,\mathrm{},\widehat{\stackrel{~}{p}}_j,\mathrm{},\stackrel{~}{p}_n)`$ if $`ij`$. The quasihomogeneous degree of $`\stackrel{~}{g}(x_1,\mathrm{},x_n)`$ is denoted by $`\stackrel{~}{d}`$.

###### Definition 1.1.

 In the above–mentioned notations the quasihomogeneous singularity is called well–formed singularity if one of the following equivalent conditions holds:

1. $`\mathrm{Diff}_{E/(𝐩)}(0)=0`$;
2. One have $`q_{ij}|\stackrel{~}{d}`$ for all $`ij`$.

In opposite case our singularity is called not well–formed singularity. Let us remark that if we choose the other quasihomogeneous weights then it can happen that the condition of well-formedness is not fulfilled. For the weakly exceptional and consequently for the exceptional singularity the primitive vector $`𝐩`$ is unique \[4, 2.13\] and therefore the notion of well-formedness is uniquely defined.

Put $`I_i=\left\{j\right|1jn,ji,q_{ij}\stackrel{~}{d}\}.`$

###### Proposition 1.2.

$$\mathrm{Diff}_E(0)=\underset{i=1}{\overset{n}{}}\left(1\frac{1}{q_i}\right)\left(C_i\underset{jI_i}{}C_{ij}\right)+\underset{i<j;jI_i}{}\left(1\frac{1}{q_{ij}q_iq_j}\right)C_{ij},$$

where the divisors $`C_i`$ and $`C_{ij}`$ on $`E`$ are defined by the equations $`x_i=0`$ and $`x_i=x_j=0`$ respectively. Moreover $`C_i^{}\stackrel{\mathrm{def}}{=}C_i_{jI_i}C_{ij}`$= $`\left\{x_i=\left(_{jI_i}x_j\right)^1g(x_1,\mathrm{},x_{i1},0,x_{i+1},\mathrm{},x_n)=0\right\}`$, where $`q_i1`$.

###### Proof.

To calculate the different we have to determine the codimension 2 components of $`\mathrm{Sing}X(𝐩)`$ in $`X(𝐩)`$ which are lying in $`E`$. The variety $`X(𝐩)`$ in the affine piece $`U_k=_{x_1,\mathrm{},x_n}^n/_{p_k}(p_1,\mathrm{},p_{k1},1,p_{k+1},\mathrm{},p_n)`$ of $`𝐩`$-blow-up is given by a polynomial $`g^{}=g(x_1,\mathrm{},x_{k1},1,x_{k+1},\mathrm{},x_n)=0`$. It is clear that the singularity defined by this polynomial $`g^{}`$ has the codimension not less then 3 in $`^{n1}\{x_k=0\}`$ and hence doesn’t influence on $`\mathrm{Diff}_E(0)`$ (recall that $`EU_k`$ is given by a polynomial $`g^{}`$ in $`\{x_k=0\}/_{p_k}`$). Thus $`\mathrm{Diff}_E(0)`$ appears only because of the action of cyclic group. Since $`q_{ij}q_iq_j=(p_1,\mathrm{},\widehat{p}_i,\mathrm{},\widehat{p}_j,\mathrm{},p_n)`$ then it remains to show the following statement: if $`q_{ij}\stackrel{~}{d}`$ then $`\{x_i=x_j=0\}`$ defines a divisor on $`E`$. Actually, this condition means that $`x_i`$ or $`x_j`$ divides any monomial of $`g`$. ∎

###### Remark 1.3.

If $`q_{ij}\stackrel{~}{d}`$ then $`g=x_ig_1(x_1,\mathrm{},x_n)+x_jg_2(x_1,\mathrm{},x_n)`$. Moreover, if $`q_i>1`$ and $`q_{ij}\stackrel{~}{d}`$ then $`q_j=1`$. Indeed, if $`q_j>1`$ then $`g=x_i^{q_i}g_1(x_1,\mathrm{},x_n)+x_j^{q_j}g_2(x_1,\mathrm{},x_n)`$. Hence $`g`$ defines a nonnormal singularity.

The next proposition is proved in the same way as the previous proposition 1.2.

###### Proposition 1.4.

$$\mathrm{Diff}_{E/(𝐩)}(0)=\underset{i<j;jI_i}{}\left(1\frac{1}{q_{ij}}\right)C_{ij}.$$

###### Corollary 1.5.

$`(`$cf. \[2, 2.10\]$`)`$ Let $`(X,0)(^n,0)`$ be a quasihomogeneous strictly log canonical, canonical outside 0 singularity with weight p. Then it is well-formed (by theorem 2.13  it is exceptional).

###### Proof.

Let $`\psi :(X(𝐩),E)(X,0)`$ be a plt blow-up. Since $`K_{X(𝐩)}+E=\psi ^{}K_X`$ is Cartier divisor then $`\mathrm{Diff}_E(0)=0`$. Hence by propositions 1.2 and 1.4 $`\mathrm{Diff}_{E/(𝐩)}(0)=0`$. In particular $`E`$ is $`K3`$ surface with Du Val singularities if $`n=4`$ \[1, 1.4.1\]. ∎

###### Example 1.6.

Let $`(X,0)(^3,0)`$ be a weakly exceptional Du Val singularity defined by a polynomial $`h`$. Then $`h`$ is one of the following polynomials: $`x_1^2+x_2^2x_3+x_3^{n1}`$$`n4`$ (type $`𝔻_n`$); $`x_1^2+x_2^3+x_3^4`$ (type $`𝔼_6`$); $`x_1^2+x_2^3+x_3^3x_2`$ (type $`𝔼_7`$); $`x_1^2+x_2^3+x_3^5`$ (type $`𝔼_8`$) (see \[5, 4.7\]). Consider $`𝐩`$-blow-ups corresponding to their quasihomogeneous weights. Then the singularities $`𝔻_n`$, $`n5`$ and $`𝔼_7`$ are not well-formed.

| Singularity | $`E^1`$ | $`\mathrm{Diff}_{E/(𝐩)}(0)`$ | $`\mathrm{Diff}_E(0)`$ |
| --- | --- | --- | --- |
| $`D_{2n},n3`$ | $`x_1+x_2^2x_3+x_3^{2n1}`$ | $`\frac{n2}{n1}P_3`$ | $`\frac{1}{2}P_1+\frac{1}{2}P_2+`$ |
|  | $`(2n1,n1,1)`$ |  | $`+\frac{2n3}{2n2}P_3`$ |
| $`D_{2n+1},n2`$ | $`x_1^2+x_2x_3+x_3^{2n}`$ | $`\frac{2n2}{2n1}P_3`$ | $`\frac{1}{2}P_1+\frac{1}{2}P_2+`$ |
|  | $`(n,2n1,1)`$ |  | $`+\frac{2n2}{2n1}P_3`$ |
| $`E_7`$ | $`x_1+x_2^3+x_3x_2`$ | $`\frac{1}{2}P_3`$ | $`\frac{1}{2}P_1+\frac{2}{3}P_2+`$ |
|  | $`(3,1,2)`$ |  | $`+\frac{3}{4}P_3`$ |

###### Definition 1.7.

Let us define a $``$-divisor $`\widehat{D}=_{i=1}^n(1\frac{1}{q_i})\{x_i=0\}`$ on $`(𝐩)`$ and a $``$-divisor $`D=_{i=1}^n(1\frac{1}{q_i})C_i`$ on $`E`$. Recall that, if $`(X,0)`$ is a well-formed singularity then $`D=\widehat{D}|_E=\mathrm{Diff}_E(0)`$ (cf. \[4, 3.12\]). Nevertheless the following proposition takes place in the general case.

###### Proposition 1.8.

$`K_E+\mathrm{Diff}_E(0)=(K_{(𝐩)}+E+\widehat{D})|_E`$.

###### Proof.

By proposition 1.4, remark 1.3 and corollary \[6, 3.10\] $`(K_{(𝐩)}+E+\widehat{D})|_E=K_E+\mathrm{Diff}_{E/(𝐩)}(\widehat{D})=K_E+_{i<j;jI_i}\left[(1\frac{1}{q_{ij}})C_{ij}+\frac{1}{q_{ij}}(1\frac{1}{q_iq_j})C_{ij}\right]+_{i=1}^n(1\frac{1}{q_i})(C_i_{jI_i}C_{ij})`$. By proposition 1.2 the last expression is equal to $`K_E+\mathrm{Diff}_E(0)`$. ∎

###### Example 1.9.

Let us return to the example 1.6. Then

| Singularity | $`D`$ | $`\mathrm{Diff}_E(0)=\mathrm{Diff}_{E/(𝐩)}(\widehat{D})`$ |
| --- | --- | --- |
| $`D_{2n},n3`$ | $`\frac{1}{2}P_1+\frac{1}{2}P_2+\frac{1}{2}P_3`$ | $`\frac{1}{2}P_1+\frac{1}{2}P_2+`$ |
|  |  | $`+\left(\frac{n2}{n1}+\frac{1}{n1}\frac{1}{2}\right)P_3`$ |
| $`D_{2n+1},n2`$ | $`\frac{1}{2}P_1+\frac{1}{2}P_2`$ | $`\frac{1}{2}P_1+\frac{1}{2}P_2+`$ |
|  |  | $`+\left(\frac{2n2}{2n1}+\frac{1}{2n1}0\right)P_3`$ |
| $`E_7`$ | $`\frac{1}{2}P_1+\frac{2}{3}P_2+\frac{1}{2}P_3`$ | $`\frac{1}{2}P_1+\frac{2}{3}P_2+`$ |
|  |  | $`+\left(\frac{1}{2}+\frac{1}{2}\frac{1}{2}\right)P_3`$ |

## 2. Investigation of three-dimensional not well-formed singularities on exceptionality

The next theorem is proved by exhaustion of all cases. This sorting is described in .

###### Theorem 2.1.

All three-dimensional not well-formed exceptional canonical quasihomogeneous singularities can be obtained by the rotations of Newton’s line passing through the following parts of type $`_2`$:

$`(1)`$ $`t^2+z^3+x^m`$, where $`x^m`$$`zx^5`$, $`zx^5y`$, $`zx^7`$, $`zx^6y`$, $`zx^5y^2`$.

$`(2)`$ $`t^2+z^4+x^m`$, where $`x^m`$$`zx^4`$, $`zx^5`$, $`zx^4y`$.

$`(3)`$ all parts of type $`\mathrm{{\rm Y}}_2^{(2)}`$ $`(`$they have the equation of the form $`t^2+z^3x+h(z,x,y))`$.

$`(4)`$ $`t^2+z^3y+x^m`$, where $`x^m`$$`x^7y`$, $`x^8y`$, $`zx^5`$, $`zx^6`$, $`zx^5y`$.

$`(5)`$ $`t^2+z^5+zx^5`$; $`t^2+z^4x+x^m`$, where $`x^m`$$`x^6`$, $`x^5y`$, $`x^4y^2`$, $`zx^5`$, $`zx^4y`$, $`zx^3y^2`$; $`t^2+z^4y+x^m`$, where $`x^m`$$`x^5y`$, $`zx^5`$, $`zx^4y`$; $`t^2+z^3x^2+x^m`$, where $`x^m`$$`zx^5`$, $`zx^4y`$, $`zx^3y^2`$; $`t^2+z^3y^2+zx^5`$; $`t^2+z^3xy+x^m`$, where $`x^m`$$`x^6`$, $`x^5y`$, $`x^7`$, $`x^6y`$, $`x^5y^2`$, $`zx^5`$, $`zx^4y`$; $`t^2+z^5+zx^3y`$; $`t^2+z^4x+zx^3y`$; $`t^2+z^4x+x^4y`$; $`t^2+z^3x^2+zx^3y`$.

$`(6)`$ $`t^3+z^2xy+x^5`$; $`t^3+z^2xy+x^4y`$; $`t^3+z^4+zx^2y`$; $`t^3+z^4+tx^3`$; $`t^3+z^3x+x^m`$, where $`x^m`$$`zx^2y`$, $`x^3y`$, $`tx^3`$, $`tx^2y`$; $`t^3+z^3y+tx^3`$; $`t^3+tz^3+x^m`$, where $`x^m`$$`z^2x^2`$, $`zx^3`$, $`zx^2y`$; $`t^3+tz^2x+zx^2y`$; $`t^3+tz^2x+x^2y^2`$; $`g_3(t,z)+zx^3`$, where $`g_3`$ is a binary form of degree 3; $`t^3+z^2x+x^m`$, where $`x^m`$$`x^4`$, $`x^5`$, $`x^4y`$, $`x^6+tx^4`$, $`x^5y`$, $`x^4y^2`$, $`x^7`$, $`x^6y`$, $`x^5y^2`$, $`x^4y^3`$, $`x^8`$, $`x^7y`$, $`x^6y^2`$, $`x^5y^3`$, $`x^4y^4`$, $`tx^3`$, $`tx^3y`$, $`tx^5`$, $`tx^4y`$, $`tx^3y^2`$; $`t^3+z^2y+x^7y`$; $`t^3+z^2y+tx^5`$.

$`(7)`$ $`t^2z+g_4(z,x)`$, where $`g_4`$ is a binary form of degree 4; $`t^2z+z^4+x^m`$, where $`x^m`$$`x^3y`$, $`x^5`$, $`x^4y`$, $`x^3y^2`$, $`tx^3`$, $`tx^2y`$;$`t^2z+z^3x+x^m`$, where $`x^m`$$`x^3y`$, $`x^5`$, $`x^4y`$, $`x^3y^2`$, $`tx^3`$, $`tx^2y`$; $`t^2z+z^3y+x^m`$, where $`x^m`$$`x^5`$, $`x^4y`$, $`zx^3`$, $`z^2x^2`$, $`tx^3`$; $`t^2z+z^2x^2+x^m`$, where $`x^m`$$`x^5`$, $`x^4y`$, $`x^3y^2`$, $`tx^2y`$; $`t^2z+z^2xy+x^5`$; $`t^2z+z^2xy+x^4y`$; $`t^2z+z^2x+x^m`$, where $`x^m`$$`x^4`$, $`zx^3+x^5`$, $`x^4y`$, $`x^6`$, $`x^5y`$, $`x^4y^2`$, $`tx^3`$, $`tx^4`$, $`tx^3y`$; $`t^2z+z^2y+x^m`$, where $`x^m`$$`x^5y`$, $`zx^3`$, $`tx^4`$.

$`(8)`$ $`t^2x+z^4+x^m`$, where $`x^m`$$`x^3y`$, $`x^5`$, $`x^4y`$, $`zx^2y`$, $`zx^4`$, $`zx^3y`$, $`zx^2y^2`$; $`t^2x+z^3x+x^m`$, where $`x^m`$$`x^2y^2`$, $`x^5`$, $`x^4y`$, $`x^3y^2`$, $`x^2y^3`$, $`x^6`$, $`x^5y`$, $`x^4y^2`$, $`x^3y^3`$, $`x^2y^4`$, $`zx^2y`$, $`zx^4`$, $`zx^3y`$, $`zx^2y^2`$; $`t^2x+z^3y+x^m`$, where $`x^m`$$`x^5`$, $`x^4y`$, $`x^6`$, $`x^5y`$, $`x^4y^2`$, $`zx^4`$, $`zx^3y`$, $`z^2x^2`$.

$`(9)`$ $`t^2y+z^4+zx^4`$; $`t^2y+z^3x+x^m`$, where $`x^m`$$`x^5`$, $`x^6`$, $`x^5y`$, $`zx^4`$, $`tx^3`$.

In investigating given singularities on the exceptionality it is convenient to use the proposition \[4, 3.13\] and corollary \[4, 3.14\]. In general case corollary \[4, 3.15\] is formulated in the following way.

###### Proposition 2.2.

Let $`(X,0)(_{x_1,\mathrm{},x_4}^4,0)`$ be a three-dimensional quasihomogeneous log canonical singularity with the weight p. Assume that $`(X,H)`$ is not lc for any hyperplane section $`H`$. Let $`D=d_iD_i`$ (see definition 1.7). If $`d_k\frac{6}{7}`$ for some $`k`$ then $`(X,0)`$ is exceptional.

###### Proof.

Suppose that $`(X,0)`$ is not exceptional. Let $`H_X=\{x_k=0\}|_X`$ be a hyperplane section and let $`H_{X(𝐩)}`$ be a proper transform of $`H_X`$ on $`X(𝐩)`$. Then $`K_{X(𝐩)}+E+H_{X(𝐩)}=\psi ^{}(K_X+H_X)+aE`$. Let $`C_k=_iC_{ki}`$, where $`C_{ki}`$ are the irreducible curves (see proposition 1.2). One have $`\mathrm{Diff}_E(0)=_id_{ki}C_{ki}+\mathrm{\Xi }`$. The coefficients $`d_{ki}6/7`$ for all $`i`$ since $`d_{ki}=\frac{m_i1}{m_i}+\frac{1}{m_i}d_kd_k`$ \[6, 3.10\]. By proposition \[4, 3.14\] there exists 1,2,3,4 or 6-complement $`D^+`$. By the definition of complement we have $`D^+_iC_{ki}+\mathrm{\Xi }`$. Therefore $`(K_E+_iC_{ki}+\mathrm{\Xi })`$ is nef and lc. Since $`\mathrm{Diff}_E(H_{X(𝐩)})=\mathrm{Diff}_{E/(𝐩)}\left(\widehat{D}+\{x_k=0\}\right)`$, where $`\{x_k=0\}`$ is the corresponding divisor on $`(𝐩)`$ then $`\mathrm{Diff}_E(H_{X(𝐩)})=_iC_{ki}+\mathrm{\Xi }`$. Hence $`a0`$ and $`(X,H_X)`$ is lc. This contradiction concludes the proof. ∎

###### Remark 2.3.

Proposition 2.2 is hypothetically true in any dimension (see \[4, 3.16\]). In the two-dimensional case we must require that $`d_k\frac{2}{3}`$ for some $`k`$.

###### Example 2.4.

In proposition 2.2 we can’t change $`D`$ on $`\mathrm{Diff}_E(0)`$. Let us show it in the following two examples.

(1) Let $`(X,0)(^3,0)`$ be a Du Val singularity of type $`𝔻_n`$ defined by a polynomial $`x_1^2+x_2^2x_3+x_3^{n1}`$, where $`n4`$. The pair $`(X,H)`$ is not lc for any hyperplane section $`H`$. If $`n5`$ then one of the coefficients of $`\mathrm{Diff}_E(0)`$ is not less then $`2/3`$ and all coefficients of $`D`$ are equal to $`1/2`$. This singularity is not exceptional by the definition since $`K_X+\{x_1=x_3=0\}_{}0`$ is 2-complement (see also example \[4, 1.8\]).

(2) Let $`(X,0)(_{t,z,x,y}^4,0)`$ be a canonical singularity defined by a polynomial $`t^3+z^2x+x^4+xy^5`$. The pair $`(X,H)`$ is not lc for any hyperplane section $`H`$. The variety $`X`$ has $`cA_2`$ singularity along the curve $`C=\{t=x=z^2+y^5=0\}=\mathrm{Sing}X`$ (outside 0). The divisor $`T=\{t=x=0\}`$ on $`X`$ is $``$-Cartier since $`3T\{x=0\}`$. Therefore the divisor $`(K_X+T)|_T=K_{_{z,y}^2}+\frac{2}{3}C`$ is klt. Thus $`(X,0)`$ is not weakly exceptional singularity \[5, 4.6\].

Now consider $`(40,45,30,18)`$-blow-up of $`^4`$. It induces a plt blow-up of $`(X,0)`$. Then $`(E,\mathrm{Diff}_E(0))=\left(t+zx+x^4+xy(4,3,1,3),\frac{2}{3}(\{t=0\}\backslash \{t=x=0\})+\frac{1}{2}\{z=0\}+\frac{4}{5}\{y=0\}+\frac{8}{9}\{t=x=0\}\right)=(^2,\frac{2}{3}l_1+\frac{1}{2}l_2+\frac{4}{5}l_3+\frac{8}{9}l_4)`$, where $`l_i`$ are the straight lines in general position. The divisor $`K_^2+\frac{2}{3}l_1+\frac{1}{2}l_2+\frac{5}{6}l_3+l_4_{}0`$ is klt 6-complement of minimal index. Since the divisor $`K_^2+\mathrm{Diff}_E(0)+\frac{1}{9}l_4`$ is anti-ample and non-klt then by criterion \[3, 2.1\] $`(X,0)`$ is not weakly exceptional singularity again. In this example all coefficients of $`D`$ are not more then $`\frac{4}{5}`$.

It is very interesting to compare the following example 2.5 with the previous example 2.4(2).

###### Example 2.5.

Let $`(X,0)(_{t,z,x,y}^4,0)`$ be a canonical singularity defined by a polynomial $`t^3+z^2x+tx^3+ty^5`$. The variety $`X`$ has $`cA_1`$ singularity along the curve $`C=\{t=z=x^3+y^5=0\}=\mathrm{Sing}X`$ (outside 0). Consider the divisor $`T=\{t=z=0\}`$ on $`X`$. If $`T`$ had been $``$-Cartier divisor then the divisor $`(K_X+T)|_T=K_{_{x,y}^2}+\frac{1}{2}C`$ would have been klt and consequently $`(X,0)`$ would have been not weakly exceptional singularity. Now we prove that $`T`$ is not $``$-Cartier divisor and $`(X,0)`$ is an exceptional singularity.

Consider $`𝐩=(30,35,20,12)`$-blow-up of $`^4`$. It induces a plt blow-up $`\psi :(X(𝐩),E)(X,0)`$. Then $`E=\left(t^3+zx+tx^3+ty(3,7,2,6)\right)`$, $`\mathrm{Diff}_E(0)=\frac{1}{2}\{z=t^2+x^3+y=0\}+\frac{3}{4}\{t=z=0\}+\frac{4}{5}\{y=0\}=\frac{1}{2}\mathrm{\Delta }+\frac{3}{4}\mathrm{\Gamma }_2+\frac{4}{5}L`$. In the next picture we illustrate the situation of these curves and the singularities of surface $`E`$.

Assume that $`T`$ is $``$-Cartier divisor. Then $`K_{X(𝐩)}+E+T_Y=\psi ^{}(K_X+T)+aE`$ is anti-ample divisor over $`X`$ since $`a>0`$. Hence $`K_E+\mathrm{Diff}_E(T_Y)=K_E+\frac{1}{2}\mathrm{\Delta }+\mathrm{\Gamma }_2+\frac{4}{5}L`$ is anti-ample divisor. It is not true since $`(K_E+\mathrm{Diff}_E(T_Y))\mathrm{\Delta }=(K_E+\frac{1}{2}\mathrm{\Delta }+\frac{3}{4}\mathrm{\Gamma }_2+\frac{11}{12}L)\mathrm{\Delta }+(\frac{1}{4}\mathrm{\Gamma }_2+(\frac{4}{5}\frac{11}{12})L)\mathrm{\Delta }=(K_{(𝐩)}+E+\frac{1}{2}\{z=0\}+\frac{11}{12}\{y=0\})\mathrm{\Delta }+\frac{1}{4}\frac{7}{60}=0+\frac{1}{4}\frac{7}{60}>0`$. Here $`\{z=0\}`$ and $`\{y=0\}`$ are denoted the corresponding divisors on $`(𝐩)`$.

The divisor $`K_E+\mathrm{Diff}_{E/(𝐩)}\left(\frac{3}{5}\{z=0\}+\frac{4}{5}\{y=0\}\right)=K_E+\frac{3}{5}\mathrm{\Delta }+\frac{4}{5}\mathrm{\Gamma }_2+\frac{4}{5}L`$ is 5-complement of minimal index. It is easy to show that there is only the following 6-complement $`K_E+\mathrm{Diff}_{E/(𝐩)}\left(\frac{1}{6}\{t=0\}+\frac{1}{2}\{z=0\}+\frac{5}{6}\{y=0\}\right)=K_E+\frac{1}{2}\mathrm{\Delta }+\frac{5}{6}\mathrm{\Gamma }_2+\frac{5}{6}L+\frac{1}{6}\{t=x=0\}_{}0`$ among 1,2,3,4 and 6-complements. It is clear that this divisor is klt. Therefore $`(X,0)`$ is an exceptional singularity.

Finally let us consider one more example showing the investigation process of given singularities on the exceptionality.

###### Example 2.6.

Let $`(X,0)(_{t,z,x,y}^4,0)`$ be a canonical singularity defined by a polynomial $`f=t^2x+z^3x+zx^2y+atz^2y+bz^2y^3+cxy^4`$. Consider $`𝐩=(6,4,5,3)`$-blow-up of $`^4`$. Then $`(E,\mathrm{Diff}_E(0))=(f(6,4,5,3),\frac{2}{3}\mathrm{\Delta }_3+\frac{1}{2}\mathrm{{\rm Y}}_4)`$, where $`\mathrm{\Delta }_3=\{z=x=0\}`$ and $`\mathrm{{\rm Y}}_4=\{x=y=0\}`$. Let us prove that the divisor $`(K_{(𝐩)}+E+\frac{1}{4}\{z=0\})|_E=K_E+\frac{3}{4}\mathrm{\Delta }_3+\frac{1}{2}\mathrm{{\rm Y}}_4+\frac{1}{2}\{t=z=0\}`$ is non-klt 4-complement if $`c=0`$ (then $`b0`$). It implies that $`(X,0)`$ is not exceptional singularity if $`c=0`$. Since the klt and lc properties of pair remain true by the finite dominant morphisms then it is enough to prove that $`(Z,\frac{1}{4}T)`$ is lc, but not klt pair, where $`Z=\{t^2x+z^3x+zx^2+atz^2+z^2=0\}_{t,z,x}^3`$ and $`T=\{z=0\}|_Z`$. Note that $`(Z,0)`$ is the singularity of type $`𝔻_5`$. The blow-up of $`^3`$ with the weights $`(3,4,2)`$ induces a blow-up $`\psi :(\stackrel{~}{Z},C)Z`$. Besides $`K_{\stackrel{~}{Z}}+C+\frac{1}{4}T_{\stackrel{~}{Z}}=\psi ^{}(K_Z+\frac{1}{4}T)`$. The exceptional divisor $`C`$ is not well-formed hypersurface $`tx+zx^2+z^2(3,2,1)`$. The divisor $`K_C+\mathrm{Diff}_C(\frac{1}{4}T_{\stackrel{~}{Z}})=K_C+\frac{1}{2}P_1+(\frac{1}{2}+\frac{1}{2}\frac{1}{2})P_2+(\frac{2}{3}+\frac{1}{3}\frac{1}{4})P_3`$ is klt. By Inversion of Adjunction we complete the proof. If $`c0`$ then by the same argument one can prove that $`(X,0)`$ is an exceptional singularity.

## 3. Classification of three-dimensional exceptional not well-formed singularities

The comments to the tables. All agreements of §4  about the notations and records in the tables remain valid. The underlying monomials (see point (I) of §4 ) are written without the coefficients before them.

If log Del Pezzo surface $`E=(\stackrel{~}{f}(\stackrel{~}{p}_1,\mathrm{},\stackrel{~}{p}_4)`$ is a linear cone, i.e. $`\stackrel{~}{d}=\stackrel{~}{p}_k`$ for some $`k`$ then we write the simplified form of $`E`$

$$E=(\frac{\stackrel{~}{p}_1}{s_1},\mathrm{},\frac{\stackrel{~}{p}_{k1}}{s_{k1}},\frac{\stackrel{~}{p}_{k+1}}{s_{k+1}},\mathrm{},\frac{\stackrel{~}{p}_4}{s_4}),\text{where}s_m=\underset{\begin{array}{c}q_{ij}\stackrel{~}{d}\\ i<j;i,jm\end{array}}{}q_{ij}$$

For the record of different $`\mathrm{Diff}_E(0)`$ one use the following notations. The curves $`C_i^{}`$, where $`i=1,\mathrm{},4`$ (see proposition 1.2) are denoted by $`\mathrm{\Gamma }`$, $`\mathrm{\Delta }`$, $`\mathrm{{\rm Y}}`$, $`\mathrm{\Omega }`$ respectively. The curves $`\{t=z=0\}`$, $`\{t=x=0\}`$, $`\{t=y=0\}`$, $`\{z=x=0\}`$, $`\{z=y=0\}`$, $`\{x=y=0\}`$ are denoted by $`\mathrm{\Gamma }_2`$, $`\mathrm{\Gamma }_3`$, $`\mathrm{\Gamma }_4`$, $`\mathrm{\Delta }_3`$, $`\mathrm{\Delta }_4`$, $`\mathrm{{\rm Y}}_4`$ respectively.

If there is a symbol * in the description of different cases then this marked case has the other different. Then the different differs by the addition of other curve marked by * or the different is completely another one marked by * also.

1. Singularity – $`t^2+z^3+g(z,x,y)`$

2. Singularity – $`t^2+z^4+g(z,x,y)`$

3. Singularity – $`t^2+z^3x+g(z,x,y)`$

4. Singularity – $`t^2+z^3y+g(z,x,y)`$

5. Singularity – $`t^2+g_5+g(z,x,y)`$.

6. Singularity – $`t^3+g(z,x,y)`$.

7. Singularity – $`t^2z+g(z,x,y)`$.

8. Singularity – $`t^2x+g(z,x,y)`$.

9. Singularity – $`t^2y+g(z,x,y)`$.


text:
# 1 Introduction

## 1 Introduction

Fermions at finite density or finite chemical potential is a subject of a wide range of interest. It is relevant to condensed matter physics, such as the Hubbard model away from half-filling. The research about nuclei and neutron stars at low and high nucleon density is actively pursued in nuclear physics and astrophysics. The subject of quark gluon plasma is important for understanding the early universe and is being sought for in relativistic heavy-ion collisions in the laboratories. Furthermore, speculation about color superconducting phase has been proposed recently for quantum chromodynamics (QCD) at very high quark density .

Although there are models, e.g. chiral models and Nambu Jona-Lasinio model which have been used to study QCD at finite quark density, the only way to study QCD at finite density and temperature reliably and systematically is via lattice gauge calculations. There have been extensive lattice calculations of QCD at finite temperature . On the contrary, the calculation at finite density is hampered by the lack of a viable algorithm.

In this talk, I shall first review the difficulties associated with the finite density algorithm with chemical potentials in Sec. 2. I will then outline in Sec. 3 a proposal for a finite density algorithm in the canonical ensemble which projects out the nonzero baryon number sector from the fermion determinant. In Sec. 4, a newly developed Noisy Monte Carlo algorithm which admits unbiased estimate of the probability is described. Its application to the fermion determinant is outlined in Sec. 5. I will discuss an efficient way, the Padé-Z2 method, to estimate the $`Tr\mathrm{log}`$ of the fermion matrix in Sec. 6. The recent progress on the implementation of the Kentucky Noisy Monte Carlo algorithm to dynamical fermions is presented in Sec. 7. Finally, a summary is given in Sec. 8.

## 2 Finite Chemical Potential

The usual approach to the finite density in the Euclidean path-integral formalism of lattice QCD is to consider the grand canonical ensemble with the partition function

$$Z_{GC}(\mu )=\underset{N}{}Z_Ne^{\mu N}=𝒟U𝑑etM[U,\mu ]e^{S_g[U]},$$
(1)
where the fermion fields with fermion matrix $`M`$ has been integrated to give the determinant. $`U`$ is the gauge link variable and $`S_g`$ is the gauge action. The chemical potential is introduced to the quark action with the $`e^{\mu a}`$ factor in the time-forward hopping term and $`e^{\mu a}`$ in the time-backward hopping term. Here $`a`$ is the lattice spacing. However, this causes the fermion action to be non-Hermitian, i.e. $`\gamma _5M\gamma _5M`$. As a result, the fermion determinant $`detM[U]`$ is complex and this leads to the infamous sign problem.

There are several approaches to avoid the sign problem:

### 2.1 Fugacity Expansion

It was proposed by the Glasgow group  that the sign problem can be circumvented based on the expansion of the grand canonical partition function in powers of the fugacity variable $`e^{\mu /T}`$,

$$Z_{GC}(\mu /T,T,V)=\underset{B=3V}{\overset{B=3V}{}}e^{\mu /TB}Z_B(T,V),$$
(2)
where $`Z_B`$ is the canonical partition function for the baryon sector with baryon number $`B`$. $`Z_{GC}`$ is calculated with reweighting of the fermion determinant

$$Z_{GC}(\mu )=\frac{detM[U,\mu ]}{detM[U,0]}_{\mu =0}.$$
(3)
Since the reweighting is based on the gauge configuration with $`\mu =0`$, it avoids the sign problem. However, this does not work, except perhaps at small $`\mu `$ or near the finite temperature phase transition. We will dwell on this later in Sec. 3. This is caused by the ‘overlap problem’  where the important samples of configurations in the $`\mu =0`$ simulation has exponentially small overlap with those relevant for the finite density. As a result, the onset of baryon begins at $`\mu m_\pi /2`$ instead of the expected $`M_N/3`$ which resembles the situation of the quenched approximation.

### 2.2 Imaginary Chemical Potential

In this approach, the chemical potential is taking an imaginary value $`\mu =i\nu `$. The fermion determinant is real in this case and one can avoid the sign problem . The partition function is

$$Z_{GC}(i\nu /T,T,V)=Tre^{\widehat{H}/T}e^{i\nu \widehat{B}/T},$$
(4)
which is periodic with respect to $`\nu `$ with a period of $`2\pi T`$. Comparing with Eq. (2), one can in principle obtain canonical partition function $`Z_B`$ from the Fourier transform

$$Z_B(T,V)=\frac{1}{2\pi T}_0^{2\pi T}𝑑\nu Z_{GC}(i\nu /T,T,V)e^{i\nu \widehat{B}/T}.$$
(5)
In this approach, one needs to integrate over the whole range of $`\nu `$ from $`0`$ to $`2\pi T`$ after one obtains the Monte Carlo configurations of $`Z_{GC}(i\nu /T,T,V)`$ at different $`\nu `$. In practice, it is proposed to calculate the following ratio in the two-dimensional Hubbard model ,

$$\frac{Z_{GC}(i\nu /T,T,V)}{Z_{GC}(i\nu _0/T,T,V)}=𝒟\varphi e^{S_{bos}}𝑑etM(i\nu _0)\frac{detM(i\nu )}{detM(i\nu _0)},$$
(6)
with a reference value $`\nu _0`$. Several patches each centered around a different reference point $`\nu _0`$ are used to cover the range of $`\nu `$. This was successful for the two-dimensional Hubbard model with a $`4^2\times 10`$ lattice up to B = 6 where the determinant was calculated exactly. While this works for a small lattice in the Hubbard model, it would not work for reasonably large lattices in QCD. This is because the direct calculation of the determinant is a $`V^3`$ (or $`V^2`$ for a sparse matrix) operation which is an impracticable task for the quark matrix which is typically of the dimension $`10^6\times 10^6`$. Any stochastic estimation of the determinant will inevitably introduce systematic error. Furthermore, this will also suffer from the ‘overlap’ problem discussed above. Any Monte Carlo simulation at a reference point $`\nu _0`$ will have exponentially small overlap with those configurations important to a nonzero baryon density.

### 2.3 Overlap Ensuring Multi-parameter Reweighting

To alleviate the sign problem with the real chemical potential and the overlap problem due to reweighting, it is proposed  to do the reweighting in the multiple parameter space. The generic partition function $`Z_{GC}`$ in Eq. (1) is parametrized by a set of parameters $`\alpha `$, such as the chemical potential $`\mu `$, the gauge coupling $`\beta `$, the quark mass $`m_q`$, etc. The partition function can be written to facilitate reweighting

$$Z_{GC}(\alpha )=𝒟U𝑑etM[U,\alpha _0]e^{S_g[U,\alpha _0]}\{e^{S_g[U,\alpha ]+S_g[U,\alpha _0]}\frac{detM[U,\alpha ]}{detM[U,\alpha _0]}\},$$
(7)
where the Monte Carlo simulation is carried out with the $`\alpha _0`$ set of parameters and the terms in the curly bracket are treated as observables. This is applied to study the end point in the T-$`\mu `$ phase diagram. In this case, the Monte Carlo simulation is carried out where the parameters in $`\alpha _0`$ include $`\mu =0`$ and $`\beta _c`$ which corresponds to the phase transition at temperature $`T_c`$. The parameter set $`\alpha `$ in the reweighted measure include $`mu0`$ and an adjusted $`\beta `$ in the gauge action. The new $`\beta `$ is determined from the Lee-Yang zeros so that one is following the transition line in the T-$`\mu `$ plane and the large change in the determinant ratio in the reweighting is compensated by the change in the gauge action to ensure reasonable overlap. This is shown to work to locate the transition line from $`\mu =0`$ and $`T=T_c`$ down to the critical point on the $`4^4`$ and $`6^3\times 4`$ lattices with staggered fermions .

While the multi-parameter reweighting is successful near the transition line, it is not clear how to extend it beyond this region, particularly the $`T=0`$ case where one wants to keep the $`\beta `$ and quark mass fixed while changing the $`\mu `$. One still expects to face the overlap problem in the latter case. For large volumes, calculating the determinant ratio will be subjected to the same practical difficulty as discussed in the previous section 2.2.

## 3 Finite Baryon Density – A Canonical Ensemble Approach

We would like to propose an algorithm to overcome the overlap problem at zero temperature which is based on the canonical ensemble approach. To avoid the overlap problem, one needs to lock in a definite nonzero baryon sector so that the exponentially large contamination from the zero-baryon sector is excluded. To see this, we first note that the fermion determinant is a superposition of multiple quark loops of all sizes and shapes. This can be easily seen from the property of the determinant

$$detM=e^{Tr\mathrm{log}M}=1+\underset{n=1}{}\frac{(Tr\mathrm{log}M)^n}{n!}.$$
(8)
Upon a hopping expansion of $`\mathrm{log}M`$, $`Tr\mathrm{log}M`$ represents a sum of single loops with all sizes and shapes. The determinant is then the sum of all multiple loops. The fermion loops can be separated into two classes. One is those which do not go across the time boundary and represent virtual quark-antiquark pairs; the other includes those which wraps around the time boundary which represent external quarks and antiquarks. The configuration with a baryon number one which entails three quark loops wrapping around the time boundary will have an energy $`M_B`$ higher than that with zero baryon number. Thus, it is weighted with the probability $`e^{M_BN_ta_t}`$ compared with the one with no net baryons. We see from the above discussion that the fermion determinant contains a superposition of sectors of all baryon numbers, positive, negative and zero. At zero temperature where $`M_BN_ta_t1`$, the zero baryon sector dominates and all the other baryon sectors are exponentially suppressed. It is obvious that to avoid the overlap problem, one needs to select a definite nonzero baryon number sector and stay in it throughout the Markov chain of updating configurations. To select a particular baryon sector from the determinant can be achieved by the following procedure : first, assign an $`U(1)`$ phase factor $`e^{i\varphi }`$ to the links between the time slices $`t`$ and $`t+1`$ so that the link $`U/U^{}`$ is multiplied by $`e^{i\varphi }/e^{i\varphi }`$; then the particle number projection can be carried out through the Fourier transformation of the fermion determinant like in the BCS theory

$$P_N=\frac{1}{2\pi }_0^{2\pi }𝑑\varphi e^{i\varphi N}𝑑etM[\varphi ]$$
(9)
where $`N`$ is the net particle number, i.e. particle minus antiparticle. Note that all the virtual quark loops which do not reach the time boundary will have a net phase factor of unity; only those with a net N quark loops across the time boundary will have a phase factor $`e^{i\varphi N}`$ which can contribute to the integral in Eq. (9). Since QCD in canonical formulation does not break $`Z(3)`$ symmetry, it is essential to take care that the ensemble is canonical with respect to triality. To this end, we shall consider the triality projection  to the zero triality sector

$$det_0M=\frac{1}{3}\underset{k=0,\pm 1}{}detM[\varphi +k2\pi /3].$$
(10)
This amounts to limiting the quark number N to a multiple of 3. Thus the triality zero sector corresponds to baryon sectors with integral baryon numbers.

Another essential ingredient to circumvent the overlap problem is to stay in the chosen nonzero baryon sector so as to avoid mixing with the zero baryon sector with exponentially large weight. This can be achieved by preforming the baryon number projection as described above before the accept/reject step in the Monte Carlo updating of the gauge configuration. If this is not done, the accepted gauge configuration will be biased toward the zero baryon sector and it is very difficult to project out the nonzero baryon sector afterwords. This is analogous to the situation in the nuclear many-body theory where it is known  that the variation after projection (Zeh-Rouhaninejad-Yoccoz method ) is superior than the variation before projection (Peierls-Yoccoz method ). The former gives the correct nuclear mass in the case of translation and yields much improved wave functions in mildly deformed nuclei than the latter.

To illustrate the overlap problem, we plot in Fig.1 $`Tr\mathrm{log}M[\varphi ]`$ for a configuration of the $`8^3\times 12`$ lattice with the Wilson action with $`\beta =6.0`$ and $`\kappa =0.150`$ which is obtained with 500 $`Z_2`$ noises. We see that the it is rather flat in $`\varphi `$ indicating that the Fourier transform in Eq. (9) will mainly favor the zero baryon sector. On the other hand, at finite temperature, it is relatively easier for the quarks to be excited so that the zero baryon sector does not necessarily dominate other baryon sectors. Another way of seeing this is that the relative weighting factor $`e^{M_BN_ta_t}`$ can be $`O(1)`$ at finite temperature. Thus, it should be easier to project out the nonzero baryon sector from the determinant. We plot in Fig. 2 a similarly obtained $`Tr\mathrm{log}M[\varphi ]`$ for a configuration of the $`8\times 20^2\times 4`$ lattice with $`\beta =4.9`$ and $`\kappa =0.182`$. We see from the figure that there is quite a bit of wiggling in this case as compared to that in Fig. 1 indicating that it is easier to project out a nonzero baryon sector through the Fourier transform at finite temperature.

We should mention that while we think we can overcome the overlap problem and the determinant $`detM[\varphi ]`$ is real in this approach, nevertheless in view of the fact that the Fourier transform in Eq. (9) involves the quark number $`N`$ the canonical approach may still have the sign problem at the thermodynamic limit when $`N`$ and $`V`$ are very large. However, we think it might work for small $`N`$ such as 3 or 6 for one or two baryons in a finite $`V`$. This should be a reasonable start for practical purposes.

While it is clear what the algorithm in the canonical approach entails, there are additional practical requirements for the algorithm to work. These include an unbiased estimation of the huge determinant in lattice QCD and, moreover, a Monte Carlo algorithm which accommodates the unbiased estimate of the probability. We shall discuss them in the following sections.

## 4 A Noisy Monte Carlo Algorithm

There are problems in physics which involve extensive quantities such as the fermion determinant which require $`V^3`$ steps to compute exactly. Problems of this kind with large volumes are not numerically applicable with the usual Monte Carlo algorithm which require an exact evaluation of the probability ratios in the accept/reject step. To address this problem, Kennedy and Kuti  proposed a Monte Carlo algorithm which admits stochastically estimated transition probabilities as long as they are unbiased. But there is a drawback. The probability could lie outside the interval between $`0`$ and $`1`$ since it is estimated stochastically. This probability bound violation will destroy detailed balance and lead to systematic bias. To control the probability violation with a large noise ensemble can be costly.

We propose a noisy Monte Carlo algorithm which avoids this difficulty with two Metropolis accept/reject steps. Let us consider a model with Hamiltonian $`H(U)`$ where $`U`$ collectively denotes the dynamical variables of the system. The major ingredient of the new approach is to transform the noise for the stochastic estimator into stochastic variables. The partition function of the model can be written as

$`Z`$ $`=`$ $`{\displaystyle [DU]e^{H(U)}}`$ (11)
$`=`$ $`{\displaystyle [DU][D\xi ]P_\xi (\xi )f(U,\xi )}.`$
where $`f(U,\xi )`$ is an unbiased estimator of $`e^{H(U)}`$ from the stochastic variable $`\xi `$ and $`P_\xi `$ is the probability distribution for $`\xi `$.

The next step is to address the lower probability-bound violation. One first notes that one can write the expectation value of the observable $`O`$ as

$$O=[DU][D\xi ]P_\xi (\xi )O(U)\mathrm{sign}(f)|f(U,\xi )|/Z,$$
(12)
where $`sign(f)`$ is the sign of the function $`f`$. After redefining the partition function to be

$$Z=[DU][D\xi ]P_\xi (\xi )|f(U,\xi )|,$$
(13)
which is semi-positive definite, the expectation of $`O`$ in Eq. (12) can be rewritten as

$$O=O(U)\mathrm{sign}(f)/\mathrm{sign}(f).$$
(14)
As we see, the sign of $`f(U,\xi )`$ is not a part of the probability any more but a part in the observable. Notice that this reinterpretation is possible because the sign of $`f(U,\xi )`$ is a state function which depends on the configuration of $`U`$ and $`\xi `$.

It is clear then, to avoid the problem of lower probability-bound violation, the accept/reject criterion has to be factorizable into a ratio of the new and old probabilities so that the sign of the estimated $`f(U,\xi )`$ can be absorbed into the observable. This leads us to the Metropolis accept/reject criterion which incidentally cures the problem of upper probability-bound violation at the same time. It turns out two accept/reject steps are needed in general. The first one is to propose updating of $`U`$ via some procedure while keeping the stochastic variables $`\xi `$ fixed. The acceptance probability $`P_a`$ is

$$P_a(U_1,\xi U_2,\xi )=\mathrm{min}(1,\frac{|f(U_2,\xi )|}{|f(U_1,\xi )|}).$$
(15)
The second accept/reject step involves the refreshing of the stochastic variables $`\xi `$ according to the probability distribution $`P_\xi (\xi )`$ while keeping $`U`$ fixed. The acceptance probability is

$$P_a(U,\xi _1U,\xi _2)=\mathrm{min}(1,\frac{|f(U,\xi _2)|}{|f(U,\xi _1)|}).$$
(16)
It is obvious that there is neither lower nor upper probability-bound violation in either of these two Metropolis accept/reject steps. Furthermore, it involves the ratios of separate state functions so that the sign of the stochastically estimated probability $`f(U,\xi )`$ can be absorbed into the observable as in Eq. (14).

Detailed balance can be proven to be satisfied and it is unbiased . Therefore, this is an exact algorithm.

## 5 Noisy Monte Carlo with Fermion Determinant

One immediate application of NMC is lattice QCD with dynamical fermions. The action is composed of two parts – the pure gauge action $`S_g(U)`$ and a fermion action $`S_F(U)=Tr\mathrm{ln}M(U)`$. Both are functionals of the gauge link variables $`U`$.

To find out the explicit form of $`f(U,\xi )`$, we note that the fermion determinant can be calculated stochastically as a random walk process

$$e^{Tr\mathrm{ln}M}=1+Tr\mathrm{ln}M(1+\frac{Tr\mathrm{ln}M}{2}(1+\frac{Tr\mathrm{ln}M}{3}(\mathrm{}))).$$
(17)
This can be expressed in the following integral

$`e^{Tr\mathrm{ln}M}={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}d\eta _iP_\eta (\eta _i)_0^1\underset{n=2}{\overset{\mathrm{}}{}}d\rho _n}`$
$`[1+\eta _1^{}\mathrm{ln}M\eta _1(1+\theta (\rho _2{\displaystyle \frac{1}{2}})\eta _2^{}\mathrm{ln}M\eta _2(1+\theta (\rho _3{\displaystyle \frac{2}{3}})\eta _3^{}\mathrm{ln}M\eta _3(\mathrm{}],`$ (18)
where $`P_\eta (\eta _i)`$ is the probability distribution for the stochastic variable $`\eta _i`$. It can be the Gaussian noise or the $`Z_2`$ noise ($`P_\eta (\eta _i)=\delta (|\eta _i|1)`$ in this case). The latter is preferred since it has the minimum variance . $`\rho _n`$ is a stochastic variable with uniform distribution between 0 and 1. This sequence terminates stochastically in finite time and only the seeds from the pseudo-random number generator need to be stored in practice. The function $`f(U,\eta ,\rho )`$ ( $`\xi `$ in Eq. (11) is represented by two stochastic variables $`\eta `$ and $`\rho `$ here) is represented by the part of the integrand between the the square brackets in Eq. (5). One can then use the efficient Padé-Z2 algorithm  to calculate the $`\eta _i\mathrm{ln}M\eta _i`$ in Eq. (5). We shall discuss this in the next section.

Finally, there is a practical concern that $`Tr\mathrm{ln}M`$ can be large so that it takes a large statistics to have a reliable estimate of $`e^{Tr\mathrm{ln}M}`$ from the series expansion in Eq. (5). In general, for the Taylor expansion $`e^x=x^n/n!`$, the series will start to converge when $`x^n/n!>x^{n+1}/(n+1)!`$. This happens at $`n=x`$. For the case $`x=100`$, this implies that one needs to have more than 100! stochastic configurations in the Monte Carlo integration in Eq. (5) in order to have a convergent estimate. Even then, the error bar will be very large. To avoid this difficulty, one can implement the following strategy. First, one notes that since the Metropolis accept/reject involves the ratio of exponentials, one can subtract a universal number $`x_0`$ from the exponent $`x`$ in the Taylor expansion without affecting the ratio. Second, one can use a specific form of the exponential to diminish the value of the exponent. In other words, one can replace $`e^x`$ with $`(e^{(xx_0)/N})^N`$ to satisfy $`|xx_0|/N<1`$. The best choice for $`x_0`$ is $`\overline{x}`$, the mean of $`x`$. In this case, the variance of $`e^x`$ becomes $`e^{\delta ^2/N}1`$.

## 6 The Padé – Z2 Method of Estimating Determinants

Now we shall discuss a very efficient way of estimating the fermion determinant stochastically .

### 6.1 Padé approximation

The starting point for the method is the Padé approximation of the logarithm function. The Padé approximant to $`\mathrm{log}(z)`$ of order $`[K,K]`$ at $`z_0`$ is a rational function $`N(z)/D(z)`$ where deg $`N(z)`$ = deg $`D(z)`$ = $`K`$, whose value and first $`2K`$ derivatives agree with $`\mathrm{log}z`$ at the specified point $`z_0`$. When the Padé approximant $`N(z)/D(z)`$ is expressed in partial fractions, we obtain

$$\mathrm{log}zb_0+\underset{k=1}{\overset{K}{}}\left(\frac{b_k}{z+c_k}\right),$$
(19)
whence it follows

$$\text{log det }𝐌=\text{Tr log}𝐌b_0Tr𝐈+\underset{k=1}{\overset{K}{}}b_kTr(𝐌+c_k𝐈)^1.$$
(20)
The Padé approximation is not limited to the real axis. As long as the function is in the analytic domain, i. e. away from the cut of the log, say along the negative real axis, the Padé approximation can be made arbitrarily accurate by going to a higher order $`[K,K]`$ and a judicious expansion point to cover the eigenvalue domain of the problem.

### 6.2 Complex Z2 noise trace estimation

Exact computation of the trace inverse for $`N\times N`$ matrices is very time consuming for matrices of size $`N10^6`$. However, the complex Z2 noise method has been shown to provide an efficient stochastic estimation of the trace . In fact, it has been proved to be an optimal choice for the noise, producing a minimum variance .

The complex Z2 noise estimator can be briefly described as follows . We construct L noise vectors $`\eta ^1,\eta ^2,\mathrm{},\eta ^L`$ where $`\eta ^j=\{\eta _1^j,\eta _2^j,\eta _3^j,\mathrm{},\eta _N^j\}^T`$, as follows. Each element $`\eta _n^j`$ takes one of the four values $`\{\pm 1,\pm ı\}`$ chosen independently with equal probability. It follows from the statistics of $`\eta _n^j`$ that

$`E[<\eta _n>]E[{\displaystyle \frac{1}{L}}{\displaystyle \underset{j=1}{\overset{L}{}}}\eta _n^j]=0,E[<\eta _m^{}\eta _n>]E[{\displaystyle \frac{1}{L}}{\displaystyle \underset{j=1}{\overset{L}{}}}\eta _m^j\eta _n^j]=\delta _{mn}.`$ (21)
The vectors can be used to construct an unbiased estimator for the trace inverse of a given matrix $`M`$ as follows:

$`E[<\eta ^{}𝐌^1\eta >]`$ $``$ $`E[{\displaystyle \frac{1}{L}}{\displaystyle \underset{j=1}{\overset{L}{}}}{\displaystyle \underset{m,n=1}{\overset{N}{}}}\eta _m^jM_{m,n}^1\eta _n^j]`$
$`=`$ $`{\displaystyle \underset{n}{\overset{N}{}}}M_{n,n}^1+({\displaystyle \underset{mn}{\overset{N}{}}}M_{m,n}^1)[{\displaystyle \frac{1}{L}}{\displaystyle \underset{j}{\overset{L}{}}}\eta _m^j\eta _n^j]`$
$`=`$ $`\text{Tr }𝐌^1.`$
The variance of the estimator is shown to be

$`\sigma _M^2`$ $``$ $`\mathrm{Var}[<\eta ^{}𝐌^1\eta >]=\mathrm{E}\left[|<\eta ^{}𝐌^1\eta >\text{Tr }𝐌^1|^2\right]`$
$`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{mn}{\overset{N}{}}}M_{m,n}^1(M_{m,n}^1)^{}={\displaystyle \frac{1}{L}}{\displaystyle \underset{mn}{\overset{N}{}}}|M_{m,n}^1|^2.`$
The stochastic error of the complex Z2 noise estimate results only from the off-diagonal entries of the inverse matrix (the same is true for Zn noise for any n). However, other noises (such as Gaussian) have additional errors arising from diagonal entries. This is why the Z2 noise has minimum variance. For example, it has been demonstrated on a $`16^3\times 24`$ lattice with $`\beta =6.0`$ and $`\kappa =0.148`$ for the Wilson action that the Z2 noise standard deviation is smaller than that of the Gaussian noise by a factor of 1.54 .

Applying the complex Z2 estimator to the expression for the $`Trlog𝐌`$ in Eq. (20), we find

$`{\displaystyle \underset{k}{}}b_k\text{Tr}(M+c_k)^1`$
$`{\displaystyle \frac{1}{L}}{\displaystyle \underset{k}{\overset{K}{}}}{\displaystyle \underset{j}{\overset{L}{}}}b_k\eta ^j(M+c_k)^1\eta ^j`$
$`={\displaystyle \frac{1}{L}}{\displaystyle \underset{j}{\overset{L}{}}}{\displaystyle \underset{k=1}{\overset{K}{}}}b_k\eta ^j\xi ^{k,j},`$ (22)
where $`\xi ^{k,j}=(M+c_k𝐈)^1\eta ^j`$ are the solutions of

$$(M+c_k𝐈)\xi ^{k,j}=\eta ^j,$$
(23)
Since $`M+c_k𝐈`$ are shifted matrices with constant diagonal matrix elements, Eq. (23) can be solved collectively for all values of $`c_k`$ within one iterative process by several algorithms, including the Quasi-Minimum Residual (QMR) , Multiple-Mass Minimum Residual (M3 R) , and GMRES. We have adopted the M3 R algorithm, which has been shown to be about 2 times faster than the conjugate gradient algorithm, and the overhead for the multiple $`c_k`$ is only 8% . The only price to pay is memory: for each $`c_k`$, a vector of the solution needs to be stored. Furthermore, one observes that $`c_k>0`$. This improves the conditioning of $`(𝐌+c_k𝐈)`$ since the eigenvalues of $`𝐌`$ have positive real parts. Hence, we expect faster convergence for column inversions for Eq. (23).

In the next section, we describe a method which significantly reduces the stochastic error.

### 6.3 Improved PZ estimation with unbiased subtraction

In order to reduce the variance of the estimate, we introduce a suitably chosen set of traceless $`N\times N`$ matrices $`𝐐^{(p)}`$, i.e. which satisfy $`_{n=1}^N𝐐_{n,n}^{(p)}=0,p=1\mathrm{}P`$. The expected value and variance for the modified estimator $`<\eta ^{}(𝐌^1_{p=1}^P\lambda _p𝐐^{(p)})\eta >`$ are given by

$`E[<\eta ^{}(𝐌^1{\displaystyle \underset{p=1}{\overset{P}{}}}\lambda _p𝐐^{(p)}))\eta >]`$ $`=`$ $`\text{Tr }𝐌^1,`$ (24)
$`\mathrm{\Delta }_M(\lambda )=\mathrm{Var}[<\eta ^{}(𝐌^1{\displaystyle \underset{\mathrm{p}=1}{\overset{\mathrm{P}}{}}}\lambda _\mathrm{p}𝐐^{(\mathrm{p})})\eta >]`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{mn}{}}|𝐌_{m,n}^1{\displaystyle \underset{p=1}{\overset{P}{}}}\lambda _p𝐐_{m,n}^{(p)})|^2,`$ (25)
for any values of the real parameters $`\lambda _p`$. In other words, introducing the matrices $`𝐐^{(p)}`$ into the estimator produces no bias, but may reduce the error bars if the $`𝐐^{(p)}`$ are chosen judiciously. Further, $`\lambda _p`$ may be varied at will to achieve a minimum variance estimate: this corresponds to a least-squares fit to the function $`\eta ^{}𝐌^1\eta `$ sampled at points $`\eta _j,j=1\mathrm{}L`$, using the fitting functions $`\{1,\eta ^{}𝐐^{(p)}\eta \},p=1\mathrm{}P`$.

We now turn to the question of choosing suitable traceless matrices $`𝐐^{(p)}`$ to use in the modified estimator. One possibility for the Wilson fermion matrix $`𝐌=𝐈\kappa 𝐃`$ is suggested by the hopping parameter — $`\kappa `$ expansion of the inverse matrix,

$`(𝐌+c_k𝐈)^1`$ $`=`$ $`{\displaystyle \frac{1}{𝐌+c_k𝐈}}={\displaystyle \frac{1}{(1+c_k)(𝐈\frac{\kappa }{(1+c_k)}𝐃)}}`$ (26)
$`=`$ $`{\displaystyle \frac{𝐈}{1+c_k}}+{\displaystyle \frac{\kappa }{(1+c_k)^2}}𝐃+{\displaystyle \frac{\kappa ^2}{(1+c_k)^3}}𝐃^2+{\displaystyle \frac{\kappa ^3}{(1+c_k)^4}}𝐃^3+\mathrm{}.`$
This suggests choosing the matrices $`𝐐^{(p)}`$ from among those matrices in the hopping parameter expansion which are traceless:

$`𝐐^{(1)}`$ $`=`$ $`{\displaystyle \frac{\kappa }{(1+c_k)^2}}𝐃,`$
$`𝐐^{(2)}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{(1+c_k)^3}}𝐃^2,`$
$`𝐐^{(3)}`$ $`=`$ $`{\displaystyle \frac{\kappa ^3}{(1+c_k)^4}}𝐃^3,`$
$`𝐐^{(4)}`$ $`=`$ $`{\displaystyle \frac{\kappa ^4}{(1+c_k)^5}}(𝐃^4\mathrm{Tr}𝐃^4),`$
$`𝐐^{(5)}`$ $`=`$ $`{\displaystyle \frac{\kappa ^5}{(1+c_k)^6}}𝐃^5,`$
$`𝐐^{(6)}`$ $`=`$ $`{\displaystyle \frac{\kappa ^6}{(1+c_k)^7}}(𝐃^6\mathrm{Tr}𝐃^6),`$
$`𝐐^{(2r+1)}`$ $`=`$ $`{\displaystyle \frac{\kappa ^{2r+1}}{(1+c_k)^{2r+2}}}𝐃^{2r+1},r=3,4,5,\mathrm{}.`$
It may be verified that all of these matrices are traceless. In principle, one can include all the even powers which entails the explicit calculation of all the allowed loops in $`TrD^{2r}`$. In this manuscript we have only included $`𝐐^{(4)}`$, $`𝐐^{(6)}`$, and $`𝐐^{(2r+1)}`$.

### 6.4 Computation of $`TrlogM`$

Our numerical computations were carried out with the Wilson action on the $`8^3\times 12`$ ($`N`$ = 73728) lattice with $`\beta =5.6`$. We use the HMC with pseudofermions to generate gauge configurations. With a cold start, we obtain the fermion matrix $`𝐌_1`$ after the plaquette becomes stable. The trajectories are traced with $`\tau =0.01`$ and 30 molecular dynamics steps using $`\kappa =0.150`$. $`𝐌_2`$ is then obtained from $`𝐌_1`$ by an accepted trajectory run. Hence $`𝐌_1`$ and $`𝐌_2`$ differ by a continuum perturbation, and $`\mathrm{log}[det𝐌_1/det𝐌_2]𝒪(1)`$.

We first calculate $`\mathrm{log}det𝐌_1`$ with different orders of Padé expansion around $`z_0=0.1`$ and $`z_0=1.0`$. We see from Table 1 that the 5th order Padé does not give the same answer for two different expansion points, suggesting that its accuracy is not sufficient for the range of eigenvalues of $`𝐌_1`$. Whereas, the 11th order Padé gives the same answer within errors. Thus, we shall choose P$`(z)`$ with $`z_0=0.1`$ to perform the calculations from this point on.

In Table 2, we give the results of improved estimations for $`\text{Tr}\mathrm{log}𝐌_1`$. We see that the variational technique described above can reduce the data fluctuations by more than an order of magnitude. For example, the unimproved error $`\delta _0=5.54`$ in Table 2 for 400 Z2 noises is reduced to $`\delta _{11}=0.15`$ for the subtraction which includes up to the $`𝐐^{11}`$ matrix. This is 37 times smaller. Comparing the central values in the last row (i.e. the $`11^{th}`$ order improved) with that of unimproved estimate with 10,000 Z2 noises, we see that they are the same within errors. This verifies that the variational subtraction scheme that we employed does not introduce biased errors. The improved estimates of $`\text{Tr}\mathrm{log}𝐌_1`$ from 50 Z2 noises with errors $`\delta _r`$ from Table 2 are plotted in comparison with the central value of the unimproved estimate from 10,000 noises in Fig. 3.

## 7 Implementation of the Kentucky Noisy Monte Carlo Algorithm

We have recently implemented the above Noisy Monte Carlo algorithm to the simulation of lattice QCD with dynamical fermions by incorporating the full determinant directly . Our algorithm uses pure gauge field updating with a shifted gauge coupling to minimize fluctuations in the trace log is the Wilson Dirac matrix. It gives the correct results as compared to the standard Hybrid Monte Carlo simulation. However, the present simulation has a low acceptance rate due to the pure gauge update and results in long autocorrelations. We are in the process of working out an alternative updating scheme with molecular dynamics trajectory to include the feedback of the determinantal effects on the gauge field which should be more efficient than the pure gauge update.

## 8 Summary

After reviewing the finite density algorithm for QCD with the chemical potential, we propose a canonical approach by projecting out the definite baryon number sector from the fermion determinant and stay in the sector throughout the Monte Carlo updating. This should circumvent the overlap problem. In order to make the algorithm practical, one needs an efficient way to estimate the huge fermion determinant and a Monte Carlo algorithm which admits unbiased estimates of the probability without upper unitarity bound violations. These are achieved with the Padé-Z2 estimate of the determinant and the Noisy Monte Carlo algorithm. So far, we have implemented the Kentucky Noisy Monte Carlo algorithm to incorporate dynamical fermions in QCD on a relatively small lattice and medium heavy quark based on pure gauge updating. As a next step, we will work on a more efficient updating algorithm and project out the baryon sector to see if the finite density algorithm proposed here will live up to its promise.

## Acknowledgments

This work is partially supported by the U.S. DOE grant DE-FG05-84ER40154. The author wishes to thank M. Faber for introducing the subject of the finite density to him and Shao-Jing Dong for providing some unpublished figures. Fruitful discussions with M. Alford, I.Barbour, U.-J. Wiese, and R. Sugar are acknowledged. He would also thank the organizer, Ge Mo-lin for the invitation to attend the conference and his hospitality.