text:
# The Effect of the (2+1)-Dimensional Gravity in the Spiral Galaxy

## I Introduction

It is known that our Milky Way looks like a thin-disk of shining stars with a diameter of roughly $`8.5`$ to $`10kpc`$ across . The thickness of the Milky Way is roughly of the order of $`100`$ to $`1000pc`$. In fact, the diameter of the shining part for a typical spiral galaxy is roughly $`10^6`$ light year across. This thin plane structure has a few problems needed to be resolved. One of the issue is the possible existence of the dark matter  in order to account for its nearly asymptotically constant velocity along its spiral arms.

Indeed, the eccentric force of the farthermost stars in the spiral arms will require a huge interior mass to keep them on track in the family of the galaxy. It is known that the $`1/r^2`$ gravitational force will require a total mass $`M(r)=v^2r/G`$ interior to the test star (with mass $`m`$) at radial distance $`r`$ in order to balance the eccentric force on the test star along its spiral arms. Here $`r`$ is the radial coordinate in the spherical coordinates.

In fact, when the estimate was extended to distances beyond the point where the light from a galaxy ceases, it was found that $`M(r)`$ continues to increase. On the other hand, the observed mass only provides roughly one percent or less . This result is the conclusion based on a $`1/r^2`$ gravitational force. Hence it is apparent that most of the mass source in our galaxy is dark in a $`(3+1)`$-dimensional spacetime. This is also true for the entire universe, namely, the observed universe is closed to $`\mathrm{\Omega }=1`$ Friedmann-Robertson-Walker space with most of its mass provider $`dark`$ throughout the universe . For a brief introduction of the dark matter please see, for example, chapter one of Ref. . There are a number of candidates for the prescribed dark matter including the recently discovered massive neutrino detected at SuperK .

One notes that the effect of the $`(3+1)`$-dimensional dark matter in the spiral galaxy can be simulated by a $`(2+1)`$-dimensional fluid model. This is because that the gravitational force for a point mass is proportional to $`1/\rho `$ in a space with two spatial dimensions. Here $`\rho `$ denotes the radial coordinate in polar coordinates. In other words, dark matter may not be needed if there is a $`\mathrm{ln}\rho `$ gravitational potential in effect in the scale of the order of the size of a typical galaxy.

Note that the $`(2+1)`$-dimensional fluid model will be introduced in section II. It will be shown that the pressure term present in the fluid model provides an expelling force that counteracts the increasing gravitational force due to the increasing mass along the spiral arms. In addition, isolated stars beyond the typical scale of the galaxy do not act as fluid any more. Hence there does not exist a pressure term beyond that scale any more. Fortunately, a $`2`$-dimensional $`\mathrm{ln}\rho `$ potential is enough to hold remote stars on track without requiring any additional interior mass. Hence this model also explains the reason why there is no need for the dark matter beyond the shining region of the spiral galaxy.

In the following section, we will also show that spiral arm is a stable configuration for the $`(2+1)`$-dimensional fluid model without dark matter similar to the $`(3+1)`$-dimensional fluid model with additional dark matter for a wide range of pressure and rotational velocity. Another way to interpret this result is that the $`2`$ spatial dimensional effect is induced by the $`3`$ spatial dimensional dark matter in a nontrivial way. Some conclusions will be drawn at the end of this paper.

## II Stability of a $`(2+1)`$-dimensional spiral galaxy

Note that it was assumed that the mass density takes the form $`\rho (x)=\sigma (\rho ,\phi )\delta (z)`$ for a $`(3+1)`$-dimensional fluid model in Ref. . We will assume a mass density of the form $`\sigma (x)=\sigma (\rho ,\phi )`$ for a $`(2+1)`$-dimensional fluid model in this paper. The Poisson equation reads

$$\frac{1}{\rho }_\rho (\rho _\rho V)+\frac{_\phi ^2}{\rho ^2}V=\frac{4\pi }{3}G\sigma .$$
(1)
Furthermore, one has

$`_t\sigma +{\displaystyle \frac{_\rho }{\rho }}(\rho \sigma u)+{\displaystyle \frac{_\phi }{\rho ^2}}(\sigma j)`$ $`=`$ $`0,`$ (2)
$`\sigma (_tu+u_\rho u+{\displaystyle \frac{j_\phi }{\rho ^2}}u{\displaystyle \frac{j^2}{\rho ^3}})`$ $`=`$ $`_\rho \mathrm{\Pi }\sigma _\rho V,`$ (3)
$`\sigma (_tj+u_\rho j+{\displaystyle \frac{j_\phi }{\rho ^2}}j)`$ $`=`$ $`_\phi \mathrm{\Pi }\sigma _\phi V`$ (4)
from the conservation of energy and momentum in this $`(2+1)`$-dimensional fluid model. Here the pressure term is parameterized as $`\mathrm{\Pi }\sigma v_t^2`$ with $`v_t`$ the averaged thermal speed. It is known that the pressure term has to be imposed as an external source . This term will play an important role for the stability of the $`(2+1)`$-dimensional fluid model for the spiral galaxy.

Assuming that $`\mathrm{\Pi }=\mathrm{\Pi }(\sigma )`$, one can expands the pressure term perturbatively following the expansion of $`\sigma =\sigma _0+\sigma _1`$. In addition, one will write $`v_t(\sigma _0)=v_0`$. Here $`|\sigma _1||\sigma _0|`$ denotes a perturbation of density around the background mass density $`\sigma _0`$. Furthermore, $`ud\rho /dt`$ and $`j\rho ^2d\phi /dt=\rho ^2\omega (\rho )`$ denote the radial velocity and angular momentum per unit mass respectively.

First assuming that there exists a stable background configuration $`(\sigma ,u,j,V)=(\sigma _0(\rho ),0,\rho ^2\omega (\rho ),V_0(\rho ))`$, one can perform a small perturbation around this background. One will also assume that the perturbed fields take the following form: $`(\sigma _{}(\rho ),u_{}(\rho ),j_{}(\rho ),V_{}(\rho ))e^{i\chi (\rho ,\phi ,t)}.`$ Here $`\chi (\rho ,\phi ,t)=\mathrm{\Phi }(\rho )+\omega _0tm\phi `$ such that the existence of such kind of stable solution represents the existence of a stable spiral galaxy with $`m`$ arms.

Note however that $`\omega _0`$ is a complex angular velocity such with an imaginary part representing a decreasing or increasing factor of the perturbed field. In addition, the WKBJ approximation $`|\mathrm{\Phi }(\rho )^{}||(\mathrm{ln}V_{})^{}|`$ will be our standard approach . This simply means that the variational dependence of $`\rho `$ in $`V_{}`$ came mostly from the phase factor $`\mathrm{\Phi }(\rho )`$.

One also assumes that  $`O(j_{}/\rho )O(u_{})`$, $`O(\sigma _0/\rho )O(_\rho \sigma _0)`$ and $`O(_\rho \sigma _{})O(\sigma _{}/\rho )`$ in addition to the WKBJ condition $`|\mathrm{\Phi }(\rho )^{}||(\mathrm{ln}V_{})^{}|`$ in solving the equations of motion. The consistency of the system requires that the following equation

$$(\omega _0m\omega )^2=k^2v_0^2+\mathrm{\Omega }^24\pi G\sigma _0/3$$
(5)
remains valid from solving the system equations. Here $`kd\mathrm{\Phi }/d\rho `$ in above expression.

Therefore the local stability of the spiral solution requires that

$$Q^21\frac{Z^2}{v_0^2}$$
(6)
from the consistency of the system equations. Here $`Q^2k^2/k_J^2`$, $`Z^2\mathrm{\Omega }^2/k_J^2`$ with $`k_J^24\pi G\sigma _0/(3v_0^2)`$. Therefore local stability is established whenever the constraint (6) is obeyed. Hence the stability of this $`(2+1)`$-dimensional fluid model can be easily achieved similar to the $`(3+1)`$-dimensional fluid model.

## III Conclusion

We have shown that the dark matter in the Milky Way may not be required at all if a galaxy exhibits a $`\mathrm{ln}\rho `$ potential. One hence goes further to show that the $`(2+1)`$-dimensional fluid model (without dark matter) does support stable spiral arms similar to the $`(3+1)`$-dimensional fluid model (with dark matter). One also points out that the effective $`(2+1)`$-dimensional galaxy is probably nothing more than an effective theory induced by the effect of the dark matter in the $`(3+1)`$-dimensional fluid model.

Note however that most spiral galaxies are not exactly $`(2+1)`$-dimensional objects. Hence the $`(2+1)`$-dimensional toy fluid model introduced here is only a simplified version. One should consider a more realistic $`(2^++1)`$-dimensional fluid model in order to account for the formation of the physical galaxies. The simplified model serves, however, to indicate that the dark matter can as well be replaced by other alternatives. Possible applications to the physics in the scale of the entire universe may also reveal more information indicated in this paper. In addition, if the dark matter can not be found, one should pay more attention to the questions such as: what the large-scale physical dimensions are, and what dictates the physical process of the dimensional-reduction mechanism.

###### Acknowledgements.

This work is supported in part by the National Science Council under the contract number NSC88-2112-M009-001.


format: text

text:
# An exact solution for 2+1 dimensional critical collapse

## I Introduction

Critical gravitational collapse at the threshold of black hole formation, as first found by Choptuik, has been studied in many systems. With the exception of a study of vacuum, axisymmetric collapse, the systems studied are spherically symmetric. Because of this symmetry, the equations describing the collapse are PDEs for functions of two variables. The critical solution is often discretely self-similar (DSS) or continuously self-similar (CSS). In the CSS case one can study the critical solution itself by assuming the CSS symmetry and thus reducing the collapse equations to a set of ODEs. In general, the equations (both the PDEs describing collapse and the ODEs describing a CSS critical solution) are sufficiently complicated that a numerical treatment is needed.

The Einstein equations in 2+1 dimensions are simpler than in 3+1 dimensions. Though to form a black hole in 2+1 dimensions, one must add a cosmological constant, this added feature still allows analytic treatment of some collapse situations. One might then hope that critical collapse in 2+1 dimensions would be more tractable. Indeed, for the collapse of thin dust rings or the collision of point particles the collapse can be treated analytically. (This is essentially because the spacetime has constant curvature outside of the zero thickness sources).

Recently, Pretorius and Choptuik performed numerical simulations of the collapse of a massless, minimally coupled scalar field with a cosmological constant in 2+1 dimensions. They find that the critical solution is CSS.

In this paper, we find the critical solution of reference in closed form. In section 2 we write the Einstein-scalar equations in an appropriately chosen coordinate system. In section 3 we make a CSS ansatz and find the solution. This solution is compared in section 4 to the numerical results of reference and perturbations of the solution are considered in section 5.

## II Field equations

The Einstein-scalar equations with cosmological constant are

$$G_{ab}+\mathrm{\Lambda }g_{ab}=\kappa T_{ab}$$
(1)
Here, $`G_{ab}`$ is the Einstein tensor, $`T_{ab}`$ is the stress energy of the scalar field $`\varphi `$

$$T_{ab}=_a\varphi _b\varphi \frac{1}{2}^c\varphi _c\varphi g_{ab}$$
(2)
and $`\kappa `$ is a constant. Following the conventions of reference we chose units such that $`\kappa =4\pi `$.

