text:
# Fano’s inequality is false for a simple Cremona transformation of five-dimensional projective space

## Introduction

It was shown in  and  that Fano’s inequalities are false for tree-dimensional projective space and for three-dimensional quadric. Both the birational transformations were of degree 13. Here we construct a simple five-dimensional Cremona transformation of degree 7 and not satisfying the inequalities of Fano .

Let $`X`$ be a non-singular $`n`$-dimensional variety such that Pic$`(X)=,`$ and the anticanonical class $`(K_X)`$ is ample. If $`(K_X)=rH`$ for a generator $`H`$ of the Picard group, then $`X`$ is said to be Fano variety of index $`r`$ ( and of the first kind ). According to the texts mentioned in survey , Fano’s inequality is the statement:

> For any birational transformation,
>
> $$f:XX$$
>
> defined by a linear system of degree $`d>1`$ (the degree is the number defined by $`f(mH)=dmH`$ ) there exists an irreducible subvariety $`YX`$ such that
>
> $$0dimYdimX2,$$
>
> $$mult_Y(f(|H|))>\frac{(dimXdimY1)d}{r}.$$

Five-dimensional projective space $`^5`$ has index 6, therefore one can reformulate Fano’s inequality:

> For any Cremona transformation,
>
> $$f:^5^5$$
>
> defined by six homogeneous polynomials of the same degree $`d>1`$ and without a common nonconstant factor,
>
> $$x_i^{}=f_i(x_0,x_1,x_2,x_3,x_4,x_5),i=0,\mathrm{},5,$$
>
> there exists an irreducible subvariety $`Y^5`$ such that
>
> $$0dimY3,$$
>
> $$mult_Y(f_i)>\frac{(4dimY)d}{6}.$$
>
> for every $`i=0,\mathrm{},5`$.

The goal of my article is to show that these inequalities do not take place for a certain Cremona transformation of degree 7. That is, I write down the formulas for a Cremona transformation of degree 7 such that for the forms $`f_0,\mathrm{},f_5`$ defining the transformation, for any surface $`S^5,`$ for any curve $`C^5`$ and for any point $`P^5,`$ one can see that

$$\underset{i}{\mathrm{min}}(mult_S(f_i))2,\underset{i}{\mathrm{min}}(mult_C(f_i))3,\underset{i}{\mathrm{min}}(mult_P(f_i))4.$$

## Construction of the example

Let us consider six homogeneous coordinates for $`^5`$ as the normalized coefficients $`x_{00},x_{01},x_{11},x_{02},x_{12},x_{22}`$ of a ternary quadratic form $`F(T_0,T_1,T_2),`$

$$F(T_0,T_1,T_2)=x_{00}T_0^2+2x_{01}T_0T_1+x_{11}T_1^2$$

$$+2x_{02}T_0T_2+2x_{12}T_1T_2+x_{22}T_2^2.$$

Let

$$D=D(x_{00},x_{01},x_{11},x_{02},x_{12},x_{22})$$

be the discriminant of the ternary quadric,

$$D=\left|\begin{array}{ccc}x_{00}& x_{01}& x_{02}\\ x_{01}& x_{11}& x_{12}\\ x_{02}& x_{12}& x_{22}\end{array}\right|.$$

The set of double points of the cubic discriminant hypersurface consists of the points of the Veronese surface , these points correspond to the ternary quadrics which are perfect squares of linear forms.

We fix parameters $`s,t`$ and consider six following forms of degree 7.

$`(f_{s,t})_{00}`$ $`=`$ $`x_{00}D^2,`$
$`(f_{s,t})_{01}`$ $`=`$ $`x_{01}D^2,`$
$`(f_{s,t})_{11}`$ $`=`$ $`x_{11}D^2,`$
$`(f_{s,t})_{02}`$ $`=`$ $`x_{02}D^2+x_{01}x_{11}^3Ds+x_{00}^4Dt,`$
$`(f_{s,t})_{12}`$ $`=`$ $`x_{12}D^2+x_{01}x_{00}^3Dt+x_{11}^4Ds,`$
$`(f_{s,t})_{22}`$ $`=`$ $`x_{22}D^2+2x_{12}x_{11}^3Ds+2x_{02}x_{00}^3Dt+x_{11}^7s^2+2x_{01}x_{11}^3x_{00}^3st+x_{00}^7t^2.`$
These six forms define a two-parameter family of rational maps

$$g_{s,t}:^5^5.$$

If $`s=t=0`$, all four forms have a common factor $`D^2`$. After cancelling this, we see that $`g_{0,0}`$ is the identity transformation. For our example, we need nonzero values of $`s,t`$. If one of the parameters $`s,t`$ is not zero, then it is clear that the six forms have no nonconstant common factor. Moreover,

$$D((f_{s,t})_{00},(f_{s,t})_{01},(f_{s,t})_{11},(f_{s,t})_{02},(f_{s,t})_{12},(f_{s,t})_{22})=D^7.$$

In fact, this identity expresses the invariance of the discriminant under triangular transformation of variables $`T_0,T_1,T_2`$. Using the latter identity, it is not hard to see that

$$(f_{s,t})_{ij}((f_{s,t})_{00},(f_{s,t})_{01},(f_{s,t})_{11},(f_{s,t})_{02},(f_{s,t})_{12},(f_{s,t})_{22})=x_{ij}D^{16},$$

that is,

$$g_{s,t}g_{s,t}=\text{the identity transformation}.$$

Thus $`g_{s,t}`$ is rationally invertible and is a Cremona transformation. More generally,

$$g_{s,t}g_{p,q}=g_{s+p,t+q},$$

and we get a two-parameter group of Cremona transformations. These transformations induce biregular automorphisms of an affine open subset of the five-dimensional projective space, the complement of the discriminant cubic hypersurface $`D=0`$. Indeed, the above formula of the discriminant transformation proves this (moreover, one can see below the exact calculation of the fundamental points of such a transformation).

Remark. The formulas for $`g_{s,t}`$ can be generalized for obtaining an infinite dimensional family of automorphisms of the complement to the discriminant hypersurface. The construction of general formulas resembles the trick used on page 8 of the Max-Planck Institute preprint . If $`U(x,y,z)`$ is a ternary cubic form, $`\varphi _m(u,d),\psi _m(u,d)`$ are binary forms of degree $`m`$,

$$\varphi =\varphi _m(U(x_{00},x_{01},x_{11}),D),$$

$$\psi =\psi _m(U(x_{00},x_{01},x_{11}),D),$$

then following transformation

$`x_{00}^{}`$ $`=`$ $`x_{00}D^{2m},`$
$`x_{01}^{}`$ $`=`$ $`x_{01}D^{2m},`$
$`x_{11}^{}`$ $`=`$ $`x_{11}D^{2m},`$
$`x_{02}^{}`$ $`=`$ $`x_{02}D^{2m}+x_{01}D^m\psi +x_{00}D^m\varphi ,`$
$`x_{12}^{}`$ $`=`$ $`x_{12}D^{2m}+x_{01}D^m\varphi +x_{11}D^m\varphi ,`$
$`x_{22}^{}`$ $`=`$ $`x_{22}D^{2m}+2x_{12}D^m\varphi +2x_{02}D^m\psi +x_{11}\varphi ^2+2x_{01}\varphi \psi +x_{00}\psi ^2.`$
is a Cremona transformation inducing a biregular automorphism of the mentioned affine subset.

Let us return to our two-parameter family. We fix nonzero values of the parameters $`s,t`$, for example, $`s=t=1`$, and consider the corresponding Cremona transformation

$`x_{00}^{}`$ $`=`$ $`x_{00}D^2,`$
$`x_{01}^{}`$ $`=`$ $`x_{01}D^2,`$
$`x_{11}^{}`$ $`=`$ $`x_{11}D^2,`$
$`x_{02}^{}`$ $`=`$ $`x_{02}D^2+(x_{01}x_{11}^3+x_{00}^4)D,`$
$`x_{12}^{}`$ $`=`$ $`x_{12}D^2+(x_{01}x_{00}^3+x_{11}^4)D,`$
$`x_{22}^{}`$ $`=`$ $`x_{22}D^2+2(x_{12}x_{11}^3+x_{02}x_{00}^3)D+(x_{11}^7+2x_{01}x_{11}^3x_{00}^3+x_{00}^7).`$
First of all, we find the points $`P^5`$ where each form on the right hand side has positive multiplicity (that is, the set of all common zeros of these right hand sides, or, in other words, the fundamental points of the transformation). The first right hand side vanishes if

* either $`x_{00}=0`$,
* or $`D=0`$,
* or simultaneously $`x_{00}=D=0`$.

If $`x_{00}=0`$ but $`D0`$ then using the other five formulas, one sees that for other five coordinates of a fundamental point, the equalities $`x_{01}=x_{11}=x_{02}=x_{12}=x_{22}=0`$ take place. This case is not a point of $`^5`$.

The case $`D=0`$ but $`x_{00}0`$ is realizable. The points satisfying

$$x_{00}0D=0,x_{11}^7+2x_{01}x_{11}^3x_{00}^3+x_{00}^7=0$$

are fundamental, but they are of multiplicity 1 for one of three forms $`x_{02}^{},x_{12}^{},x_{22}^{}`$ at least. Indeed, if such a point is double on the discriminant hypersurface, then then it is non-singular on the hypersurface of seventh degree. If both the expressions $`(x_{01}x_{11}^3+x_{00}^4)`$ , $`(x_{01}x_{00}^3+x_{11}^4)`$ in $`x_{02}^{},x_{12}^{}`$ vanish , and the point is fundamental , then the multiplicity of $`x_{22}^{}`$ at the point is equal to 1.

If $`x_{00}=D=0`$ and the point is fundamental, then also $`x_{11}=0`$ and either $`x_{01}=0`$ or $`(2x_{02}x_{12}x_{01}x_{22})=0,`$ or both the expressions vanish . If the point is simple on the discriminant cubic and has multiplicity more than 2 on every $`x_{ij}^{},`$ then all the homogeneous coordinates vanish. If the point is double on $`D=0`$, then $`x_{01}=x_{02}=x_{12}=0,`$ that is the only non-zero coordinate of the point is $`x_{22}.`$ The multiplicity of $`x_{22}^{}`$ at the point is 4.


format: text

text:
# Spectral fluctuations effects on conductance peak height statistics in quantum dots

## Abstract

Within random matrix theory for quantum dots, both the dot’s one-particle eigenlevels and the dot-lead couplings are statistically distributed. While the effect of the latter on the conductance is obvious and has been taken into account in the literature, the statistical distribution of the one-particle eigenlevels is generally replaced by a picket-fence spectrum. Here we take the random matrix theory eigenlevel distribution explicitly into account and observe significant deviations in the conductance distribution and magnetoconductance of closed quantum dots at experimentally relevant temperatures.

