text:
# Exact Asymptotics for One-Dimensional Diffusion with Mobile Traps

## Abstract

We consider a diffusing particle, with diffusion constant $`D^{}`$, moving in one dimension in an infinite sea of noninteracting mobile traps with diffusion constant $`D`$ and density $`\rho `$. We show that the asymptotic behavior of the survival probability, $`P(t)`$, satisfies $`lim_t\mathrm{}[\mathrm{ln}P(t)]/\sqrt{\rho ^2Dt}=4/\sqrt{\pi }`$, independent of $`D^{}`$. The result comes from obtaining upper and lower bounds on $`P(t)`$, and showing that they coincide asymptotically. We also obtain exact results for $`P(t)`$ to first order in $`D^{}/D`$ for an arbitrary finite number of traps.

PACS numbers: 05.40.-a, 02.50.Ey, 82.20.-w

The asymptotics of the survival probability, $`P(t)`$, for a particle diffusing among mobile traps is a longstanding challenge. Over ten years ago Bramson and Lebowitz  showed by rigorous arguments that in one dimension $`P(t)`$ has a stretched exponential decay for large $`t`$,

$$P(t)\mathrm{exp}(\lambda t^{1/2}).$$
(1)
The computation of the constant $`\lambda `$, however, has thus far proved intractable, both analytically and numerically. This longstanding problem is resolved in this Letter.

The problem addressed by Bramson and Lebowitz was first posed almost 20 years ago in the seminal paper of Toussaint and Wilczek (TW) . These authors introduced the two-species annihilation process, $`A+B0`$, as a model of monopole-antimonopole annihilation in the early universe, though applications to chemical kinetics and condensed matter physics are more numerous . TW showed that if the initial densities, $`\rho _A(0)`$ and $`\rho _B(0)`$, of the $`A`$ and $`B`$ particles are equal (and the particles are randomly distributed in space), the densities decay asymptotically as $`t^{d/4}`$ in space dimensions $`d<4`$. If the initial densities are different, however, the density of the minority species ($`A`$, say) decays much more rapidly.

The connection with the trapping problem is as follows. At late times $`\rho _A(t)\rho _B(t)`$ and the $`A`$ particles can be regarded as independently diffusing in a background of the majority $`B`$ particles, which act as traps for the $`A`$ particles by virtue of the annihilation reaction. An equivalent problem, therefore, is to consider a single $`A`$ particle moving among $`B`$ particles (which do not interact with each other) and ask for the probability, $`P(t)`$, that the $`A`$ particle survives to time $`t`$. In the context of the original $`A+B0`$ process (or $`A+BB`$ , which has the same asymptotics ) the “particle” $`A`$ and “traps” $`B`$ are taken to have the same diffusion constants, but for generality we will take them to be different, $`D^{}`$ and $`D`$ respectively.

To our knowledge, no analytical result has been obtained up to now for the constant $`\lambda `$ in Eq. (1). Furthermore, attempts to determine $`\lambda `$ by numerical simulations (or even to confirm the predicted stretched-exponential decay with exponent $`1/2`$) are severely hampered by large, slowly decaying transients . In ref. , a sophisticated numerical approach enabled data to be obtained down to $`P(t)10^{35}`$, but still the asymptotic time dependence could not be unambiguously established.

For this trapping problem, the traps are infinite in number, and distributed randomly on the interval ($`\mathrm{}`$, $`\mathrm{}`$) with density $`\rho `$. By contrast, if the number of traps is finite, the problem is equivalent to the much studied predator-prey problem , where the traps are predators and the particle is the prey. In this case the prey survival probability decays as a power law,

$$P(t)t^{\theta (N_L,N_R,D^{}/D)},$$
(2)
where $`N_L`$ ($`N_R`$) is the number of predators initially to the left (right) of the prey. To our knowledge the only exactly-solved examples are for $`N_L+N_R2`$. Obtaining analytical results for more than two predators is another longstanding challenge.

In this Letter we obtain two exact analytical results. First, we finally resolve the question of the asymptotics of $`P(t)`$ for the trapping problem. We verify the asymptotic form (1), and determine exactly the value of the constant $`\lambda `$, namely $`\lambda =4\rho (D/\pi )^{1/2}`$, i.e.

$$P(t)\mathrm{exp}[4\rho (Dt/\pi )^{1/2}].$$
(3)
Note that this asymptotic result depends only on the density, $`\rho `$, and diffusion constant, $`D`$, of the traps and is independent of the diffusion constant $`D^{}`$ of the particle. The value of $`D^{}`$, however, does affect the rate of approach to asymptopia. Eq. (3) is derived by obtaining upper and lower bounds for $`P(t)`$, and showing that they converge for large $`t`$. Furthermore, the form of the lower bound suggests why large corrections to the asymptotic behavior might be expected.

Secondly, we outline  a perturbation theory in $`D^{}/D`$, based on a path-integral formulation, for the computation of the exponent $`\theta `$ in Eq. (2). For small $`D^{}/D`$ and arbitrary values $`N_L`$ and $`N_R`$, we find

$$\theta =\frac{N}{2}+\frac{1}{\pi }(N\mathrm{\Delta }^2)\frac{D^{}}{D}+O\left(\frac{D^2}{D^2}\right),$$
(4)
where $`N=N_L+N_R`$ is the total number of predators and $`\mathrm{\Delta }=N_LN_R`$ measures their left-right asymmetry with respect to the prey.

We will first describe the treatment of the trapping problem with infinitely many traps. The bounds that lead to Eq. (3) are obtained as follows.

1. Upper Bound. An obvious upper bound, $`P_U(t)`$, on $`P(t)`$, for any $`D^{}`$, is provided by the problem with $`D^{}=0`$, in which the particle stays at its initial position, which we call $`x=0`$. Although we have as yet been unable to make this bound rigorous , it is intuitively clear that when, as here, the traps are (statistically) symmetrically placed with respect to the particle, the particle will on average survive longer if it stays still than if it diffuses. This assertion has been checked, using the algorithm outlined in ref., for all (lattice) walks up to time $`t=28`$. It is also supported by Eq. (4), which shows that for any symmetric case ($`\mathrm{\Delta }=0`$), the decay of $`P(t)`$ is faster for small $`D^{}`$ than for $`D^{}=0`$, for any finite $`N`$.

For $`D^{}=0`$ (sometimes called the “scavenger model”  in predator-prey terminology), $`P(t)`$ is just the probability that none of the moving traps has reached the origin up to time $`t`$. This problem is exactly soluble . Since similar techniques will be needed to derive the lower bound, we outline the solution here.

The $`N`$ traps move independently according to the Langevin equations

$$\dot{x_i}=\eta _i(t),i=1,\mathrm{},N,$$
(5)
where $`\eta _i(t)`$ is Gaussian white noise with mean zero and correlator

$$\eta _i(t)\eta _j(t^{})=2D\delta _{ij}\delta (tt^{}).$$
(6)
The quantity $`P(t)`$ is just the product of the individual trap probabilities. For a given trap starting at $`x_i`$, the required probability is  $`P_1(x_i,t)=\mathrm{erf}(|x_i|/\sqrt{4Dt})`$. So our upper bound is

$$P_U(t)=\underset{i=1}{\overset{N}{}}\mathrm{erf}\left(|x_i|/\sqrt{4Dt}\right),$$
(7)
where $`\mathrm{}`$ means an average over the initial positions of the traps. Since the $`x_i`$ are also independent, the latter average also factors. Using $`N=\rho L`$, where $`L`$ is the length of the system, and each $`x_i`$ is uniformly distributed in $`(L/2,L/2)`$, gives

$`P_U(t)`$ $`=`$ $`\left[11/L{\displaystyle _{L/2}^{L/2}}𝑑x\mathrm{erfc}(|x|/\sqrt{4Dt})\right]^{\rho L}`$ (8)
$``$ $`\mathrm{exp}[4\rho (Dt/\pi )^{1/2}],`$ (9)
where the final result follows on taking the limit $`L\mathrm{}`$.

2. Lower Bound. Consider the same system as before but with a pair of absorbing boundaries at $`x=\pm l/2`$. We consider the subset of initial conditions in which all the traps lie outside the interval $`(l/2,l/2)`$ (and the particle is at $`x=0`$), and trajectories in which neither the particle nor any of the traps has crossed a boundary up to time $`t`$. We calculate the probability, $`P_L(t)`$, of such an occurrence over the ensemble of all initial conditions and trajectories. These restricted initial conditions and trajectories are a subset of all the possible initial conditions and trajectories in which the particle never meets a trap. It follows that $`P(t)P_L(t)`$, i.e. $`P_L(t)`$ is a lower bound.

The probability that there are no traps in the interval $`(l/2,l/2)`$ at $`t=0`$ is $`\mathrm{exp}(\rho l)`$. Given that there are no traps in this interval at $`t=0`$, the probability that no traps enter the interval up to time $`t`$ is given by the same result as in the derivation of the lower bound, namely $`\mathrm{exp}[4\rho (Dt/\pi )^{1/2}]`$. Finally, the probability that the particle, starting at $`x=0`$, has not left the interval $`(l/2,l/2)`$ up to time $`t`$ is given, for times $`tl^2/D^{}`$, by  $`(4/\pi )\mathrm{exp}(\pi ^2D^{}t/l^2)`$. Assembling these contributions gives

$$P(t)(4/\pi )\mathrm{exp}[4\rho (Dt/\pi )^{1/2}(\rho l+\pi ^2D^{}t/l^2)].$$
(10)
Since this inequality holds for all $`l`$, the best lower bound is obtained by maximizing with respect to $`l`$. The optimum value is $`l=(2\pi ^2D^{}t/\rho )^{1/3}`$, and the best lower bound is

$$P_L(t)=\frac{4}{\pi }\mathrm{exp}[4\rho (Dt/\pi )^{1/2}3(\pi ^2\rho ^2D^{}t/4)^{1/3}].$$
(11)
Since the second term in the exponent is negligible compared to the first as $`t\mathrm{}`$, the two bounds converge to yield the asymptotic form, Eq. (3), for $`P(t)`$. More precisely, we can take the logarithm of $`P(t)`$ and divide out the leading $`\sqrt{t}`$ dependence to get

$$\frac{4}{\sqrt{\pi }}\frac{\mathrm{ln}P(t)}{(\rho ^2Dt)^{1/2}}\frac{4}{\sqrt{\pi }}+3\left(\frac{\pi }{2}\right)^{2/3}\frac{(D^{}/D)^{1/3}}{(\rho ^2Dt)^{1/6}},$$
(12)
giving $`lim_t\mathrm{}[\mathrm{ln}P(t)]/(\rho ^2Dt)^{1/2}=4/\sqrt{\pi }`$.