We now consider an appropriate choice of coordinate system for the metric. Since we want to study the CSS critical solution, we want a coordinate system in which the solution appears manifestly CSS. This is not the case for the coordinates of reference . Instead we use the method of Christodoulou to choose a coordinate system where the coordinates are geometric quantities. We choose as a radial coordinate $`\overline{r}`$ such that $`2\pi \overline{r}`$ is the length of the circles of symmetry. This is the analog of the usual area coordinate used in spherical symmetry. We choose as our “time” coordinate, the null coordinate $`u`$ defined as follows: $`u`$ is constant on outgoing light rays, and on the world line of the central observer, $`u`$ is equal to the proper time of that observer. Finally, we choose a coordinate $`\theta `$ so that $`/\theta `$ is the Killing vector. The metric then takes the form

$$ds^2=e^{2\nu }du^22e^{\nu +\lambda }dud\overline{r}+\overline{r}^2d\theta ^2$$
(3)
where $`\nu `$ and $`\lambda `$ are functions of $`u`$ and $`\overline{r}`$. In reference the Einstein-scalar equations were found for spherical symmetry in any number of dimensions. For our purposes, we specialize the results of reference to 2+1 dimensions, generalize them to add a cosmological constant and change the convention to $`\kappa =4\pi `$. We begin by introducing null vectors $`l^a`$ and $`n^a`$ defined by

$$l^a=e^\lambda \left(\frac{}{\overline{r}}\right)^a$$
(4)
$$n^a=e^\nu \left(\frac{}{u}\right)^a\frac{1}{2}e^\lambda \left(\frac{}{\overline{r}}\right)^a$$
(5)
Then the Einstein-scalar equations with cosmological constant are satisfied provided that the scalar field satisfies the wave equation and that the following components of Einstein’s equation are satisfied

$$G_{ab}l^al^b=4\pi \left(l^a_a\varphi \right)^2$$
(6)
$$G_{ab}l^an^b=\mathrm{\Lambda }$$
(7)
Equations (6) and (7) become

$$\frac{1}{\overline{r}}e^{2\lambda }\frac{}{\overline{r}}(\lambda +\nu )=4\pi e^{2\lambda }\left(\frac{\varphi }{\overline{r}}\right)^2$$
(8)
$$\frac{1}{2\overline{r}}e^{2\lambda }\frac{}{\overline{r}}(\nu \lambda )=\mathrm{\Lambda }$$
(9)
The solution of equations (8) and (9) is most easily expressed by defining the quantities $`ge^{\nu +\lambda }`$ and $`\overline{g}e^{\nu \lambda }`$. Then we have

$$g=\mathrm{exp}\left[4\pi _0^{\overline{r}}\overline{r}\left(\frac{\varphi }{\overline{r}}\right)^2𝑑\overline{r}\right]$$
(10)
$$\overline{g}=12\mathrm{\Lambda }_0^{\overline{r}}\overline{r}g𝑑\overline{r}$$
(11)
The wave equation for $`\varphi `$ becomes

$$2\frac{^2\varphi }{u\overline{r}}+\frac{1}{\overline{r}}\frac{\varphi }{u}\frac{1}{\overline{r}}\frac{}{\overline{r}}\left(\overline{r}\overline{g}\frac{\varphi }{\overline{r}}\right)=0$$
(12)
## III Critical solution

We now make the ansatz that the scalar field is CSS and use this to solve the field equations. Choose the origin of $`u`$ to be at the singularity, and define two new coordinates $`T`$ and $`R`$ by

$$ue^T$$
(13)
$$\overline{r}e^TR$$
(14)
Then demand that $`\varphi `$ take the form

$$\varphi =cT+\psi (R)$$
(15)
where $`c`$ is a constant. This ansatz requires that we neglect the cosmological constant, which in turn means that $`\overline{g}=1`$ and thus reduces equation (12) to the flat space wave equation. Putting the ansatz in equation (15) into equation (12) we obtain

$$R(12R)\psi ^{\prime \prime }+(13R)\psi ^{}c=0$$
(16)
where a prime denotes derivative with respect to $`R`$. The solution of this equation that is regular at the origin is

$$\psi =2c\mathrm{ln}\left[\frac{1}{2}\left(1+\sqrt{12R}\right)\right]$$
(17)
which leads to a scalar field given by

$$\varphi =c\left(T2\mathrm{ln}\left[\frac{1}{2}\left(1+\sqrt{12R}\right)\right]\right)$$
(18)
Then using equation (10) we find that the metric function $`g`$ is

$`g=\mathrm{exp}\left[4\pi {\displaystyle _0^R}R\left({\displaystyle \frac{\varphi }{R}}\right)^2𝑑R\right]`$ (19)
$`=\left[{\displaystyle \frac{(1+\sqrt{12R})^2}{4\sqrt{12R}}}\right]^{8\pi c^2}`$ (20)
The metric of the CSS critical solution must be singular at $`u=0`$. However, our critical solution (equations (18) and (19)) appears to have an additional singularity at $`R=1/2`$ which is the past light cone of the $`u=0`$ singularity. We now consider whether the apparent singularity at $`R=1/2`$ is a real singularity or a coordinate singularity. Note that from equation (19) it follows that $`R=1/2`$ (for any value of $`T`$) is a marginally outer trapped surface and that the Christodoulou coordinates go bad at just such surfaces. Now define a new coordinate $`v`$ by $`v=(u+2\overline{r})`$. Then the metric is

$`ds^2=\left[{\displaystyle \frac{u}{16}}\left(1+\sqrt{{\displaystyle \frac{v}{u}}}\right)^4\right]^{4\pi c^2}v^{4\pi c^2}dvdu`$ (21)
$`+{\displaystyle \frac{1}{4}}(u+v)^2d\theta ^2`$ (22)
Now define the number $`q`$ and the coordinate $`w`$ by $`1/(2q)=14\pi c^2`$ and $`w^{2q}=v`$. Then the metric is

$`ds^2=\left[{\displaystyle \frac{u}{16}}\left(1+{\displaystyle \frac{1}{\sqrt{u}}}w^q\right)^4\right]^{1(2q)^1}2qdudw`$ (23)
$`+{\displaystyle \frac{1}{4}}\left(u+w^{2q}\right)^2d\theta ^2`$ (24)
This metric is smooth at $`w=0`$ provided that $`q=n`$ where $`n`$ is a positive integer. That is, the metric is smooth for values of $`c`$ given by

$$c=\pm \sqrt{\frac{1}{4\pi }\left(1\frac{1}{2n}\right)}$$
(25)
For $`n=1`$ the spacetime is the 2+1 dimensional Robertson-Walker metric.

We now consider the question of whether it is physically necessary to impose the condition that the metric be smooth at $`w=0`$. If $`q=n`$, then one can show using equation (23) that the spheres of symmetry are trapped surfaces for $`w<0`$. However, the critical solution cannot have trapped surfaces, since it forms the boundary between those evolutions that result in trapped surfaces and those that do not. Therefore, one should expect the numerical critical solution to approach our CSS solution only inside the past light cone of the singularity (that is for $`w>0`$). It is therefore not physically necessary that the CSS solution be smoothly extendible past $`w=0`$ since such an extension cannot correspond to the behavior of the numerical critical solution.

## IV Comparison with numerical results

In a near-critical collapse, the evolution at first approaches the critical solution, and then diverges from it as the single unstable mode grows. Therefore, in comparing a numerical simulation of near-critical collapse to a proposed critical solution, one should make the comparison at an intermediate time: late enough for the critical solution to be approached, but early enough so that the unstable mode does not have appreciable amplitude.

To compare the analytic and numerical results, one must express both in the same coordinate system. In the coordinates of reference the metric takes the form

$$ds^2=\frac{e^{2A}}{\mathrm{cos}^2(r/l)}(dr^2dt^2)+l^2\mathrm{tan}^2(r/l)e^{2B}d\theta ^2$$
(26)
where $`l=\sqrt{1/\mathrm{\Lambda }}`$ and $`A`$ and $`B`$ are functions of $`r`$ and $`t`$. The coordinates of reference are related to ours by $`\overline{r}=l\mathrm{tan}(r/l)e^B`$ and $`u`$ is that function of $`tr`$ that at $`r=0`$ is equal to the proper time of the central observer. Therefore, it is fairly straightforward to take scalars and tensors in the coordinates of reference and express them in our coordinate system.

Our solution has an unknown parameter $`c`$, which should be chosen for best fit with the numerical results. To make this choice, we note that the analytical solution (equation (18)) satisfies $`c=\varphi /T|_{R=0}`$. Thus the values of $`\varphi /T|_{R=0}`$ for the numerical solution should allow us to determine the value of $`c`$. Figure 1 shows a plot of $`\varphi /T|_{R=0}`$ vs $`T`$ for the numerical solution. Note that there is a range of intermediate times for which this quantity is approximately constant. The line in the solution corresponds to the value of $`c`$ given by equation (25) with $`n=4`$. While it is not necessary that $`c`$ be given by equation (25), we find that this particular value of $`c`$ gives excellent agreement with the numerical data. From now on, we will assume that $`c`$ has this value.

Figure 2 shows a comparison between the numerical and analytic results for the scalar field $`\varphi `$. Here the dots are the numerical results and the curve is the analytic one. The comparison is made at $`T=9`$. The freedom to add a constant to the scalar field is used to set the value of the scalar field to zero at the origin. Note that there is excellent agreement between analytic and numerical results. Though this figure shows the comparison at only one time, the agreement persists for a large range of intermediate times.

Due to the symmetry of the spactime, the metric is determined entirely by the matter. This is made explicit in equations (10) and (11). Therefore, since the scalar fields of the numerical and analytic solutions agree, the metrics must agree also. Nonetheless, for illustrative purposes we present a comparison of the metrics of the numerical and analytic solutions. Figure 3 shows a comparison of the quantity $`M=(\overline{r}^2/l^2)(\overline{g}/g)`$ with the numerical solution given by the dots and the analytic solution given by the curve. Note the excellent agreement between the two solutions. Here the comparison is made at $`T=9`$, but the agreement persists for a large range of intermediate times. The metric is determined by $`M`$ and $`\overline{g}`$ which is 1 for the analytic solution and very close to 1 for the numerical solution in the intermediate range of times. Thus, in the intermediate range of times, the scalar field and metric of the numerical and analytic solutions agree.

In figures 2 and 3, we have performed the comparison of analytic and numerical solutions only up to the past light cone of the singularity. However, the numerical solution certainly continues beyond this light cone, and since $`q=4`$, the analytic solution also extends. Do the two solutions still agree in this region? I argued at the end of the previous section that they cannot agree, since the analytic solution contains a closed trapped surface beyond the light cone. Nonetheless, it would be useful to compare the two solutions in this larger region to see the extent and nature of their disagreement.

To make this comparison, a different coordinate from $`R`$ is needed, since the solution is not smooth in $`R`$ across the light cone. Here we use as our new coordinate the affine parameter $`\lambda `$ along an outgoing radial null geodesic. The affine parameter is made scale invariant by choosing the null geodesic to have inner product $`u`$ with the central observer. The analytic solution can be expressed parametrically in terms of $`\lambda `$ as follows: (here we specialize to the case $`q=4`$). Introduce the variable $`x(u)^{1/8}w`$. Then $`x=1`$ at the origin and $`x=0`$ at the past light cone of the singularity. We have