The universal statistical fluctuations observed at the low-energy part of the spectrum of quantum systems whose associated classical dynamics are chaotic can be described by random matrix theory (RMT). This type of description can be justified for diffusive quantum dots and quantum dots with irregular shapes Alhassid01  which makes quantum dots a particular example for the study of RMT fluctuations. While in open quantum dots (which have a strong dot-lead coupling) the electron-electron interaction effects are mostly neglected, in closed quantum dots this interaction leads to the Coulomb blockade (see CBreview  for a review): The low-temperature conductance is heavily suppressed due to the large charging energy for adding an electron to the quantum dot, except for the Coulomb blockade peaks at which the potential of the quantum dot is adjusted such that $`N`$ and $`N+1`$ electrons in the dot have the same energy. The RMT approach Jalabert92 ; Alhassid01  successfully described the mesoscopic fluctuations of these Coulomb blockade peaks, i.e., the statistical distribution of their height $`P(G^{\mathrm{max}})`$ and its dependence upon magnetic field ChangFolk96 . On the other hand, recent improved experiments show significant deviations from the RMT prediction, suggesting that interaction effects beyond charging should be considered as well. In particular, dephasing of the single-particle states due to interactions modifies the conductance peak height statistics (see Alhassid01 ; Aleiner01  and references therein). In a recent experiment, Patel et al. Patel98  found that the statistical distribution has a smaller ratio of standard deviation to mean peak height $`\sigma (G^{\mathrm{max}})/G^{\mathrm{max}}`$ than predicted by RMT Alhassid98  and attributed this to dephasing effects. In another experiment, Folk et al. Folk00  measured the change of the conductance in a magnetic field $`B`$

$$\alpha =\frac{G^{\mathrm{max}}_{B0}G^{\mathrm{max}}_{B=0}}{G^{\mathrm{max}}_{B0}},$$
(1)
as a probe of dephasing times. This is the closed dot analog of the weak localization magnetoconductance which had proven to be an effective measure for open dots dephasing times Huibers98 . It was pointed out that $`\alpha =1/4`$ as long as the transport is dominated by elastic scattering Alhassid98 ; Beenakker00 . Therefore, any deviation of the measured $`\alpha `$ from $`1/4`$ was considered an indication for dephasing.

In this paper we discuss the effects of spectral fluctuations of the RMT one-particle eigenlevels on the statistical distribution $`P(G^{\mathrm{max}})`$ and the weak localization correction $`\alpha `$. Previous works Alhassid01 ; Jalabert92 ; Alhassid98 ; Beenakker00 ; Eisenberg01  have generally considered a picket fence spectrum, i.e., a rigid level spacing between successive eigenlevels in the quantum dot, for the calculation of the conductance. This ignores the effect of spectral eigenlevel fluctuations. The picket fence spectrum is a good approximation for both very high temperatures and very low temperatures Alhassid01 , and a comparison of $`P(G^{\mathrm{max}})`$ with full RMT statistics and a picket fence spectrum without spin-degeneracy at three temperatures showed only minor deviations Alhassid98b .

In the present paper, we study the full RMT statistics in detail with and without spin-degeneracy, and find significant differences compared to the picket fence spectrum, in particular in an experimentally relevant regime $`k_BT\mathrm{\Delta }`$. The spectral fluctuations lead to lower values of $`\alpha `$ than $`1/4`$ such that this value is not universal, even in the absence of any dephasing mechanism. One therefore has to be careful while using $`\alpha `$ as a probe for dephasing in this temperature regime.

Within the constant interaction model, the conductance of a quantum dot is given by the formula Beenakker91

$$G=\frac{e^2}{kT}\underset{p=1}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }_i^L\mathrm{\Gamma }_i^R}{\mathrm{\Gamma }_i^L+\mathrm{\Gamma }_i^R}P_{\mathrm{eq}}(N)P(E_i|N)[1f(E_i\mu )]$$
(2)
where $`\mathrm{\Gamma }_i^{L(R)}`$ is the tunneling rate between the $`i`$th one-particle eigenlevel of the dot and the left (right) lead, $`P_{\mathrm{eq}}(N)`$ is the equilibrium probability to find $`N`$ electrons in the dot with the Coulomb blockade allowing for $`N`$ and $`N+1`$ electrons, $`P(E_i|N)`$ is the canonical probability to have the level $`i`$ occupied given the presence of $`N`$ electrons in the dot, and $`f(E)`$ is the Fermi function at the effective chemical potential $`\mu `$, which includes the charging energy. In a typical experimental situation, the charging energy is much large than temperature, and thus only one term contributes to the sum over $`N`$. In Eq. (2), $`\mathrm{\Gamma }_i^{L(R)}`$ is Porter-Thomas distributed in the Gaussian orthogonal ensemble (GOE) and Gaussian unitary ensemble (GUE) without and with a magnetic field, respectively, and the eigenlevel energies $`E_i`$ obey the respective RMT distribution Alhassid01 . In contrast, the picket fence spectrum has $`E_{2i}=E_{2i1}=i\mathrm{\Delta }`$ in the case of spin-degeneracy and $`E_i=i\mathrm{\Delta }/2`$ without spin-degeneracy. The first term in the sum $`\mathrm{\Gamma }_i^L\mathrm{\Gamma }_i^R/(\mathrm{\Gamma }_i^L+\mathrm{\Gamma }_i^R)`$ depends only on the eigenfunctions of the dot, and thus is uncorrelated with the spectrum within the RMT approach. The ensemble average of this term in the absence (GOE) or presence (GUE) of a magnetic field is