As an aside we note that, while the left-hand inequality in (12) holds for all $`t`$, since Eq. (9) does, the right-hand inequality is strictly a large $`t`$ result. This is because the factor $`(4/\pi )\mathrm{exp}(\pi ^2D^{}t/l^2)`$ in Eq. (10) comes from the lowest mode in the Fourier decomposition of the survival probability of the particle in the interval $`(l/2,l/2)`$. This mode dominates for $`D^{}tl^2`$, which requires $`\rho ^2D^{}t1`$. A lower bound on $`P(t)`$ valid for all $`t`$ can easily be written down by including all Fourier modes, but the large-$`t`$ form (11) is sufficient for present purposes.

In Figure 1, the left and right hand sides of Eq. (12), representing the two asymptotic bounds, are plotted and compared to the numerical data of ref.. The data are plotted using the dimensionless time $`\rho ^2Dt`$, and the axes are chosen to test the asymptotic form (3). The data lie between the bounds (except at early times where lattice effects are important ). The bounds displayed are for the case $`D^{}=D`$. The very slow convergence of the bounds is consistent with the observed trend in the data.

Some comments on Eq. (12) are in order. First, the right-hand inequality already suggests the possibility that the approach to asymptopia could be very slow, as observed in simulations, since the second term vanishes only as $`t^{1/6}`$. For the case $`D^{}=D`$ that is usually simulated, the asymptotic result (3) gives, for a dimensionless time $`\rho ^2Dt=10^4`$ say, a survival probability of order $`10^{98}`$. This is far smaller than can be reached in simulations (or experiments!) even with sophisticated methods . Yet even at this large time, the second term on the right-hand side of (12) is still 39% as big as the first! This term would give a further factor of $`10^{38}`$ in the survival probability, i.e. a total probability of order $`10^{136}`$! Of course, this second term is only a bound, and the true correction to the asymptotic form could be smaller. The form of the bound is nonetheless suggestive and illustrates how large subdominant terms could arise. A clearer understanding of these subdominant terms is necessary before a detailed comparison with numerical or experimental data can be made.

A second interesting point is that the mean-square displacement of the particle, averaged over the surviving trajectories used to derive the lower bound, grows as $`l^2(D^{}t/\rho )^\nu `$ with $`\nu =2/3`$. This can be compared with the estimate $`\nu =0.5`$ to 0.6 obtained from recent simulations . Note, however, that these simulations were not in the asymptotic regime, so the numerical estimates should be treated with caution. Furthermore, the value $`\nu =2/3`$ obtained from the lower bound for $`P(t)`$ does not necessarily represent any kind of bound for $`\nu `$. Nevertheless, the rough agreement between the measured value and that obtained from our very simple arguments is again suggestive.

The model can be generalized to include $`n`$ diffusing particles starting from the same point . Our bounds can easily be generalized to this case . The probability that all $`n`$ particles survive at time $`t`$ satisfies the same inequalities (12) except that the final term on the right-hand side acquires an additional factor $`n^{1/3}`$. So the bounds converge for $`t\mathrm{}`$ and the same asymptotic form, Eq. (3), is obtained for all $`n`$.

Our approach is also readily generalized to dimensions $`d>1`$ . For $`d<2`$ we find $`P(t)\mathrm{exp}[a_d\rho (Dt)^{d/2}]`$, where $`a_d=(2/\pi d)(4\pi )^{d/2}\mathrm{sin}(\pi d/2)`$ while, for $`d=2`$, $`P(t)\mathrm{exp}[4\pi \rho Dt/\mathrm{ln}t]`$. The latter agrees with the functional form obtained in , but with a precise value for the constant. For $`d>2`$, simple exponential decay is obtained, in agreement with , but the bounds no longer converge and the decay constant cannot be determined.

In the remainder of this Letter we sketch  the derivation of the perturbation theory, for a finite number of traps, that leads to Eq. (4). We recall that the $`N`$ traps move independently according to Langevin equations (5) with noise correlator (6). We call the particle coordinate $`x_0(t)`$. It obeys the Langevin equation $`\dot{x}_0=\eta _0(t)`$, where $`\eta _0(t)`$ is independent of the other noise terms, with correlator $`\eta _0(t)\eta _0(t^{})=2D^{}\delta (tt)`$, corresponding to a particle diffusion constant $`D^{}`$. It is convenient to introduce the relative coordinates $`y_i=x_ix_0`$, $`i=1,\mathrm{},N`$. The corresponding Langevin equations are $`\dot{y}_i=\xi _i(t)`$ where the noise terms, $`\xi _i(t)=\eta _i(t)\eta _0(t)`$, have correlators

$$\xi _i(t)\xi _j(t^{})=2D_{ij}\delta (tt^{}),$$
(13)
where $`D_{ij}=D\delta _{ij}+D^{}`$. The equivalent Fokker-Planck equation, $`P/t=_{i,j=1}^ND_{ij}^2P/y_iy_j`$, has to be solved subject to absorbing boundary conditions at $`y_i=0`$ for any $`i`$ (i.e. when any trap meets the particle). In principle this equation can be solved by a coordinate transformation that diagonalizes the matrix $`D`$ to give isotropic diffusion. This transformation, however, also rotates the edges of the absorbing region so that they are no longer mutually orthogonal but form the edges of an $`N`$-dimensional hyperwedge. For $`N=2`$ this becomes a two-dimensional wedge, and the problem can be solved exactly , but to our knowledge there are no exact solutions for $`N3`$.

An alternative starting point is the path-integral representation for the survival probability of the particle,

$$P(t)=\frac{_RD\stackrel{}{y}(t)\mathrm{exp}(S[\stackrel{}{y}])}{D\stackrel{}{y}(t)\mathrm{exp}(S[\stackrel{}{y}])},$$
(14)
where $`\stackrel{}{y}`$ is a shorthand for $`(y_1,\mathrm{},y_n)`$, $`S[\stackrel{}{y}]=(1/4)_{i,j=1}^N(D^1)_{ij}_0^t𝑑t^{}\dot{y}_i(t^{})\dot{y}_j(t^{})`$, and the subscript $`R`$ indicates that the path integral is restricted to “surviving” paths, in which none of the $`y_i`$ has changed sign up to time $`t`$. The matrix $`D`$ is easily inverted to give $`(D^1)_{ij}=(1/D)(\delta _{ij}\lambda )`$, where $`\lambda =D^{}/(D+ND^{})`$ will be our expansion parameter. Thus

$`S[\stackrel{}{y}]`$ $`=`$ $`S_0[\stackrel{}{y}]S_1[\stackrel{}{y}]`$ (15)
$`S_0[\stackrel{}{y}]`$ $`=`$ $`{\displaystyle \frac{1\lambda }{4D}}{\displaystyle \underset{i}{}}{\displaystyle _0^t}𝑑t^{}[\dot{y}_i(t^{})]^2`$ (16)
$`S_1[\stackrel{}{y}]`$ $`=`$ $`{\displaystyle \frac{\lambda }{2D}}{\displaystyle \underset{i<j}{}}{\displaystyle _0^t}𝑑t^{}\dot{y}_i(t^{})\dot{y}_j(t^{}),`$ (17)
where the diagonal terms, $`i=j`$, have been absorbed into $`S_0`$.

A convenient normalization of the path integral is provided by the path integral for the survival probability, $`P_0(t)t^{N/2}`$, of the same problem with $`D^{}=0`$ and $`DD_0=D/(1\lambda )`$. Then

$$\frac{P(t)}{P_0(t)}=\frac{\mathrm{exp}(S_1[\stackrel{}{y}])_R}{\mathrm{exp}(S_1[\stackrel{}{y}])},$$
(18)
where the average in both cases is over paths weighted by $`\mathrm{exp}(S_0[\stackrel{}{y}])`$.

The final step is to expand numerator and denominator using the cumulant expansion. To first order in $`\lambda `$ (i.e. first order in $`S_1`$) one has

$$P(t)/P_0(t)=\mathrm{exp}(S_1_RS_1).$$
(19)
Since the different values of $`i`$ decouple in $`S_0`$ we have

$$S_1_R=(\lambda /2D)\underset{i<j}{}_0^t𝑑t^{}\dot{y}_i(t^{})_R\dot{y}_j(t^{})_R.$$
(20)
The corresponding expression for the unrestricted paths vanishes, since $`\dot{y}_i(t^{})=0`$ by symmetry. The quantity $`\dot{y}_i(t^{})_R`$ is independent of $`i`$ apart from the sign: $`\dot{y}_i(t^{})_R>0`$ ($`<0`$) for traps that stay to the right (left) of the particle. The calculation of $`\dot{y}_i(t^{})_R`$ is straightforward. We only require the result in the regime $`t_0t^{}t`$, where $`t_0=y_{i0}^2/D`$ and $`y_{i0}`$ is the initial value of trap $`i`$. In this regime one finds (for traps which start on the right of the particle)

$$\dot{y}_i(t^{})_R=2(D_0/\pi t^{})^{1/2},t_0t^{}t.$$
(21)
We now insert this result into Eq. (20). Let there be $`N_L`$ ($`N_R`$) traps to the left (right) of the particle. Then the $`[N_L(N_L1)/2+N_R(N_R1)/2]`$ pairs $`(i,j)`$ of traps whose members are on the same side of the particle contribute with positive weight to (20) while the $`N_LN_R`$ pairs whose members are on opposite sides enter with negative weight. Defining $`\mathrm{\Delta }=N_LN_R`$ (and $`N=N_L+N_R`$) this gives

$$S_1_R=\frac{\lambda D_0}{\pi D}(\mathrm{\Delta }^2N)\mathrm{ln}t$$
(22)
to leading logarithmic accuracy for large $`t`$. The factor $`\mathrm{ln}t`$ comes from the integral $`_{t_0}^t(dt^{}/t^{})`$ that appears when Eq. (21) is substituted into Eq. (20). Using the form (21) for all $`t^{}`$, with $`t`$ and $`t_0`$ as upper and lower cut-offs, is correct to leading logarithmic accuracy . Note that the contributions from the lower cut-off, which depend on the initial positions of the traps, do not contribute to the leading logarithm. Using (22) in (19), recalling that $`S_1=0`$ and $`P_0(t)t^{N/2}`$, and noting that, to leading order in $`D^{}`$, $`\lambda =D^{}/D`$ and $`D_0=D`$, one finds $`P(t)t^\theta `$ with $`\theta `$ given by Eq. (4). The power-law decays hold when $`ty_{0i}^2/D`$ for all $`i`$. The magnitudes of the initial coordinates $`y_{0i}`$ determine the amplitude of the power-law, while their signs determine the exponent through the value of $`\mathrm{\Delta }`$.