$$\lambda =\frac{1}{2\sqrt{2}}_x^1\left(1+x^4\right)^{7/2}𝑑x$$
(27)
$$\varphi =\sqrt{\frac{7}{8\pi }}\mathrm{ln}\left(\frac{1+x^4}{2}\right)$$
(28)
$$M=x^7\left(\frac{2}{1+x^4}\right)^{7/2}$$
(29)
From equation (27) it follows that the past light cone is at $`\lambda 0.838`$. Using equation (26), the affine parameter $`\lambda `$ can be found numerically from the numerical data. Note, however that the numerical simulations have only a certain range in $`\lambda `$ since a critical numerical solution must stop before the time $`t`$ of singularity formation. In our case, the maximum value of $`\lambda `$ for the numerical data is approximately 1.

Figures 4 and 5 show a comparison between the numerical and analytic solutions for the range $`0\lambda 1`$. Here the values of $`\varphi `$ are compared in figure 4, while the mass aspect $`M`$ is compared in figure 5. As in figures 2 and 3, the comparison is made at $`T=9`$. From the figures it is not clear whether there is a disagreement between the two solutions in the range $`\lambda >0.838`$ that is beyond the past light cone. To see the disagreement clearly, one would need a larger range of $`\lambda `$. This could perhaps be provided by a numerical simulation that used double null coordinates.

## V Perturbations

We now turn to a treatment of perturbations of the critical solution. This treatment should in principle allow us to determine analytically (1) the correct value of the parameter $`c`$ and (2) the value of the exponent in scaling relations for near-critical collapse. Unfortunately, as we will see below, it is unclear what boundary conditions to impose on the perturbations.

The critical solution, when perturbed, has one unstable mode that grows as $`e^{kT}`$ for some constant $`k`$. Therefore the correct value of $`c`$ is the one for which there is exactly one unstable mode. The quantity $`k`$ is related to scaling laws for near-critical collapse. In the collapse of a one parameter family of initial data, a quantity $`Q`$ with dimension $`(\mathrm{length})^s`$ obeys a scaling relation $`Q|pp^{}|^{s/k}`$ where $`p`$ is the parameter and $`p^{}`$ is its critical value. Pretorius and Choptuik examine scaling in maximum scalar curvature for subcritical collapse, a quantity that has dimension $`(\mathrm{length})^2`$. They find $`k0.83`$.

In the approximation that $`\overline{g}=1`$, the scalar field satisfies the flat space wave equation. Therefore, the perturbed scalar field $`\delta \varphi `$ also satisfies this equation. Then making the ansatz $`\delta \varphi =e^{kT}S(R)`$ and using equation (12) we obtain

$$R(12R)S^{\prime \prime }+(1[3+2k]R)S^{}kS=0$$
(30)
The solution is

$$S=F(k,1/2,1,2R)$$
(31)
where $`F`$ is a hypergeometric function. Writing the hypergeometric function in integral form, we have

$$S=\frac{2}{\pi }_0^{\pi /2}\left(12R\mathrm{sin}^2\chi \right)^k𝑑\chi $$
(32)
In principle, one should now impose boundary conditions on the behavior of the perturbation at $`R=1/2`$ and these boundary conditions would then determine the allowed values of $`k`$. Unfortunately, it is not clear what boundary conditions are physically reasonable, since the critical numerical solution should match the analytical solution only for $`R<1/2`$. Therefore, though we have the solution of the perturbation equation, we cannot use this solution to determine the critical exponent.

## VI Acknowledgements

I would like to thank Matt Choptuik and Frans Pretorius for many helpful discussions and for providing the comparison between their numerical work and the exact solution of this paper. This work was partially supported by NSF grant PHY-9722039 to Oakland University.


text:
# Quantum Hall Transition in the Classical Limit

## I Introduction

Integer quantum Hall transition is a localization-delocalization transition experienced by a quantum mechanical particle moving in a two dimensional plane in a magnetic field perpendicular to the plane and in a random scalar potential $`V(x)`$. As the energy of the particle approaches its critical value, the particle becomes more and more delocalized. In spite of the fact that the problem has been studied numerically and experimentally for 20 years, it has resisted exact theoretical descriptions.

In this paper we examine a new method to probe the properties of the quantum Hall transitions recently introduced by Girvin et al . We show that this method provides a natural framework in which the classical limit of quantum Hall transition (classical percolation) can be studied. We show, both numerically and analytically, that even in this classical limit the percolating particle exhibits normal diffusion, just as the quantum critical particle would. The critical exponent $`\nu `$ turns out to be equal to the percolation exponent $`\frac{4}{3}`$. In this framework, the renormalization group flow from percolation to the quantum Hall transition can be studied.

The method introduced in  is based on the introduction of the density-density correlation function

$$G(x,t)=\underset{n}{}n|\rho (x,t)\rho (0)|n$$
(1)
where $`\rho `$ is the density operator and the sum goes over the Hilbert space of states. This correlation function has a simple meaning. It is the probability that a quantum particle, localized at the origin at the initial moment in time, is found at the position $`x`$ at a time $`t`$. Indeed, set up a wave function of a particle localized at the origin,

$$\psi (x)=\underset{n}{}\psi _n^{}(0)\psi _n(x),$$
(2)
where $`\psi _n`$ are the eigenstates of the Hamiltonian of the problem. At time $`t`$ the wave function will be (we set $`\mathrm{}=1`$)

$$\psi (x,t)=\underset{n}{}\psi _n^{}(0)\psi _n(x)e^{iE_nt},$$
(3)
where $`E_n`$ are the eigenenergies. The probability is equal to

$$\left|\psi (x,t)\right|^2=\underset{m,n}{}\psi _n^{}(0)\psi _n(x)\psi _m^{}(x)\psi _m(0)e^{i\left(E_nE_m\right)t}.$$
(4)
The matrix elements of the density operator are

$$\rho _{m,n}(x)=\psi _m^{}(x)\psi _n(x).$$
(5)
It is therefore obvious that the correlation function

$$G(x,t)=\underset{m,n}{}\rho _{m,n}(x,t)\rho _{n,m}(0)$$
(6)
is equal to the probability (4).

It sounds reasonable that such a correlation function, which tests how far a particle can go in a given amount of time, is useful when studying localization. In the localization literature, however, a different function is usually studied. It is called the spectral function and is given by the square of the single particle Green’s function

$$S(x,\omega ;E)=\underset{m,n}{}\frac{\psi _m^{}(x)\psi _m(0)\psi _n(x)\psi _n^{}(0)}{\left(E\frac{\omega }{2}iϵE_n\right)\left(E+\frac{\omega }{2}+i\delta E_m\right)},$$
(7)
with $`ϵ`$ and $`\delta `$ infinitesimally small numbers. As is well known, the conductance of the system is simply related to the spectral function $`S(x,\omega ;E)`$. This function contains more information than the correlation $`G`$. It tells us how far a particle can go if it propagates along the state with a given energy $`E`$. The correlation function $`G`$ is related to the integral of $`S`$ over $`E`$,

$$2\pi G(x,\omega +iϵ+i\delta )=_{\mathrm{}}^{\mathrm{}}𝑑ES(x,\omega ;E),$$
(8)
where $`G(x,\omega )`$ is a Fourier transform of $`G(x,t)`$. The last formula can be obtained by directly integrating (4) and (7).

One could argue that the density-density correlator $`G(x,t)`$ we just introduced would not give much information about the localization properties. Indeed, in a transition such as Anderson transition in 3D, this function will always be long ranged due to the contributions of the delocalized states to (8). In contrast, the spectral function $`S(x,\omega ;E)`$ tests only the states at the energy $`E`$.

The insight of Girvin at al was that in quantum Hall transition only the centers of the Landau bands are believed to be delocalized. Moreover, the physics of quantum Hall transitions allows us to restrict the particle to the lowest Landau level. In practical terms, that means restricting the sum in (6) to go only over the states in the lowest Landau level. Then there is only one delocalized state. In (8) only the states whose localization length is larger than $`x`$ contribute to the integral over all the energies. Consequently the density-density correlation function includes in itself some useful information about quantum Hall transitions.

It follows from the scaling theory of phase transitions that the localization length $`\xi (E)`$ is divergent as we approach the localization transition with the exponent $`\nu `$,

$$\xi (E)|EE_c|^\nu $$
(9)
where $`E_c`$ is the critical energy. It has also been argued  that a quantum particle in 2D at the localization-delocalization transition undergoes normal diffusion, which in particular results in a finite critical conductivity. Heuristically people sometimes estimate the Fourier-transformed density-density correlator by an interpolating method based on (8). Until the particle hits the localization length, we assume that

$$S(k,t;E)\mathrm{exp}\left(Dt\frac{k^2}{4}\right),\mathrm{if}Dt\xi ^2(E),$$
(10)
which describes diffusion with the diffusion constant $`D`$. This assumption is true only if $`k0`$ . The diffusion at small $`k`$ reflects the existence of critical conductivity in quantum Hall effect. On the other hand, at larger $`k`$ the motion of the particle becomes more complicated. That reflects the fact that critical conductance is not Ohmic and in fact dependent on the shape of the sample. Obviously, (10) will not work if the particle is deeply in the localized regime, when $`|EE_c|E_{\mathrm{max}}`$, with $`E_{\mathrm{max}}`$ describing the scale where (10) and (9) break down.

If we wait long enough,

$$S(k,t;E)\mathrm{exp}\left(\frac{k^2}{4}\xi ^2(E)\right),\mathrm{if}Dt\xi ^2(E).$$
(11)
Roughly speaking, $`G(x,t)`$ can be estimated as

$$G(k,t)_{E_cE_{\mathrm{max}}}^{E_c+E_{\mathrm{max}}}𝑑E\mathrm{exp}\left(\frac{k^2}{4\left(\frac{1}{Dt}+\frac{1}{\xi ^2(E)}\right)}\right),$$
(12)
where the function being integrated interpolates between the limiting behavior (10) and (11). Using (9) and changing the variables $`E=E_c+zt^{\frac{1}{2\nu }}`$ , we can estimate this integral to go as at

$`G(k,t)=t^{\frac{1}{2\nu }}{\displaystyle _0^{t^{2\nu }E_{\mathrm{max}}}}𝑑z\mathrm{exp}\left({\displaystyle \frac{k^2t}{4\left(D^1+z^{2\nu }\right)}}\right)=`$ (13)
$`E_{\mathrm{max}}t^{\frac{1}{2\nu }}{\displaystyle _0^{t^{2\nu }E_{\mathrm{max}}}}𝑑z\left[1\mathrm{exp}\left({\displaystyle \frac{k^2t}{4\left(D^1+z^{2\nu }\right)}}\right)\right]1t^{\frac{1}{2\nu }}g(k^2t),`$ (14)
where $`g`$ is some complicated function. We used that the integral on the second line can be extended to infinity. $`G(x,t)`$ is diffusive at small $`k`$ (a function of $`\frac{x^2}{t}`$) with a prefactor which is $`t`$ dependent. Knowing the prefactor allows us to determine the critical exponent $`\nu `$. This formula has been tested numerically in  and was found to give values of $`\nu `$ for the quantum Hall transition around the value accepted by the community.

Imprecise as the derivation of this formula may seem, we will show later in this paper how it can be given a precise meaning in the classical limit.

## II Equation of Motion

Now we would like to demonstrate that, as was first shown in , the density-density correlation function $`G(x,t)`$, when projected to the lowest Landau level and in the presence of a scalar potential $`V(x)`$, satisfies a simple equation.

Choose the basis of states in the lowest Landau level as in (see  and references therein)