$$\frac{\mathrm{\Gamma }_i^L\mathrm{\Gamma }_i^R}{\mathrm{\Gamma }_i^L+\mathrm{\Gamma }_i^R}=\{\begin{array}{cc}1/4\hfill & \text{GOE}\hfill \\ 1/3\hfill & \text{GUE}\hfill \end{array}.$$
(3)
This yields the value $`\alpha =1/4`$ if the weights $`P(E_i|N)`$ are the same for both ensembles. This should be the case in the low temperature regime $`k_BT\mathrm{\Delta }`$ since only one level $`E_0`$ contributes with maximal weight, $`P(E_i|N)\delta _{i0}`$. In general, the main contribution to the sum comes from $`O(k_BT/\mathrm{\Delta })`$ levels around the Fermi energy which gives the same contribution at large temperatures $`k_BT\mathrm{\Delta }`$ for the GOE and GUE, $`\alpha =1/4`$ in this regime as well.

However, for $`k_BT\mathrm{\Delta }`$, the probability to have more than one level in an energy window $`k_BT`$ around the Fermi energy is increased for the RMT eigenlevel distribution compared to the picket-fence spectrum. These additional levels enhance the conductance. Since there are more close-by levels for the GOE case, due to the weaker level repulsion, the GOE conductance is enhanced more, and $`\alpha `$ is suppressed.

A second important effect is the optimization of the chemical potential for the Coulomb blockade peak. This effect was generally ignored, as it is technically cumbersome to consider, and is not significant for both very low and very high temperatures. Disregarding this effect means that a theorist optimized the chemical potential w.r.t. the averaged conductance, instead of optimizing for every realization as in the experiment. Whenever there is a close-by level, the position of the peak is shifted to optimize the contribution from both levels. Typically, a level with very low tunneling rates (and, thus, suppressed conductance peak) would get enhanced significantly by contributions from its neighbors. If the tunneling rate of a neighboring level is much higher, the peak position $`\mu ^{\mathrm{max}}`$ is shifted towards it. As the distribution of level spacings is different depending on the existence of magnetic field, this enhancement mechanism is again more effective in the absence of magnetic field (GOE), where probabilities of small spacing and of small conductances are higher. Thus, this effect which was neglected in Alhassid98b  suppresses $`\alpha `$ even further.

We evaluated the sum (2) numerically by drawing $`\mathrm{\Gamma }_i^{L(R)}`$ from the Porter-Thomas distribution and $`E_i`$ according to the Wigner-Dyson distribution. Levels within a window of $`\pm 4k_BT`$ around the Fermi energy have been taken into account and the Fermi energy $`\mu `$ in Eq. (2) has been adjusted to yield $`G^{\mathrm{max}}`$ for every realization.

Figure 1 compares the probability distribution $`P(G^{\mathrm{max}})`$ for a picket fence spectrum vs. the full RMT level statistics. As explained above, RMT spectral fluctuations enhance the conductance. In particular, the probability to have a very low $`G^{\mathrm{max}}`$ is reduced and the probability to have an intermediate $`G^{\mathrm{max}}`$ is enhanced. The reason for the reduction is that a very low $`G^{\mathrm{max}}`$ requires $`\mathrm{\Gamma }^L`$ or $`\mathrm{\Gamma }^R`$ in Eq. (2) to be low. RMT spectral fluctuations enhance the contributions from close-by levels, which typically do not have a low value of $`\mathrm{\Gamma }^{L(R)}`$ at the same time. Thus, the peak position of $`\mu `$ is shifted towards a close-by level and the conductance occurs through both levels. Notably, the effect of phase-breaking inelastic scattering processes leads to similar changes Eisenberg01 .

Deviations of $`\alpha `$ from the “universal” value $`1/4`$ have been interpreted as being a result of dephasing. While dephasing would certainly suppress $`\alpha `$, we note here that in the regime $`k_BT\mathrm{\Delta }`$, the spectral fluctuation effects discussed above, lead to a similar effect. In Figure 2 we present the results for $`\alpha `$ as a function of the scaled temperature $`k_BT/\mathrm{\Delta }`$, for both spin-degenerate spectrum and the case of broken symmetry. While the effect seems to be small, one should keep in mind that in the low temperature regime, even very strong dephasing does not suppress $`\alpha `$ to zero Beenakker00 , and thus the correction due to spectral fluctuations is comparable with or even larger than the effect of dephasing Beenakker00 ; Eisenberg01 . One should therefore cautiously use $`\alpha `$ as a probe of dephasing in this regime.

In conclusion, we have shown that RMT spectral fluctuation effect the probability distribution function $`P(G^{\mathrm{max}})`$, leading to non-negligible deviations in measurable quantities in the regime $`0.1\mathrm{\Delta }<k_BT<0.8\mathrm{\Delta }`$. In particular, the weak localization correction $`\alpha `$ which was recently used as a probe of dephasing in closed quantum dots is affected: $`\alpha `$ is different from $`1/4`$ and, moreover, turns out to be temperature dependent, even in the absence of dephasing. At low temperatures, $`\alpha `$ is reduced down to $`\alpha 0.2`$, which can be below the lower limit of a picket fence model with dephasing. This should also be taken into account while analyzing the ongoing experiments aimed at measuring dephasing times in closed dots in the low temperature regime $`k_BT\mathrm{\Delta }`$.

Finally, we would like to note that during the completion of the present paper some of the results have been independently arrived at in Alhassid02 .


text:
# Superoutbursts, superhumps and the tidal-thermal instability model

## 1 Introduction

Dwarf novae are cataclysmic variables (CVs) which exhibit 2-6 mag recurrent outbursts (see Warner w95 (1995) for an encyclopaedic review). It is generally accepted that these outbursts are due to a thermal-viscous instability that occurs in an accretion disc in which the viscosity is given by the $`\alpha `$-prescription (Shakura & Sunyaev s73 (1973)). The instability arises when hydrogen becomes partially ionized and the opacities vary steeply with temperature; this occurs when the disc material reaches temperatures of order of 8000 K (see Lasota l01 (2001) for a review of the model). A number of characteristics of many light curves (characteristic time scales, amplitudes, …) are reproduced by the disc instability model (DIM) provided the viscosity parameter $`\alpha `$ is different in quiescence and in outburst. However several subclasses of dwarf novae show specific features that are not easy to explain with the standard version of the DIM.

SU UMa stars are defined as dwarf novae in which normal outbursts cycles are regularly interrupted by anomalously bright and long outbursts called superoutbursts. These outbursts are 0.7 magnitudes brighter and last 5 to 10 times longer than the normal ones. The mean time $`T_s`$ between two consecutive superoutbursts defines the period of a supercycle $`T_s`$; in most SU UMa stars, $`130<T_s<350`$ days. There are however two extreme exceptions: ER UMa stars have anomalously short supercycles (between 19 and 45 days), while WZ Sge stars show only superoutbursts separated by very long quiescence periods that can reach 30 years. In addition, the light curve is modulated during superoutbursts at a period slightly longer than the orbital period. These so-called superhumps usually disappear shortly after the decline of the superoutburst. Whitehurst (w88 (1988)) explained these superhumps by the precession of a distorted disc. The name ‘superoutburst’ is reserved for dwarf nova eruptions in which a superhump is present. This property differentiates superoutbursts from wide outbursts observed in many dwarf novae; van Paradijs (vanP83 (1983)) suggested that superoutbursts are just wider outbursts in the bimodal outburst-length distribution. In addition, a 1985 outburst of U Gem had all the properties of a superoutburst (of WZ Sge stars) except for the superhump (see Kuulkers et al. ku99 (1999) and references therein). This could suggest that the superoutburst phenomenon is independent of the disc distortion. We discuss this possibility in the Section 4.

As almost all SU UMa stars have periods below the period gap of CVs, they all have low secondary to primary mass ratio $`q`$, and the primary Roche lobe is large. The disc can therefore extend to relatively large radii, and this led Osaki (o89 (1989)) to propose the thermal-tidal instability (TTI), in which superoutbursts are caused by the disc reaching the 3:1 resonance radius, as a model for SU UMa stars. SPH simulations by Whitehurst (w88 (1988)) showed that the disc becomes eccentric and precesses when its radius reaches the 3:1 resonance radius, but it must be kept in mind that the hydrodynamical simulations of Stehle (s99 (1999)) and Kornet & Różyczka (kr00 (2000)) did not show such an effect.

Detailed observations of SU UMa stars point out several limitations of the TTI model. The mass transfer rate from the secondary is observed to be enhanced during outburst; this increase, possibly due to irradiation, is not taken into account, although it can drastically alter the light curves (Hameury et al. hlw00 (2000)).

The TTI model cannot easily reproduce the properties of ER UMa and WZ Sge subclasses. In order to obtain the very short ER UMa supercycles (19 to 45 days), one must assume that the tidal instability shuts off prematurely (Osaki 1995a ) for very small mass ratios. The very long quiescence time of WZ Sge systems would on the other hand require a very small viscosity in quiescence (Smak s93 (1993)).

Occasionally, some systems show echo outbursts at the end of a superoutburst: several consecutive small outbursts with very short recurrence time are triggered before the disc returns to quiescence. During these echo-outbursts, superhumps are still present. These so-called late superhumps (LS) have also been observed in ER UMa stars. However, the TTI model simulations predict that the end of a superoutburst should coincide with the cessation of the tidal instability.

These different problems have lead several authors to investigate alternative possibilities. Hameury et al. (hlw00 (2000)) have shown that if one takes into account the disc and secondary illumination in the standard DIM, one is able to reproduce light curves reminiscent of those of systems such as RZ LMi or EG Cnc. Recently Hellier (h01 (2001)) suggested that some difficulties could be solved in the framework of a slightly modified TTI model; he proposed that the superoutburst must end before the end of eccentricity in low mass ratio systems, thus giving an explanation for LS.

We have included in our DIM model (Hameury et al. hmd98 (1998); Buat-Ménard et al. b01 (2001), hereafter Paper I) an enhanced tidal torque when the disc radius exceeds the 3:1 resonance radius. We do reproduce the earlier results of Osaki (o89 (1989)), obtained using a simplified model. We show that in a number of cases, a cooling wave can end a superoutburst while the disc is still eccentric, thereby confirming the proposal of Hellier (h01 (2001)). However, this by itself is not sufficient to reproduce the properties of ER UMa or WZ Sge systems, and additional effects such as illumination and mass transfer variations must be included.

## 2 The TTI model

The equation of angular momentum conservation in a disc is (Hameury et al. hmd98 (1998)):

$`j{\displaystyle \frac{\mathrm{\Sigma }}{t}}={\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}\left(r\mathrm{\Sigma }jv_\mathrm{r}\right)`$ $`+`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}\left({\displaystyle \frac{3}{2}}r^2\mathrm{\Sigma }\nu \mathrm{\Omega }_\mathrm{K}\right)`$ (1)
$`+`$ $`{\displaystyle \frac{j_2}{2\pi r}}{\displaystyle \frac{\dot{M}_{\mathrm{tr}}}{r}}{\displaystyle \frac{1}{2\pi r}}T_{\mathrm{tid}}(r)`$
where $`\mathrm{\Sigma }`$ is the surface column density, $`\dot{M}_{\mathrm{tr}}`$ is the rate at which mass is incorporated into the disc at radius $`r`$, $`v_\mathrm{r}`$ is the radial velocity in the disc, $`j=(GM_1r)^{1/2}`$ is the specific angular momentum of material at radius $`r`$ in the disc, $`\mathrm{\Omega }_\mathrm{K}=(GM_1/r^3)^{1/2}`$ is the Keplerian angular velocity (where $`M_1`$ is the primary mass), $`\nu `$ is the kinematic viscosity coefficient, and $`j_2`$ the specific angular momentum of the material transferred from the secondary. $`T_{\mathrm{tid}}`$ is the torque due to the tidal forces, for which we use the formula of Smak (s84 (1984)), derived from the determination of tidal torques by Papaloizou & Pringle (pp77 (1977)):

$`T_{\mathrm{tid}}=c\omega r\nu \mathrm{\Sigma }\left({\displaystyle \frac{r}{a}}\right)^5`$ (2)
where $`\omega `$ is the angular velocity of the binary orbital motion, $`a`$ the binary orbital separation and $`c`$ a numerical constant taken so as to give an average disc radius equal to some chosen value.

In the TTI model, $`c`$ is no longer a constant, but depends on whether the disc is eccentric or not. Ichikawa et al. (iho93 (1993)) reproduced supercycles assuming that the disc becomes eccentric when its radius reaches $`r_{3:1}0.46a`$; $`c`$ is then increased by a large factor. Afterward the disc starts shrinking and can no longer maintain its eccentricity when its radius becomes smaller than $`r_{\mathrm{crit0}}0.35a`$, at which point $`c`$ returns to its normal value. We choose to take a value $`c_0`$ when the disc is axisymmetric and a larger value $`c_1`$ when the disc is eccentric. This is our Prescription $`I`$ for $`c`$, quite similar to that of Ichikawa et al. (iho93 (1993)) who extended the work of Osaki (o89 (1989)) by taking into account a non zero tidal torque during quiescence.

We have included this variable tidal torque in the disc instability model used in paper I, which is a modified version of the DIM described in Hameury et al. (hmd98 (1998)) that includes heating by the stream impact and the tidal torque (We have corrected Eq. (1) of paper I, in which the heating rate was overestimated by a factor $``$ 2). In the following, we also assume that $`\alpha =\alpha _{\mathrm{cold}}=0.04`$ in quiescence and $`\alpha =\alpha _{\mathrm{hot}}=0.2`$ in outburst. The inner radius is fixed at $`r_{\mathrm{in}}=10^9`$ cm.

We obtain the cycle shown in Fig. 1. After a superoutburst, the disc mass and radius are small. The disc mass slowly increases, and several normal outbursts are triggered due to the thermal instability. During each outburst, only a small amount of mass is accreted by the white dwarf, so that the average disc mass increases; a significant disc expansion occurs during the rise to maximum, when the heat front reaches the outer edge of the disc. The disk then shrinks when a cooling front brings the disc back to a cool state, causing a significant reduction of the outward angular momentum flow. After the end of an outburst, the disc radius very slowly decreases as a result of mass addition and viscous diffusion, but the net balance is an increase of the radius from one outburst to the next one. Eventually, $`r_{\mathrm{out}}`$ becomes larger than $`r_{3:1}`$ during an outburst and the tidal torque is increased, as well as the associated energy dissipation. The disc is then maintained in the hot state until most of the mass has been accreted. During this phase the radius decreases until $`r<r_{\mathrm{crit0}}`$. The disc returns to a circular shape, the anomalous tidal dissipation stops, and a cooling wave starts from the disc outer edge, bringing the whole disc in a cool state.

### 2.1 Radial dependence of $`c`$

There is, however, no compelling reason to assume that the whole disc is affected by the development of the tidal instability. We therefore use another prescription for $`c`$ (hereafter Prescription $`II`$). We assume that when the disc is circular, $`c=c_0`$; when eccentric, $`c(r)=c_1`$ for $`r>0.9r_{\mathrm{crit0}}`$ and $`c=c_0`$ for $`r<0.8r_{\mathrm{crit0}}`$. A linear interpolation is assumed between $`0.8r_{\mathrm{crit0}}`$ and $`0.9r_{\mathrm{crit0}}`$. We compare and discuss the results obtained using both prescriptions in 3.1.

### 2.2 Time delays

Temporal variations of $`c(r)`$ are not instantaneous. Observations show that normal superhumps appear within a few of days after the onset of a superoutburst, and that they develop within a day, so that it does not take long for the disc to become eccentric (Semeniuk s80 (1980)). However, it takes the disc at least a dynamical time to change its geometrical shape; we assumed that, when the tidal instability sets in, $`c`$ increases linearly on a timescale $`t_\mathrm{r}=\mathrm{2\hspace{0.33em}10}^4`$ s. Note that in WZ Sge systems, the early superhumps that first appear seem to be of a different nature from the normal superhumps; Osaki & Meyer (om02 (2002)) suggested that they could be related to the 2:1 resonance.

Similarly, the disc must dissipate its excess energy in order to return to a circular shape. The relative energy difference between a circular and an elliptic orbit with eccentricity $`e`$ with the same angular momentum is $`\mathrm{\Delta }E/E_{\mathrm{circ}}=e^2`$. As the dissipation at radius $`r`$ is $`DGM\mathrm{\Sigma }\nu /r^3`$, the time needed for the disc to return to a normal state is:

$$t_{\mathrm{e}\mathrm{c}}=\frac{\mathrm{\Delta }E}{D}=e^2\frac{r^2}{2\nu }=e^2t_\nu $$
(3)
where $`t_\nu `$ is the viscous time in the outer part of the disc. When the disc reverts to its normal state, the cooling wave has already started or will start with the end of tidal instability (see section 3). $`t_\nu `$ therefore refers to the cool state, and is of the order of the recurrence time, since the outbursts obtained here are of the inside-out type. Typical values are tens of days. Numerical simulations of the development of the tidal instability show that $`e`$ is of order of 0.1–0.3 (Murray m96 (1996)), which is confirmed by observations (Patterson et al. p00 (2000)), and $`t_{\mathrm{e}\mathrm{c}}`$ should be a few percent of $`t_\nu `$. We then assume that $`c`$ linearly decreases when the disc returns to a circular shape on a time $`t_\mathrm{f}=\mathrm{2\hspace{0.33em}10}^5`$ s.

We use three different sets of parameters that are representative for SU UMa, a typical ER UMa star and a typical WZ Sge star respectively, shown in Table 1. $`r_{\mathrm{circ}}`$ is the circulation radius at which a particle leaving the Lagrangian point would stay in circular orbit if there were no accretion disc (see Paper I for details). One should note that $`r_{\mathrm{circ}}`$ is calculated taking into account the gravitational potential of the secondary, and therefore differs slightly from the Keplerian value. In the following section, we will use the term ultra-low mass-ratio (ULMR) stars when referring to simulations with ER UMa parameters.

## 3 Decoupling the thermal and tidal instabilities

### 3.1  A property of ultra-low mass ratio systems ?

Hellier (h01 (2001)) suggested that the ultra-low mass ratio of ER UMa and WZ Sge stars is responsible for their peculiarities. In their case, the Roche lobe extends far beyond the 3:1 resonance radius, and the disc could in principle be larger than $`r_{3:1}`$. Hellier proposed that regions outside the resonance radius suffer a lower tidal dissipation than the portions in the 3:1 resonance zone. Then, during a superoutburst, heating by the tidal torque will not be as efficient as in other SU UMa stars, and a cooling wave could start while the disc is still eccentric, and will remain so several days after the end of the superoutburst.

This would then account for late superhumps. In addition, Hellier suggested that this model could also explain the echo outbursts in WZ Sge stars and the short ER UMa supercycle: at the end of a superoutburst, the disc would be still eccentric, and the tidal dissipation sufficiently large to trigger short outbursts during which the heating wave does not propagate to the outer edge (one should add that such outbursts are not only short but also have amplitudes lower than outbursts in which the heating front reaches the outer disc’s edge). ER UMa stars would then be systems always in an eccentric state, while WZ Sge stars would finally return to complete quiescence with circular discs. The difference between the supercycle duration of ER UMa and WZ Sge stars would be due to higher mass-transfer rates in ER UMa’s.

Fig 1 shows the simulated light curve for the parameters of SU UMa, assuming a factor 20 between $`c_1`$ and $`c_0`$. The superoutburst ends when the disc radius becomes lower than $`r_{\mathrm{crit0}}`$, i.e. the end of the tidal instability causes the end of the superoutburst. Fig 2 shows the simulated light curve for the parameters of ER UMa, with the same $`c_1/c_0`$. The mass transfer rate of $`10^{16}`$ g s-1 chosen here is close to the critical value above which the system is stable (within a factor 1.3). In this case, as predicted by Hellier, the superoutburst ends because the thermal instability stops, while the tidal instability is still effective. However, contrary to Hellier’s hypothesis (h01 (2001)), this phase is very short because of the rapid shrinking of the disc. Therefore, the duration of the supercycle is much larger than observed (there is practically no quiescence phase in real ER UMa stars): short outbursts are due to the usual thermal-viscous instability. As a consequence, contrary to Hellier’s assumption, the heating front in short outbursts always reaches the outer disc’s edge. Mass transfer rates closer to the stability limit do not lead to shorter supercycles, but rather to longer superoutbursts; this is actually the reason for which it is impossible to obtain very short supercycles in the standard TTI model without changing the tidal instability condition (Osaki 1995a ).

We have investigated the influence of the ratio $`c_1/c_0`$ and of the Roche size on our results. Fig 3 compares the superoutburst behavior for $`c_1/c_0`$ factors of 20, 50 and 80 for ULMR stars with prescription I. As can be seen, the width of the superoutburst strongly depends on $`c_1`$. Very small and very large values of $`c_1`$ both give short superoutburst, either because the tidal heating is too small (small $`c_1`$), or because the disc contraction under the effect of the tidal torque is too large (large $`c_1`$). Typical $`c_1/c_0`$ should be in the range $``$ 20 – 50 in order to reproduce the observed superoutburst durations.

The decoupling of the tidal and thermal instability was found for values of $`c_1/c_0`$ lower than 80, for both prescription I and II for ER UMa parameters, thereby confirming the first part of Hellier (h01 (2001)) proposal. On the other hand, for the parameters of SU UMa, the instability stopped because the disc shrinks below $`0.35a`$, except for $`c_1/c_0<10`$; this limit is increased to $``$ 20 when one uses prescription II.

These results strongly support the idea that the decoupling of thermal and tidal instability is a property of the ultra low mass ratio systems. One should however keep in mind that this conclusion depends somewhat on the assumed radial dependence of $`c`$, and on the assumed $`c_1/c_0`$; extreme values (less than 10 or more than 50) cause too short superoutburst to be triggered but ER UMa stars sometimes have short superoutbursts.

### 3.2 The superoutburst cooling wave

As shown in the previous section, decoupling of tidal and thermal instabilities is expected in ultra low mass ratio SU UMa stars. However, this is not sufficient to explain the late superhumps or echo outbursts observed in ER UMa and WZ Sge systems, as suggested by Hellier (h01 (2001)). When a cooling front starts to propagate, the radius decreases rapidly (in less than a few days) to its minimum value in the supercycle. If this minimum is less than $`r_{\mathrm{crit0}}`$, then the disc will return to a circular shape in a time scale $`t_\mathrm{f}`$, and one does not expect to observe late superhumps. If on the other hand the minimum radius is larger than $`r_{\mathrm{crit0}}`$, the disc always remains eccentric, the superhumps are permanent, but one then loses the very mechanism that was supposed to trigger superoutbursts.

The disc shrinkage at the end of an outburst or a superoutburst is therefore a key question. The radius evolves under the influence of three main effects, as shown in Eq. (1): (i) the addition of mass with low specific angular momentum at the outer edge of the disc and (ii) the tidal torque that both act to contract the disc, and (iii) the viscous outward transport of angular momentum that tends to make the disc larger. These effects are parameterized by $`\dot{M}_{\mathrm{tr}}`$, $`c`$ and $`\alpha `$ respectively. In order to determine their influence on the disc radius variations, we have changed their values by a factor 2 (linear variation on a time scale of 104s) for a steady disc in a hot state (Fig 4). The disc radius in the new steady state is unchanged, except in the case where $`c`$ is modified; the short term influence of these variations are summarized in table 2 which shows the value of $`d\mathrm{log}r/d\mathrm{log}x`$, $`10^4`$ s after the beginning of the variation, $`x`$ being one of the three parameters. As can be seen, $`d\mathrm{log}r/d\mathrm{log}c`$ is 5 to 10 times smaller than the two other quantities. This explains why $`c_1/c_0`$ has to be larger than 10 in order to obtain a superoutburst, and why this ratio needs to be larger than 20 in order for the radius to decrease fast enough during a superoutburst and reach $`r_{\mathrm{crit0}}`$ before the cooling wave starts.

In principle, the disc shrinkage could be slower, or possibly delayed if the mass transfer rate from the secondary were reduced, possibly as a result of a decrease of the illumination of the secondary. However, when the cooling wave starts, $`\alpha `$ decreases from $`\alpha _{\mathrm{hot}}`$ to $`\alpha _{\mathrm{low}}`$ in the outer disc regions, i.e. by a factor 5. $`c`$ has to be kept fixed, since one requires the disc to remain eccentric. Because of the large change in viscosity when a cooling front starts to propagate, $`\dot{M}_{\mathrm{tr}}`$ would then have to be reduced by an extremely large factor; we made several numerical experiments with various prescriptions for a change in $`\dot{M}_{\mathrm{tr}}`$ (amplitude and possible delay with respect to the formation of a cooling front), and we were unable to find a case in which the minimum disc radius was not obtained shortly after the end of a superoutburst.

## 4 Possible solutions

### 4.1 ER UMa supercycles

The standard TTI model cannot reproduce the very short supercycles observed in ultra-low mass-ratio systems. Osaki (1995a ) reproduced RZ LMI’s supercycle by ending prematurely the superoutburst assuming $`r_{\mathrm{crit0}}=0.42a`$ instead of the usual value of $`0.35a`$ (note that in the TTI framework, one assumes $`r_{\mathrm{crit0}}=0.35a`$ when modelling WZ Sge, even though its mass ratio is presumably smaller than in most SU UMa stars). In addition, the presence of superhumps after the end of the superoutburst implies that the eccentricity stops much later than the end of a superoutburst. The only ways for the disc to remain eccentric after the end of a superoutburst are then: (i) $`t_\mathrm{f}`$ is very long, or (ii) the accretion disc is always eccentric in ER UMa systems.

In the first case, the transition time $`t_\mathrm{f}`$ must be much larger than the one used here. However, as shown in Fig. 3, the interval between a superoutburst and the next normal outburst is far too long compared to observations. The second solution implies that the tidal torque is no longer responsible for superoutbursts, since $`c(r)=c_1`$ is a constant. However, the presence of large outbursts is still possible, as the (modified) DIM predicts the alternation of narrow and wide outbursts (see paper I) for large enough mass transfer rates. If the mass transfer rate is constant, the light curve consists of one or two small outbursts surrounded by large ones, but small variations of the mass transfer rate can easily lead to ER UMa type light curves, provided that these variations show some regularity. Irradiation of the disc and of the secondary can also account for the presence of long and short outbursts; Hameury et al. (hlw00 (2000)) have included these effects in the standard disc instability model and produced light curves typical of systems such as RZ LMi.

We therefore conclude that ER UMa stars should be dwarf novae with a permanently eccentric accretion disc, thereby accounting for superhumps, and where the illumination of the disc and the secondary star plays an important role. We thus predict that superhumps should exist at all phases of the supercycle of ER UMa stars; this apparently agrees with observations (Gao et al. g99 (1999)).

### 4.2  Echo outbursts in WZ Sge stars

WZ Sge stars have very long supercycles and superoutbursts. In addition, no normal outbursts between two consecutive superoutbursts are observed. Some WZ Sge stars also show echo outbursts at the end of the superoutburst: several small outbursts, spaced every tens of days, during which superhumps are still present.

Smak (s93 (1993)) deduces the mass transfer rate of WZ Sge from the luminosity of the hot spot. It is observed to increase at least by a factor 10 during superoutburst, and it decreases afterwards, remaining larger than the quiescent value during several tens of days. Such mass transfer rate fluctuation could result from irradiation of the secondary star by the white dwarf. Hameury et al. (hlw00 (2000)) have shown that irradiation (including disc irradiation) could indeed account for some peculiarities of WZ Sge stars: for example, they reproduced the echo outbursts phenomenon without including the tidal instability in the DIM; however, they did not reproduce a full cycle with long recurrence times, and the echo outbursts they obtained were slightly too dim.

The long recurrence times can be due to low alpha value (Smak s93 (1993); Osaki 1995b ), possibly due to a decay of the MHD turbulence that would lead to a time-dependant $`\alpha `$prescription (Osaki et al. o01 (2001)), thereby explaining the echo outbursts. A low viscosity could result from the secondary being a brown dwarf (Meyer & Meyer-Hofmeister mm99 (1999)), but one would have to explain why the viscosity is so much lower in these systems as compared to other SU UMa systems which have comparable or even shorter orbital periods. Another possibility (Lasota et al. lhh95 (1995), lkc99 (1999); Hameury et al. hlh97 (1997)) is that WZ Sge stars are in a stable low state between superoutbursts, thus explaining the absence of normal outbursts. This requires a hole in the central regions of the disc, as a result of either a moderate magnetic field, or of evaporation. The superoutbursts would then be triggered by an externally imposed increase of the mass transfer rate, the long duration of the outburst and the large mass accreted onto the white dwarf being due to the irradiation of the secondary. They did not include the tidal instability and, as in the case of ER UMa stars, the presence of superhumps and late superhumps is explained if the accretion disc is always eccentric. If one combines these results with those of Hameury et al. (hlw00 (2000)) on the echo outburst, one should be able to reproduce a WZ Sge star light curve with echo outbursts, which would however be slightly different than the observed ones.

We have used our TTI model including irradiation and the presence of an inner hole for the accretion disc as determined by equation (7) of Hameury et al. (hlw00 (2000)). We use the parameters of EG Cnc. The mass transfer rate from the secondary is assumed to be affected by irradiation according to:

$$\dot{M}_{\mathrm{tr}}=\mathrm{max}(\dot{M}_2,\gamma <\dot{M}_{\mathrm{acc}}>)$$
(4)
where $`\dot{M}_2`$ is the mass transfer rate that would be obtained in the absence of illumination, $`\gamma `$ is a constant in the range \[0-1\], and $`<\dot{M}_{\mathrm{acc}}>`$ is some average of the mass accretion onto the white dwarf, defined as:

$$<\dot{M}_{\mathrm{acc}}>=_{\mathrm{}}^{t_0}\dot{M}_{\mathrm{acc}}e^{(t_0t)/\mathrm{\Delta }t}𝑑t$$
(5)
$`t_0`$ being the actual time. This prescription differs slightly from that of Hameury et al. (hlh97 (1997)), in that the average of the mass accretion rate is used instead of its value at time $`t_0`$; this accounts for the fact that a fraction of the accretion luminosity is released at later times when the white dwarf cools. In the following, we have taken $`\mathrm{\Delta }t`$ = 14 days. We start with a stable system in the low state with a mass transfer rate of 1.5 $`10^{15}`$ g.s-1, in agreement with Smak (s93 (1993)) estimate for WZ Sge ($`\mathrm{2.\hspace{0.33em}10}^{15}`$g s-1). At $`t=10`$ days, the mass transfer rate $`\dot{M}_2`$ is increased to $`\mathrm{3\hspace{0.33em}10}^{16}`$g.s-1; 20 days later, $`\dot{M}_2`$ is reduced to its initial value. Our results are insensitive on how $`\dot{M}_2`$ returns to its normal value, because illumination dominates in Eq. (4).

This model predicts that the luminosity should increase just before the triggering of a superoutburst. The increase need not be large – one simply requires the mass transfer rate to exceed the critical rate for stability; this does therefore not contradict the observations of Ishioka et al. (i02 (2002)) that no strong orbital hump due to a high mass transfer rate was detected during the early phase of WZ Sge 2001 outburst.

We obtain a light curve with a very long superoutburst followed by 3 normal outbursts that last $`3`$ days and occur within two months of the main outburst (Fig. 5). The disc radius remains always larger than 0.35 $`a`$, and is always in an eccentric state; permanent superhumps are thus expected. However, the interval between normal outbursts is twice longer than for EG Cnc echo outbursts, which are brighter than observed. We did not obtain a better fit with ER UMa light curves; in view of the very crude assumptions made to derive the mass transfer rate from the secondary under the influence of illumination, this is not really a surprise.

## 5 Conclusion

Hellier (h01 (2001)) suggested that several limitations of the TTI model could be overcome if, in ultra-low mass ratio systems, the thermal instability ends while the disc is still eccentric. This would explain the late superhumps, and Hellier (h01 (2001)) also suggested that this would account for the very short ER UMa supercycles and the WZ Sge stars echo outbursts.

We do find that the decoupling of the two instabilities can happen much more easily in ULMR systems than in other SU UMa stars. However, the tidal instability stops almost immediately after the end of the thermal-viscous instability (delay less than a day), because the disc shrinks very rapidly after the departure of a cooling wave. We estimate that it takes only a few days for the disc to return to a circular shape after the cessation of the tidal instability. This is still insufficient to explain late superhumps that are found tens of days after the superoutburst in some systems.

We propose instead that the discs in ER UMa systems are permanently eccentric, and that they should therefore have permanent superhumps; there are observational indications of this (Gao et al. g99 (1999)). The short supercycle is then due to either an alternation of narrow and wide outbursts as in SS Cyg (see Buat-Ménard et al. b01 (2001)), or to the effects of irradiation (see Hameury et al. hlw00 (2000)). Further observations of late superhumps should confirm their permanence in ER UMa stars.

As in Hameury et al. (hlh97 (1997)), we propose that WZ Sge stars are in a stable cold state during quiescence. The superoutburst is then triggered by an enhancement of mass transfer rate, and echo outbursts are triggered while the mass transfer rate is still high after the superoutburst. However, simulated light curves are still somewhat different from observations; this might be due to our crude treatment of the secondary illumination.

###### Acknowledgements.

This work was supported in part by a grant from Programme National de Physique Stellaire of the CNRS. We thank J.P. Lasota for a very careful reading of the manuscript and for helpful comments. We are grateful to the referee, Prof. Y. Osaki, for pointing us an error in one of the sets of parameters used.


text:
# Critical Percolation in High Dimensions

## Abstract

We present Monte Carlo estimates for site and bond percolation thresholds in simple hypercubic lattices with 4 to 13 dimensions. For $`d<6`$ they are preliminary, for $`d6`$ they are between $`20`$ to $`10^4`$ times more precise than the best previous estimates. This was achieved by three ingredients: (i) simple and fast hashing which allowed us to simulate clusters of millions of sites on computers with less than 500 MB memory; (ii) a histogram method which allowed us to obtain information for several $`p`$ values from a single simulation; and (iii) a variance reduction technique which is especially efficient at high dimensions where it reduces error bars by a factor up to $`30`$ and more. Based on these data we propose a new scaling law for finite cluster size corrections.

In spite of decades of intensive studies sa94 , percolation remains an active subject of research. While there has been enormous progress in understanding percolation in 2 dimensions cardy , mainly because of conformal invariance, progress in high dimensions has been much slower. It has been known since long time that $`d=6`$ is the upper critical dimension sa94 , and expansions of $`p_c`$ in $`1/(2d1)`$ have been given already more than 20 years ago. But up to now there exists no detailed numerical study of logarithmic corrections in $`d=6`$, finite size corrections are not yet understood for $`d>6`$, and even numerical estimates of $`p_c`$ in $`6`$ dimensions are very poor. One reason for this is obviously that straightforward simulations of large lattices in $`6`$ dimensions require huge amounts of fast memory. This lack of stimulus by numerical verification certainly was part of the reason for the slow analytical progress.

It is the purpose of the present note to improve on this situation by presenting precise numerical estimates of $`p_c`$ (and of finite cluster size corrections) for site and bond percolation on simple hypercubic lattices with $`d=6`$ to $`d=13`$.

Our main results are summarized in Table 1, where we also include preliminary results for $`d=4`$ and $`d=5`$. We also give the best previous estimates for $`p_c`$ and expansions in $`1/(2d1)`$. We shall discuss them later in more detail, but here we just point out that our new estimates are vastly better than all previous ones. They were possible, with rather modest effort (we used only fast PCs and Alpha work stations, with altogether ca. $`10^3`$ CPU hours), due to several important ingredients:

1) We used as basic routine a standard breadth-first version of Leath’s algorithm which simulates single clusters. We do not use the popular Hoshen-Kopelman method since that would require prohibitively large memory if we want to simulate large clusters. In Leath’s method, one writes the coordinates of each cluster site (which consist of a single integer – see item 2 below) into a first-in-first-out queue Q, where each new entry represents a newly wetted neighbor of the oldest entry in the queue.

2) We used a simple but very efficient form of hashing s90  for storing the information whether a site has already been wetted or not. On Compaq Alpha work stations with 64 bit long integers, we labelled lattice sites by a single long integer. Using as lattice size $`L`$ an odd number slightly smaller than $`2^{64/d}`$, we label the neighbors of site $`i`$ as $`i\pm 1,i\pm L,\mathrm{}i\pm L^{d1}`$. If we want to simulate bond percolation clusters with roughly $`N`$ sites, we find first the power of 2 nearest to $`N`$, $`2^kN`$, and use it to obtain for each site $`i`$ its key $`m_i=(i\mathrm{mod}\mathrm{\hspace{0.33em}2}^k)`$ (notice that this is done most efficiently by bitwise AND). Assume now that site $`i`$ with key $`m_i`$ is the $`n`$-th site wetted. Then an entry is written into the $`m_i`$-th element of an array of pointers S of size $`2^k`$. This element points to the $`n`$-th element of a structure (L,Q) where Q is the above queue and L is a linked list. In Q, the coordinate $`i`$ is stored. The element of L remains empty, if the key $`m_i`$ had not been encountered before. Otherwise, if some other site $`j`$ with the same key $`m_j=m_i`$ had been wetted in an earlier step $`n^{}<n`$, the old element of S (which had pointed to $`n^{}`$) is written in the $`n`$-th element of L. In this way we can deal with virtual lattices of $`2^{64}`$ sites, using $`2^k+2N_{\mathrm{max}}`$ storage places, where $`N_{\mathrm{max}}`$ is an upper bound on the size of clusters to be simulated. The algorithm is slightly different for site percolation where a tested site has to be excluded from further growth even if it is not wetted, in contrast to bond percolation. It also has to be modified on machines with only 32 bit long integers where a pair of numbers replaces $`i`$ and a pair of co-prime odd numbers $`L_1`$ and $`L_2`$, both slightly smaller than $`2^{32/d}`$, replaces $`L`$.

This is not as storage efficient as the recent algorithm of pzs01 . But it works with usual (pseudo-)random number generators (we used the 4-tap generator with period $`2^{9689}1`$ of ziff ), while the algorithm of pzs01  needs a random function generator. The most easily available random function generator today is the Data Encryption Standard press  which is rather slow when implemented in software and of unproven quality for this application (it was developed for entirely different purposes, and lacks any published theoretical justification).

3) In order to estimate cluster statistics for several values of $`p`$ from a single run at nominal value $`p_0`$, we use a trick similar to the histogram methods used by Dickman dickman  for the contact process (see also balle ). If a cluster with $`n`$ wetted sites and $`b`$ non-wetted boundary sites was generated with nominal value $`p_0`$, it contributes to the ensemble with $`p_0`$ replaced by $`p`$ with weight

$$W=(p/p_0)^n((1p)/(1p_0))^b.$$
(1)
Instead of collecting histograms for cluster numbers with fixed $`n`$ and $`b`$ (which would have led to excessively large arrays) we calculated on the fly three distributions: One for the nominal $`p_0`$ (which was chosen close to $`p_c`$ as estimated from short test runs and from Eqs.(7) resp. (8)) and two for neighbouring values $`p_\pm =p_0\pm \delta p`$, using Eq.(1) for the latter. Observables at $`p`$-values in between (including $`p_c`$) were obtained by geometric (i.e. linear in logarithm) interpolation. Having three values of $`p`$ instead of just two allowed us to check that the error due to the interpolation was negligible.

4) Our main observable will be the number $`M(t)`$ of wetted sites with “chemical distance” $`t`$ from the seed of the cluster (i.e. the number of sites infected at time $`t`$, if cluster growth is interpreted as spreading of an epidemic). For $`d>6`$ we expect its average $`M(t)`$ to become a constant at the critical point, since the process is basically a branching process with small corrections. But instead of using $`M(t)`$ itself, we obtain a less noisy signal by the following trick which would give the exact ensemble average of $`M(t)`$ if the cluster growth indeed were a branching process footnote .

Let us assume we have a (still growing) cluster $`C`$ with $`M(t)`$ sites wetted at step $`t`$, and denote by $`M^+(t)`$ the number of free neighbors, i.e. the number of sites which can be wetted at step $`t+1`$. The actual number wetted will fluctuate, but the expected average number, conditioned on $`C`$ and thus also on $`M(t)`$, is exactly given by

$$𝐄[M(t+1)|C]=pM^+(t).$$
(2)
Thus the expected geometric increase of the number of wetted sites, still conditioned on $`C`$, is

$$𝐄[M(t+1)/M(t)|C]=pM^+(t)/M(t)$$
(3)
and its weighted sample average over all clusters is

$$r(t)\frac{\underset{C}{}M(t)𝐄[M(t+1)/M(t)|C]}{_CM(t)}=\frac{pM^+(t)}{M(t)}.$$
(4)
Our estimate for the true ensemble average of $`M(t)`$ is then finally

$$\widehat{M(t)}=\underset{t^{}=0}{\overset{t1}{}}r(t^{}).$$
(5)
Since we measured also the direct estimate $`M(t)`$ and the (co-)variances of both estimates, we can also compute the variance of any linear combination of both. For $`d=4`$ and $`d=5`$, where both variances are comparable and the covariance is negative, a substantial achievement is obtained by taking as the final estimate the linear combination with the smallest variance.

An expression similar to Eq.(5) can be obtained also for the rms. radius, if we replace the ratios in Eqs.(3,4) by differences and the product in Eq.(5) by a sum. This would also be exact and non-fluctuating if the cluster growth were a branching process with translation invariance.

The variance reduction due to Eq.(5) is largest for small $`t`$. Yet, for bond percolation in $`d=11`$, it gave even for the largest $`t`$ (=200) a factor $`1/1000`$ over using just $`M(t)`$. For $`d=6`$ and $`t=2000`$, the reduction was still by a factor $`140`$. Indeed, there were substantial improvements even for $`d=4`$ and 5, while the improvement in $`d=3`$ was marginal. For site percolation the improvements were similar but somewhat less dramatic.

We should note that we calculated also $`P(t)`$, the probability that a cluster survives at least $`t`$ steps (i.e. has “chemical radius” $`t`$), the cluster size distribution $`P(n)`$, and the spatial extent of clusters with $`n`$ sites. All of them gave vastly more noisy signals (since we could not use a similar variance reduction trick as for $`M(t)`$) and were not used in estimating the critical point.

Results for $`\widehat{M(t)}`$ are shown in Figs.1 to 3 for $`d=6,7`$, and 11. In all these figures we show results for bond percolation. Results for site percolation are similar albeit somewhat more noisy. In the first two cases we checked explicitly that no cluster was larger than the virtual lattice size $`L`$ (which was $`>500`$ in both cases), so there are strictly no finite lattice size effects. For $`d9`$ this was no longer possible for the cluster sizes used here (typically up to $`10^410^6`$ sites), but we can easily convince ourselves that also there finite size effects are negligible.

In each of the 3 figures the critical point $`p_c`$ is characterized by $`d\widehat{M(t)}/dt0`$ for $`t\mathrm{}`$. For $`d>6`$ we also have $`\widehat{M(t)}const`$ for $`p=p_c`$, while we see a logarithmic divergence in $`d=6`$ as predicted by the renormalization group egg78  (see Fig.1). Unfortunately, the detailed behaviour of $`M(t)`$ in $`d=6`$ has not yet been calculated, though the results of egg78  and the fact that $`\nu _t`$ (the exponent controlling the correlation time) is 1, suggest $`M(t)[\mathrm{log}t]^{2/7}`$ to leading order. Therefore, and since it is notoriously difficult to verify logarithmic terms (see, e.g., g97 ; ghs94 ; po01 ), we have not attempted any detailed analysis.

From Figs.2 and 3 we also see that corrections to scaling decrease strongly with dimension for $`d>6`$. In Fig.2 we see a straight line for $`p=p_c`$ when plotting $`\widehat{M(t)}`$ against $`1/\sqrt{t}`$, showing that the leading correction term is $`t^{0.5}`$ in $`d=7`$. Similarly, a straight line is obtained for $`d=11`$ when using $`t^{2.18}`$ (Fig.3). All these (and similar results for other values of $`d>6`$ and for site percolation, not shown here) strongly suggest anomalous scaling

$$M(t)=M_{\mathrm{}}const/t^{\omega (d)},$$
(6)
similar to the scaling for self avoiding walks in $`d>d_c`$ found in po01 . But while the exponents were simply $`(dd_c)/2`$ in po01 , they seem to depend less trivially on $`d`$ in the present case – although we cannot exclude the possibility that $`\omega (d)=(dd_c)/2`$ also here, and the observed deviations are due to higher order corrections. The latter is indeed suggested by the results of aah84 .

The constants $`M_{\mathrm{}}`$ defined in Eq.(6) are plotted in Fig.4 against $`d6`$ on doubly logarithmic scale. They seem to fall on parallel straight lines, suggesting a universal law $`M_{\mathrm{}}1(d6)^a`$ with $`a=0.73\pm 0.03`$. But a closer look reveals that deviations from this are significant (although they are small), suggesting that it holds neither for $`d\mathrm{}`$ nor for $`d6`$ exactly.

Let us finally discuss the $`p_c`$ values given in Table 1. They should be compared to the predictions gr78

$$p_{c,\mathrm{bond}}=s+5s^3/2+15s^4/2+57s^5+\mathrm{}$$
(7)
and gsr76

$$p_{c,\mathrm{site}}=s+3s^2/2+15s^3/4+83s^4/4+\mathrm{}$$
(8)
with $`s=1/(2d1)`$. The dots in these equations stand for higher powers of $`s`$. It was suggested in gsr76  that they can be approximated, for site percolation at least, by adding 2/3 of the last term,

$$p_{c,\mathrm{site}}s+3s^2/2+15s^3/4+415s^4/12.$$
(9)
The full series are presumably only asymptotic. It is thus a priori not clear whether any of these equations should be good approximations to the present data. From Table 1 we see that Eq.(9) is excellent in the range studied here, but it has wrong asymptotic behaviour and should be worse than Eq.(8) for $`d15`$. As seen from Fig.5, the agreement with Eqs.(7) and (8) is indeed better than could have been expected: For bond percolation the difference decreases roughly as $`s^{7.1}`$ (instead of $`s^6`$), while for site percolation it decreases as $`s^{5.7}`$ instead of $`s^5`$. Obviously the next terms in Eqs.(7) and (8) would be needed for a more detailed comparison.

Finally we should remind of several heuristic formulas for $`p_c`$ values on various lattices. All early ansatzes of this type were already refuted in gb84  because they contradicted Eqs.(7) or (8). More recently, such heuristics have been discussed again in gm01  and in the papers quoted there. We have not attempted any detailed comparison in view of their complete lack of theoretical basis.

In summary, we have presented vastly improved estimates for percolation thresholds on high-dimensional hypercubic lattices. They should be compared to improved series expansions and/or rigorous bounds. At present such results are not available, partly because it had seemed that they could not be compared to any numerical estimates. Apart from this, the methods used in the present paper should also be of use in other similar problems. These include simulations of percolation backbones, conductivity exponents, percolation on more exotic lattices, directed percolation in high dimensions, and self avoiding walks. In all these cases both the hashing and the variance reduction should be of help in simulating larger systems with higher precision.

I thank Walter Nadler and Hsiao-Ping Hsu for discussions and for carefully reading the manuscript.


text:
# Lapse function for maximally sliced Brill-Lindquist initial data

## Abstract

For binary black holes the lapse function corresponding to the Brill-Lindquist initial value solution for uncharged black holes is given in analytic form under the maximal slicing condition. In the limiting case of very small ratio of mass to separation between the black holes the surface defined by the zero value of the lapse function coincides with the minimal surfaces around the singularities.

PACS number(s): 04.70.Bw, 04.20.Ex, 04.20.Fy

In a recent paper , the authors re-obtained the Brill-Lindquist initial value solution for two electrically uncharged black holes  using a point-mass model for the black holes followed by Hadamard’s “partie finie” regularization procedure in 3-dimensional space. (It can be shown that the same result appears by applying the dimensional regularization of Ref. .) In the present paper we supplement this result by the explicit calculation of the initial lapse function. By the time-symmetry assumption in the Brill-Lindquist initial value solution, the vanishing of the initial shift function automatically occurs. Therefore, all initial value metric coefficients will be known. Our result, for the first time, allows the determination of the motion of test fields and particles near the initial hypersurface for time scales which are small against the time scale of the motion of the two black holes. Similarly, the result will allow discussions of quantum gravitational problems near the initial hypersurface. Furthermore, the knowledge of the initial 4-metric is a step towards an (approximate) analytic treatment of the motion of two black holes starting from Brill-Lindquist initial value data. Finally, numerical head-on-collision simulations starting from those initial data may compare the numerically obtained initial lapse function with the analytic one of the present paper to check the accuracy of the numerical code. In our brief report, we additionally will address the question of how the surface of vanishing initial lapse function, which is the surface of infinite redshift for observers momentarily at rest in the initial hypersurface, is related in some limiting case to the minimal 2-surfaces associated with the Einstein-Rosen bridges of the Brill-Lindquist 3-geometry (see Ref. ), knowing that for the Schwarzschild black hole geometry the both types of surfaces coincide (see, e.g., Ref. ).

We use units in which $`16\pi G=c=1`$, where $`G`$ is the Newtonian gravitational constant and $`c`$ is the velocity of light. In the paper all vectors and their lengths are defined in the 3-dimensional Euclidean space endowed with a standard Euclidean metric; $`𝐱`$ is the position of an arbitrary point in this space.

In the generalized isotropic ADM gauge \[see Eqs. (7-4.22) and (7-4.23) in Ref. \] the 3-metric takes the form

$$g_{ij}\left(1+\frac{1}{8}\varphi \right)^4\delta _{ij}+h_{ij}^{\text{TT}}$$
(1)
and for the canonical conjugate $`\pi ^{ij}`$ of the 3-metric, which in terms of the extrinsic curvature $`K_{ij}`$ of the spacelike hypersurfaces $`t=\text{const}`$ is given by (see, e.g., Ref. )

$$\pi ^{ij}=g^{1/2}(g^{il}g^{jm}g^{ij}g^{lm})K_{lm},$$
(2)
where $`g^{ij}`$ is the inverse of $`g_{ij}`$ and $`gdet(g_{ij})`$, the condition

$$\pi ^{ii}0$$
(3)
holds.

In case of time-symmetric initial conditions the black-hole linear momenta $`p_{ai}`$ ($`a=1,2`$) and field canonical conjugate $`\pi ^{ij}`$, and thus the shift function, $`N_i=g_{0i}`$, too, are zero

$$p_{ai}=0,\pi ^{ij}=0,N_i=0.$$
(4)
Conformal flatness, i.e.

$$g_{ij}=\left(1+\frac{1}{8}\varphi \right)^4\delta _{ij},$$
(5)
additionally implies the vanishing of the transverse-traceless part of the 3-metric,

$$h_{ij}^{\text{TT}}=0.$$
(6)
In this case, the coordinate condition (3) coincides with the maximal slicing condition, i.e., $`g^{ij}K_{ij}=0`$. Under the above conditions, the lapse function $`N_0`$ ($`N_0^2=g_{00}+g^{ij}N_iN_j=g_{00}`$, as $`N_i=0`$) fulfills the equation

$$g^{1/2}g^{ij}D_iD_jN_0=\frac{1}{4}N_0g^{1/2}R,$$
(7)
where $`D_i`$ means the covariant derivative with respect to $`g_{ij}`$, and $`R`$ is the intrinsic curvature of the hypersurface $`t=\text{const}`$. The Eq. (7) is easily obtained from the field equation for $`\pi ^{ii}`$, i.e., from $`\pi ^{ii}/t0`$ \[see, e.g., Eq. (2.7a) in Ref. \].

After substituting the ansatz

$$N_0=\frac{1\psi /8}{1+\varphi /8}$$
(8)
into the Eq. (7), we obtain the following equation for the function $`\psi `$:

$$\left(1+\frac{1}{8}\varphi \right)\mathrm{\Delta }\psi =\left(1\frac{1}{8}\psi \right)\mathrm{\Delta }\varphi .$$
(9)
To derive Eq. (9) we have employed the relation, valid for application on scalar functions,

$$g^{1/2}g^{ij}D_iD_j=\frac{}{x^i}\left(g^{1/2}g^{ij}\frac{}{x^j}\right)$$

and the following formula, valid when $`h_{ij}^{\text{TT}}=0`$ (see, e.g., Ref. ):

$$g^{1/2}R=\left(1+\frac{1}{8}\varphi \right)\mathrm{\Delta }\varphi .$$

The Brill-Lindquist solution reads, ,

$$\varphi =8\left(\frac{\alpha _1}{r_1}+\frac{\alpha _2}{r_2}\right),$$
(10)
where the parameters $`\alpha _1`$ and $`\alpha _2`$ are constants and where $`r_1`$ and $`r_2`$ denote the coordinate distances between the field point and the black hole positions in the conformally related flat space. The parameters $`\alpha _1`$, $`\alpha _2`$ are related to the bare masses $`m_1`$, $`m_2`$ of the two black holes through the formulas (see Appendix A in Ref. )

$`\alpha _1`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(2r_{12}+{\displaystyle \frac{m_2m_1}{16\pi }}\right)+{\displaystyle \frac{1}{4}}r_{12}\sqrt{4+{\displaystyle \frac{m_1+m_2}{4\pi r_{12}}}+\left({\displaystyle \frac{m_1m_2}{16\pi r_{12}}}\right)^2},`$ (12)
$`\alpha _2`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(2r_{12}+{\displaystyle \frac{m_1m_2}{16\pi }}\right)+{\displaystyle \frac{1}{4}}r_{12}\sqrt{4+{\displaystyle \frac{m_1+m_2}{4\pi r_{12}}}+\left({\displaystyle \frac{m_1m_2}{16\pi r_{12}}}\right)^2},`$ (13)
where $`r_{12}`$ denotes the coordinate distance of the two black holes.

Equation (9) is invariant against $`\psi \varphi `$, $`\varphi \psi `$, which obviously suggests identical structures of the functions $`\varphi `$ and $`\psi `$ involved in this equation. We therefore search for a solution of the function $`\psi `$ of the form

$$\psi =8\left(\frac{\beta _1}{r_1}+\frac{\beta _2}{r_2}\right),$$
(14)
where $`\beta _1`$ and $`\beta _2`$ are some new constants. After substitution Eqs. (10) and (14) into Eq. (9) we make use of the regularization procedure of Ref.  (or Ref. ) which relies here on replacing $`f(𝐱)\delta (𝐱𝐱_a)`$ by $`f_{\mathrm{reg}}(𝐱_a)\delta (𝐱𝐱_a)`$, where the regularized value $`f_{\mathrm{reg}}(𝐱_a)`$ of the function $`f`$ at its (possibly) singular point $`𝐱=𝐱_a`$ we define by means of the Hadamard’s “partie finie” procedure. We expand $`f\left(𝐱_a+\epsilon 𝐧\right)`$ (where $`𝐧`$ is a unit vector and $`\epsilon >0`$ is a number) into a Laurent series around $`\epsilon =0`$ and as the regularized value of the function $`f`$ at $`𝐱_a`$ we take the zero-order coefficient of the series averaged over all directions $`𝐧`$:

$$f\left(𝐱_a+\epsilon 𝐧\right)=\underset{m=N}{\overset{\mathrm{}}{}}a_m(𝐧)\epsilon ^m,f_{\mathrm{reg}}\left(𝐱_a\right)\frac{1}{4\pi }𝑑\mathrm{\Omega }a_0(𝐧).$$
(15)
The usage of the definition (15) boils down in our computation to simple formulas

$`{\displaystyle \frac{\delta (𝐱𝐱_1)}{|𝐱𝐱_1|}}=0,{\displaystyle \frac{\delta (𝐱𝐱_1)}{|𝐱𝐱_2|}}={\displaystyle \frac{\delta (𝐱𝐱_1)}{|𝐱_1𝐱_2|}},`$ (17)
$`{\displaystyle \frac{\delta (𝐱𝐱_2)}{|𝐱𝐱_2|}}=0,{\displaystyle \frac{\delta (𝐱𝐱_2)}{|𝐱𝐱_1|}}={\displaystyle \frac{\delta (𝐱𝐱_2)}{|𝐱_1𝐱_2|}}.`$ (18)
Making use of Eqs. (15) we have obtained, by comparing the coefficients at Dirac delta functions, the following algebraic equations for the constants $`\beta _1`$ and $`\beta _2`$:

$`\left(1+{\displaystyle \frac{\alpha _2}{r_{12}}}\right)\beta _1`$ $`=`$ $`\left(1{\displaystyle \frac{\beta _2}{r_{12}}}\right)\alpha _1,`$ (20)
$`\left(1+{\displaystyle \frac{\alpha _1}{r_{12}}}\right)\beta _2`$ $`=`$ $`\left(1{\displaystyle \frac{\beta _1}{r_{12}}}\right)\alpha _2.`$ (21)
The unique solution of the both equations is easily obtained. It reads

$`\beta _1`$ $`=`$ $`\alpha _1{\displaystyle \frac{r_{12}+\alpha _1\alpha _2}{r_{12}+\alpha _1+\alpha _2}},`$ (23)
$`\beta _2`$ $`=`$ $`\alpha _2{\displaystyle \frac{r_{12}+\alpha _2\alpha _1}{r_{12}+\alpha _1+\alpha _2}}.`$ (24)
In the limit $`r_{12}\mathrm{}`$, where two completely separated Schwarzschild black holes have to arise, it holds: $`\beta _1=\alpha _1`$ and $`\beta _2=\alpha _2`$. If $`\alpha _1=0`$, then $`\beta _1=0`$ and $`\beta _2=\alpha _2`$ (analogously, if $`\alpha _2=0`$, then $`\beta _2=0`$ and $`\beta _1=\alpha _1`$). Also, if $`\alpha _1=\alpha _2`$, then $`\beta _1=\beta _2`$. Thus, the solution for $`\psi `$ does have correct limiting values. Furthermore, the post-Newtonian expansion of our solution for $`\psi `$ coincides with the corresponding part of the lapse function obtained in Ref.  through second post-Newtonian order.

For completeness, we give the initial 4-dimensional line element. It takes the form

$$ds^2=\left(\frac{1\psi /8}{1+\varphi /8}\right)^2dt^2+(1+\varphi /8)^4(dx^2+dy^2+dz^2),$$
(25)
where $`x`$, $`y`$, $`z`$ denote Cartesian coordinates and where the Eqs. (10), (14), and (Lapse function for maximally sliced Brill-Lindquist initial data) \[and also Eqs. (10) if the dependence in terms of bare masses is needed\] have to be substituted.

In the case of Schwarzschild black holes (there, $`\psi =\varphi `$), the surface $`N_0=0`$ (being there the cross-section of a $`t=\text{const}`$ slice with the event horizon) coincides with the throat of the Einstein-Rosen bridge which is a minimal 2-surface in 3-space. For binary black-hole systems, the surfaces $`N_0=0`$ need not to be identical with minimal 2-surfaces, respectively throats, in 3-space. However, an explicit calculation has shown that the initial value surface $`N_0=0`$ (or, $`\psi =8`$), coincides with the minimal surfaces around the singularities in the limiting case of very large separation between the black holes, i.e., we exactly have reproduced the Eq. (27) in Ref.  for uncharged black holes.

P.J. thanks the Theoretisch-Physikalisches Institut of the FSU Jena for hospitality during the realization of this work. The work of P.J. was supported by the EU Programme “Improving the Human Research Potential and the Socio-Economic Knowledge Base” (Research Training Network Contract HPRN-CT-2000-00137) and by the Polish KBN Grant No. 5 P03B 034 20.


text:
# Fractal Plate Tectonics

## Introduction

According to the theory of plate tectonics, the outermost $`100`$ km of the solid Earth is composed of a relatively small number of internally rigid plates (lithosphere) that move over a weak substrate (asthenosphere). The plate sizes and positions change over time and are driven by the internal mantle convection. DeMets et al.  performed a global inversion to determine the relative rotation rates of the 12 largest plates (the NUVEL-1 model) later refined into the NUVEL-1A solution \[DeMets et al., 1994\]. The edges of these plates, where they move against each other, are sites of intense geologic activity, such as earthquakes, volcanoes, and mountain building. These plate edges are not sharp narrow boundaries but are often constituted of complex systems of competing faults and other geological structures of width extending over several hundreds of kilometers for transform and subduction boundaries and up to thousands kilometers for continental collisions. It is now common lore to view such tectonic deformation as possessing some kind of self-similarity or fractal properties \[Scholz and Mandelbrot, 1989\], or better a hierarchical structure \[Ouillon et al., 1996\]. In addition, several researchers have repeatedly proposed to describe this multi-scale organization of faulting by a hierarchy of blocks of multiple sizes \[Gelfand et al., 1976; Gorshkov et al., 2001\]. In these models, blocks slide against each other, rotate, lock at nodes or triple junctions which may represent the loci of major earthquakes, in a way similar to (but at a reduced scale compared to) the relative motion of the major tectonic plates. An example of a complex network of faults is observed in the broad San Andreas transform boundary between the Pacific and the North American plates in California, for which Bird and Rosenstock  have suggested $`21`$ possible microplates within southern California alone.

Keeping in mind these ingredients of a few major plates at large scales on one hand and a hierarchical self-similar organization of blocks at the boundary scale on the other hand, the recent reassessment of present plate boundaries on the Earth by Bird  is particularly inspiring: taking into account relative plate velocities from magnetic anomalies, moment tensor solutions, and/or geodesy, to the 14 large plates whose motion was described by the NUVEL-1A poles, model PB2001 \[Bird, 2002\] includes 28 additional small plates, for a total of 42 plates. Bird  suggests that the cumulative-number/area distribution for his model follows a power-law for plates of less than 1 steradian of area.

## Statistical analysis of the size distribution of the $`42`$ plates

### Analysis with rank-order statistics and the Pareto law

The statistical analysis of such a data set is very difficult due to its small size $`N=42`$ but is not impossible. Figure 1 shows the complementary cumulative number $`N(A)`$ of plates as a function of area $`A`$ in steradians, i.e., the number of plates with an area equal to or larger than $`A`$. Most of the data except for a few largest plates follow a linear dependence in the double log-scale of the figure. The slope of this straight line is close to $`0.25`$. We first compare this sample with a pure power law (the Pareto distribution). For this purpose we use the rank-ordering method (see \[Sornette et al., 1996\]). We put the sample into descending order $`A_1A_2\mathrm{}A_N`$. The probability density function (PDF) of the $`n`$-th rank, denoted $`\varphi _{n,N}(x)`$, is well-known

$$\varphi _{n,N}(x)=(Nn+1)\left(\genfrac{}{}{0pt}{}{N}{n}\right)F^{Nn}(x)(1F(x))^{n1}f(x),$$
(1)
where $`F(x),f(x)`$ are the distribution function (DF) and PDF of the random values in question. Putting the Pareto law $`F(x)=1(a/x)^\mu `$, $`xa`$, into (1), we get

$$\varphi _{n,N}(x)\left(1(a/x)^\mu \right)^{n1}x^{\mu (Nn+1)1}.$$
(2)
The mode $`M_{n,N}`$ of the PDF (2) (i.e., the maximum of the PDF) is the most probable value of the random variable $`A_n`$:

$$M_{n,N}=a\left(\frac{N\mu +1}{n\mu +1}\right)^{1/\mu }.$$
(3)
Besides, an interval around the mode containing some prescribed probability (say, $`90\%`$) can be derived from the density $`\varphi _{n,N}(x)`$. The dependence (3) with $`\mu =1/4`$ is shown as the long-dashed line in Fig. 1. The two short-dashed lines represent the upper/lower limits of the $`90\%`$-confidence interval. The data are well accounted for by the power law, except, perhaps, for the three smallest ranks, i.e. the three largest plates, the Pacific, the Antarctica, and the Africa plates, which fall outside the confidence interval.

### “Finite-size” constraint on the PDF of plate sizes

From a visual inspection of figure 1, it might be argued that the deviation from the power law prediction (3) occurs somewhat earlier, say, at rank $`n=7`$, i.e., the seven largest plates with area more than $`1`$ steradian belong to a different population than the rest of the plates. The simplest possibility is that the largest plates are sensitive to a “finite-size” effect, namely that the sum of areas over all plates must sum up to $`4\pi `$. The largest continent and ocean plates are commensurable with the size of the convection cells in the upper mantle and they would thus feel the constraints due to the finite size of the Earth surface. In contrast, one could argue that the smaller plates appeared as a result of the interaction and collision between the larger plates and do not feel any superficial restriction. More precisely, the smaller plates may be the result of fragmentation/creation processes in the neighborhood of the borders and especially near triple points of the seven largest plates.

The distribution $`g_{N,C}(x)`$ of sample values $`x_1,\mathrm{},x_N`$ conditioned by the constraint

$$S_N=x_1+\mathrm{}+x_N=C,$$
(4)
where $`C`$ is a constant ($`4\pi `$ for the plates) is modified from its unconditional Pareto density expression $`f(x)=\mu a^\mu /x^{1+\mu }`$ for $`xa`$ and $`f(x)=0`$ for $`x<a`$. The lower threshold (the minimum value) for the plate data is $`a=0.00178`$ steradian and corresponds to the smallest documented plate in \[Bird, 2002\]. Denoting the unconditional density of the sum $`S_k=x_1+\mathrm{}+x_k,k=N1,N`$ by $`s_k(x)`$, we have $`g_{N,C}(x)=s_{N1}(Cx)f(x)/s_N(C)`$, for $`axC`$. Thus, the constraint (4) decreases the unconditional Pareto density $`f(x)`$ by a factor $`s_{N1}(Cx)/s_N(C)`$ which acts as a “taper.”

In order to use the Maximum Likelihood (ML) method for the estimation of the exponent $`\mu `$, we need the vectorial distribution of conditional sample to take into account the interdependence between the different variables (areas of the plates) $`x_1,\mathrm{},x_{42}`$ induced by the constraint (4). The corresponding likelihood function is therefore

$$\varphi (x_1,\mathrm{},x_{42}|\mu )=\frac{\delta (x_1+\mathrm{}+x_{42}4\pi )}{s_{42}(4\pi )}f(x_1|\mu )\mathrm{}f(x_{42}|\mu ).$$
(5)
The resulting ML estimate is $`\mu =0.25\pm 0.05`$. With this value, we generate an artificial sample of 1045 conditional 42-dimensional (42D) vectors with the condition (4) with $`C=4\pi `$. Rank-ordering each of these 1045 vectors, we determine their sample medians $`M_1,\mathrm{},M_{42}`$, where $`M_j`$ is the median of the $`j`$-th rank. These conditional medians are slightly smaller than given by (3) for the unconditional Pareto distribution. The conditional distribution (5) allows us to construct a confidence domain for the 42D random vectors, defined as a “corridor” of the form $`[cM_j;(1/c)M_j]`$, the constant $`c=0.244`$ being chosen such that $`95\%`$ of vectors fall within this corridor. The medians $`M_1,\mathrm{},M_{42}`$ and their corridor are shown in Fig. 2: all samples of the tectonic plates falls within the $`95\%`$ confidence corridor, showing that (in contrast with the “pure” Pareto used in Fig.1) the Pareto model together with the total area constraint (4) accounts satisfactorily for all the data, including the largest plates.

We have performed many other tests which all confirm the goodness of fit of the constrained Pareto distribution (or of its approximation by the truncated Pareto distribution). Using Wilks’ theorem \[Rao, 1965\], we used the so-called Gamma-distribution with complementary cumulative distribution function $`G(x)`$, truncated from both sides $`G(x|a,b)=C(a,b)\frac{e^{ax}}{x^{1+b}}`$ for $`u_0xu_1`$ and find that (1) the introduction of the additional parameter $`a`$ in the exponential does not increase significantly (probability of rejection $`38\%`$) the likelihood as compared with the truncated Pareto and (2) the truncated Pareto distribution is much preferable compared with the truncated exponential distribution (significance $`110^8`$). This conclusion is strongly confirmed by tests using so-called sufficient statistics for both the power law and exponential models. We also tested the popular lognormal DF, which is a natural model for fragmentation since it is the limit distribution for processes of subdivision of rock under rather general assumptions (basically multiplicative processes) \[Kolmogorov, 1941\]. Using the Kolmogorov test and the bookstrap method, we estimate that the data deviates significantly from the lognormal DF at the $`95\%`$ significance level.

## Discussion

Our main result is that Bird’s suggestion  of a power-like distribution of plates areas is confirmed up to the largest plates, when taking into account the very natural constraint that the plates areas must sum up to $`4\pi `$ steradians.

Given this strength of the validity of the Pareto distribution for the $`42`$ plates documented in \[Bird, 2002\], one can expect that it can be extrapolated beyond this range to smaller yet undetected plates. Using the complementary Pareto distribution $`(a/x)^\mu `$ with $`a=0.00178`$ steradian and $`\mu =0.25`$, this extrapolation predicts $`100`$ plates larger than $`1.3510^4`$ steradian ($`5,500`$km2 or $`74\times 74`$ km) \[Bird, 2002\]. The total area needed to define another 58 plates ranging from this size up to the smallest size in PB2001 would be only about $`0.028`$ steradians, which could be taken from large plates like EU and NA without materially affecting their areas. As discussed by Bird , the places where additional small plates are most likely to be recognized is within the zones of distributed deformation identified in \[Bird, 2002\], which have total area of $`0.838`$ steradians (6.7% of Earth).

The main suggestion inspired by Bird’s discovery is that the separation of scales between a plate tectonics involving a few major plates at large scales and tectonic deformation at smaller scales may be an artifact that resulted historically from insufficient statistics and more recently from a not quite correct interpretation of “visual” deviations from the pure power law. With the improved extended compilation of \[Bird, 2002\] which will probably still improve in the years to come, we conjecture that there is a hierarchy of plates and block sizes from the largest ones to the smallest scales, which may be the result of the self-organization of the plate motions, creation and destruction over millions of years. In other words, we suggest that our view of plate tectonics should be revised in terms of a dynamical model of plates with creation, fragmentation and destruction that may occur at all scales. This view may provide a link between the large organization of major plates down to the block structures at scales of tens of kilometers or less discussed in the introduction.

What could be the mechanism(s) at the origin of the observed power law distribution of plate areas? We propose a fragmentation model with source terms, which can be progressively enriched by taking into account the kinematic compatibility constraints between plates. There is a large class of such fragmentation models that can be solved explicitely within the scaling theory of \[Cheng and Redner, 1988\]. In a nutshell, their approach is to write down a linear integro-differential coupled equation for the time evolution of fragment populations. The fragment distribution $`p(l)`$ is characterized by the moments $`m_\alpha =_0^{\mathrm{}}l^\alpha p(l)𝑑l`$. Note that calculating a moment corresponds to taking the Mellin transform of $`p(l)`$. The moments $`m_\alpha `$ give access to $`p(l)`$ by taken the inverse Mellin transform. The equations of evolution for the fragment population then transform into the following recurrence equation for the moments: $`m_{\alpha +\beta }=\omega \frac{1\alpha }{L_\alpha 1}m_\alpha `$, where $`\beta `$ is the exponent describing the overall rate of breakup (assumed to be proportional to $`l^\beta `$), $`\omega `$ is a normalizing constant and $`L_\alpha =_0^1x^\alpha b(x)𝑑x`$, where $`b(x)`$ is the probability that the fragment ratio be $`x`$. From the factor $`\frac{1}{L_\alpha 1}`$, if there is a value $`\alpha ^{}`$ such that $`L_\alpha ^{}=1`$, then all moments with $`\alpha >\alpha ^{}`$ will become infinite. Provided reasonable analyticity conditions hold, it follows that the value $`\alpha ^{}`$ is a pole of $`m_\alpha `$. Taking the inverse Mellin transform of $`m_\alpha `$ then allows us to get $`p(l)`$ and, using the existence of the pole at $`\alpha ^{}`$, this immediately predict that $`p(l)l^{(1+\beta +\alpha ^{})}`$ is a power law distribution. This model was used by Ouillon et al.  to account for the power law distribution of joint spacing in a granitic massif in Saudi Arabia, with exponent $`\mu _s=1/2`$. Note that block sizes are proportional to spacing taken to the second power, which then predicts $`\mu =\mu _s/2=1/4`$, which is in agreement with our best estimate $`\mu =0.25`$ for the plate area distribution. This agreement between two very different systems may be fortuitous but is worth keeping in mind, especially in view of our claim of a continuity of the physics of fragmention from the large plates to the blocks. We note also that the fragmentation model of \[Cheng and Redner, 1988\] predicts generically log-periodic corrections to the main power law since, for a large variety of $`b(x)`$, there are solutions of $`L_\alpha ^{}=1`$ with complex exponents $`\alpha ^{}`$. This has been shown to signal the existence of discrete scale invariance \[Ouillon et al.. 1996\]. It is tempting to interpret in this way the three approximately regular oscillations observed in figure 1, with a prefered scaling ratio approximately equal to $`10`$. This value is in good agreement with the square of the prefered scaling ratio $`3.5`$ found in the discrete hierarchy of linear sizes in fracture of rocks by Sadovskiy et al. .

Our present conclusion on a universal fractal character of plate fragmentation does not prevent the existence of some underlying distinctions. For instance, P. Bird (private communication and \[Bird,2001\]) proposes three main tectonic origins for plates: plume-assisted rifting for the larger plates (ranks 1-8 in Figs. 1 and 2), continental collision stress for the intermediate plate sizes (ranks 8-20), and back-arc spreading for the smallest plates (ranks 21-42). In Figs. 1 and 2, one can discern slightly different patterns of the tail behavior of these three subgroups. However, any formal statistical analysis of these possible distinctions would be, to our opinion, excessive due to the extreme smallness of the data set.

Acknowledgments. We acknowledge stimulating discussions with P. Bird and Y.Y. Kagan and thank P. Bird warmly for sharing his data.