Eq. (3) can be compared with the exactly solved cases $`N_L=N_R=1`$ (“surrounded prey”) and $`N_L=0,N_R=2`$ (“chased prey”) , for which $`\theta =\pi /2\alpha `$ and $`[\pi /2(\pi \alpha )]`$ respectively, where $`\alpha =\mathrm{cos}^1[D^{}/(D+D^{})]`$. Expanding these to first order in $`D^{}`$ gives agreement with Eq. (4).

In summary, we have obtained the exact asymptotic form of the survival probability of a particle moving in an infinite background of mobile traps. The asymptotic form is independent of the diffusion constant of the particle. A perturbation expansion in the ratio of the particle and trap diffusion constants has been developed for the case of a finite number of traps.

We thank V. Mehra and P. Grassberger for supplying their raw data, and for a helpful correspondence. This work was supported by EPSRC grant GR/R53197.


format: text

text:
# Untitled Document

INTERACTING LINEAR POLYMERS ON THREE–DIMENSIONAL SIERPINSKI FRACTALS

Jelena Maričić 1 and Sunčica Elezović–Hadžić 2
1 Physics Department, University of Hawaii at Manoa

2 Faculty of Physics, University of Belgrade

Abstract

Using self–avoiding walk model on three–dimensional Sierpinski fractals (3d SF) we have studied critical properties of self–interacting linear polymers in porous environment, via exact real–space renormalization group (RG) method. We have found that RG equations for 3d SF with base $`b=4`$ are much more complicated than for the previously studied $`b=2`$ and $`b=3`$ 3d SFs. Numerical analysis of these equations shows that for all considered cases there are three fixed points, corresponding to the high–temperature extended polymer state, collapse transition, and the low–temperature state, which is compact or semi–compact, depending on the value of the fractal base $`b`$. We discuss the reasons for such different low–temperature behavior, as well as the possibility of establishing the RG equations beyond $`b=4`$.

Introduction

At low temperatures and in poor solvents linear polymers are in compact globule state, while at high temperatures and in good solvents they have extended coil configuration . Self–interacting self–avoiding walk (SISAW) on a lattice is a good model for linear polymer in solvent. If $`\epsilon <0`$ is the attractive energy per pair of nearest–neighbors monomers (not directly connected) and $`w=`$exp$`(\epsilon /k_BT)`$ the corresponding Boltzman factor, then for $`N`$–step SISAW and $`N1`$ the following behavior of the mean-squared end-to-end distance $`R_N^2`$ is expected:

$$R_N^2N^{2\nu }.$$
(1)
There is a critical value of interaction $`w=w_c`$, such that for any $`w<w_c`$ critical exponent $`\nu =\nu _{\mathrm{S}AW}`$, which corresponds to extended coil phase, whereas for any $`w>w_c`$ exponent $`\nu `$ has smaller value, equal to $`1/d`$ for $`d`$–dimensional homogeneous lattices. The value $`w=w_c`$ corresponds to the collapse transition, with $`\nu =\nu _c`$, between the low– and high–temperature $`\nu `$ values: $`1/d<\nu _c<\nu _{\mathrm{S}AW}`$.

On non–homogeneous lattices collapse transition also occurs, with the different values of the exponent $`\nu `$. Three–dimensional Sierpinski fractals (3d SFs) are nonhomogeneous lattices where exact real space renormalization group (RG) method can be applied for studying asymptotic properties of SISAW .

Renormalization group scheme

Each member of the 3d SF family is labelled by $`b=2,3,4,\mathrm{}`$ and can be constructed recursively starting with a tetrahedron of base $`b`$, containing $`b(b+1)(b+2)/6`$ unit tetrahedrons – so called generator $`G_1(b)`$ (see fig.1).

Figure 1: First two stages of the construction of 3d SF with $`b=4`$

The subsequent fractal stages are constructed self–similarly, by replacing each unit tetrahedron of the initial generator by a new generator. To obtain the $`l`$th–stage fractal lattice $`G_l(b)`$, which we shall call the $`l`$th order generator, this process of construction has to be repeated $`l1`$ times, and the complete fractal is obtained in the limit $`l\mathrm{}`$. Fractal dimension $`d_f`$ of 3d SF is equal to

$$d_f=\mathrm{ln}[(b+2)(b+1)b/6]/\mathrm{ln}b.$$
(2)
In this case RG parameters are restricted generating functions, which are for the $`G_l(b)`$ equal to

$$A^{(l)}(x,w)=\underset{N,P}{}A_{N,P}x^Nw^P,B^{(l)}(x,w)=\underset{N,P}{}B_{N,P}x^Nw^P,$$
(3)
where $`x`$ is weight of any step of a SISAW (fugacity), and $`A_{N,P}`$ ($`B_{N,P}`$) is the number of $`N`$–step SISAWs with $`P`$ nearest–neighbor contacts, traversing the $`G_l(b)`$ once (twice) – see fig.2. For elementary tetrahedrons $`G_0(b)`$ one has

$$A^{(0)}(x,w)=x+2x^2w+2x^3w^3,B^{(0)}(x,w)=x^2w^4.$$
(4)
For any $`b`$ RG equations have the form

$`A^{(l+1)}`$ $`=`$ $`{\displaystyle \underset{n,m}{}}𝒜_{n,m}\left(A^{(l)}\right)^n\left(B^{(l)}\right)^m,`$
$`B^{(l+1)}`$ $`=`$ $`{\displaystyle \underset{n,m}{}}_{n,m}\left(A^{(l)}\right)^n\left(B^{(l)}\right)^m,`$ (5)
where $`𝒜_{n,m}`$ ($`_{n,m}`$) is the number of SISAW configurations consisting of $`n`$ branches of ”$`A^{(l)}`$”—type and $`m`$ branches of ”$`B^{(l)}`$”—type, traversing $`G_{l+1}(b)`$ once (twice).

Figure 2: Examples of $`A`$ and $`B`$ type of SISAW on $`G_1(2)`$

Analysis of the $`b=4`$ case

As $`b`$ increases the number of possible SISAW configurations quickly grows, which is the reason why we don’t quote the corresponding extremely cumbersome RG equations  for the studied $`b=4`$ case here. Numerical analysis of these equations shows that there are three physically relevant fixed points.

(1) The high–temperature fixed point ($`A^{}=0.289912,B^{}=0.0122393\mathrm{}`$) is reached for any $`w<w_c=5.028\mathrm{}`$ and corresponding critical fugacity $`x=x^{}(w)`$. Linearizing RG equations around this fixed point one gets only one relevant eigenvalue $`\lambda _1=8.6924\mathrm{}`$, and the critical exponent $`\nu `$ is equal to $`\nu _{\mathrm{S}AW}=\mathrm{ln}b/\mathrm{ln}\lambda _1=0.641076\mathrm{}`$.

(2) The low–temperature fixed point ($`A^{}=0,B^{}=0.356883\mathrm{}=22^{1/3}`$) is reached when starting with $`w>w_c`$. There is one relevant eigenvalue $`\lambda _1=16`$, and the critical exponent is $`\nu =1/2`$.

(3) The tricritical fixed point ($`A^{}=0.19292\mathrm{},B^{}=0.3388\mathrm{}`$) corresponds to the collapse transition, when $`w=w_c`$. In this case there are two relevant eigenvalues: $`\lambda _1=15.4194\mathrm{}`$ and $`\lambda _2=5.52427\mathrm{}`$. The critical exponent $`\nu `$ is equal to $`\nu _c=\mathrm{ln}b/\mathrm{ln}\lambda _1=0.506756\mathrm{}`$.

Discussion

Values of $`\nu `$ for SISAW on all so far studied 3d SF \[2-4\] are depicted on fig. 3. One can see that for any fractal $`\nu `$ has its largest value in the high–temperature phase, and that it decreases as $`b`$ grows, approaching the corresponding Euclidean value $`\nu _{\mathrm{SAW}}^E0.588`$ . Respective $`\nu `$ values are smaller at the collapse transition, also decreasing with $`b`$, being close to the collapse Euclidean value $`\nu _c^E=0.5`$ . Low–temperature behavior is more interesting. In contrast to the $`b=2`$ case for which $`\nu =1/d_f`$ (meaning that polymer densely fills the underlying lattice, as on homogeneous lattices), critical exponent $`\nu `$ for $`b=3`$ and $`b=4`$ is larger than $`1/d_f`$, indicating that polymer is in semi–compact, rather than compact phase. This is due to the topological frustration, which is not the same for $`b=3`$ and $`b=4`$, being also the reason for strange behavior of the low-temperature RG fixed point: $`B^{}`$ is equal to $`22^{1/3}`$ for both $`b=2`$ and $`b=4`$, whereas for $`b=3`$ it is $`B^{}=\mathrm{}`$ (of course, $`A^{}=0`$ for any $`b`$). It would be interesting to find out what happens beyond $`b=4`$, but unfortunately the number of SISAW configurations grows extremely fast with $`b`$. Estimated computer time required for enumeration and sorting SISAWs for $`b=5`$ is not accessible, so one should search other methods for calculating $`\nu `$.

Figure 3: Critical exponent $`\nu `$ for 3d SFs with $`b=2,3`$, and 4

References

 P.G. de Gennes, ”Scaling Concepts in Polymer Physics” (Itacha, NY: Cornell University Press, 1979)

 D. Dhar and J. Vannimenus, J.Phys. A 20 (1987) 199

 M. Knežević and J. Vannimenus, J.Phys. A 20 (1987) L969

 J. Maričić, ”Random Walks on Three–dimensional Sierpinski Fractals” – diploma work (Faculty of Physics, University of Belgrade, 1998)

 D. Lee et al., J.Stat.Phys. 80 (1995) 661; D. MacDonald et al., J.Phys. A 33 (2000) 5973


text:
# EROS 2 proper motion survey: Constraints on the halo white dwarfs Based on observations made with the Marly telescope, located at the European Southern Observatory, La Silla, Chile.

## 1 Introduction

The rotation of the Milky Way is believed to be sustained by a massive dark halo, at least beyond one disk scale length. Old, cool white dwarfs (WDs) are one of the proposed constituents of that dark matter. A number of constraints have been placed on the total amount of halo white dwarfs (HWDs). One is based on the light emitted by WD progenitors (Graff et al. 1999). Extragalactic metallicity measurements also place strong constraints through the metal and helium enrichment due to the HWDs (Freese 2000), although mechanisms have been proposed to evade those constraints.

Over the past three years a great excitement has come from proper motion surveys used to find fast WDs in the Solar neighbourhood. The surveys conducted by Luyten (1979a) and Knox et al. (1999) are jointly sensitive to $`1500\mathrm{pc}^3`$ for $`M_V=17.5`$ HWDs (Flynn et al. 2001). They find no halo WD (Liebert et al. 1988), although this is disputed in the case of Luyten’s survey, whose efficiency is still a matter of strong debate, 23 years after publication. On the other hand, several collaborations have reported proper motion detection of fast, cool WDs some of them being subsequently confirmed spectroscopically (Ibata et al. 2000; Oppenheimer et al. 2001a). Nelson et al. (2002) also reported the discovery of new HWD candidates in the HST Groth–Westphal strip, although using only two epochs. Finally, Monet et al. (2000) discovered new high velocity WDs but the small sample of new objects confirms the high efficiency of Luyten (1979b), thus decreasing the maximum density of halo WDs.