$$\psi _n(z,z^{})=\frac{z^n}{\sqrt{\left(2l^2\right)^nn!\pi }}\mathrm{exp}\left(\frac{zz^{}}{4l^2}\right),$$
(15)
where, as always, $`z=x+iy`$, complex coordinates on the 2D plane and $`l`$ is the so called magnetic length. Consider the Fourier transformed density operator

$$\rho _k=\mathrm{exp}\left[\frac{i}{2}\left(kz^{}+k^{}z\right)\right],$$
(16)
where similarly to the complex coordinates $`k=k_x+ik_y`$. Let us project this operator to the lowest Landau level, following the methods developed in . The operator consists of the product of two factors. One of them, $`\mathrm{exp}\left[\frac{i}{2}k^{}z\right]`$ is analytic in $`z`$ and so acts entirely within the lowest Landau level. As such, it does not have to be projected. Indeed, if we multiply one of the wave functions (15) of the lowest Landau level by this factor, the result can be reexpressed as a linear combination of other functions of (15). The other factor involves $`z^{}`$ and therefore does not act within the lowest Landau level. It has to be projected. We note, however, that the operator

$$\widehat{z}=2l^2\frac{}{z}+\frac{1}{2}z^{}$$
(17)
acts within the lowest Landau level,

$$\widehat{z}\psi _n(z,z^{})=\sqrt{2l^2n}\psi _{n1}$$
(18)
Furthermore,

$$d^2z\left(z^{}\right)^n\left(z^{}\right)^sz^m\mathrm{exp}\left(\frac{zz^{}}{2l^2}\right)=d^2z\left(z^{}\right)^n\mathrm{exp}\left(\frac{zz^{}}{4l^2}\right)\left(2l^2\frac{}{z}+\frac{1}{2}z^{}\right)^s\left[z^m\mathrm{exp}\left(\frac{zz^{}}{4l^2}\right)\right].$$
(19)
In other words,

$$n|\left(z^{}\right)^s|m=n\left|\left(2l^2\frac{}{z}+\frac{1}{2}z\right)^s\right|m$$
(20)
This allows us to write down the projected density operator as in

$$\rho _k=\mathrm{exp}\left[\frac{i}{2}k\widehat{z}\right]\mathrm{exp}\left[\frac{i}{2}k^{}z\right].$$
(21)
Normally the density operators commute with each other, but the projected density operators do not. Multiplying two density operators together and using (17) we obtain

$$\rho _k\rho _q=\mathrm{exp}\left[\frac{l^2}{2}\left(kq+ik\times q\right)\right]\rho _{k+q},$$
(22)
where $`kq`$ is the dot and $`k\times q`$ is the cross products of two vectors $`(k_x,k_y)`$ and $`(q_x,q_y)`$.

Finally we note that the Hamiltonian of a particle moving in a lowest Landau level in the presence of the scalar potential $`V(x)`$ is just that potential itself projected to the lowest Landau level,

$$H=𝑑xV(x)\rho (x)=𝑑kV(k)\rho _k.$$
(23)
Therefore, the equation of motion for the operator $`\rho `$ is

$$i\mathrm{}\frac{}{t}\rho _k=[H,\rho _k]=d^2q2i\mathrm{sin}\left(\frac{l^2}{2}k\times q\right)V(kq)\mathrm{exp}\left[\frac{l^2}{2}\left(k^2kq\right)\right]\rho _q,$$
(24)
where we restored the $`\mathrm{}`$ factor. Clearly $`G(k,t)=\mathrm{tr}(\rho _k(t)\rho _k(0))`$ satisfies

$$i\mathrm{}\frac{}{t}G(k,t)=d^2q2i\mathrm{sin}\left(\frac{l^2}{2}k\times q\right)V(kq)\mathrm{exp}\left[\frac{l^2}{2}\left(k^2kq\right)\right]G(q,t),$$
(25)
with the boundary condition

$$G(k,0)=\mathrm{tr}\left(\rho (k)\rho (k)\right)=N\mathrm{exp}\left(\frac{l^2}{2}k^2\right),$$
(26)
N being the number of states in the lowest Landau level. Solving this equation allows us in principle to determine $`G(k,t)`$. Therefore, it encodes all the properties of quantum Hall transitions. Unfortunately, while this equation is simple to write down, solving it is not a trivial task.

## III Classical Limit

Let us consider the classical limit of the quantum Hall transition. The classical limit $`\mathrm{}0`$ is equivalent here to the limit $`l0`$ since $`l^2\mathrm{}`$. In this limit, the equation (25) becomes

$$\frac{G(x,t)}{t}=\frac{l^2}{\mathrm{}}ϵ_{ij}_iV_jG(x,t).$$
(27)
Since $`l^2\mathrm{}`$, the planck constant completely drops out. In what follows, we choose the units where $`\mathrm{}/l^2=1`$.

The formal solution to this equation is straightforward to obtain:

$$G(x,t)=\delta (x(t)x),$$
(28)
where $`x(t)`$ is a trajectory satisfying the following equation

$$\frac{dx_i(t)}{dt}=ϵ_{ij}_jV(x(t)).$$
(29)
Studying the solutions to the equation (27) reduces to studying the trajectories satisfying (29).

The tangent to the trajectory $`x(t)`$ is always perpendicular to the gradient of the potential $`V(x(t))`$. Therefore, $`V(x(t))`$ is conserved along the trajectory. In other words, the particle moves along the equipotential lines, as shown on Fig. 1. The value of the potential along the trajectory becomes an important parameter $`V_0=V(x(t))`$ which can be assigned to each trajectory satisfying (29). In a way, the potential $`V(x)`$ can be thought of as a mountainous landscape with peaks and valleys. As water fills the valleys of the landscape to the height $`V_0`$, the equipotential lines of Fig 1 are the shapes of the lakes formed in this way.

In this form the problem we study becomes equivalent to classical percolation . Indeed, the potential $`V(x)`$ is random and independent in different points in space. We can consider the value of the potential at a point $`x`$ as a random variable. Take a particular trajectory $`x(t)`$ with a particular value of $`V(x(t))=V_0`$. Let us say that the point $`x`$ is “on” if $`V(x)>V_0`$ and “off” otherwise. Then the trajectory surrounds the clusters of the “on” sites. The trajectory becomes the so-called “hull” of the percolation cluster, and the area it surrounds is the percolation cluster itself (or the islands in the landscape analogy, for $`V_0>0`$). In the percolation problem, the site is “on” with a probability $`p`$. In our problem it can be expressed as in

$$p=_{V_0}^{\mathrm{}}P(V)𝑑V$$
(30)
where $`P(V)`$ is the probability for the values of the potential at different points in space (sites). In this way, we recover the classical percolation picture of the quantum Hall transitions . The electron density percolates along the equipotential lines of the scalar potential. The trajectories are closed. However, as $`V_0\left(p\frac{1}{2}\right)`$ approaches zero, the trajectories become more and more extended in space. Denoting by $`\xi `$ the root mean square distance between points on the trajectory, we say that the percolation hulls exhibit the delocalization transition, with $`\xi V_0^\nu `$.

Define the area of a cluster as

$$A=d^2x\theta (V(x)V_0),$$
(31)
$`\theta (x)=1`$ when $`x>0`$ and $`\theta (x)=0`$ when $`x<0`$.

It turns out that the time it takes to go along the boundary of a cluster is just $`dA/dV_0`$, or

$$T=d^2x\delta (V_0V(x)).$$
(32)
This can be derived directly from (29). Alternatively, one can argue that (29) are the Hamiltonian equations, with $`V`$ being the Hamiltonian. Then $`S`$ is the “adiabatic invariant” and its derivative is obviously time.

From here one can conclude that, since the area of a cluster goes as $`AV_0^\gamma `$, the time goes as $`TV_0^{\gamma 1}`$, where $`\gamma =\frac{43}{18}`$ is another well known percolation exponent . However, as we will see below, $`T`$ is not equal to the time it takes for a particle to go around the cluster. Rather it is much bigger having to do with the “swiss cheese” shape of a typical cluster. As a result, $`T`$ sums the time it takes go around the cluster, and around each hole in the cluster.

The classical percolation picture of the quantum Hall transitions was suggested a long time ago. It was subsequently demonstrated that it does not capture the physics of quantum Hall transitions. Indeed, the exponent $`\nu `$ computed in this way coincides with the its value for percolation $`\nu =\frac{4}{3}`$ , while the true value for $`\nu `$ is believed to be around $`\nu =2.35\frac{7}{3}`$. We can say that the classical density-density correlation function $`G(x)`$ forgets about the phase a quantum particle carries with itself. It was also known that the classical percolation picture can be derived as a classical limit of the Chalker-Coddington model , which is believed to lie in the same universality class as quantum Hall transitions.

Some time ago, F. Evers undertook a thorough numerical study of the solutions to the equations (29) at zero energy $`V_0`$, motived by classical limit of quantum Hall transition, in the paper . Its results will be useful for us below.

Now we are going to study percolation as a classical limit of quantum Hall transition. The approach we are using allows us to give a precise definition of what percolation is for quantum Hall transition, which quantum Hall quantities percolate in the classical picture and finally how they change as we allow for finite $`l`$ to account for the quantum effect.

## IV Percolation

To study the percolation limit, we suggest a lattice model of the equation (29). As a possible lattice model, we could have taken the classical limit Chalker-Coddington model . This is the version of the model where a particle turns left or right at each node with a certain probability, instead of a certain quantum probability amplitude as in the original work . However, to make a closer connection with percolation, we would like to consider a somewhat different model.

Consider a triangular lattice. Associate with each of its sites a random number $`V`$ drawn from a probability distribution $`P(V)`$. Now associate the points of the dual lattice with the path $`x(t)`$. Each path is assigned a number $`V_0`$, the value of the energy on this particular path. At each step the path moves in such a way that the potential of the left node of the bond the particle crosses is greater than $`V_0`$ which in turn is greater than the potential of the right node. For example, the particle on Fig 2 enters the triangle through the horizontal bond, as shown by the solid line. Therefore, $`a>V_0>b`$. If $`c>V_0`$, it would leave the triangle through the right bond, as shown by the dashed line. If $`c<V_0`$, then it would leave the triangle through the left bond.

The trajectory then becomes the hull of percolation clusters. A typical hull is shown on Fig 3.

The equation (27) also contains dynamical information. The velocity of the particle which follows (29) is proportional to the gradient of the potential. We model it by saying that the time it takes for a particle to traverse a bond of the triangular lattice is inversely proportional to the potential difference at the ends of the bond. For example, the time it takes to traverse the bond of the Fig. 2 (the process denoted by the dashed line) is $`1/(cb)`$.

Using this model, we performed extensive numerical simulations. In the simulations we assumed that the values of the potential at the nodes of the triangular lattice are random numbers uniformly distributed on the interval $`[1,1]`$. The lattice size was 8500 by 8500, and we discarded the trajectories which went outside the lattice. The typical trajectories were, however, of much smaller size, about 500. We discovered that the time it takes for a particle to go along a trajectory is essentially proportional to the trajectory’s length. This is because the velocity of the trajectory remains roughly the same as $`V_0`$ is varied, having no singularity or zero as $`V_0`$ passes through 0. Indeed, it is easy to estimate the average velocity by assuming that all the potentials at the corners of a triangle are random independent numbers. Then all we need to do is to average $`1/(V_1V_2)`$ where $`V_1`$ and $`V_2`$ are the values of the potential at the ends of the bond the trajectory is traversing. For the potential uniformly distributed on the interval $`[1,1]`$ we obtain that the average velocity is very weakly $`V_0`$ dependent, with its value being $`\mathrm{log}(4)`$ when $`V_0=0`$. This was confirmed numerically. In fact, that coincides with the earlier observation in .

Interestingly, that allows us to find the time $`T`$ it takes for a particle to go around the cluster. The length of a trajectory is known to go as

$$LV_0^\gamma ^{},$$
(33)
where $`\gamma ^{}=2`$ . Since the average velocity is only weakly $`V_0`$ dependent, we obtain that $`TV_0^2`$. That should be contrasted with the time it takes to go around the cluster and around each hole in the cluster which we found in the previous section.

Replacing time with length we studied $`\overline{x^2(t)}`$ where the average means averaging over different trajectories with the same value of $`V_0`$. In our simulations we averaged over about $`10^6`$ trajectories. Fig. 4 summarizes the results of the numerics.

It is clear that $`\overline{x^2(t)}`$ is a straight line

$$\overline{x^2(t)}=Dt,tt^{}$$
(34)
until a certain time $`t^{}`$ which varies with $`V_0`$. After that,

$$\overline{x^2(t)}=\xi ^2,tt^{},$$
(35)
with the parameter $`\xi `$ depending on $`V_0`$ but not on time.

The length $`\xi `$ is none other than the localization (correlation) length of the percolation problem we introduced above. Indeed, if we wait long enough, the particle satisfying (29) starts winding around the closed trajectories such as on Fig. 1. Then $`\overline{x^2(t)}`$ measures the mean square distance between two points on the percolation cluster’s hull which is the trajectory we study. So we can write down

$$\xi V_0^{\frac{4}{3}},$$
(36)
which we also confirmed by our numerical simulations.

The numerical result that $`\overline{x^2(t)}t`$ for $`t<t^{}`$ shows that the particle undergoes a normal diffusion. While this is what is expected of a quantum particle in a quantum Hall transition , it is far less obvious that a classical particle should do the same. It is also not obvious how to explain this from a purely percolation point of view.

This result was previously obtained in  where $`\overline{x^2(t)}`$ for the trajectories at $`V_0=0`$ was studied. There it was shown that the diffusion was due to a “mysterious” cancellation of exponents in percolation. Perhaps, a certain scaling relation between those exists which leads to the cancellation, but at the moment we are not aware of it.

Now we would like to note that a convenient way to compute $`\nu `$ would be just summing all the trajectories of Fig. 4 together. As above, we expect this to be

$$\underset{V_0}{}\overline{x^2(t)}=𝑑V_0P(V_0)\frac{1}{\frac{1}{Dt}+V_0^{2\nu }}t^{1\frac{1}{2\nu }}$$
(37)
Numerically this exponent is around $`0.62`$ which corresponds to $`\nu =1.30`$, very close to the exact value. This method proves to be very convenient when studying quantum Hall transition in the quantum regime as well .

Knowing $`x(t)`$ it is not difficult to compute $`G(x,t)`$ according to (28).

$$G(k,t)=\underset{V_0}{}\overline{\mathrm{exp}\left(ikx(t)\right)},$$
(38)
so that the average is taken over all the trajectories. Compare this formula with (8). We are led to conjecture that

$$S(x,t;V_0)=2\pi P(V_0)\overline{\mathrm{exp}\left(ikx(t)\right)},$$
(39)
where the average, as always, is taken over the trajectories with the same $`V_0`$. Indeed, not only (8) is automatically satisfied, but also this form of $`S`$ reflects the fact that the particle propagates with a given energy.

With this definition of $`S`$ and $`G`$ it is easy to check that the formulae (10)-(13) are indeed true and precisely defined. Indeed, their applicability was at small $`k`$ where

$$\overline{\mathrm{exp}\left(ikx(t)\right)}1\frac{k^2}{4}\overline{x^2(t)}\mathrm{exp}\left(\frac{k^2}{4}\overline{x^2(t)}\right).$$
(40)
From this point on, (10)-(13) can be carried on without any further obstacles. The only difference between classical and quantum case will be the value of $`\nu `$, which in the classical case is equal to $`4/3`$.

It is interesting to note that at large $`k`$ the functions $`S`$ and $`G`$ will be very different in quantum and classical case. In the quantum case, $`S`$ is a function of the ratio $`q^2/\omega `$. In the classical case this question was investigated in . There it was shown that $`S(x,t;0)`$, computed at zero energy, was a function of $`q^{\frac{7}{4}}/\omega `$, with $`7/4`$ being the fractal dimension of the infinite percolating cluster. We refer the reader to  for further information on this subject. In the meanwhile, the form of $`G`$ or even $`S`$ at generic values of energy in the classical case remains unknown.

Finally, we would like to comment on the interesting feature of the numerical data of Fig. 4. The diffusion constant, when read off this data (given in units where the bonds of the triangular lattice of Fig. 3 have unit length), is equal to

$$D=1.008\pm 0.014$$
(41)
This number is remarkably close to $`1`$

Assuming that $`D=1`$ is indeed an exact relation, we can demonstrate that the critical conductivity is equal to

$$\sigma _{xx}=\frac{\sqrt{3}}{4}\frac{e^2}{h},$$
(42)
an exact result recently derived in  by a very different approach. To do that, we need to use the relationship between the diffusion constant and the critical conductivity proposed by F. Evers and W. Brenig in .

According to the Einstein’s relation, the conductivity is given by

$$\sigma _{xx}=e^2\rho \frac{D}{4},$$
(43)
where $`e`$ is the electron charge and $`\rho `$ is the density of states per unit area per unit energy at the transition. The problem is to find $`\rho `$ in the lattice model introduced here.

We can use a standard quasiclassical argument to find the density of states for a particle moving in a magnetic field and a potential $`V(x)`$. The total number of states whose energy is below $`E`$ is given, in the quasiclassical limit, by the area of the region where $`V(x)<E`$, divided by $`2\pi l^2`$,

$$A(E)=\frac{1}{2\pi l^2}d^2x\theta \left(EV(x)\right).$$
(44)
The density of states per unit area is

$$\rho (E)=\frac{1}{2\pi l^2A}d^2x\delta \left(EV(x)\right),$$
(45)
where $`A`$ is the total area.

Comparing $`\rho (E)`$ with (32) we see that $`\rho (E)`$ coincides with the time $`T`$ it takes for a particle to go over all the trajectories with energy $`E`$. Restoring $`\mathrm{}`$ and $`l`$ dependence in (32) we find

$$\rho (E)=\frac{T}{hA}$$
(46)
Therefore, the conductivity reduces to

$$\sigma =\frac{e^2}{h}\frac{DT}{4A}$$
(47)
It is now clear that the overall answer is independent of the units in which we measure $`D`$, as long as they are the same as the ones for $`A`$ and $`T`$.

In the lattice model of Fig. 3 the total time is equal to the total number of links available for a particle with zero energy, the energy of the transition. Assume that the total number of triangles is equal to $`N`$. The area $`A`$ of the lattice is equal to $`N\times \frac{\sqrt{3}}{4}`$, where $`\frac{\sqrt{3}}{4}`$ is the area of a single triangle. The total number of links is $`L=\frac{3N}{4}`$. We took into account that each triangle has 2 links with a probability $`3/4`$. This is because with a probability of $`1/4`$ all three energies at the vertices of a triangle are either less or greater than 0, in which case there are no trajectories in the triangle at all. But each link belongs to two triangles, hence the overcounting and the answer of $`\frac{3N}{4}`$. Substituting this into (47) and taking into account that $`D1`$, we arrive at

$$\sigma _{xx}=D\frac{\sqrt{3}}{4}\frac{e^2}{h}\frac{\sqrt{3}}{4}\frac{e^2}{h},$$
(48)
in very good agreement with .

This result is also in agreement with the one of  where the critical conductivity was estimated at $`\sigma _{xx}0.45\frac{e^2}{h}`$.

Notice that the conductance of a finite sample would not be Ohmic since the particle does not undergo a simple random walk, as reflected in a complicated structure of $`G`$ or $`S`$ at finite momenta $`k`$. Only for very thin samples, when it takes little time for a particle to reach one end of a sample from another, can the Ohm’s law be applicable, in agreement with .

## V Conclusions

To summarize, we showed that in the classical limit the density of a particle undergoing the quantum Hall transition moves along the percolation hull trajectories. At small times this particle exhibits normal diffusion, while at large times it localizes with the localization length exponent given by the percolation $`\nu =\frac{4}{3}`$. Our paper sets a new framework in which the classical percolation to quantum Hall effect renormalization group flow can be studied. Finally, we confirm numerically with a high degree of accuracy the result $`\sigma _{xx}=\sqrt{3}/4`$ for the critical conductivity of a classical quantum Hall transition.

## VI Acknowledgments

The authors are grateful to S. Boldyrev for many useful discussions and comments, to F. Evers who painstakingly explained his method  of relating the lattice diffusion constant to conductivity to us, to J. Cardy for interesting comments and to to E. H. Rezayi for making his computer facilities available to us.

This work was supported by the NSF grant PHY 94-07194


text:
# Predictions about the behaviour of diamond, silicon, SiC and some 𝐴^{𝐼𝐼𝐼}𝐵^𝑉 semiconductor materials in hadron fields

## 1 Introduction

The crystalline materials for semiconductor devices used in high fluences of particles are strongly affected by the effects of radiation. After the interaction between the incoming particle and the target, mainly two classes of degradation effects are observed: surface and bulk material damage, the last due to the displacement of atoms from their sites in the lattice. After lepton irradiation, the effects are dominantly at the surface, while heavy particles (hadrons and ions) produce both types of damages.

Up to now, in spite of the experimental and theoretical efforts, the understanding of the behaviour of semiconductor materials in radiation fields, the identification of the induced defects and their characterisation, as well as the explanation of the degradation mechanisms are still open problems.

The utilisation of semiconductor materials as detectors and devices operating in high radiation environments, at the future particle colliders, in space applications, in medicine and industry, makes necessary to obtain radiation harder materials.

Diamond, SiC and different $`A^{III}B^V`$ compounds (GaAs, GaP, InP, InAs, InSb) are in principle, possible competitors for silicon in the realisation of different electronic devices.

All analysed materials have a zinc-blend crystalline structure, with the exception of SiC, that presents the property of polytypism . The polytypism refers to one-dimensional polymorphism, i.e. the existence of different stackings of the basic structural elements along one direction. More than 200 polytypes have been reported in literature , but only few of them have practical importance. These include the cubic form $`3C(\beta )`$, and the $`4H`$ and $`6H`$ hexagonal forms. For the cubic polytype, the symmetry group is $`T_d^2`$, while for the hexagonal ones this is $`C_{6v}^4`$. Silicon is at the base of electronic industry, diamond and the $`A^{III}B^V`$ compounds present attractive electrical and/or luminescence properties, of interest for different applications, while the utilisation of SiC as a radiation detector, both in high energy physics and in the field of X-ray astronomy is now under extensive investigation more conferences in the field having sections dedicated to SiC.

The diamond has the reputation of being a radiation hard material and it is considered as a good competitor to silicon, but non all its properties as a radiation hard material have been proved experimentally.