On the microlensing side, Lasserre et al. (2000) put a 30% upper limit (95% C.L.) on the contribution of HWDs of any age – or any compact object of similar mass – to a standard halo, while Alcock et al. (2000) interpret their 13 microlensing candidates towards the Large Magellanic Cloud as 0.15 – $`0.9\mathrm{M}_{\mathrm{}}`$ objects contributing $`20_6^{+11}`$% to the standard halo.

The interesting controversy developed around the Oppenheimer et al. (2001b) results addresses the question of the nature of their findings: thick disk or halo objects. It is in any case agreed that these relatively bright $`M_I14`$, young WDs contribute very little to the Galactic dark matter, a few percent at most. It has been suggested however that these could be the bright tail of a fainter, older population (Hansen 2001).

Our proper motion survey began in 1996 using the EROS 2 cameras. It was designed to study the contribution of older, fainter WDs to the halo, and also addresses the question raised by Hansen (2001). Like other proper motion searches, ours takes advantage of the large expected velocity dispersion of halo objects, and it is not intended to address the question of the low velocity, or brighter, WDs that most recent surveys probe (Oppenheimer et al. 2001a; Majewski & Siegel 2002). Here we report our results based on the analysis of $`250^2`$ that were observed during at least three epochs. These yield a sensitivity about three times that of Luyten (1979b), depending on HWD colours and luminosities.

## 2 The EROS 2 proper motion survey

### 2.1 Instrument and observations

The Eros 2 wide-field imager (Bauer & de Kat 1998) was designed to search for microlensing effects towards the Galactic bulge and disk and the Magellanic Clouds. Its two $`1^2`$ CCD cameras are mounted at the Cassegrain focus of the 1-m Marly telescope at La Silla (Chile), with a pixel size of $`0\mathrm{}.6`$. These two $`8k\times 4k`$ mosaics are illuminated through a dichroic beam splitter, which, together with the CCD efficiency, defines the bandpasses. The colour transformations between the Eros system ($`𝒱_\mathrm{E}`$, $`_\mathrm{E}`$) and the standard Johnson-Cousins $`(V,I)`$ system are determined by observing Landolt (1992) standards and Ogle (Paczynski et al. 1999) secondary standards in the Baade window: $`𝒱_\mathrm{E}=V(0.32\pm 0.03)(VI)`$ and $`_\mathrm{E}=I+(0.06\pm 0.03)(VI)`$. We estimate the precision of the colour transformations to be $`0.1^{\mathrm{mag}}`$ for $`VI=`$0 – 2, the range of expected colours for HWDs. Comparison between our $`_\mathrm{E}`$ photometry and I as measured by Denis (Fouqué et al. 2000) reveals a small offset, $`I_{\mathrm{Eros}}=I_{\mathrm{Denis}}0.06`$, while a comparison between our photometry and the Photometric GSC (Lasker et al. 1988) reveals no significant offset, with a $`0.05^{\mathrm{mag}}`$ uncertainty.

Proper motion observations were performed one to two hours per dark night, within 90 minutes of the meridian to minimize atmospheric refraction. Our limiting magnitudes are about $`V=21.5`$ and $`I=20.5`$. We observed $`190^2`$ in the Southern Galactic Hemisphere ($`79^{}<b_{\mathrm{gal}}<48^{}`$), in the following strips along the $`\alpha `$ coordinate, $`\mathrm{\Delta }\delta =1.4\mathrm{deg}`$ wide: 22 h 16 min$`<\alpha <3`$h 44 min at $`\delta =44^{}\mathrm{\hspace{0.17em}45}`$min, 23 h 31 min$`<\alpha <1`$h 34 min at $`\delta =40^{}\mathrm{\hspace{0.17em}09}`$min and 22 h 24 min$`<\alpha 3`$h 28 min at $`\delta =38^{}\mathrm{\hspace{0.17em}45}`$min, and $`187^2`$ in the Northern Hemisphere ($`41^{}<b_{\mathrm{gal}}<59^{}`$) with 10 h 57 min$`<\alpha 13`$h 23 min at $`\delta =12^{}`$ and 10 h 57 min$`<\alpha <12`$h 53 min at $`\delta =4^{}\mathrm{\hspace{0.17em}36}`$min.

Only those fields ($`250^2`$) with three epochs separated by one-year intervals (two-year total baseline), in one or two bands, are taken into account in this letter.

### 2.2 Proper motion catalogue

The reduction software for source detection, classification and catalogue matching was written in the framework of the Eros Peida++ package (Ansari 1996). Since photon noise dominates the astrometric errors for most of the search volume, we use a simple two-dimensional Gaussian PSF-fitting algorithm to determine stellar positions. A rough star/galaxy classification is performed to limit galaxy contamination, with cuts chosen so that few stars are misclassified. The catalogues of the three epochs are geometrically aligned, with the deepest one taken as reference, using a linear transformation fitted to the 40 brightest stars, and matched within a search radius corresponding to a $`6^{\prime \prime }\mathrm{yr}^1`$ proper motion. The RMS distance between matched stars provides an upper limit to the total astrometric error, which is 25 mas (1$`\sigma `$) for bright objects, degrading to 150 mas for $`V=21`$ or $`I=20`$.

In order to minimize contamination by spurious candidates we require that all objects be observed at least three times, separated by one-year intervals. The proper motions are derived from linear motion fits. The $`\chi ^2`$ distributions deviate from the expected ones only for confidence levels lower than 0.5%, indicating the onset of systematic effects. The 4% of objects so affected, mostly poorly measured galaxies, stars close to CCD defects or visual binaries, are removed from the analysis.

## 3 Search for halo white dwarfs

### 3.1 Selection of halo candidates

The choice of the criteria used to distinguish halo objects from the much more numerous disk and thick disk stars, and from extragalactic objects, is of prime importance. We apply two cuts: We first require that the halo candidates detected in the $`𝒱_\mathrm{E}`$ band exhibit a reduced proper motion, $`H_V=V+5\mathrm{log}\frac{\mu }{1^{\prime \prime }\mathrm{yr}^1}+5=M_V+5\mathrm{log}\frac{v_{}}{1\mathrm{km}\mathrm{s}^1}3.38>22.5`$, where $`\mu `$ is the proper motion and $`v_{}`$ the tangential velocity, and that those detected in $`_\mathrm{E}`$ have $`H_I>21.5`$. This cut select candidates that are both fast and faint (see Fig.1).

We then require a proper motion $`\mu >0\mathrm{}.7\mathrm{yr}^1`$. This cut removes objects whose proper motion reality cannot be confirmed by examining our images. It affects our sensitivity, depending on the mean WD magnitude: $`\mu _{\mathrm{min}}=0\mathrm{}.7\mathrm{yr}^1`$ corresponds to a transverse velocity of $`v_{}=\frac{60\mathrm{pc}}{d}\times 200\mathrm{km}\mathrm{s}^1`$ at distance $`d`$. A typical $`200\mathrm{km}\mathrm{s}^1`$ transverse velocity is expected for halo objects, and $`d=60`$pc is the typical distance of WDs in our effective volume for $`M_V=17`$ stars. In fact, a detailed simulation shows that this cut reduces our sensitivity by 20% (resp. 45%) for $`M_V=18`$ (resp. 17) HWDs.

The proper motion cut was set to $`\mu >0\mathrm{}.7\mathrm{yr}^1`$ as a result of the following optimization: on the one hand this value is low enough to ensure a high detection efficiency of HWDs and, on the other hand, it is sufficiently high to eliminate candidates which are likely to be spurious because of a poor measurement of their proper motion. Indeed, when the proper motion cut is lowered to $`0\mathrm{}.6\mathrm{yr}^1`$, 10 candidates pass the cut, but for those whose proper motions were measured in both bands (6 among the 10 candidates) we find that for all of them the measured values of their visible and red proper motion are incompatible.

No candidates survive the two cuts: $`H_V>22.5`$ or $`H_I>21.5`$, and $`\mu >0\mathrm{}.7\mathrm{yr}^1`$. Had we required $`\mu >0\mathrm{}.8\mathrm{yr}^1`$, our upper limits on the HWD contribution to the halo given below would move up by less than 1% of the standard dark halo.

### 3.2 Comparison with Luyten’s catalogue

We cross match our proper motion catalogue with the NLTT catalogue (Luyten 1979b), by searching for bright stars in our catalogue with proper motion higher than $`0\mathrm{}.1\mathrm{yr}^1`$, within 1’ of Luyten positions. It turns out that 85% of Luyten’s stars are recovered with compatible proper motion and proper motion angle, which demonstrates our high detection efficiency for stars similar to Luyten’s stars. The 15% loss is partly due to the expected loss due to defects in the CCDs (4%) and to pointing dispersion between epochs (4%). Additionally, our proper motion measurements for Luyten’s stars are significant at the $`3\sigma `$ level in each band, and are compatible at the $`2\sigma `$ level. This gives us confidence in our proper motion accuracy, as most Luyten stars have proper motion much lower than our $`0\mathrm{}.7\mathrm{yr}^1`$ halo selection cut. However, this only applies to bright stars, as there are no Luyten stars close to our photometric detection limit. Unfortunately we can draw no conclusion regarding Luyten’s completeness based on this study, as we lack a reliable estimate of the contamination of our disk-like proper motion catalogue. This study will be reported in more detail in a forthcoming article.

## 4 Discussion

### 4.1 Halo model predictions

In order to check our sensitivity to HWDs, we perform a full simulation of our observations. We measure the detection efficiency of a star on each frame by adding simulated stars to real images. Using the Besançon model of the Galaxy (Robin et al. 2000), to which we add one HWD per field, we then construct a simulate catalogue for each of our fields and each of our images, and process these catalogues in the same way as our observations. The Besançon model is a consistent hydrodynamical description of Galactic star counts and proper motion, based on initial mass function and stellar evolution.

We use various halo models with different kinematics. We use isotropic velocity dispersions of $`\sigma _{1\mathrm{D}}=100`$, 130 or $`156\mathrm{km}\mathrm{s}^1`$, or anisotropic, spheroid-like dispersions, with $`\sigma _{1\mathrm{D}}=85`$, 105 or $`125\mathrm{km}\mathrm{s}^1`$ along the $`V`$ (Galactic rotation) direction. We suppose halo rotations of $`50`$, 0 or $`50\mathrm{km}\mathrm{s}^1`$. The expected number of detections depend on these parameters at the 5 to 10% level. For a given kinematics, flattened halos correspond to higher local densities, and thus to proportionately more detections.

We consider HWDs of $`16.5M_V18`$ and $`0.5VI1.5`$, in steps of 0.5 mag. We also generate HWDs with hydrogen atmosphere according to the Chabrier (1999) luminosity functions for halo ages of 14 and 15 Gyr. The predicted number of detections, in either the $`𝒱_\mathrm{E}`$ or the $`_\mathrm{E}`$ band, are reported in Table 1, for the standard isotropic, isothermal halo, which has a local density of $`8.10^3\mathrm{M}_{\mathrm{}}\mathrm{pc}^3`$. Statistical errors are 5% while systematic errors are estimated to be less than 25%, the quadrature sum of errors in calibration (15%), efficiencies (10%) and parameters of the halo kinematics (15%).

We have not taken binarity effects into account. Double stars of similar magnitude would appear brighter than individual WDs, so that objects at the brighter, slower end of the distribution could be lost by the reduced proper motion cut. Furthermore, the number of systems to which we are sensitive would be smaller than the actual number of stars; our constraints would have to be relaxed accordingly.

### 4.2 Constraints on the halo WD fraction and age

The constraints on the contribution of WD to the halo are obtained from the above predictions and the lack of halo candidates, and shown in Fig. 2. We find that a 14 Gyr halo cannot be made of more than 5% 0.6 or $`1.0\mathrm{M}_{\mathrm{}}`$ HWDs, at the 95% C.L., while a 15 Gyr halo cannot be made of more than 15% of $`0.6\mathrm{M}_{\mathrm{}}`$ HWDs, or 45% of $`1.0\mathrm{M}_{\mathrm{}}`$ HWDs.

### 4.3 Comparison with other surveys

Our results, based on a larger volume than the study by Flynn et al. (2001), confirm their conclusion. They are marginally compatible with results of Ibata et al. (2000), depending on the exact WD colours and absolute magnitudes.

Regarding Oppenheimer et al. (2001a), we do not exclude a 1 to 2% contribution by bright $`M_I`$=13 – 15 HWDs. We stress that our selection cut, $`\mu >0\mathrm{}.7\mathrm{yr}^1`$, is more severe than theirs, while preserving most of the halo sensitivity for the faint magnitudes of interest here. Reid et al. (2001), Reylé et al. (2001) and Koopmans & Blandford (2002) also discuss this survey. We do constrain the proposition of (Hansen 2001) that the WDs detected by Oppenheimer et al. (2001a) be the bright tail of a fainter population. If this population has a halo-like kinematics, then the older part of the luminosity function, older by several Gyrs, cannot make up more than $`7\%`$ of the halo.

Finally, our constraints are fully compatible with the Eros microlensing results. The constraints from the present study are even stronger than those from microlensing, for HWDs of absolute magnitude $`M_V<18`$ and also if the halo age is smaller than 15 Gyr (using the Chabrier (1999) luminosity function for DAs).

## 5 Conclusion

In this letter we have presented the analysis of $`250^2`$ of the Eros proper motion survey, dedicated to the search of old, faint white dwarfs. We found no halo white dwarf candidates. We exclude a white dwarf contribution to the halo above 5% for $`M_V=17`$, $`1VI1.5`$, at the 95% C.L. Assuming a Chabrier (1999) cooling model, a 14 Gyr halo cannot be made of more than 5% of WDs with hydrogen atmosphere. Our limits degrade to 20% for a 15 Gyr halo of $`1\mathrm{M}_{\mathrm{}}`$ white dwarfs or for $`M_V=18`$, $`0VI1.5`$ white dwarfs.

###### Acknowledgements.

We thank Annie Robin and her collaborators for allowing us to use the latest version of the Besançon model prior to publication. We thank also Nicolas Epchtein and the Denis collaboration for providing infrared unpublished photometry of our sources in order to check our $`_\mathrm{E}`$ calibration, and Gilles Chabrier for his comments. This research has made use of the Simbad database, operated at CDS, Strasbourg, France. We are grateful to Daniel Lacroix and the technical staff at the Observatoire de Haute Provence and to André Baranne for their help in refurbishing the MARLY telescope and remounting it in La Silla. We also thank the technical staff of ESO La Silla for their support of Eros. We thank Jean-François Lecointe and Adelino Gomes for assistance with the online computing. Work by A.G. was supported by NSF grant AST 97-27520 and by a grant from le Centre français pour l’accueil et les échanges internationaux.


text:
# Partition function zeros of the 𝑄-state Potts model on the simple-cubic lattice

## I introduction

The $`Q`$-state Potts modelwu  in two and three dimensions exhibits a rich variety of critical behavior and is very fertile ground for the analytical and numerical investigation of first- and second-order phase transitions. With the exception of the two-dimensional $`Q=2`$ Potts (Ising) model in the absence of an external magnetic field, exact solutions for arbitrary $`Q`$ are not known. However, some exact results at the critical temperature have been established for the two-dimensional $`Q`$-state Potts model for no magnetic field. From the duality relation the ferromagnetic critical temperature is known to be $`T_c=J/k_B\mathrm{ln}(1+\sqrt{Q})`$ for the isotropic square lattice. The ferromagnetic Potts model has a first-order phase transition for $`Q>4`$ in two dimensionsbaxter73 , and has been the most important means to understand first-order transitions which are omnipresent in nature, not only in condensed-matter physics but also in subatomic physics, for example, the deconfining or chiral symmetry breaking transition in quantum chromodynamicsbrown ; alves90a  and phase transitions in the very early Universekolb . In three dimensions the ferromagnetic $`Q`$-state Potts model for $`Q>Q_c2.45(1)`$jlee  has a first-order transition, and the three-state Potts model is related to the order parameter (the expectation value of the Wilson line) for the deconfinement transition in quantum chromodynamicssvetitsky .

On the other hand, the antiferromagnetic Potts model is much less well understood than the ferromagnetic model. One of the most interesting properties of the antiferromagnetic model is that for $`Q>2`$ the ground-state is highly degenerate and the ground-state entropy is nonzero. The Baxter conjecturebaxter82  for the critical temperature of the Potts antiferromagnet ($`J<0`$) on the square lattice, $`T_c=J/k_B\mathrm{ln}(\sqrt{4Q}1)`$, gives the known exact value for $`Q=2`$, a critical point at zero temperature for $`Q=3`$, and no critical point for $`Q>3`$. In three dimensions the Potts model on the simple-cubic lattice has the antiferromagnetic ordered phase for $`Q4`$banavar . It is now commonly accepted that the order-disorder phase transition of the three-state Potts antiferromagnet on the simple-cubic lattice belongs to the universality class of the three-dimensional $`XY`$ modelbanavar ; wang . Recently Rosengren et al.rosengren93 ; kundrotas  claimed the existence of two additional phase transitions below the order-disorder phase transition for the simple-cubic three-state Potts antiferromagnet. However, Kolesik and Suzukikolesik , Heilmann et al.heilmann , and Oshikawaoshikawa  reported different results in disagreement with observations by Rosengren’s group.

By introducing the concept of the zeros of the partition function in the complex magnetic-field plane, Yang and Leeyang  proposed a mechanism for the occurrence of phase transitions in the thermodynamic limit and yielded a new insight into the unsolved problem of the Ising model in an arbitrary nonzero external magnetic field. It has been shownyang ; tdlee ; creswick97 ; creswick98 ; creswick99 ; janke01  that the distribution of the zeros of a model determines its critical behavior. Lee and Yangtdlee  also formulated the celebrated circle theorem which states that the partition function zeros of the Ising ferromagnet lie on the unit circle in the complex magnetic-field ($`X=e^{\beta H}`$) plane. In three dimensions the properties of the zeros in the complex $`X`$ plane of the Ising model have been studied by exact enumerationsuzuki  and series expansionkortman . However, for the $`Q`$-state Potts model with $`Q>2`$ the partition function zeros in the complex $`X`$ plane lie close to, but not on, the unit circle with the two exceptions of the critical point $`X=1`$ ($`H=0`$) itself and the zeros in the limit $`T=0`$kim98a ; kim99 ; kim00a .

Fisherfisher65  emphasized that the partition function zeros in the complex temperature ($`Y=e^{\beta J}`$) plane are also very useful in understanding phase transitions, and showed that for the square lattice Ising model in the absence of an external magnetic field the zeros in the complex $`Y`$ plane lie on two circles (the ferromagnetic circle $`Y_{FM}=1+\sqrt{2}e^{i\theta }`$ and the antiferromagnetic circle $`Y_{AFM}=1+\sqrt{2}e^{i\theta }`$) in the thermodynamic limit. In particular, using the Fisher zeros both the ferromagnetic phase and the antiferromagnetic phase can be considered at the same time. Fisher also showed that the logarithmically infinite specific heat singularity of the Ising model results from the properties of the density of zeros. The zeros in the complex $`Y`$ plane of the simple-cubic Ising model have been studied extensively to understand the analytical properties and critical behavior of the model by exact enumerationono ; pearson ; itzykson ; martin83a ; bhanot90a  and to determine the critical point and exponents by Monte Carlo simulationmarinari ; bhanot87 ; alves90b ; alves00 ; janke01 . Recently it has been shown that for self-dual boundary conditions near the ferromagnetic critical point $`Y_c=1/(1+\sqrt{Q})`$ the zeros of the $`Q`$-state Potts model on a finite square lattice in the absence of an external magnetic field lie on the circle with center $`1/(Q1)`$ and radius $`\sqrt{Q}/(Q1)`$ in the complex $`Y`$ plane, while the antiferromagnetic circle of the Ising model completely disappears for $`Q>2`$martin86 ; cnchen ; kim98b ; kim00b ; kim01 . Using the distribution of the zeros in the complex temperature plane of the square-lattice Potts model, the critical behavior in both the ferromagnetic and antiferromagnetic phases has been studiedkim98b ; kim00b ; kim01 , and the unknown properties of the model for nonzero external magnetic field have been invetigatedkim98b . The partition function zeros in the complex $`Y`$ plane of the three-state Potts model on the simple-cubic lattice have been studied to understand the ferromagnetic critical behavior by exact enumeration for the $`3\times 3\times 9`$ simple-cubic latticemartin83b  and by Monte Carlo simulationalves91 ; janke01 .

In this paper we study the partition function zeros of the $`Q`$-state Potts model for $`2Q7`$ on the simple-cubic lattice to unveil some of the rich structures of the model. The partition function zeros of the three-dimensional Potts model have never been obtained for $`Q>3`$ in the literature. In the next section we describe briefly the microcanonical transfer matrix to evaluate the density of states, from which the exact partition function of the simple-cubic Potts model is obtained. In Section 3 we discuss the ferromagnetic and antiferromagnetic phase transitions and the complex-temperature singularity (nonphysical critical point) of the simple-cubic Potts model using the partition function zeros in the absence of an external magnetic field. We also discuss the properties of the complex-temperature singularity using the low-temperature series expansion. In section 4 we study the unknown properties of the simple-cubic Potts model for nonzero magnetic field using the partition function zeros and the high-field, low-temperature series expansion. We discuss the Yang-Lee edge singularity and the Fisher edge singularity, which are nonphysical critical points for nonzero magnetic field, and relationship between the Yang-Lee and Fisher edge singularities.