The main goal of this paper is to give some theoretical predictions about the behaviour of different semiconductors in hadron fields (pions, protons), these materials representing potential candidates for detectors and electronic devices working in hostile environments.

The treatment of the interaction between the incident particle and the solid can be performed from the point of view of the projectile or of the target. In the first case, the relevant quantity is the energy loss (or equivalently the stopping power) and in the second situation the effects of the interactions are described by different physical quantities characterising material degradation. There is no a physical quantity dedicated to the global characterisation of the effects of radiation in the semiconductor material. A possible choice is the concentration of primary radiation induced defects on the unit particle fluence (CPD), introduced by . It permits the correlation of damages produced in different materials at the same kinetic energy of the incident hadron. For the comparison of the effects of different particles in the same semiconductor material, the non ionising energy loss (NIEL) is useful.

As a measure of the degradation to radiation, in the present paper the energy lost by the incident particle in the nuclear interaction and the concentration of primary defects induced in semiconductor bulk are calculated. If the energy loss of the incident projectile is, in principle, a measurable physical quantity, the concentration of primary defects is not directly observable and measurable and can be put in evidence only indirectly, e.g. from the variation of macroscopical parameters of the material or/and of electronic devices. It is to be noted that there exists also a kinetics of the defects induced by irradiation giving rise to the annealing process. A general treatment of the evolution processes is not possible, so, only some particular models exist in literature, see for example  and the references cited therein in the case of silicon. In these circumstances, in the present paper only the primary process of defect generation is modelled .This way, the energy range of incident hadrons for which the concentration of primary defects do not affect irreversibly the device properties could be established (theoretically predicted).

## 2 Model of the degradation

### 2.1 Energy loss

At the passage of the incident charged particle in the semiconductor material same of its energy is deposited into the target. The charged particles interact with both atomic and electronic systems in a solid. The total rate of energy loss, could, in general, be divided artificially into two components, the nuclear and the electronic part.

The energy lost due to interactions with the electrons of the target gives rise to material ionisation, while the energy lost in interactions with nuclei is at the origin of defect creation.

A comprehensive theoretical treatment of electronic stopping, which covers all energies of interest, cannot be formulated simply because of different approximations concerning both the scattering and contribution of different electrons in the solid. For fast particles with velocities higher then the orbital velocities of electrons, the Bethe-Bloch formula is to be used . At lower velocities, inner electrons have velocities greater than particle velocity, and therefore do not contribute to the energy loss. This regime has been modelled for the general case by Lindhard and Scharff and particular cases have been treated, e.g. in reference . If the particle has a positive charge, and a velocity close to the orbital velocity of its outer electrons, it has a high probability of capturing an electron from one of the atoms of the medium through which it passes. This process contributes to the total inelastic energy loss since the moving ion has to expend energy in the removal of the electrons which it captures.

The nuclear stopping depends on the detailed nature of the atomic scattering, and this in turn depends intimately on the form of the interaction potential. At low energies, a realistic potential based on the Thomas-Fermi approximation has been used in the literature  and at higher energies, where scattering results from the interaction of unscreened nuclei, a Rutherford collision model is to be used.

### 2.2 Bulk defect production

The mechanism considered in the study of the interaction between the incoming particle and the solid, by which bulk defects are produced, is the following: the particle, heavier than the electron, with electrical charge or not, interacts with the electrons and with the nuclei of the crystalline lattice. The nuclear interaction produces bulk defects. As a result of the interaction, depending on the energy and on the nature of the incident particle, one or more light particles are produced, and usually one or more heavy recoil nuclei. These nuclei have charge and mass numbers lower or at least equal to those of the medium. After the interaction process, the recoil nucleus or nuclei, if they have sufficient energy, are displaced from the lattice positions into interstitials. Then, the primary knock-on nucleus, if its energy is large enough, can produce the displacement of a new nucleus, and the process could continue as a cascade, until the energy of the nucleus becomes lower than the threshold for atomic displacements.

The concentration of the primary radiation induced defects on unit fluence has been calculated starting from the following equation:

$$CPD\left(E\right)=\frac{1}{2E_d}\underset{\text{i}}{}\frac{d\sigma _i}{d\mathrm{\Omega }}L\left(E_{Ri}\right)d\mathrm{\Omega }$$

where $`E`$ is the kinetic energy of the incident particle, $`E_d`$ the threshold energy for displacements in the lattice, $`E_{Ri}`$ the recoil energy of the residual nucleus, $`L\left(E_{Ri}\right)`$ the Lindhard factor describing the partition between ionisation and displacements and $`d\sigma _i/d\mathrm{\Omega }`$ the differential cross section for the process responsible in defect production. In the concrete calculations, all nuclear processes, and all mechanisms inside each process are included in the summation over index $`i`$. Because of the regular nature of the crystalline lattice, the displacement energy is anisotropic.

In the concrete evaluation of defect production, the nuclear interactions must be modelled, see for example references . The primary interaction between the hadron and the nucleus of the lattice presents characteristics reflecting the peculiarities of the hadron, especially at relatively low energies. If the inelastic process is initiated by nucleons, the identity of the incoming projectile is lost, and the creation of secondary particles is associated with energy exchanges which are of the order of MeV or larger. For pion nucleus processes, the absorption, the process by which the pion disappears as a real particle, is also possible.

The energy dependence of cross sections, for proton and pion interaction with the nucleus, presents very different behaviours: the proton-nucleus cross sections decrease with the increase of the projectile energy, then have a minimum at relatively low energies, followed by a smooth increase, while the pion nucleus cross sections present for all processes a large maximum, at about 160 MeV, reflecting the resonant structure of interaction (the $`\mathrm{\Delta }_{33}`$ resonance production), followed by other resonances, at higher energies, but with much less importance. Due to the multitude of open channels in these processes, some simplifying hypothesis have been done .

The process of partitioning the energy of the recoil nuclei (produced due the interaction of the incident particle with the nucleus, placed in its lattice site) by new interaction processes, between electrons (ionisation) and atomic motion (displacements) is considered in the frame of the Lindhard theory .

The factor characterising recoil energy partition between ionisation and displacements has been calculated analytically, solving the general equations of the Lindhard theory in some physical approximations. Details about the hypothesis used could be found in reference . All curves start, at low energies, from the same curve; they have at low energies identical values of the energy spent into displacements, independent on the charge and mass number of the recoil. At higher energies, the curves start to detach from this main branch. This happens at lower energies if their charge and mass numbers are smaller. The maximum energy transferred into displacements corresponds to recoils of maximum possible charge and mass numbers. The curves present then a smooth increase with the energy. For the energy range considered here, the asymptotic limit of the displacement energy is not reached.

For binary compounds, the Lindhard curves have been calculated separately for each component of the material, and the average weight Bragg additivity has been used. In this case, a threshold for atomic displacements must be considered for each atomic species and for each direction in the crystal. In the concrete calculations, a weighted value, independent on the crystalline direction has been used.

## 3 Results, discussions and some possible applications

The nuclear stopping power presents an energy dependence with a pronounced maximum. It is greater for heavier incident particles: protons compared to pions. In a given medium, the position of the maximum is the same for all particles with the same charge. In figure 1, the nuclear energy loss in diamond, silicon, silicon carbide, GaP, GaAs, InP, InAs and InSb is represented for protons and pions respectively, as a function of their kinetic energy. In the same medium, the position of the maximum is the same for pions and protons.

The behaviour of these materials in proton and pion fields is characterised by the CPD. In Figure 2, the dependence of the CPD as a function of the protons kinetic energy and medium mass number is presented for diamond, silicon, SiC GaAs and InP - see reference  and references cited therein. The values for diamond degradation are from reference , the corresponding ones for silicon are averaged values from references  and , SiC - from reference  and those for GaAs and InP are from reference . Low kinetic energy protons produce higher degradation in all materials. The discontinuity in the surface is related to differences in the behaviour of the CPD for monoatomic materials (or binary ones with close elements), and binary ones with remote elements in the periodic table.

For pion induced degradation, the energy dependence of CPD (as well as of the NIEL) presents two maxima, the relative importance of which depends on the target mass number: one in the region of the $`\mathrm{\Delta }_{33}`$ resonance, more pronounced for light elements and compounds containing light elements, and another one around 1 GeV kinetic energy, more pronounced for heavy elements. At higher energies, an weak energy dependence is observed, and a general $`A_{average}^{3/2}`$ dependence of the NIEL can be approximated . In Figure 3, the CPD for all analysed materials (diamond, Si, SiC, GaP, GaAs, InP, InAs, InSb) is represented as a function of the pion kinetic energy and of material average mass number. The differences in the behaviour of these materials are clearly suggested by the discontinuity in the mesh surfaces.

In the energy range considered in the paper, it could be observed that the CPD produced by pions and protons, and characterising the bulk degradation, are very different and reflect the peculiarities of the interactions of the two particles with the semiconductors. For pions, there are two maxima, one in the region of the $`\mathrm{\Delta }_{33}`$ resonance, corresponding to about 140 - 160 MeV kinetic energy, and the other at higher energies, around 1 GeV. The relative importance of these maxima depends on the mass number of the material. In comparison with this behaviour, the CPD produced by protons decreases abruptly with the increase of energy at low energies, followed by a smooth and slow increase at higher energies.

In relation to their behaviour in pion fields, these materials could be separated into two classes, the first with monoatomic materials or materials with relatively close mass numbers (diamond, silicon, GaAs and InSb), and the second comprising binary materials with remote mass numbers of the elements (SiC, GaP, InP, InAs) with similar behaviours inside each group. The diamond is the hardest material from all considered here. A slow variation of the primary defect concentration has been found for pion irradiation of diamond, silicon, SiC, GaP and GaAs, in the whole energy range of interest, with less than 2 displacements/cm/unit of fluence. In contrast to these materials, there are others, characterised by a low CPD in the energy range up to 200 MeV, (which represents this way the upper limit of the energy range where their utilisation in pion field is recommended), followed by a pronounced increase of displacement concentration with energy to more than 8 displacements/cm/unit fluence for InSb.

It is to be mentioned that, in the present model hypothesis, for SiC, negligible differences have been found between different polytypes in what regards the effects of pion and proton degradation , conclusions in accord with the experimental results . The behaviour of SiC in radiation fields is between the corresponding one of diamond and silicon.

As it is well known, the analysed semiconductors are possible materials for detectors and electronic devices which have to work long time in particle physics experiments, space applications, etc., in intense fields of hadrons, and in experimental configurations which impose high reliability of devices, and must present a controlled degradation of their parameters. As possible applications we will analyse two hypothetical cases: the utilisation of diamond, silicon, SiC or GaAs as detectors at the Large Hadron Collider (LHC) at CERN, and the long time exposure of the electronic devices in the field produced by cosmic rays.

For the LHC, the standard physics programme is based on the study of proton - proton interactions, at about 7 TeV beam energy, on an integrated luminosity of $`5x10^5`$ pb-1 which corresponds to 9 year of operation, for an annual operation time of 1.9x107 s. The irradiation background is continuous. The charged hadrons are produced in the primary interactions, while the neutrons are albedo particles. The charged pions are the dominant particles, followed by protons, antiprotons and kaons. As an illustration of the above calculations for the degradation of different semiconductors in proton and pion fields, the results of the simulation of Gorfine and Taylor  have been chosen, for pion and proton fluxes in the inner detector assembly region, parallel to the beam axis. Both protons and pions are transported down to thermal energies in the detector, by nuclear interactions. The particle flux energy spectra have been simulated for a first layer of Si detectors (situated at 11.5 cm), and with complete moderator. The obtained spectra have been found to be slowly dependent on the material of the inner detector and so, in the present paper the same hadron spectra for diamond, silicon carbide, silicon and gallium arsenide have been utilised.