## II density of states

The $`Q`$-state Potts model in an external magnetic field $`H_q`$ on a lattice $`G`$ with $`N_s`$ sites and $`N_b`$ bonds is defined by the Hamiltonian

$$_Q=J\underset{i,j}{}\delta (\sigma _i,\sigma _j)H_q\underset{k=1}{\overset{N_s}{}}\delta (\sigma _k,q),$$
(1)
where $`J`$ is the coupling constant (ferromagnetic model for $`J>0`$ and antiferromagnetic model for $`J<0`$), $`i,j`$ indicates a sum over nearest-neighbor pairs, $`\delta `$ is the Kronecker delta, $`\sigma _i=1,2,\mathrm{},Q`$, and $`q`$ is a fixed integer between 1 and $`Q`$. The partition function of the model is

$$Z_Q=\underset{\{\sigma _n\}}{}e^{\beta _Q},$$
(2)
where $`\{\sigma _n\}`$ denotes a sum over $`Q^{N_s}`$ possible spin configurations and $`\beta =(k_BT)^1`$. If we define the restricted density of states with energy $`0EN_b`$ and magnetization $`0MN_s`$ by

$$\mathrm{\Omega }_Q(E,M)=\underset{\{\sigma _n\}}{}\delta (E\underset{i,j}{}[1\delta (\sigma _i,\sigma _j)])\delta (M\underset{k}{}\delta (\sigma _k,q)),$$
(3)
which takes on only integer values, then the partition function can be written as

$$Z_Q(Y,X)=Y^{N_b}\underset{E=0}{\overset{N_b}{}}\underset{M=0}{\overset{N_s}{}}\mathrm{\Omega }_Q(E,M)Y^EX^M,$$
(4)
where $`Y=e^{\beta J}`$ and $`X=e^{\beta H_q}`$, and states with $`E=0`$ ($`E=N_b`$) correspond to the ferromagnetic (antiferromagnetic) ground states. From Eq. (4) it is clear that $`Z_Q(Y,X)`$ is simply a polynomial in $`Y`$ and $`X`$. If we sum the restricted density of states over $`M`$, we obtain the density of states

$$\mathrm{\Omega }_Q(E)=\underset{M=0}{\overset{N_s}{}}\mathrm{\Omega }_Q(E,M),$$
(5)
whose sum over $`E`$ is equal to the $`Q^{N_s}`$ spin configurations,

$$\underset{E=0}{\overset{N_b}{}}\mathrm{\Omega }_Q(E)=Q^{N_s}.$$
(6)
In the absence of an external magnetic field, the partition function of the model is given by

$$Z_Q(Y)=Y^{N_b}\underset{E=0}{\overset{N_b}{}}\mathrm{\Omega }_Q(E)Y^E,$$
(7)
which is again a polynomial in $`Y`$.

Here we describe briefly the microcanonical transfer matrix ($`\mu `$TM)bhanot90a ; bhanot90b ; stosic ; stodolsky ; creswick95 ; kim01  to obtain exact integer values for the density of states $`\mathrm{\Omega }_Q(E)`$ of the $`Q`$-state Potts model on an $`L\times L\times N`$ simple-cubic lattice with periodic boundary conditions in the $`x`$-$`y`$ plane (horizontal area of $`L\times L`$) and free boundaries in the $`z`$ direction (height of $`N`$). First, an array, $`\omega ^{(1)}`$, which is indexed by the energy $`E`$ and variables $`\sigma _i`$, $`1iL^2`$ for the first $`x`$-$`y`$ plane of sites is initialized as

$$\omega ^{(1)}(E;\sigma _1,\sigma _2,\mathrm{},\sigma _{L\times L})=\delta (E\underset{i=1}{\overset{L\times L}{}}[1\delta (\sigma _i,\sigma _{i+1})]).$$
(8)
Now each spin in the plane is traced over in turn, introducing a new spin variable from the next $`x`$-$`y`$ plane,

$$\stackrel{~}{\omega }^{(2)}(E;\sigma _1^{},\sigma _2,\mathrm{},\sigma _{L\times L})=\underset{\sigma _1}{}\omega ^{(1)}(E[1\delta (\sigma _1^{},\sigma _1)];\sigma _1,\sigma _2,\mathrm{},\sigma _{L\times L}).$$
(9)
This procedure is repeated until all the spins in the first plane have been traced over, leaving a new function of the $`L\times L`$ spins in the second plane. The bonds connecting the spins in the second $`x`$-$`y`$ plane are then taken into account by shifting the energy,

$$\omega ^{(2)}(E;\sigma _1^{},\sigma _2^{},\mathrm{},\sigma _{L\times L}^{})=\stackrel{~}{\omega }^{(2)}(E\underset{i=1}{\overset{L\times L}{}}[1\delta (\sigma _i^{},\sigma _{i+1}^{})];\sigma _1^{},\sigma _2^{},\mathrm{},\sigma _{L\times L}^{}).$$
(10)
This procedure is then applied to each $`x`$-$`y`$ plane in turn until the final ($`N`$th) plane is reached. The density of states is then given by

$$\mathrm{\Omega }_Q(E)=\underset{\sigma _1^{}}{}\underset{\sigma _2^{}}{}\mathrm{}\underset{\sigma _{L\times L}^{}}{}\omega ^{(N)}(E;\sigma _1^{},\sigma _2^{},\mathrm{},\sigma _{L\times L}^{}).$$
(11)
The permutation symmetry of the $`Q`$-state Potts model allows us to freeze the last spin $`\sigma _{L\times L}=1`$ of each $`x`$-$`y`$ plane. Now we need to consider only $`Q^{L\times L1}`$ possible spin configurations in each plane instead of $`Q^{L\times L}`$ configurations, and we save a great amount of memory and CPU time. For example, Table I shows the density of states $`\mathrm{\Omega }_{Q=7}(E)`$ of the seven-state Potts model on the $`3\times 3\times 3`$ simple-cubic lattice.

The above procedure can be modified in a straightforward way to calculate exact integer values for the restricted density of states $`\mathrm{\Omega }(E,M)`$stosic ; creswick95 . For the first $`x`$-$`y`$ plane one has, instead of Eq. (8),

$$\omega ^{(1)}(E,M;\sigma _1,\sigma _2,\mathrm{},\sigma _{L\times L})=\delta (E\underset{i=1}{\overset{L\times L}{}}[1\delta (\sigma _i,\sigma _{i+1})])\delta (M\underset{k=1}{\overset{L\times L}{}}\delta (\sigma _k,q)).$$
(12)
The spins in the first plane are traced over in sequence just as in Eq. (9) leading to a new function of the spins in the second $`x`$-$`y`$ plane;

$$\stackrel{~}{\omega }^{(2)}(E,M;\sigma _1^{},\sigma _2,\mathrm{},\sigma _{L\times L})=\underset{\sigma _1}{}\omega ^{(1)}(E[1\delta (\sigma _1^{},\sigma _1)],M\delta (\sigma _1^{},q);\sigma _1,\sigma _2,\mathrm{},\sigma _{L\times L}).$$
(13)
Once all the spins in an $`x`$-$`y`$ plane have been traced over, the function for the next plane is completed as in Eq. (10)

$$\omega ^{(2)}(E,M;\sigma _1^{},\sigma _2^{},\mathrm{},\sigma _{L\times L}^{})=\stackrel{~}{\omega }^{(2)}(E\underset{i=1}{\overset{L\times L}{}}[1\delta (\sigma _i^{},\sigma _{i+1}^{})],M;\sigma _1^{},\sigma _2^{},\mathrm{},\sigma _{L\times L}^{}),$$
(14)
and the restricted density of states is given by

$$\mathrm{\Omega }_Q(E,M)=\underset{\sigma _1^{}}{}\underset{\sigma _2^{}}{}\mathrm{}\underset{\sigma _{L\times L}^{}}{}\omega ^{(N)}(E,M;\sigma _1^{},\sigma _2^{},\mathrm{},\sigma _{L\times L}^{}).$$
(15)
## III partition function zeros in the complex temperature plane

In this section we discuss the properties of the $`Q`$-state Potts models on the simple-cubic lattice using the partition function zeros $`\{Y_0\}`$ or $`\{A_0\}`$ of the models in the complex temperature ($`Y=e^{\beta J}`$ or $`A=e^{\beta J}`$) plane which are obtained by solving the equation

$$Z_Q(Y_0)=Y_0^{N_b}\underset{E=0}{\overset{N_b}{}}\mathrm{\Omega }_Q(E)Y_0^E=0$$
(16)
for the exact partition functions on finite lattices. Figure 1 shows the partition function zeros in the complex $`Y=e^{\beta J}`$ plane of the $`Q`$-state Potts models for $`Q`$=2, 4, and 7 on the $`3\times 3\times 3`$ simple-cubic lattice with periodic boundary conditions in the $`x`$-$`y`$ plane and free boundary conditions in the $`z`$ direction. In Figure 1, near the positive real axis, the partition function zeros show one-dimensional distributions for all the Potts models which become denser as $`Q`$ increases. Since the value of latent heat of a system is proportional to the density of zeros in the complex temperature planecreswick98 , for $`Q3`$ the latent heat increases at the transition temperature as Q increases, and the first-order phase transition becomes stronger. Moreover, the first zeros of these linear distributions will cut the real axis and be the ferromagnetic (FM) critical points in the thermodynamic limit. In Figure 1, as $`Q`$ increases, the first zeros become closer to the real axis and to the values of the FM critical temperatures estimated in the literature (Table II). This fact suggests that the first zero obtained for a smaller size of simple-cubic lattice approach very closely the real critical point in the limit $`Q\mathrm{}`$ where the mean-field theory is known to be exactmittag ; pearce .

Even for the three-dimensional Ising model the critical temperature has never been known exactly except for some conjectures, for example, the Rosengren conjecturerosengren1 . However, for the simple-cubic lattice Hajduković conjecturedhajdukovic  the FM critical temperature

$$Y_c^{FM}(Q)=\frac{1}{(Q+\sqrt{2Q1})^{1/3}}$$
(17)
of the $`Q`$-state Potts model based on two assumptions. The first assumption is that $`Y_c^{FM}(Q)`$ reduces to the mean-field result $`Y_c^{FM}=Q^{1/3}`$ in the limit $`Q\mathrm{}`$mittag ; pearce . And the second assumption is that in $`d`$ dimensions $`Y_c^{FM}(Q)`$ is the solution of an equation