The convolution of the pion and proton spectra with the energy dependence of the CPD has been done in the energy range 50 MeV - 10 GeV, and 10 MeV - 10 GeV for pions and protons respectively. Below 50 MeV, a realistic estimation of materials degradation to pions is very difficult due to the lack of experimental data on pion - nucleus interaction and also to the increase of the weight of Coulomb interaction.

The results of these calculations are summarised in Fig. 4, for the diamond, silicon carbide, silicon and GaAs options, both for pion and proton degradation. In the analysed case, diamond and SiC are the hardest materials in both pion and proton fields; the diamond is harder to pions than to protons. The behaviour of silicon is similar in both particle fields. The GaAs option is not recommended because of an order of magnitude higher degradation in comparation with all other considered materials.

Another possible utilisation of semiconductor devices in radiation fields is related to space applications. In the primary cosmic radiation, the most abundant particles are the protons . Other charged particles (for example $`\pi ^{+/}`$, $`e^{+/}`$, $`\mu ^{+/}`$, $`\nu _{mu}`$, etc.) are produced in the interaction of the primary cosmic rays in air. The damage induced in diamond, Si, SiC, GaAs and InP by the primary cosmic field has been estimated for protons in the energy range 10 MeV - 10 GeV, and the results are presented in Figure 5. The devices have been supposed to be exposed directly to the cosmic field. In this case too, diamond has been found to be the hardest material. GaAs and InP suffer a degradation of a factor of about 50 times higher in comparation with diamond and this behaviour can affect irreversibly the properties of these materials for long time operation.

The degradation produced by the particle field at LHC (pions and protons) and by the free protons from the primary cosmic rays respectively are of the same order of magnitude for each of the materials investigated.

## 4 Summary

A systematic theoretical study has been performed, investigating the interaction of charged hadrons with semiconductor materials and the mechanisms of defect creation by irradiation.

The nuclear stopping power has been found to be greater for heavier incident particles (protons compared to pions), and for lighter media. The position of its maximum is the same for protons and pions in the same medium.

The mechanisms of the primary interaction of the hadron with the nucleus (nuclei) of the semiconductor lattice have been explicitly modelled and the Lindhard theory of the partition between ionisation and displacements has been applied.

For protons, the low kinetic energy particles produce higher degradation in all materials.

For pions, the energy dependence of CPD presents two maxima, the relative importance of which depends on the target mass number: one in the region of the $`\mathrm{\Delta }_{33}`$ resonance, more pronounced for light elements and compounds containing light elements, and another one around 1 GeV kinetic energy, more pronounced for heavy elements. At higher energies, an weak energy dependence is observed. A slow variation of the primary defect concentration has been found for pion irradiation of diamond, silicon, GaP and GaAs, in the whole energy range of interest, with less than 2 displacements/cm/unit of fluence. In contrast to this situation, for the other semiconductor materials analysed, a low CPD is estimated in the energy range up to 200 MeV (which represent the energy range up to their utilisation in pion field is recommended), followed by a pronounced increase of displacement concentration to more than to more than 8 displacements/cm/unit of fluence at high energies.

The behaviour of this semiconductor materials has been analysed comparatively both in relation to particle physics experiments (inner part of the detection system at LHC) and to space applications (the devices being considered to be exposed directly to the cosmic ray field.

Figure captions

Figure 1: The nuclear energy loss in diamond, silicon, silicon carbide, GaP, GaAs, InP, InAs and InSb as a function of the kinetic energy of the incident particle: protons (up) and pions (down) respectively.

Figure 2a: The concentration of primary defects on unit fluence (CPD) in diamond, silicon, SiC, GaAs and InP induced by protons, as a function of the kinetic energy and average mass number of the semiconductor material. The mesh surfaces are drawn only to guide the eye.

Figure 2b: The dependence of CPD as a function of proton kinetic energy, for the same semiconductors.

Figure 3a: The energy and material dependence of the CPD on unit pion fluence for diamond, Si, SiC, GaP, GaAs, InP, InAs and InSb. The mesh surfaces are drawn to guide the eyes.

Figure 3b: The CPD as a function of the kinetic energy of incident pions, for the same semiconductors.

Figure 4: Estimated CPD on unit fluence induced in diamond, silicon, SiC and GaAs, by the simulated flux energy spectra of pions and protons in the inner detector at LHC.

Figure 5: Estimated CPD on unit fluence induced by the primary cosmic ray flux energy spectra in diamond, Si, SiC, GaAs and InP (only the effects produced by protons are considered), in the hypothesis that these semiconductor materials are exposed directly in the radiation field.


text:
# Convergence of Monte Carlo Simulations to Equilibrium

\[

## Abstract

We give two direct, elementary proofs that a Monte Carlo simulation converges to equilibrium provided that appropriate conditions are satisfied. The first proof requires detailed balance while the second is quite general.

\]

Monte Carlo simulations are widely used in statistical physics. If the algorithm satisfies the detailed balance condition, it is easy to show that the desired distribution is a stationary distribution, i.e. if the system is, by some means, put into the desired distribution, it will subsequently stay in this distribution. It is harder, but nonetheless crucial, to also show that, starting from a general distribution, the algorithm will converge to the desired distribution. Although this can be proved without too much difficulty for a system with a finite number of states , this proof is unfamiliar to most physicists. The argument in the original paper of Metropolis et al. makes convergence plausible but is not a proof. Most physics texts on Monte Carlo methods either do not give a proof of convergence or refer the reader to rather abstract derivations in the mathematics literature, or rely on the Frobenius-Perron theorem which is unfamiliar to most physicists.

In this paper we present two proofs of convergence which are self contained and use only elementary methods. We feel that it is useful to present these derivations here because (i) it is not widely known in the physics community that it is not difficult to prove convergence, at least for systems with a finite number of states, and (ii) our proofs are (to the best of our knowledge) different from and as simple as existing proofs in the mathematical literature. Even the proof of Ref., which is of comparable simplicity, is hard to understand physically; this is discussed more fully near the end of the paper.

The first proof relies on the Monte Carlo algorithm satisfying the condition of detailed balance. Although this is true in essentially all Monte Carlo simulations, it is not strictly required for the algorithm to converge to the equilibrium distribution. In the second half of this paper, we shall also present a more general proof which relaxes this condition, and which is very different from the proof assuming detailed balance.

Throughout this paper, we shall assume that the system being considered has a finite number of states. Most systems in physics can be approximated as such by discretizing any continuous variables sufficiently finely; for instance, in molecular simulations, it should be permissible to limit the phase space for any particle to a sufficiently large region, and then discretize it in small intervals. A rigorous proof of convergence for systems with an infinite number of states is much more complicated .

In the simplest case, the desired distribution that the Monte Carlo method seeks to simulate is the Boltzmann distribution at some temperature. However, in dealing with glassy systems, it is sometimes more efficient to simulate a different distribution, as in the multicanonical ensemble, the $`1/k`$ ensemble, and parallel tempering. The results presented in this paper are valid whenever the desired distribution is a stationary distribution of the algorithm, and detailed balance — or, more generally, ergodicity — is satisfied. This includes the cases mentioned above.

The essential ingredients of the Monte Carlo method in statistical physics are the (non-negative) “transition rates”, $`w_{lm}`$ , defined to be the probability that, given the system is in state $`l`$ at “time” $`t`$, then it will be in state $`m(l)`$ at time $`t+1`$. We define time to be incremented by one every Monte Carlo move (not sweep) and assume initially that all moves are equivalent, so the $`w_{lm}`$ do not depend on time. An example would be flipping a single spin chosen at random. The important case of sequential updating will be discussed later.

The probability that the system is in state $`l`$ at time $`t`$ is defined to be $`P_l(t)`$. The evolution of these probabilities is governed by the “master equation”,

$$P_l(t+1)P_l(t)=\underset{ml}{}\left[P_m(t)w_{ml}P_l(t)w_{lm}\right].$$
(1)
The first term on the right hand side describes transitions into state $`l`$ from $`m`$ (which therefore increases $`P_l`$ and so has a plus sign) while the second term describes transitions out of state $`l`$, which decreases $`P(l)`$. Note that only terms with $`ml`$ contribute. We can also define $`w_{ll}`$ to be the probability that the system stays in state $`l`$, i.e. $`w_{ll}=1_{ml}w_{lm}`$, or equivalently,

$$\underset{m}{}w_{lm}=1.$$
(2)
Eq. (2) implies that the master equation can be written

$$P_l(t+1)=\underset{m}{}P_m(t)w_{ml},$$
(3)
where the term $`m=l`$ is now included.

A necessary condition for the method to work is that the desired distribution, $`P^{\mathrm{eq}}`$, is stationary, i.e. if $`P_l(t)=P_l^{\mathrm{eq}}`$ for all $`l`$ then $`P_l(t+1)=P_l^{\mathrm{eq}}`$. This means that the right hand side of Eq. (1) must vanish for $`P=P^{\mathrm{eq}}`$. The condition of detailed balance consists of the assumption that each term on the right hand side of Eq. (1) separately vanishes for $`P=P^{\mathrm{eq}}`$, i.e.

$$P_l^{\mathrm{eq}}w_{lm}=P_m^{\mathrm{eq}}w_{ml}.$$
(4)
In the first part of this paper, we shall assume that Eq. (4) is satisfied.

We start with the following quantity, which is a measure of the deviation from equilibrium,

$$G=\underset{l}{}\frac{1}{P_l^{\mathrm{eq}}}\left(P_lP_l^{\mathrm{eq}}\right)^2=\underset{l}{}\left(\frac{P_l^2}{P_l^{\mathrm{eq}}}\right)1$$
(5)
evaluated at time $`t`$, where the last expression follows because $`P`$ and $`P^{\mathrm{eq}}`$ are normalized.

At time $`t+1`$ we indicate (for compactness of notation) the probabilities by $`P_l^{}`$ and the corresponding value of $`G`$ by $`G^{}`$. We will show that $`G`$ monotonically decreases, i.e.

$$\mathrm{\Delta }GG^{}G0,$$
(6)
where the equality only holds if $`G`$ and $`G^{}`$ both vanish, so the system is in equilibrium. This shows that the system will eventually approach arbitrarily close to the equilibrium distribution.

Using Eqs. (3) and (5), $`\mathrm{\Delta }G`$ can be written as

$$\mathrm{\Delta }G=\underset{l,m,n}{}\left[w_{ml}w_{nl}\frac{P_mP_n}{P_l^{\mathrm{eq}}}\right]\underset{l}{}\frac{P_l^2}{P_l^{\mathrm{eq}}}.$$
(7)
In the first term on the right hand side of Eq. (7) we use the detailed balance condition, Eq. (4), to replace $`w_{ml}`$ by $`w_{lm}P_l^{\mathrm{eq}}/P_m^{\mathrm{eq}}`$, and in the second term we can use Eq. (2) to insert a factor of $`_mw_{lm}`$ (and interchange the indices $`l`$ and $`m`$). This gives

$$\mathrm{\Delta }G=\underset{l,m,n}{}\left[w_{lm}w_{ln}P_l^{\mathrm{eq}}\frac{P_mP_n}{P_m^{\mathrm{eq}}P_n^{\mathrm{eq}}}\right]\underset{l,m}{}w_{ml}\frac{P_m^2}{P_m^{\mathrm{eq}}}$$
(8)
Applying the detailed balance relation again and incorporating a factor of $`_nw_{ln}`$, the last term in the above equation can be written as

$$\underset{l,m,n}{}w_{lm}w_{ln}P_l^{\mathrm{eq}}\left(\frac{P_m}{P_m^{\mathrm{eq}}}\right)^2.$$
(9)
Taking the half the sum of this and the same expression with $`m`$ replaced by $`n`$, we finally get

$$\mathrm{\Delta }G=\frac{1}{2}\underset{l,m,n}{}w_{lm}w_{ln}P_l^{\mathrm{eq}}\left(\frac{P_m}{P_m^{\mathrm{eq}}}\frac{P_n}{P_n^{\mathrm{eq}}}\right)^2,$$
(10)
where terms with $`m=l`$ and $`n=l`$ are included.

Eq. (10) is the main result of this part of the paper. It shows that $`\mathrm{\Delta }G`$ is definitely negative unless, for every state $`l`$, all states which can be reached from $`l`$ in a single move — equivalently, with detailed balance, all states from which $`l`$ can be reached in a single move — have probabilities proportional to the equilibrium probabilities. The most natural scenario is that all states satisfy this with the same proportionality constant (which must be unity) i.e. the system is in equilibrium. However, $`\mathrm{\Delta }G`$ also vanishes if $`P_m/P_m^{\mathrm{eq}}`$ assumes different values for states which have no common one-step descendants. Hence, to achieve full equilibrium, the algorithm must also be ergodic, i.e. starting from a given state, after a sufficiently long time there is non-zero probability for the system to be in any state. The condition of ergodicity is sufficient to ensure that even if $`\mathrm{\Delta }G`$ is accidentally zero at some time step, it must decrease later, since any two states must have common descendants after several time steps. If in addition $`w_{ll}0`$ for all $`l,`$ which is usually true, the one-step descendants of a set of states must include the set itself, so that it is not possible to break up all the states of the system into subsets with no common one-step descendants across two subsets. Thus if this condition is satisfied, $`\mathrm{\Delta }G`$ cannot be zero (without the system being in equilibrium) at any time step.

We will distinguish between a process which is ergodic and one which satisfies the lesser condition of being ”irreducible”. In the latter, the system will eventually sample all states starting from a given initial state, but, at a fixed later time, the probability for some of the states is zero. A familiar example which is irreducible but not ergodic is the Ising model at infinite temperature simulated using Metropolis updating, for which the probability to flip is unity in this limit. Clearly after an odd time, the number of flipped spins must be odd and vice-versa. For such a non-ergodic system, it is possible for $`P_m/P_m^{\mathrm{eq}}`$ and $`P_n/P_n^{\mathrm{eq}}`$ to be different for states which have no common descendant at any fixed later time (for the Ising model example given here, states which differ by an odd number of spin flips).

Note that Eq. (10) does not give an estimate for how fast equilibrium is reached.

For random updating considered so far the probability of making a transition is the same for every move, i.e. writing Eq. (1) as

$$P_l(t+1)=\underset{m}{}\mathrm{\Gamma }_{lm}P_m(t),$$
(11)
then $`\mathrm{\Gamma }`$, the transition matrix (related to $`w`$ by $`\mathrm{\Gamma }_{lm}=w_{ml}`$), is the same for each “time” $`t`$. However, for sequential updating, the transition matrix depends on which site is being updated, so, for a complete sweep, we have

$$\mathrm{\Gamma }=\mathrm{\Gamma }^{(1)}\mathrm{\Gamma }^{(2)}\mathrm{}\mathrm{\Gamma }^{(N)},$$
(12)
where $`\mathrm{\Gamma }^{(i)}`$ is the transition matrix for updating spin $`i`$. Although the $`\mathrm{\Gamma }^{(i)}`$ individually satisfy the detailed balance condition, the transition matrix for the whole sweep, $`\mathrm{\Gamma }`$, does not, because the probability of the reverse transition, $`ml`$ say, for a whole sweep, is related to the probability of transition $`lm`$ in the desired way only if the spins are updated in the reverse order. Despite overall lack of detailed balance, convergence to the equilibrium distribution is still obtained for sequential updating because $`G`$ decreases at each step, as long as each of the transition probabilities, $`\mathrm{\Gamma }^{(i)}`$, satisfies the detailed balance condition.

More generally, detailed balance is not necessary, see Ref. and references therein. In the rest of this paper, we give a derivation of the necessary and sufficient conditions for convergence to $`P^{\mathrm{eq}}`$:

1. the algorithm has $`P^{\mathrm{eq}}`$ as a stationary distribution
2. for any pair of states $`(i,j)`$, there exists some $`T_{ij}`$ and some state $`k_{ij}`$ such that at the $`T_{ij}`$’th time step the probabilities to have reached $`k`$ from $`i`$ and $`j`$ are both non-zero.

We shall see that (ii) (with (i)) implies ergodicity.

We first prove that conditions (i) and (ii) are sufficient. Note that Eq. (1) is linear in the probabilities $`P_l,`$ so that if we define $`\delta P_l=P_lP_l^{\mathrm{eq}},`$ then (with condition (i)) $`\delta P_l`$ also satisfies Eq. (1), with $`_l\delta P_l=0.`$

It is convenient to use a measure of the deviation from equilibrium that is different from Eq. (5):

$$L=\underset{l}{}|\delta P_l|.$$
(13)
Like $`G,`$ the quantity $`L`$ is positive unless the distribution has converged to $`P^{\mathrm{eq}}.`$ We denote by $`W_{lm}(t)`$ the probability that the system, starting out in state l, reaches state m after t time steps. This is obtained by ‘iterating’ the transition rates $`\{w\}.`$ Then

$$L(t)=\underset{m}{}|\underset{l}{}W_{lm}(t)\delta P_l(0)|\underset{m,l}{}W_{lm}(t)|\delta P_l(0)|.$$
(14)
Using the result $`_mW_{lm}(t)=1,`$ we see that $`L(t)L(0).`$ We now compare the two sides of the inequality in Eq.(14), by choosing a specific value of $`m`$ and carrying out the $`l`$-summation. If all states $`l`$ for which $`W_{lm}0`$ have the same sign for $`\delta P_l(0)`$, it is clear that the $`l`$-summation on both sides are equal. On the other hand, if some of these states have $`\delta P(0)>0`$ and others have $`\delta P(0)<0,`$ the $`l`$-summation on the left hand side has both positive and negative terms, and must be less than the corresponding sum on the right hand side. Thus so long as there is at least one state $`m`$ which receives ‘contributions’ from two states $`(i,j)`$ with opposite $`\delta P(0)`$, i.e. $`\delta P_i(0)`$ and $`\delta P_j(0)`$ have opposite signs and $`W_{im}0`$ and $`W_{jm}0,`$ we see that $`L(t)<L(0).`$ Condition (ii) ensures that for any $`(i,j)`$ this is the case for $`t=T_{ij}.`$ Thus $`L(t)`$ stays constant until $`t=\mathrm{min}[T_{ij}],`$ where the minimum is taken over all $`(ij)`$ for which $`\delta P_i`$ and $`\delta P_j`$ have the opposite sign, and then decreases at that time step. The time $`t`$ can be reinitialized to zero at this point, and the whole argument repeated again. Note that nothing in the argument requires the transition rates $`w_{lm}`$ to be time independent, so long as conditions (i) and (ii) are always satisfied. Also, as with the first approach above, no estimate has been obtained for how fast $`L(t)`$ approaches zero.

It may seem surprising that $`L(t)`$ need not decrease at every time step, whereas Eq. (10) shows that $`G`$ must. If we start out with $`\delta P`$ positive and negative on two states that are well separated from each other (in the sense that many time steps are required before the two states have common descendants), it is clear that $`L(t)`$ does not change at first. However, $`G`$ decreases because the positive and negative $`\delta P`$’s both get ‘smeared out’ over several states. Thus different measures can depict the approach to equilibrium at different rates. Since any physical observable $`O`$ is given by $`O=_lO_lP_l,`$ the error $`\delta O_l|O_l\delta P_l|L\mathrm{max}_l[|O_l|],`$ so that $`L`$ is a conservative indicator of how observables approach equilibrium.

We have proved the sufficiency of conditions (i) and (ii); the necessity can be easily demonstrated. The necessity of (i) is obvious. If (ii) is violated, starting with an initial condition $`\delta P_i(0)=\delta P_j(0)`$ and all other $`\delta P_l`$’s equal to zero, the positive and negative regions for $`\delta P`$ stay separate for all time and cannot neutralize each other.

Condition (ii) (with (i)) is equivalent to the (seemingly stronger) condition of ergodicity. If the system were not ergodic, then starting with $`P_l(0)=\delta _{il}`$ would lead to a $`P(t)`$ with $`P_j(t)=0`$ for some $`j,`$ since all states would not be accessible from the state $`i`$ at time $`t.`$ Therefore $`P(t)`$ could not have converged to $`P^{\mathrm{eq}}.`$ Since we have seen that (i) and (ii) imply convergence, they must imply ergodicity.

The proof of convergence in Ref. is comparable in simplicity to the proof given above. However, it applies the transpose $`\mathrm{\Gamma }^T`$ of the transition matrix to the initial state of the system (represented as a column vector), and shows that under iteration the initial state evolves to a column vector whose entries are all identical. Although it is possible  to deduce the desired properties of $`\mathrm{\Gamma }`$ from this, there is no physical interpretation of evolution under $`\mathrm{\Gamma }^T,`$ which does not even conserve probability. In contrast, the proof given above shows how fluctuations from the desired distribution are ”mixed” by time evolution, with positive fluctuations annihilating negative ones. It is also physically clear why ergodicity is essential for this annihilation process to proceed to completion: positive and negative fluctuations are then assured of encountering each other under time evolution.

Although we have considered the conditions for convergence to $`P^{\mathrm{eq}},`$ in practice Monte Carlo simulations are not carried out using such an ensemble average, but by simulating a single system and taking a time average of many measurements. If we restrict ourselves to the case when the transition matrix $`\mathrm{\Gamma }`$ is time-independent, this defines a homogeneous Markov chain . It can then be shown that if the Markov chain is irreducible and has $`P^{\mathrm{eq}}`$ as a stationary distribution, the time average $`P(t)_t`$ converges to $`P^{\mathrm{eq}}.`$ We indicate how to prove this result here: since $`L(t)L(0)`$ in general, none of the eigenvalues of $`\mathrm{\Gamma }`$ can have modulus greater than unity . Irreducibility ensures that the eigenvalue 1 has a unique eigenvector, $`P^{\mathrm{eq}}`$. The time averaging removes contributions from oscillatory eigenvalues of $`\mathrm{\Gamma },`$ of the form $`e^{i\theta }`$ ($`\theta 0`$). By comparison, the stronger condition of ergodicity (ii) (with (i)) is enough to rule out oscillatory eigenvalues, since no time averaging is needed for convergence to $`P^{\mathrm{eq}}.`$ It is not clear how this proof that irreducibility is sufficient would generalize to time dependent transition rates, since the evolution of $`P`$ is then not directly related to the eigenvalues of $`\mathrm{\Gamma }(t)`$.

We thank M. E. J. Newman for useful correspondence.