$$Y^dc_1(d)Y^{d/2}c_2(d)=Q2,$$
(18)
where $`c_1(d)`$ and $`c_2(d)`$ have to be determined, for example, $`c_1=2`$ and $`c_2=1`$ for $`d=2`$. As shown in Table II the Hajduković conjecture is close to numerical estimates in the literature for $`2Q7`$. However, for $`Q=1`$ the Hajduković conjecture, $`Y_c^{FM}=0.7937\mathrm{}`$, completely disagrees with the numerical estimate $`0.753`$gaunt ; stauffer . We can guess a simpler form than the Hajduković conjecture for the critical temperature,

$$Y_c^{FM}(Q)=\frac{1}{1/\pi +Q^{1/3}},$$
(19)
satisfying the first assumption. Eq. (19) reasonably agrees with numerical estimates for $`1Q7`$ as shown in Table II. Even though the Hajduković conjecture and Eq. (19) are incorrect expressions for the FM critical temperature of the $`Q`$-state Potts model on the simple-cubic lattice, we can infer that the true expression for the critical temperature may be something between the Hajduković conjecture and Eq. (19).

Recently using cluster-variation calculations Rosengren and Lapinskasrosengren93  claimed the existence of two additional phase transitions at $`T_{RL1}=0.787(J/k_B)`$ and $`T_{RL2}=0.7715(J/k_B)`$ below the well-known antiferromagnetic (AFM) critical temperature $`T_c^{AFM}=1.2259(7)(J/k_B)`$wang  for the three-state AFM Potts model ($`J<0`$) on the simple-cubic lattice. Using the Monte Carlo method Kundrotas, Rosengren, and Lapinskaskundrotas  obtained $`T_{KRL}=0.68(J/k_B)`$ and supported their earlier claim. However, Kolesik and Suzukikolesik  and Heilmann et al.heilmann  observed different results by Monte Carlo simulations which are in disagreement with observations by Rosengren’s group, and obtained $`T_{KH}=1.0(J/k_B)`$ as the intermediate transition temperature in the region below $`T_c^{AFM}`$. Furthermore, based on a renormalization-group picture Oshikawaoshikawa  argued that there exists only one transition at $`T_c^{AFM}`$ for the three-state AFM Potts model on the simple-cubic lattice. Figure 2 shows the partition function zeros in the complex $`A=e^{\beta J}`$ plane of the three-state Potts model on the $`3\times 3\times 30`$ simple-cubic lattice with periodic boundary conditions in the $`x`$-$`y`$ plane and free boundary conditions in the $`z`$ direction. The interval $`0A1`$ ($`0T\mathrm{}`$) is antiferromagnetic ($`J<0`$) while the interval $`1A\mathrm{}`$ ($`\mathrm{}T0`$) ferromagnetic ($`J>0`$). There exists a clear one-dimensional distribution of zeros near $`A_{KRL}=0.2298`$, in Figure 2, which favors the existence of additional phase transitions below $`A_c^{AFM}=0.4424`$. However, we do not know if the first zero of the distribution near $`A_{KRL}`$ cuts the real axis in the thermodynamic limit, and if there exist the true phase transitions below $`A_c^{AFM}`$.

We have also obtained two pairs of the complex-temperature (CT) singularities (Figure 2) using the low-temperature seriesvohwinkel ; bhanot ; guttmann1  for the three-state Potts model on the simple-cubic lattice. Figure 2 exhibits accumulations of partition function zeros near the CT singularities. These accumulations imply that some thermodynamic quantities may diverge at the CT singularities which are not physical. The critical exponents associated with the CT singularities, which we call the CT critical exponents, are defined in the usual way,

$$C_{CT}^0\left(1\frac{Y}{Y_{CT}}\right)^{\alpha _{CT}},$$
(20)
$$m_{CT}^0\left(1\frac{Y}{Y_{CT}}\right)^{\beta _{CT}},$$
(21)
and

$$\chi _{CT}^0\left(1\frac{Y}{Y_{CT}}\right)^{\gamma _{CT}},$$
(22)
where $`Y_{CT}`$ is the location of a CT singularity, and $`C_{CT}^0`$, $`m_{CT}^0`$ and $`\chi _{CT}^0`$ are the singular parts of the specific heat, spontaneous magnetization, and susceptibility, respectively. Here, the superscript $`0`$ represents no external magnetic field. Because the density of zeros near a CT singularity is given byfisher65

$$g(Y)\left(1\frac{Y}{Y_{CT}}\right)^{1\alpha _{CT}},$$
(23)
the density of zeros diverges at the CT singularity if $`\alpha _{CT}>1`$.

Guttmann et al.guttmann2 ; guttmann3 ; guttmann4  studied the locations of the CT singularities and the CT critical exponents for the Ising model in detail, introducing a new method of low-temperature series analysis. For the Ising model on the simple cubic, body-centered cubic, and face-centered cubic lattices there exist the CT singularities, closer to the origin than the physical critical points, which control the low-temperature series and spoil the analysis of the real critical behavior. Guttmann et al.guttmann2  reported that for the Ising model on the simple-cubic lattice the CT singularity is located at $`(Y_{CT})^2=0.2860(2)`$ ($`|(Y_{CT})^2|<(Y_c^{FM})^2=0.4120`$) with $`\alpha _{CT}1`$, $`\frac{1}{8}\beta _{CT}\frac{1}{16}`$, and $`\gamma _{CT}\frac{9}{8}`$. Itzykson et al.itzykson  studied the zeros of the exact partition function in the complex temperature plane of the Ising model on the $`4\times 4\times 4`$ simple-cubic lattice, and found that the distribution of zeros supports the real existence of the CT singularities. By the analysis of the longer series they obtained $`\alpha _{CT}=1.25(15)`$, $`\beta _{CT}=0.05`$, and $`\gamma _{CT}=1.1(1)`$ for the Ising model on the simple-cubic lattice, and argued that the values of the CT critical exponents of the Ising model are same for the simple cubic, body-centered cubic, and face-centered cubic lattices and satisfy $`\alpha _{CT}+2\beta _{CT}+\gamma _{CT}2`$. Finally, Itzykson et al. conjectured that the behavior at the CT singularities is universal. Recently Guttmann and Entingguttmann5  obtained $`(Y_{CT})^2=0.2853(3)`$ with $`\alpha _{CT}1.03`$, $`\beta _{CT}0.01`$, and $`\gamma _{CT}1.01`$ using the low-temperature series generated by the finite-lattice method of series expansion for the simple-cubic Ising model. They proposed $`\alpha _{CT}=1`$, $`\beta _{CT}=0`$, and $`\gamma _{CT}=1`$ with logarithmic corrections, satisfying the Rushbrooke scaling law $`\alpha _{CT}+2\beta _{CT}+\gamma _{CT}=2`$.

Using Dlog Padé approximantsguttmann6  to the low-temperature seriesguttmann1  of the specific heat, spontaneous magnetization, and susceptibility for the three-state Potts model on the simple-cubic lattice, we found the locations of two pairs of the CT singularities

$$Y_{CT}^{(1)}=0.0254(11)\pm 0.4948(9)i$$
(24)
and

$$Y_{CT}^{(2)}=0.4981(45)\pm 0.2453(35)i.$$
(25)
Since $`|Y_{CT}^{(1)}|=0.4955(9)`$ and $`|Y_{CT}^{(2)}|=0.5552(42)`$, these CT singularities lie closer to the origin than the estimated values of the FM critical point shown in Table II. The CT critical exponents for the simple-cubic three-state Potts model are same at $`Y_{CT}^{(1)}`$ and $`Y_{CT}^{(2)}`$ as shown in Table III, and satisfy the scaling law $`\alpha _{CT}+2\beta _{CT}+\gamma _{CT}2`$. Because in three dimensions the values of the CT critical exponents for the three-state Potts model are very close to those for the Ising model, it is difficult not to suggest that the values of the CT critical exponents are independent of $`Q`$ in three dimensions. That is, our results support the universality of the CT singularities in three dimensionsitzykson . For the two-dimensional $`Q`$-state Potts model on the square lattice it has been shownmatveev  that two pairs of the CT singularities of the model are mapped each other through the dual transformation

$$Y\frac{1Y}{1+(Q1)Y},$$
(26)
which determines the FM critical temperature

$$Y_c^{FM}(Q)=\frac{1}{1+\sqrt{Q}}.$$
(27)
Therefore, one expects the true expression for the FM critical points of the simple-cubic Potts model to provide the mapping between $`Y_{CT}^{(1)}`$ and $`Y_{CT}^{(2)}`$.

## IV partition function zeros for an external magnetic field

In this section we discuss the properties of the simple-cubic $`Q`$-state Potts models for nonzero magnetic field using the high-field, low-temperature series expansionssykes65 ; sykes73 ; bhanot92 ; enting ; straley  and the zeros of the exact partition function $`Z_Q(Y,X)`$ on a finite lattice in the complex temperature ($`Y=e^{\beta J}`$) or magnetic-field ($`X=e^{\beta H_q}`$) plane. Figure 3 shows the partition function zeros in the complex $`X`$ plane of the three-state Potts model on the $`3\times 3\times 3`$ simple-cubic lattice with periodic boundary conditions in the $`x`$-$`y`$ plane and free boundary conditions in the $`z`$ direction. The zeros lie close to the unit circle at lower temperatures, but they move away from the positive real axis and the unit circle as the temperature is increased. The zeros at $`T=0`$ are uniformly distributed on the circle with radius $`(Q1)^{1/N_s}`$ which approaches unity in the thermodynamic limit, independent of $`Q`$, whereas the zeros at $`T=\mathrm{}`$ are $`N_s`$-degenerate at $`X=1Q`$, independent of lattice sizekim98a ; kim99 ; kim00a .

The edge zero, which we call the first zero, and its complex conjugate of a circular distribution of zeros in the complex $`X`$ plane cut the positive real axis at the physical critical point $`X_c=1`$ for $`TT_c`$ in the thermodynamic limit. However, for $`T>T_c`$ the first zero does not cut the positive real axis in the thermodynamic limit, that is, there is a gap in the distribution of zeros around the positive real axis. Within this gap, the free energy is analytic and there is no phase transition. Kortman and Griffithskortman  carried out the first systematic investigation of the magnetization at the first zero, based on the high-field, high-temperature series expansion for the Ising model on the square lattice and the diamond lattice. They found that above $`T_c`$ the magnetization at the first zero diverges for the square lattice and is singular for the diamond lattice. For $`T>T_c`$ we rename the first zero as the Yang-Lee edge singularity. The divergence of the magnetization at the Yang-Lee edge singularity means the divergence of the density of zeros, which does not occur at a physical critical point. Fisherfisher78  proposed the idea that the Yang-Lee edge singularity can be thought of as a new second-order phase transition with associated critical exponents and the Yang-Lee edge singularity can be considered as a conventional critical point. The critical point of the Yang-Lee edge singularity is associated with a $`\varphi ^3`$ theory, different from the usual critical point associated with the $`\varphi ^4`$ theory. The crossover dimension of the Yang-Lee edge singularity is $`d_c=6`$. The edge critical exponent $`\sigma _e=1/\delta _e=(d2+\eta _e)/(d+2\eta _e)`$ at a Yang-Lee edge singularity $`X_e`$ for $`Y>Y_c`$ is defined by

$$m_e^X\left(1\frac{X}{X_e}\right)^{\sigma _e},$$
(28)
where $`m_e`$ is the singular part of the magnetization in the complex $`X`$ plane. The study of the Yang-Lee edge singularity has been extended to the classical $`n`$-vector modelkurtze1 , the quantum Heisenberg modelkurtze1 , the spherical modelkurtze2 , the quantum one-dimensional transverse Ising modeluzelac , the hierarchical modelbaker , the one-dimensional Potts modelmittag84 ; glumac , branched polymersparisi , fluid models with repulsive-core interactionslai ; ypark , etc. Dhardhar  calculated the edge critical exponent $`\sigma _e`$ in two dimensions by solving a particular model of three-dimensional directed animals and mapping the solution to the hard hexagon model. Using Fisher’s idea and conformal field theory, Cardy cardy  studied the Yang-Lee edge singularity for a two-dimensional $`\varphi ^3`$ theory. It is generally accepted that the value of the edge critical exponent $`\sigma _e`$ depends only on dimension. The values of $`\sigma _e`$ are known to be $`\sigma _e=\frac{1}{2}`$ in one dimension and $`\sigma _e=\frac{1}{6}`$dhar ; cardy  in two dimensions, whereas $`\sigma _e0.08`$lai  in three dimensions.

The partition function zeros in the complex $`X`$ plane for real $`Y`$ have been extensively studied and well understood. However, the partition function zeros in the complex $`Y`$ plane for real $`X`$ are much less well understood than the zeros in the complex $`X`$ plane. Figure 4 shows the partition function zeros in the complex $`Y`$ plane of the three-state Potts model for an external magnetic field $`X1`$ ($`H_q0`$) on the $`3\times 3\times 3`$ simple-cubic lattice with periodic boundary conditions in the $`x`$-$`y`$ plane and free boundary conditions in the $`z`$ direction. For $`X=1`$ some zeros distribute along the one-dimensional locus near $`Y_c`$, and the first zero of the distribution cuts the positive real axis at $`Y_c`$ in the thermodynamic limit as explained in the previous section. As $`H_q`$ increases, all the zeros move away from the origin. In the limit $`H_q\mathrm{}`$ ($`X\mathrm{}`$) the field $`H_q`$ favors the state $`q`$ for every site and the $`Q`$-state Potts model is transformed into the one-state model whose the partition function zeros are the roots of the equation $`Y^{N_b}=0`$kim98b . That is, $`|Y|`$ for all the zeros increases without bound as $`X`$ increases. Note that for $`X>1`$ there is accumulation of the zeros as we approach the first zero. Matveev and Shrockmatveev96a  studied the properties of the thermodynamic functions at the first zero in the complex temperature plane for the square-lattice Ising model with $`X>1`$. They observed that the thermodynamic functions and the density of zeros diverge at the first zero when $`X>1`$, and that the values of the critical exponents associated with the first zero for $`X>1`$ are nearly independent of $`X`$ and satisfy the Rushbrooke scaling law approximately. We call the first zero in the complex temperature plane for $`X>1`$ as the Fisher edge singularity. The edge critical exponents $`\alpha _e`$, $`\beta _e`$, and $`\gamma _e`$ at the Fisher edge singularities are given by

$$C_e\left(1\frac{Y}{Y_e}\right)^{\alpha _e},$$
(29)
$$m_e^Y\left(1\frac{Y}{Y_e}\right)^{\beta _e},$$
(30)
and

$$\chi _e\left(1\frac{Y}{Y_e}\right)^{\gamma _e},$$
(31)
where $`Y_e`$ is the location of a Fisher edge singularity, and $`C_e`$, $`m_e^Y`$ and $`\chi _e`$ are the singular parts of the specific heat, magnetization, and susceptibility, respectively, for $`X>1`$. Kim and Creswickkim98b  extended the study of the Fisher edge singularity to the $`Q`$-state Potts model ($`Q3`$) on the square lattice. They found that the values of the edge critical exponents at the Fisher edge singularities are nearly independent of $`Q`$ ($`Q2`$).

In three dimensions to study the critical behavior at the Fisher edge singularity we have used the high-field, low-temperature series expansions for the Ising modelsykes65 ; sykes73 ; bhanot92  and the three-state Potts modelenting ; straley  on the simple-cubic lattice. Table IV shows the locations ($`Y_e`$) and the edge critical exponents ($`\alpha _e`$, $`\beta _e`$, and $`\gamma _e`$) of the Fisher edge singularities estimated by Dlog Padé approximants to the specific heat, magnetization, and susceptibility series for $`x=100`$, 200, and 500. The edge critical exponents, which are independent of $`X`$, have the values of $`\alpha _e=1.01.1`$, $`\beta _e0.4`$, and $`\gamma _e1.0`$ for $`Q=2`$ and for $`Q=3`$ similar values of $`\alpha _e`$ and $`\gamma _e`$ and $`\beta _e0.6`$. So in three dimensions we seem to have a strong violation of the scaling law $`\alpha _e+2\beta _e+\gamma _e=2`$. Another interesting observation is that there exist apparently three pairs of the Fisher edge singularities in three dimensions (Figure 4 and Table IV). The estimates of the edge critical exponents at each of the three edge singularities indicate that they are equal within error bars.

On the other hand, everything can be evaluated exactly in one dimension, and there exist only one pair of the Fisher edge singularities characterized by the edge critical exponents $`\alpha _e=\frac{3}{2}`$, $`\beta _e=\frac{1}{2}`$, and $`\gamma _e=\frac{3}{2}`$, all independent of $`Q`$. These values satisfy the scaling law $`\alpha _e+2\beta _e+\gamma _e=2`$. Table V shows the locations and the edge critical exponents of the Fisher edge singularities for the $`Q`$-state Potts models on the square lattice. In two dimensions we have also found two pairs of the Fisher edge singularities which have the same values of the edge critical exponents within error bars. The edge critical exponents have the values of $`\alpha _e=1.131.24`$, $`\beta _e=0.110.20`$, and $`\gamma _e=1.111.26`$. In two dimensions the edge critical exponents also satisfy the scaling law $`\alpha _e+2\beta _e+\gamma _e=2`$. Now we can conclude that the number of pairs of the Fisher edge singularities is equal to the dimensionality.

Until now the Yang-Lee and Fisher edge singularities have been studied and understood independently. However, the Yang-Lee edge singularity seems to have a direct relationship to the Fisher edge singularity. From the edge critical exponents of the Fisher edge singularity in one dimension we get $`y_t^e=2`$ and $`y_h^e=2`$, which also follows from $`\sigma _e=\frac{1}{2}`$ at the Yang-Lee edge singularity. The equivalence of the Yang-Lee and Fisher edge singularities is not surprising in light of the transformation which maps the temperature and field variables into each other in one dimensionglumac . In two dimensions it is known that $`\sigma _e=\frac{1}{6}`$ from which it follows that $`y_h^e=\frac{12}{5}`$. If we assume $`y_t^e=y_h^e=\frac{12}{5}`$ for the Fisher edge singularity, we have $`\alpha _e=\frac{7}{6}`$, $`\beta _e=\frac{1}{6}`$, and $`\gamma _e=\frac{7}{6}`$, which are not far from the values of the exponents estimated by series expansions. In three dimensions if we take $`\alpha _e=1.0`$ then $`y_t^e=3`$, and $`\sigma _e0`$ (logarithmic singularity), which seems to be the historical trend, then $`y_h^e=3`$ as well. So in $`d=1`$, 2 and 3 we seem to find $`y_t^ey_h^e`$. If $`y_t^e=y_h^e`$, the Fisher and Yang-Lee edge singularities are not independent behaviors, and they seem to have only one relevant variable.

## V conclusion

We have investigated the interesting properties of the $`Q`$-state Potts model on the simple-cubic lattice using the zeros of the exact partition function on a finite lattice and series expansions. The critical behavior of the Potts model in the ferromagnetic and antiferromagnetic phases has been studied at the same time based on the distribution of the partition function zeros in the complex temperature plane. The distribution of the zeros in the complex temperature plane reveals the possibility of the intermediate phase transitions below the order-disorder phase transition of the Potts antiferromagnet. The behavior at the CT singularities, where some thermodynamic quantities diverge with the characteristic CT exponents $`\alpha _{CT}`$, $`\beta _{CT}`$, and $`\gamma _{CT}`$, has been discussed using the partition function zeros and the low-temperature series expansion of the Potts model in the absence of an external magnetic field. In the complex temperature plane the CT singularities coexist with the physical critical points for $`H_q=0`$.

We have also studied the veiled properties of the Potts model for nonzero magnetic field using the partition function zeros and the high-field, low-temperature series expansion. As $`H_q`$ increases, all the zeros in the complex $`Y=e^{\beta J}`$ plane move away from the origin. For $`H_q>0`$ the physical critical points $`Y_c^{FM}`$ and $`Y_c^{AFM}`$ disappear and the Fisher edge singularities $`Y_e`$ appear in the complex temperature plane. For different values of $`X=e^{\beta H_q}`$ we have calculated the locations and the edge critical exponents of the Fisher edge singularities of the Ising and three-state Potts models on the simple-cubic lattice. The values of the edge critical exponents $`\alpha _e`$, $`\beta _e`$, and $`\gamma _e`$ at a Fisher edge singularity are universal, that is, independent of $`X`$ and $`Q`$ but dependent on dimension $`d`$. We found that the number of pairs of the Fisher edge singularities is equal to $`d`$. We have considered the Yang-Lee edge singularity $`X_e`$, which appears in the complex $`X`$ plane for $`Y>Y_c`$ instead of the physical critical point $`X_c=1`$, and its edge exponent $`\sigma _e=1/\delta _e`$. In one, two, and three dimensions the Yang-Lee edge singularity is related to the Fisher edge singularity. They seem to have only one relevant variable $`y_t^ey_h^e`$.

ACKNOWLEDGMENTS

The author is indebted to Prof. R. J. Creswick for very helpful discussions. The author is grateful to Prof. A. J. Guttmann for kindly providing him with the missed coefficients in the paperguttmann1 . The author also wishes to thank Prof. M. E. Fisher for valuable comments and for kindly drawing his attention to the paperypark .