text:
# Protein Unfolding and the Diffusion Collision Model

## I. Introduction

In the diffusion-collision model (DCM) of protein folding (Bashford et al 1988, Karplus and Weaver, 1976, 1979, 1994) the protein is modeled using a collection of spheres connected by flexible strings. The spheres represent the secondary structural elements such as $`\alpha `$-helices or $`\beta `$-sheets, called microdomains, that make up the protein.

The folding process from a completely unfolded protein to the the final native state is accomplished via diffusion through the solvent, collision, and finally coalescence of the microdomains. The state of the protein is defined by the number of pairings (hydrophobic interactions) between the microdomains that are present at a given time $`t`$. The rate equations for transitions between these states can be written as

$$\frac{d𝐏(t)}{dt}=\widehat{K}𝐏(t)$$
(1)
where $`𝐏(t)`$ is the vector of states and $`\widehat{K}`$ is a matrix containing the transition rates between the different states.

As an example, consider a simple monomeric protein containing two $`\alpha `$-helices A and B joined by a connecting string. This gives us a one-pair/two-state system. We will call state $`1`$ the state when the two microdomains are not in hydrophobic contact, and state $`2`$ when they are hydrophobically docked. In this case (1) can easily be written out explicitly

$$\begin{array}{c}\dot{P}_1(t)=k_{12}P_1(t)+k_{21}P_2(t)\\ \dot{P}_2(t)=k_{12}P_1(t)k_{21}P_2(t)\end{array}$$
(2)
or in matrix form as

$$\frac{d}{dt}\left(\begin{array}{c}P_1(t)\\ P_2(t)\end{array}\right)=\left[\begin{array}{cc}k_{12}& k_{21}\\ k_{12}& k_{21}\end{array}\right]\left(\begin{array}{c}P_1(t)\\ P_2(t)\end{array}\right).$$
(3)
A more complicated protein having, say, $`n`$ microdomains would involve $`p=n(n1)/2`$ pairings, $`2^p`$ states $`P_i(t)`$ and a $`2^p\times 2^p`$ rate matrix $`\widehat{K}`$. In general the calculation of the elements of the rate matrix $`\widehat{K}`$ is somewhat involved. The forward rates are the rates of microdomain coalescence. In the diffusion-collision model the forward rates are calculated assuming the microdomains diffuse through a solvent environment, the space of which is limited by the length of the intervening strings and the van der Waals radii of the microdomains. These microdomains are assumed to be nascently formed, and their degree of formation is given by a helix-coil transition theory (Zimm and Bragg, 1959, Lifson and Roig, 1961) as in AGADIR (Lacroix et al, 1998, Munoz and Serrano, 1994 a,b,c, 1997) calculation in the case of $`\alpha `$-helices, or via a combination of theory (Munoz et al, 1998) and experiment (Dinner et al, 1999) in the case of $`\beta `$-sheets. As the microdomains undergo diffusion, they occasionally collide. When this happens the microdomains coalesce with a probability $`\gamma `$, being held together by hydrophobic interactions in the case of $`\alpha `$-helices, or a combination of hydrophobic and hydrogen bond interactions in the case of $`\beta `$-sheets. The coalescence probability $`\gamma `$ is given by the likelihood that the microdomain is in $`\alpha `$-helical or $`\beta `$-sheet form, the percentage of hydrophobic area, and the likelihood of proper geometrical orientation upon collision.

The forward folding times in the mean first passage time approximation (Weiss, 1967, Weaver, 1979, Szabo et al 1980) are given by

$$\tau _f=\frac{l^2}{D}+\frac{LV(1\gamma )}{\gamma DA}$$
(4)
where *$`V`$* is the diffusion volume available to the microdomain pair, *$`A`$* is the target area for collisions, *$`D`$* is the relative diffusion coefficient, $`\gamma `$ is the probability of coalescence upon collision and *$`l`$* and *$`L`$* are geometrical parameters calculated for diffusion in a spherical space. The inverse of the first passage time-scales $`\tau _f`$ are the forward folding rates *$`k_f`$* that are used in the rate matrix $`\widehat{K}`$*.* In the example given by (2) and (3), $`k_f`$ is $`k_{12}`$.

The microdomain pairings can also dissociate. To date, in typical diffusion-collision model calculations, the form of the backward folding, or unfolding times $`\tau _b`$ used for two microdomains A and B comes from from the Van’t Hoff-Arrhenius law (Van’t Hoff, 1884, Arrhenius, 1889) given by

$$\tau _b=\nu ^1e^{fA_{AB}/k_BT}$$
(5)
where *$`f`$* is the free energy change per unit buried hydrophobic area in the pairing (Chothia, 1974), *$`A_{AB}`$* is the buried area (Lee and Richards, 1971), *$`k_B`$* is Boltzmann’s constant, *$`T`$* is the temperature and $`\nu `$ is an attempt rate.

In diffusion collision model calculations to date the attempt rate is a parameter whose value must be set by hand, usually requiring some guess work or adjustment to obtain the desired result (see for example Burton et. al 1998). The typical values used lie in the range $`1\times 10^9s^11000\times 10^9s^1`$, obtained from estimates of covalent bond oscillation frequencies (Fersht, 1999). This procedure is not explicitly stated in most diffusion-collision model studies; it is justified, however, when the equilibrium occupation probabilities of the states are known. In fact, the rate constants can be calculated from such probabilities (Chandler, 1987). In cases where the final occupation probabilities are unknown, for instance in studies of protein mis-folding and non-native kinetic intermediates (Beck et al., 2000), such methods are clearly not possible.

In this work we show how the parameter $`\nu `$, or more generally the unfolding rates, can be determined from thermal fluctuations providing a means of avoiding the guesswork. This makes the diffusion-collision model more predictive and enables it to be used in situations where the final occupation probabilities are unknown.

## II. Calculation of the Unfolding Rates

In order to calculate the backward rate it is convenient to view the unfolding process as that of diffusion within, and escape from, an effective one-dimensional potential well. This is a good approximation if only motion perpendicular to the hydrophobic contact surface is important in microdomain pair dissociation.

Consider the pair of microdomains A and B connected by a string. Microdomain coalescence and dissociation can be approximated by diffusion in a potential like the one depicted in Figure 1. The quantity that diffuses is the separation between the microdomains and we can think of

$$dP(x)=\rho (x)dx$$
(6)
as the fraction of microdomains with a distance between them in the range $`x`$ and $`x+dx`$. In a potential $`V(x)`$ the probability density $`\rho (x)`$ satisfies the one dimensional Smoluchowski equation

$$\frac{\rho (x,t)}{t}=\frac{}{x}\left[D\left(\frac{\rho (x,t)}{x}+\rho (x,t)\frac{V(x)}{x}\right)\right].$$
(7)
The outer boundary on the far right is given by an infinite potential, arising from the two microdomains coming to the end of their connecting string. The inner boundary is also infinite, stemming from the hard-core repulsion of the van der Waals contact between the microdomains. The well depth $`E_b=fA_{AB}`$ is the free energy difference between paired and unpaired microdomains and the well width $`L`$ is set to the diameter of a water molecule. A separation larger than $`L`$ exposes the microdomain to the solvent, the free energy savings is lost, and the microdomains separate, resulting in an escape from the potential well.

The forward folding rates can be satisfactorily found by the mean first passage time using the method outlined above (4). As a consequence of this we shall forego direct analysis of the entire diffusion process and concentrate on the escape of a bound microdomain out of the potential well in order to find the backward folding rate.

If we concentrate on the part of the potential where $`x<L`$ in Figure 1, we can think of the potential as being perfectly reflecting on one side and partially permeable on the other. The permeability arises from the thermal fluctuations in the energy which allow the microdomains to escape the bounded region, provided they have sufficient energy, and it can be expressed using an appropriate choice of boundary conditions. We proceed as follows.

For the inner boundary, representing van der Waals contact, we impose perfectly reflecting boundary conditions, meaning the flux through the boundary is zero

$$J(0,t)=D\frac{d\rho (x,t)}{dx}|_{x=0}=0.$$
(8)
The flux at the permeable boundary at $`x=L`$ depends on the density of microdomains at that boundary and the probability that the microdomain energy is high enough to thermally tunnel through the boundary. The differential element of flux at the boundary $`L`$ of microdomains with relative velocity between *$`v`$* and *$`v+dv`$* is given by

$$dJ(L,t)=v\rho (L,t)dN(v)$$
(9)
where

$$dN(v)=\left(\frac{\mu }{2\pi k_BT}\right)^{1/2}e^{\mu v^2/2k_BT}dv$$
(10)
is the Maxwell-Boltzmann distribution, namely the fraction of microdomains with relative velocities between $`v`$ and $`v+dv`$ and $`\mu `$ is the reduced mass given by

$$\mu =\frac{m_Am_B}{(m_A+m_B)}$$
(11)
where $`m_A`$ and $`m_B`$ are the masses of the two microdomains. The distribution is normalised so that the bounds of the integral go from $`\mathrm{}`$ to $`\mathrm{}`$.

It is clear that in order to find the total flux through the outer boundary at *$`L`$* we must integrate over all microdomain velocities larger than $`+\sqrt{E_b/2m}`$ since the potential well height is $`E_b`$, and only microdomains with velocity higher than that can escape and thus contribute to the flux leaving the well. This yields a flux out of state $`2`$

$$J_2^{out}(L,t)=\rho (L,t)\left(\frac{k_bT}{2\pi \mu }\right)^{1/2}e^{E_b/k_BT}$$
(12)
If we assume that the probability is uniformly distributed across the well, then $`\rho (x,t)=P_2(t)/L`$ everywhere within the well. This is a reasonable assumption because the permeability of the well is so low. Furthermore it can be shown that typically the damping time (or velocity autocorrelation) time is considerably smaller than the well crossing time and the mean first passage time. Clearly then, diffusion in the interior of the well distributes the microdomain separations evenly among the available space and therefore

$$J_2^{out}(L,t)=\frac{P_2(t)}{L}\left(\frac{k_bT}{2\pi \mu }\right)^{1/2}e^{E_b/k_BT}.$$
(13)
The interesting result is that by using a stationary flux method (Kramers, 1940), or equivalently setting $`P_2(t)=1,`$ the outward flux at the boundary in one dimension $`J^{out}(L,t)`$ is the fraction of microdomain pairs crossing the boundary at $`L`$ per unit time and thus clearly identical to the unfolding rate $`k_b`$ for the $`AB`$ microdomain pairing. Below we explain this in more detail and show more generally how the rate equations of the diffusion-collision model (1, 2) can be derived from the diffusion equation.

To understand this, we remember that the Smoluchowski equation (7) is simply a detailed statement of the equation of continuity

$$\frac{\rho (𝐱,t)}{t}=𝐉(𝐱,t)$$
(14)
The quantity in brackets in (7) is the flux $`𝐉(𝐱,t)`$, which includes the ordinary $`\rho (𝐱,t)`$ term, and the second term which takes into account external forces. Integration of both sides of (14) over the bound part of diffusion volume, using the divergence theorem on the R.H.S. and assuming isotropic flow into and out of the boundary of the diffusion volume yields

$$\frac{dP_2(t)}{dt}=𝐧𝐉(t)A=\{J_2^{in}(t)J_2^{out}(t)\}A.$$
(15)
The LHS of this equation is the rate of change of probability in the bound region. On the RHS $`𝐧`$ is a unit vector normal to the boundary, $`𝐧𝐉(t)`$ is the net flux which traverses the boundary surface and $`A`$ is the area of that boundary surface. We suggestively write the net flux $`𝐧𝐉(t)`$ in terms of $`J_2^{in}(t)`$ and $`J_2^{out}(t)`$, the probability fluxes in and out of the diffusion volume in question.

We see from comparison to (2) that (15) can be written out explicitly for our case

$$\{J_2^{in}(t)J_2^{out}(t)\}A=k_{12}P_1(t)k_{21}P_2(t).$$
(16)
One can easily see that the positive terms on either side of (16) are equal to each other as are the negative terms. Referring to (12), this quantity is

$$J_2^{out}(t)A=k_{21}P_2(t)$$
(17)
making the rate

$$k_{21}=\frac{J_2^{out}(t)A}{P_2(t)}.$$
(18)
In this expression the one-dimensional “area” through which the probability is flowing is simply a point and therefore $`A=1`$ so identifying $`J_2^{out}(t)`$ with (13) we can solve for the unfolding rates in $`\widehat{K}`$, $`k_{21}`$ in (2) and (3), by dividing (13) by $`P_2(t)`$, or setting $`P_2(t)=1`$ in (18) according to the stationary flux method. We find the backward folding rate in the one-dimensional case to be

$$k_{21}=\frac{1}{L}\left(\frac{k_BT}{2\pi \mu }\right)^{1/2}e^{E_b/k_BT}.$$
(19)
The terms preceding the exponential correspond to our prediction for the Van’t Hoff-Arrhenius attempt rate $`\nu `$ in (5). As an example, the attempt rate found for a coalesced pair of 16-residue Regan-Degrado (Regan and Degrado 1988) helices with a combined hydrophobic area loss of $`600`$Å2 is $`64\times 10^9s^1`$.

We believe that the most probable dissociation of a microdomain-microdomain pairing occurs via relative motion along a vector connecting the centers-of-mass of the two microdomains. This implies that typically the one dimensional case is the best way to view the dissociation event and therefore that one should use (19) to go about calculating the rate. The one dimensional case should also be sufficient if there is an “unzipping” of paired $`\alpha `$-helices.

It is possible, however, that dissociation may also include a relative rolling motion or other motion perpendicular to the axis connecting the microdomain pairs. In this case one needs to repeat the calculation above with a few minor differences: The relative velocity distribution of the microdomains is still the one dimensional Maxwell-Boltzmann distribution because the degrees of freedom parallel to the surface through which the probability is flowing do not contribute to escape from the well, the probability in the bound region is assumed to be evenly distributed in a two-dimensional “volume”, namely $`\rho (x,t)=P_2(t)/\pi L^2`$ in the two-dimensional analogue of (12), and the flux goes through a two-dimensional “area” $`A=2\pi L`$ in (18). This calculation yields the result

$$k_{21}=\frac{2}{L}\left(\frac{k_bT}{2\pi \mu }\right)^{1/2}e^{E_b/k_BT}.$$
(20)
Due to the steric clashing of the side chains it seems rather unlikely that dissociation would include a relative sliding motion along the axes of the microdomains. For completeness, however, we include the three-dimensional result derivable from analogous considerations to the ones given above

$$k_{21}=\frac{3}{L}\left(\frac{k_bT}{2\pi \mu }\right)^{1/2}e^{E_b/k_BT}.$$
(21)
Although the three-dimensional result seems an unlikely candidate for protein unfolding it may be relevant in the context of molten globules.

This approach succeeds in removing the free parameter $`\nu `$ from the model, and allows us to find the backward rates from a simple energetic model based on diffusion in a potential with appropriate boundary conditions. Moreover our results for the one, two and three-dimensional unfolding rates have a $`\sqrt{T/\mu }`$ dependence that could be used to distinguish between this and other proposals for the mechanism of helix-helix dissociation.

The removal of the parameter $`\nu `$ is important when considering folding processes which do not involve the native state. In previous applications of the diffusion-collision model, the folding kinetics from a denatured or random coil state to the final native state were followed. In such a case, it is reasonable to set the parameter $`\nu `$ such that the native state achieves 90 or 95 percent of the probability, because we know that the final state is attained at the end of the folding process. In studying intermediate processes or more importantly, non-native intermediates (Beck et al, 2000), where the occupation probability may be completely unknown, such reasonable estimates of $`\nu `$ are not available. In this case, elimination of $`\nu `$ as a free parameter is crucial.

## III. Concluding Remarks

We have presented a calculation for the helix-helix dissociation rate using a simplified potential surface, a square potential having a depth equal to the free energy savings of hydrophobic docking and width equal to the diameter of a water molecule. We have found the unfolding rates arising from thermal fluctuations out of this potential well to be in good agreement with currently accepted values of $`\nu `$. The potential itself is due to hydrophobic forces, which to date are not well understood and for this reason the potential chosen was simple.

The initial motivation of this work was to eliminate the free parameter $`\nu `$ from the diffusion-collision model. In the context of this model our result should allow us to perform kinetics simulations than were not possible before, for instance time developments of non-native intermediates, in which the occupation probabilities are unknown, and reasonable estimates for $`\nu `$ are not possible.

The results presented here also predict a $`\nu \sqrt{T/\mu }`$ dependence in all cases which can be distinguished experimentally from other proposals such as the covalent bond model where $`\nu k_BT/\mathrm{}`$ (Fersht, 1999). Another difference is the dependence of the unfolding rates on the states, not only through the hydrophobic area, but also through the reduced mass $`\mu `$ of the microdomains or groups of microdomains undergoing dissociation. This is markedly different from typical diffusion-collision model calculations where the attempt rate $`\nu `$ is assumed to be the same for all dissociation events within the protein.

Along the way we have re-examined the idea that couching the problem of association and dissociation of microdomain pairs via diffusion over a potential surface affords us a clear and simple picture of the protein folding process. Indeed, the rate equations of the diffusion-collision model can be derived from such a picture. Interestingly, the potential energy surface is an element of other models such as folding funnels (for an overview see Bryngelson et al 1995).

Admittedly, our approximation is rather rough. The thrust of future work should therefore be in constructing more realistic potentials including activation energies of hydrophobic dockings.

Generalization to more complex proteins is straightforward because every interaction between helices or clusters of helices can be considered a two-state process similar to the one described above; generally, the folding or unfolding of a given protein involves several such processes. In this case the diffusion space for each possible pairing could be constructed, and the forward and backward rates for each transition can be found, as outlined above, to construct the transition rate matrix $`\widehat{K}`$ in order to find the time evolution of the state vector $`𝐏(t)`$. This is a much more tenable proposition than directly solving the Smoluchowski equation (7) on such a complicated potential surface, although it could be done in principle.

## Acknowledgments

We would like to thank Ken Olum, David Weaver, and Larry Ford for many useful discussions. We are further grateful to David Weaver for pointing out a mistake in our calculations and Ken Olum for helping us fix it.

## References

Arrhenius S. 1889. *Z Phys Chem 4*:226.

Bashford D, Cohen FE, Karplus M, Kuntz ID, Weaver DL 1988. Diffusion-collision model for the folding kinetics of myoglobin. *Proteins 4*:211-227.

Beck C, Siemens X, Weaver DL. 2000. Non-Native Kinetic intermediates and the Diffusion-Collision model of Protein Folding. *In preparation*.

Burton RE, Meyers JK, Oas TG. 1998. Protein folding dynamics: quantitative comparison between theory and experiment. *Biochemistry 37*:5337-5343.

Bryngelson JN, Onuchic ND, Socci ND, Wolynes PG. 1995. Funnels, pathways, and the energy landscape of protein folding: a synthesis. *Proteins: Struct Func and Gen 21*:need pages

Chandler D. 1987. Introduction to Modern Statistical Mechanics. Oxford University Press.

Chothia C. 1974. Hydrophobic bonding and accessible surface area in proteins. *Nature 248*:338-339.

Dinner AR, Lazaridis T, Karplus M. 1999. Understanding $`\beta `$-hairpin formation. *Proc Nat Acad Sci USA 96*:9068-9073.

Fersht A. Structure and mechanism in protein science. 158-159. W.H. Freeman and Co. New York. 1999.

Hanggi P, Talkner P, Borkovec M. 1990. Reaction-rate theory \- 50 years after Kramers. *Rev Mod Phys 62*:251-341.

Karplus M, Weaver DL. 1976. Protein folding dynamics. *Nature 260*:404-406.

Karplus M, Weaver DL. 1979. Diffusion-collision model for protein folding. *Biopolymers 18*:1421-1437.

Lacroix E, Viguera AR, Serrano L. 1998. Elucidating the folding problem of $`\alpha `$-helices: Local motifs, long-range electrostatics, ionic strength dependence, and prediction of of NMR parameters. *J Mol Biol 284*:173-191.

Lee B, Richards FM. 1971. The interpretation of protein structure: Estimation of static accessibility. *J Mol Biol 55*:379-400.

Munoz V, Henry ER, Hofrichter J, Eaton WA. 1998. A statistical mechanical model for $`\beta `$-hairpin kinetics. *Proc Nat Acad Sci USA 95*:5872-5879.

Munoz V, Serrano L. 1994a. Elucidating the folding problem of helical peptides using empirical parameters. *Nature: Struct Biol 1*:399-409.

Munoz V, Serrano L. 1994b. Elucidating the folding problem of helical peptides using empirical parameters II: Helix macrodipole effects and rational modification of the helical content of natural peptides. *J Mol Biol 245*:275-296.

Munoz V, Serrano L. 1994c. Elucidating the folding problem of helical peptides using empirical parameters III: Temperature and pH dependence. *J Mol Biol 245*:297-308.

Munoz V, Serrano L. 1997. Development of the multiple sequence approximation within the AGADIR model of $`\alpha `$-helix formation. Comparison with Zimm-Bragg and Lifson-Roig formalisms. *Biopolymers 41*:495-509.

Nagi AD, Anderson KS, Regan L 1999. Using loop length variant to dissect the folding pathway of a four-helix-bundle protein. *J Mol Biol 286*:257-265.

Pappu RV, Weaver DL. 1997. The early folding kinetics of apomyoglobin. *Protein Science 7*:480-490.

Regan L, Degrado WF. 1988. Characterization of a helical protein designed from first princliples. *Science 241*:976-978.

Szabo A, Schulten K, Schulten Z. 1980. First passage time approach to diffusion controlled reactions. *J Chem Phys 72*:4350.

Van’t Hoff JH. 1884. Studies in chemical dynamics. London, 1896.

Weaver DL. 1979. Some exact results for one-dimensional diffusion with absorption. *Phys Rev B 20*:2558.

Weiss GH. 1967. First passage time problems is chemical physics. *Adv Chem Phys 13*:1-18.

Zimm BH, Bragg JK. 1959. Theory of phase transition between helix and random coil in polypeptide chains. *J Chem Phys 31*:526-535.


format: text

text:
# Tree structure of the percolating Universe

## Abstract

We present a numerical study of topological descriptors of initially Gaussian and scale-free density perturbations evolving via gravitational instability in an expanding universe. We carefully evaluate and avoid numerical contamination in making accurate measurements on simulated fields on a grid in a finite box. Independent of extent of non linearity, the measured Euler number of the excursion set at the percolation threshold, $`\delta _\mathrm{c}`$, is positive and nearly equal to the number of isolated components, suggesting that these structures are trees. Our study of critical point counts reconciles the clumpy appearance of the density field at $`\delta _\mathrm{c}`$ with measured filamentary local curvature. In the Gaussian limit, we measure $`|\delta _\mathrm{c}|>\sigma `$ in contrast to widely held belief that $`|\delta _\mathrm{c}|\sigma `$, where $`\sigma ^2`$ is the variance of the density field.

The large scale structure of the universe can be characterized topologically at various scales by studying the excursion sets of the density contrast $`\delta (𝐱)`$; i.e., the sets $`E_{\delta _{\mathrm{th}}}^+=\{𝐱|\delta (𝐱)\delta _{\mathrm{th}}\}`$ or $`E_{\delta _{\mathrm{th}}}^{}=\{𝐱|\delta (𝐱)\delta _{\mathrm{th}}\}`$. An interesting choice of $`\delta _{\mathrm{th}}`$ in either cases is the percolation threshold, $`\delta _\mathrm{c}`$, for which the excursion set includes an infinite, connected structure .

Accurate determination of $`\delta _\mathrm{c}`$ is a challenge both from the theoretical point of view – no analytic calculation of $`\delta _\mathrm{c}`$ exist so far in the three-dimensional case; the experimental point of view – numerous spurious effects, such as discreteness, finiteness of the sampled volume, etc., can contaminate the measurements, and the results depend on the method employed . Here, we proceed as in Ref.  and consider site percolation on a grid, i.e., neighboring sites above/below the density threshold are connected through faces.

We analyze the topology of $`E_{\delta _{\mathrm{th}}}^\pm `$ in terms of the three dimensional local curvature of $`\delta (𝐱)`$ given by the Hessian matrix, $`^2\delta /x_ix_j`$. The number of negative eigenvalues, $`I`$, of the Hessian matrix at any point classifies the local density structure into one of four distinct classes: region with $`I=3`$ is in a “clump”(c); $`I=2`$ in a “filament” (f); $`I=1`$ in a “pancake” (p) and $`I=0`$ in a “void”(v). In particular, this classification applies to critical points of the field, where the gradient $`\delta /x_i=0`$. Connectivity of the excursion set $`E_{\delta _{\mathrm{th}}}^\pm `$ is dictated solely by the number, $`𝒞_I`$, of critical points of each class $`I`$ within $`E_{\delta _{\mathrm{th}}}^\pm `$. For clarity, we replace the numeral index $`I`$ of $`𝒞_I`$ by the alphabetic label $`I=0,1,2,3\mathrm{v},\mathrm{p},\mathrm{f},\mathrm{c}`$.

At a simple qualitative level, the role of critical points in outlining the connectivity of excursion sets can be understood quite intuitively. Connectivity at percolation is along special field lines (“ridges” and “river beds”) threading the critical points. Filament saddle points (f-saddles) are thus crucial for percolation in overdense regions ($`E_{\delta _{\mathrm{th}}}^+`$). Indeed, a large fraction of them lie along ridges connecting two local maxima. When $`\delta _{\mathrm{th}}>\delta `$ at a saddle point, the point drops out of the excursion set, thus disconnecting the corresponding clumps. Analogous reasoning can be applied to percolation in underdense regions where a large fraction of pancake saddle points (p-saddles) lie on river beds that connect local minima.

In detail, the situation can be more complicated than this simple vision, e.g., saddle-points connected through ridges or river beds. However, the work of Morse rigorously establishes a simple link between the number distribution of critical points (local property related to curvature) and global connectivity (topology) of the excursion sets of any smooth field lending credence to our above intuitive line of argument . A simple, but powerful, relation exists between the Euler characteristic, $`\chi `$ , of $`E_{\delta _{\mathrm{th}}}^\pm `$ and a linear combination of the four critical point counts, $`𝒞_I`$. It can indeed be shown that despite all the possible complexity ,

$$\chi ^\pm =\pm \left[𝒞_\mathrm{c}𝒞_\mathrm{f}+𝒞_\mathrm{p}𝒞_\mathrm{v}\right].$$
(1)
The observed structures in the galaxy distribution are believed to have evolved through gravitational instability from small initial fluctuations, usually postulated to be a Gaussian random field. We pay particular attention to the Gaussian limit since at large scales the memory of initial conditions is retained in the evolved system. At small scales, non linear clustering creates an asymmetry between overdense and underdense regions possibly changing the properties of the percolating network and the value of $`\delta _\mathrm{c}`$. To study clustering in the non linear regime we use $`N`$-body simulations.

In this paper we have two goals: (i) obtaining an accurate measurement of $`\delta _\mathrm{c}`$ as a function of smoothing scale; (ii) characterizing the topology of the percolating network through a combination of the Euler number, critical points and local curvature. We attempt to exploit eq. (1) to understand the global properties, such as percolation, of the evolving density contrast in the universe in terms of the distribution of (local) critical points.

In our numerical experiment, we consider $`N`$-body simulations (described in ) with scale-free Gaussian initial conditions, where the power-spectrum as a function of wavenumber, $`k`$, spectral index, $`n`$, and amplitude, $`A`$, is given by $`P(k)=Ak^n`$, $`n=2,1,0`$. This covers the observationally relevant range of values in the scaling regime $`1\mathrm{}100`$ Mpc. The simulations were generated with a particle-mesh code  with $`256^3`$ particles in a periodic cubic box of size $`L`$. A mesh resolution $`N_\mathrm{r}=256`$ was used to compute the forces. The density field was computed using Cloud-In-Cell interpolation  on grids of various sizes $`N_\mathrm{g}=64,128,256`$ and then convolved with Gaussian windows of scales $`\mathrm{}/L=1/N_\mathrm{g}`$, $`2/N_\mathrm{g},\mathrm{},\mathrm{\hspace{0.17em}0.04}`$, $`0.08`$ to insure sufficient differentiability . The corresponding variance of the smoothed density distribution will be denoted by $`\sigma ^2(\mathrm{})`$. To increase the dynamic range of our measurements, we take advantage of self-similar behavior for scale-free initial conditions  and combine several time snapshots of each simulation with the appropriate scaling $`(L,\mathrm{})(L/s,\mathrm{}/s)`$, where $`s`$ is the linear theory correlation length ($`\sigma _{\mathrm{linear}}^2(s)1`$). For Gaussian smoothing, $`s=[a^2A\mathrm{\Gamma }[(n+3)/2]/(4\pi ^2)]^{1/(n+3)}`$, where $`a`$ is the expansion factor of the Universe, arbitrarily fixed to unity at the beginning of the simulation.

To study the Gaussian limit, a set of $`10`$ pure Gaussian realizations of the density field, $`\delta (𝐱)`$, at each grid size and spectral index, $`(N_\mathrm{g},n)`$, were generated and smoothed at the same set of scales, $`\mathrm{}/L`$, as in the $`N`$-body samples.

The gradient and the Hessian of the density field are computed by least square fitting a second order hypersurface to the values $`\delta `$ computed at each site and its 26 closest neighbors.

We have explored the wide range of smoothing scales $`\mathrm{}/L`$ and box sizes $`N_\mathrm{g}`$ to thoroughly quantify the numerical effects contaminating topological properties. The optimal set of parameters $`(N_\mathrm{g},\mathrm{}/L)=(256,0.01)`$ is adopted to minimize the following contamination effects (partly discussed in ):

(i) The grid must be sufficiently thin compared to the smoothing scale, $`N_\mathrm{g}\mathrm{}/L2.5`$. With this constraint, our procedure for computing numerical derivatives of $`\delta `$ provides extremely accurate results.

(ii) Periodic boundaries enhance percolation in a finite box but the effect is negligible if $`\mathrm{}/L1`$. On the other hand, finite box size suppresses large scale $`[𝒪(L)]`$ power, which can be viewed as an increased effective spectral index. From our Gaussian samples, we find that combination of these two competing effects is inconsequential if $`\mathrm{}/L0.01`$.

(iii) Inherent limitations of the $`N`$-body approach: spatial resolution of the forces imposes $`N_\mathrm{r}\mathrm{}/L1`$. To evade effects of $`N`$-body relaxation and unrealistic non linear couplings between the large $`𝒪(L)`$ modes, it is required that $`2/N_\mathrm{g}^{1/3}s/L`$ and $`s/L0.1`$, respectively. Finally, discreteness effects are insignificant since an $`8`$-fold dilution of the number of particles does not change the results substantially.

First column of panels in Fig. 1 shows the measured percolation threshold, $`\delta _\mathrm{c}`$, as a function of scale, $`\mathrm{}/s`$. In agreement with previous works , the asymmetry due to gravitational clustering between overdense ($`+`$) and underdense regions ($``$), $`\delta _\mathrm{c}^+>|\delta _\mathrm{c}^{}|`$, increases with the level of non linearity and with $`n`$. It is very small for $`n=0`$. For underdense regions, the asymptotic regime $`\delta _\mathrm{c}^{}1`$ is naturally reached when $`\mathrm{}/s1`$. The value of $`\delta _\mathrm{c}^+/\sigma `$ is maximum near transition to the nonlinear regime (see also Ref. ), $`\mathrm{}/s1`$, where initial structures are optimally enhanced by large-scale motions , and is larger than the Gaussian limit, the difference increasing with $`n`$. As noted later, $`\delta _\mathrm{c}^+/\sigma `$ can become small in the highly nonlinear regime due to lost coherence.

In the Gaussian limit, we find slight increase in $`\delta _\mathrm{c}`$ with $`n`$; $`|\delta _\mathrm{c}/\sigma |=1.6\pm 0.1`$, $`1.4\pm 0.1`$, $`1.3\pm 0.1`$ for $`n=2`$,$`1`$,$`0`$, respectively. Error bars reflect the dispersion in the $`10`$ realizations and possible contaminations (i) and (ii) discussed above. These results should be contrasted with that of Ref. , where $`|\delta _\mathrm{c}/\sigma |1`$ is obtained independently of $`n`$. However, the measurements in Ref.  use a smaller grid, $`N_\mathrm{g}=64`$, and might be affected by effects (ii) mentioned above. The other difference is that the field is softened by a power-spectrum cut-off at small scales instead of our Gaussian smoothing.

Second column of panels in Fig. 1 displays the fraction $`f_I`$ of space occupied by clump, filament, pancake and void regions in $`E_{\delta _\mathrm{c}}^+`$. As shown earlier from a qualitative point of view , the simulated universes at percolation are seen to be composed mostly of filamentary and clumpy regions, in the respective fractions of $`55\%`$ and $`3540\%`$ in the Gaussian limit. The filamentary nature becomes dominant in the non linear regime and the results do not depend significantly on $`n`$.

Third column of panels in Fig. 1 shows the measured fraction of critical points of each kind. Note that maxima are dominant, explaining the rather clumpy appearance of the largely filamentary (in terms of curvature) percolating network in Fig. 2. In fact, clumpiness augments with non linearity as the fraction of maxima increases and explains the rather small $`\delta _\mathrm{c}^+/\sigma `$ despite increased filamentarity when $`\sigma 1`$.

Fourth column of panels in Fig. 1 shows the “reduced” Euler number $`\chi /𝒞_{\mathrm{tot}},𝒞_{\mathrm{tot}}=_I𝒞_I`$, as a function of density threshold measured at various stages of the simulations, corresponding to $`a=1,32,128`$ for $`n=0,1`$ and $`a=1,16,32`$ for $`n=2`$. The value of $`a`$ increases with the abscissa of the left local maximum of $``$. For $`a=1`$, the measurements agree almost perfectly with the analytic prediction in the Gaussian limit , as expected. The $``$ curves retain overall shape with time, except that they are progressively compressed near the origin and stretched at large $`\delta _{\mathrm{th}}/\sigma `$. This perhaps indicates that non linear clustering does not significantly influence the topology of the system.

There are three regimes dictated by the sign of the first derivative of $``$. With increasing $`\delta _{\mathrm{th}}/\sigma `$, (i) at small $`\delta _{\mathrm{th}}/\sigma `$, mostly local minima drop out of $`E_{\delta _{\mathrm{th}}}^+`$, creating cavities, thus increasing $``$; (ii) at intermediate $`\delta _{\mathrm{th}}/\sigma `$, p-saddles start to drop out too, and cavities connect together, thus the value of $``$ decreases; (iii) in the large $`\delta _{\mathrm{th}}/\sigma `$ region even f-saddles drop out, breaking up the ridges to create isolated clusters, thus $``$ increases again. At this last stage, the contribution of minima and p-saddles to the Euler number can be largely neglected, $`\chi 𝒞_\mathrm{c}𝒞_\mathrm{f}`$.

Fourth column of Fig. 1 embodies the main result of the paper. It shows that percolation occurs at positive values of the Euler characteristic $`\chi >0`$ in regime (iii), and specifically at $`\chi N_{\mathrm{clus}}`$, where $`N_{\mathrm{clus}}`$ is the number of “isolated clusters” (disconnected regions of $`E_{\delta _{\mathrm{th}}}^+`$). This means that at percolation threshold, isolated clusters, including the largest percolating one, have close to tree structure, i.e. they contain no or few loops. Indeed, in this case each cluster contributes unity to the total Euler number adding up to $`\chi N_{\mathrm{clus}}`$. The f-saddles, the only type available in number in the regime (iii) can produce only loops, which decrease Euler number of individual clusters, and being abundant would decrease the total $`\chi `$ much below $`N_{\mathrm{clus}}`$, which is not what is observed. Thus we conclude that the percolating network is a tree. Fig. 2 gives a visual presentation of a percolating tree.

In summary, we used very simple but novel statistics in cosmology – critical points of the density field and its local curvature to characterize the topological properties of the universe at percolation. The most important results of our analyses  are the following:

(a) At percolation, the Euler characteristic is positive, equal to the number of isolated clusters suggesting that most of these have tree structure (with almost no loops). This result is valid for overdense as well as underdense excursions. In particular, we find in the Gaussian limit that $`|\delta _\mathrm{c}|/\sigma >1`$, at variance with common belief (e.g. Ref. ) that site percolation in $`E_{\delta _{\mathrm{th}}}^\pm `$ should occur at $`|\delta _\mathrm{c}|/\sigma =1`$, where the Euler number $`\chi ^\pm `$ of $`E_{\delta _\mathrm{c}}^\pm `$ is zero.

(b) In terms of curvature, the overdense structures at percolation are mostly filamentary and clumpy in agreement with previous works . Filamentarity increases with level of non linearity. However, maxima are the dominant critical points in the excursion set. With isotropic smoothing, this explains its rather clumpy appearance as illustrated by Fig. 2.

Our efforts to identify and address sources of contamination in our numerical experiments suggest that it is actually difficult to measure the percolation threshold accurately in a real galaxy catalog, where one of the main issues is contamination by discreteness effects. These impose a large smoothing scale, a fact difficult to reconcile with the requirement $`\mathrm{}/L0.01`$, even for current or forthcoming large redshift surveys such as the Two-degree Field  or the Sloan Digital Sky Survey .

However, it is still possible to measure local curvature and critical point counts, because these statistics should not require as small values of $`\mathrm{}/L`$ to be measured with sufficient accuracy. Although, here we focus on critical point counts at the percolation threshold, these counts as functions of $`\delta _{\mathrm{TH}}`$ are of interest. In particular, the statistics $`𝒞_I`$ and $`f_I`$ can be computed explicitly as functions of $`\delta _{\mathrm{TH}}`$ in the Gaussian limit . As the genus of the boundary of $`E_{\delta _{\mathrm{th}}}^\pm `$ , they could be used in real galaxy catalogs to constrain the nature of the primordial density field. To do that properly, several issues remain to be addressed, such as the anisotropies induced by redshift distorsion in three-dimensional galaxy catalogs or the effects of biasing between the observed galaxy distribution and the underlying total matter distribution.

We thank E. Hivon for kindly providing the $`N`$-body simulations and S. F. Shandarin and J. R. Bond for fruitful discussions. TS acknowledges support through NSF CAREER grant AST-9875031. Computer time for the $`N`$-body simulations was allocated by the scientific council of IDRIS, Orsay. DP & TS thank IAP and SC thanks CITA for hospitality.


text:
# CaB6: a new semiconducting material for spin electronics.

\[

## Abstract

Ferromagnetism was recently observed at unexpectedly high temperatures in La-doped CaB6. The starting point of all theoretical proposals to explain this observation is a semimetallic electronic structure calculated for CaB6 within the local density approximation. Here we report the results of parameter-free quasiparticle calculations of the single-particle excitation spectrum which show that CaB6 is not a semimetal but a semiconductor with a band gap of 0.8 eV. Magnetism in LaxCa1-xB6 occurs just on the metallic side of a Mott transition in the La-induced impurity band.

\]

The recent observation of ferromagnetism in La doped alkaline-earth hexaboride compounds at high temperatures presents three puzzles. Firstly, ferromagnetism is usually associated with elements with a partly filled 3d, 4f or 5f shell. Secondly, ferromagnetic ordering is only observed for a narrow dopant concentration range and for a surprisingly low dopant concentration. The maximum observed moment in Ca1-xLaxB6 is 0.07 Bohr magnetons per La atom for $`x=0.005`$. Thirdly, and most surprising of all is the observation of very large Curie temperatures; for La0.01Ca0.99B6 a value close to 900 K has recently been reported. Such a large value virtually excludes the possibility of the observed magnetism being related to magnetic impurities.

In an initial comment on the experimental observations, Ceperley suggested that it might be an example of the long-predicted but never observed ferromagnetic phase of a dilute electron gas; improved calculations increase the estimated density at which this might occur. An alternative explanation proposed by Zhitomirsky was that the ferromagnetic hexaborides might be doped excitonic insulators . This explanation requires that the exciton binding energy should be comparable in size to the single-particle band gap. Jarlborg suggested that the magnetism may be conventional itinerant magnetism. Although magnetism in materials which do not have partly filled d or f shells is rare, it is not unprecedented. A handful of materials does exist, of which ZrZn2 is the best known example, in which the Fermi energy falls at or close to an exceptionally narrow peak in the electronic density of states so that the Stoner criterion for the occurence of itinerant magnetism is fulfilled. The problem posed by the hexaborides is not so much the occurence of magnetic ordering but rather the strength of the magnetic coupling as reflected by the very high Curie temperature, and this was not estimated.

All three suggestions are based on electronic structure calculations which indicate that stoichiometric CaB6 is a semimetal with a very small overlap between the conduction and valence bands  or might have a very small gap . Since the theoretical predictions rely on details of this electronic structure, it is important to examine their validity. Massidda’s and Rodriguez’ studies of the electronic structure were based on the FLAPW calculation method which is probably the most accurate available. These calculations were carried out within the framework of density functional theory (DFT) using the local density approximation (LDA). Such calculations are known to be capable of yielding total energies with a very useful accuracy. However, it is also well known that the eigenvalue spectrum which results from solving the Kohn-Sham equations of DFT cannot be interpreted unreservedly as an excitation spectrum and sometimes the results are spectacularly wrong. In particular, the band gaps of semiconductors are typically underestimated by 50% and in extreme cases such as Ge the conduction and valence bands are found to overlap resulting in metallic or semimetal character.

In order to study single-particle excitations, one should solve Dyson’s equation for the single-particle Green’s function expressed in terms of the self-energy operator $`\mathrm{\Sigma }`$. $`\mathrm{\Sigma }`$ can be expanded as a perturbation series in the Green’s function $`G`$ and the dynamically screened Coulomb interaction $`W`$. The so-called $`GW`$ approximation introduced by Hedin includes only the first term in this series. In addition, one usually assumes a quasi-particle (QP) expression for the Green’s function $`G`$. For a large number of semiconductors and insulators such calculations produce band gaps which are in very close quantitative agreement with experimental single-particle band gaps. The main purpose of the present paper is to apply these parameter-free calculations with their proven predictive capability to CaB6. From our results we shall conclude that, contrary to what is currently being assumed, the parent material is actually a semiconductor with quite a large band gap. This finding has far-reaching consequences for understanding the basic properties of the doped, ferromagnetic phase and opens up the prospect of a range of novel applications.

The starting point for our $`GW`$ study is an LDA calculation for CaB6 performed with a plane-wave basis and using norm-conserving pseudopotentials to describe the interaction between the valence electrons and the ionic cores. Using the experimentally determined structure, we reproduce essentially perfectly the energy bands calculated by Massidda using an all-electron method. The results are shown in Fig. 1(a). There are ten valence bands (bonding orbitals derived from the boron 2s and 2p levels) well-separated from the conduction bands (antibonding states) with the exception of an overlap of about 0.3 eV in a small region around the X point of the highest occupied and lowest unoccupied bands. The occurrence of such small band overlap (or the presence of a comparably small band gap) is an important criterion for the stability of an excitonic insulator phase. The effective masses in units of the free electron mass for the bands near X, $`m_e^{}=0.22`$ and $`m_e^{}=0.50`$ for the electrons and $`m_h^{}=0.20`$ and $`m_h^{}=1.8`$ for the holes, are close to values reported earlier . We use the LDA electronic wave functions and eigenvalues as input for the $`GW`$ calculations for which we adopt the space-time approach suggested by Rojas in which all operators are represented on grids in real and reciprocal space, in the time and energy domains; details of our implementation have been given by van der Horst. By varying the size of these grids we estimate that the calculated QP energies are converged within 0.1 eV. The densest grids we used were a $`(12\times 12\times 12)`$ real-space grid and a $`(6\times 6\times 6)`$ k-grid. The Green’s function is constructed including the LDA wave functions and energies of the lowest 300 bands.

The results of the $`GW`$ calculations are plotted in Fig. 1(b). Throughout the Brillouin zone, the QP corrections to particular LDA bands are rather uniform so that the dispersion of the $`GW`$ bands is very similar to the LDA band dispersion and the effective masses differ only slightly from the values obtained from the LDA calculation. They are reduced on average by less than 10%, to $`m_e^{}=0.20`$ and $`m_e^{}=0.48`$ for the electrons and $`m_h^{}=0.17`$ and $`m_h^{}=1.7`$ for the holes. The sign and size of the corrections to the energy bands depend on their energy and their wavefunction character and vary considerably from band to band. Of particular importance are the relative shifts of the bands at X near the Fermi level. The hole band and electron bands are moved, respectively, downwards and upwards in energy resulting in the opening of a sizeable band gap of about 0.8 eV. The relative shifts can be understood from the bonding and antibonding character of the wave functions in the valence and conduction bands, respectively. This is quite analogous to the situation in silicon. The largest shift in the occupied bands is calculated for the lowest valence band which is lowered in energy by about 1.9 eV.

Our finding that bulk CaB6 is a semiconductor is consistent with the large value of the low-temperature resistivity and with its increase with decreasing temperature recently reported by Ott and by Vonlanthen and also with recent results of angle-resolved photoemission experiments in which a gap of a size similar to what we calculate was found. Because of the sensitivity of photoemission experiments to surface electronic structure, the authors interpreted their results exclusively in terms of the surface electronic structure. Our result indicates that there might indeed be a slight increase in the value of the band gap at the surface, something which is known to happen in heteropolar semiconductors.

The result for the stoichiometric phase has important implications for the models proposed to explain the ferromagnetism observed in the doped phase. Using our calculated values for the effective masses and static dielectric constant of 5, we obtain a scaled density parameter $`r_\mathrm{s}`$ of 1.5, and an exciton binding energy of about 0.07 eV. The scaled density parameter is still much lower than the values of 80 and 20 reported respectively, by Ceperley and Ortiz, for which a ferromagnetic phase of the low density electron gas might be expected. The estimated excitonic binding energy is much smaller than the calculated QP band gap whereas the excitonic insulator models require them to be comparable in size. Instead we note that the ferromagnetism occurs for an impurity concentration which is just on the metallic side of a Mott transition in the impurity band: $`n^{1/3}a_H=0.4`$, where $`n`$ is the dopant concentration $`7\times 10^{19}`$ electrons cm-3 and $`a_H=10`$ Å is the Bohr radius of an isolated effective mass impurity.

Support for this possibility can be found in the very recent experimental studies on doped hexaborides reported by Terashima. He found that some of his La-doped samples were paramagnetic, exhibiting the Curie-Weiss behaviour which is to be expected if the dopant concentration were below Mott’s critical value. In other samples he found no indication of saturation of the magnetic moment even in the highest magnetic field reported (10 kOe). This behaviour might be expected for partly filled narrow bands. In addition Terashima found much larger values of the magnetic moment per dopant atom than reported by Young and Ott. He also attributed a number of reported de Haas-van Alphen oscillation frequencies to Al inclusions. We note that it is extremely important to know where the La goes and if all of it actually occupies Ca sites as usually assumed since the effective valence of an impurity atom in a semiconductor is very sensitive to the local environment. Assuming that the magnetism occurs mainly in the impurity band, we are currently attempting to estimate the Curie temperature from total energy calculations and to determine the parameters which govern the physical behaviour of the doped system such as the hopping matrix elements, the charge and spin fluctuation parameters (Hubbard U and Stoner I, respectively) and the binding energy of the impurity state. The last quantity, the position of the impurity state with respect to the bottom of the conduction band, is particularly difficult to estimate reliably. It requires knowing the central cell potential which is known to make a substantial correction to the binding energy of the lowest effective mass states.

Irrespective of the origin of the high Curie temperature in the doped material, the existence of a large gap for the undoped material means a new class of magnetic semiconductors is available with the unique characteristic of having a Curie temperature sufficiently high to allow room-temperature operation of semiconducting spin devices, something which has only very recently been achieved at low temperatures. With the hexaborides it should be possible to inject a spin-polarized current from the doped into the undoped material in the ballistic regime and study the spin dynamics as a function of temperature, and of current density without the complication of having to apply an external magnetic field or without the problems encountered when attempting to inject spins from a ferromagnetic metal into a conventional semiconductor in the diffusive regime. It should also be possible to make a field-effect device to modulate a spin-polarized current in the doped hexaboride material and achieve gain. If a $`p`$-type hexaboride with some of the Ca replaced by a monovalent ion were also found to be ferromagnetic, it should even be possible to make the spin-analogues of bipolar devices.

Both Ca and B are light elements so that the spin-orbit coupling is small. Since the lowest conduction band and highest valence band are both singly degenerate, the intrinsic spin-flip scattering should be weak with spin-flip scattering lengths comparable to what has recently been observed in carbon nanotubes. Because of the small size of the spin-orbit coupling, because the magnetic moment is so tiny, and because the hexaborides are cubic, the material should be magnetically extremely soft with a very small intrinsic magnetocrystalline anisotropy. This could be important for sensor applications.

P.J. Kelly is grateful to G.E.W. Bauer and J. van den Brink for useful discussions. P. van Gelderen and P.A. Bobbert wish to acknowledge discussions with R. Monnier and S. Massidda. This work is also part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) and Philips Research.


text:
# 1 Introduction

## 1 Introduction

It has been conjectured recently  that the Pomeron instead of being an infinite gluon ladder  may appear as a finite sum of gluon ladders corresponding to a finite sum of Regge multipoles with increasing multiplicities. The first term in the $`\mathrm{ln}s`$ series contributes to the total cross-section with a constant term and can be associated with a simple pole, the second one (double pole) goes as $`\mathrm{ln}s`$, the third one (triple pole) as $`\mathrm{ln}^2s`$, etc. All Pomeron poles have unit intercepts. A new prong opens each time the available energy exceeds the threshold value. The values of threshold parameter were found  from a fit to the experimental data. Another important set of parameters is related to the coupling of the gluons to hadrons, giving the relative weight of various diagrams. In ref.  they were fitted to the forward scattering data and subsequently they were also calculated in QCD. Below we study this problem both in the forward and non-forward directions. Because of its complexity, we do not consider here the effect of the progressively opening channels as it was done in ref. . Instead we first consider a model with a simple and dipole Pomeron \[DP\] contribution used earlier to fit elastic hadron scattering , low-$`x`$ structure functions  and photoproduction of vector mesons at HERA .

The DP ansatz reads

$$P(s,t)=isg_0\underset{i=1}{\overset{2}{}}c_iR_i^2(s)e^{R_i^2(s)t},$$
(1)
where $`c_1=1`$, $`c_2=\lambda b1=\epsilon `$; $`R_1^2`$$`(s)=\alpha ^{}(b+\mathrm{ln}(i\frac{s}{s_0}))`$; $`R_2^2`$$`(s)=\alpha ^{}\mathrm{ln}(i\frac{s}{s_0})`$ . Apart from the normalization factor $`g_0`$, this model contains 4 adjustable parameters: $`\lambda `$, $`b`$, $`s_0`$ and $`\alpha ^{}`$, moreover their number can be still reduced within the cone region (i.e. $`|t|1`$ $`GeV^2`$ by setting $`\lambda =\frac{1}{b}`$). Here a linear Pomeron trajectory, $`\alpha \left(t\right)=1+\alpha ^{^{}}t`$ is implied. Generalization to include nonlinear trajectory is straightforward. The model is strongly constrained by an integral relation between the residue at the simple and double poles , that produces the observed dip-bump structure in the differential cross-section in accordance with the experimental data.

$`\begin{array}{cc}{}_{}{}^{}\text{e-mail address:}\hfill & \text{jeno@kontr.uzhgorod.ua}\hfill \end{array}`$

$`\begin{array}{cc}{}_{}{}^{}\text{e-mail address:}\hfill & \text{sasha@len.uzhgorod.ua}\hfill \end{array}`$

$`\begin{array}{cc}{}_{}{}^{}\text{e-mail address:}\hfill & \text{iep@iep.uzhgorod.ua}\hfill \end{array}`$

The Dipole Pomeron model with a unit Pomeron intercept  was used to describe successfully $`pp`$ and $`\overline{p}p`$ elastic scattering at the ISR and the Collider energies, however the ratio $`\sigma _{el}/\sigma _{tot}`$ was found  to decrease asymptotically for all physical values of the parameters. To obtain an increasing function for this ratio the authors of  introduced a factor corresponding to the supercritical Pomeron behavior . However in this case the unitarity is broken. It was shown  that the unitarity violation occurs at energies only slightly above the Tevatron energy of 1.8 TeV, and therefore it is a problem of the present and not of the future. To avoid this problem, we consider a model containing a finite series of Pomeron terms up to $`\mathrm{ln}^2s`$, in accordance with the unitarity constrains.

## 2 The model

Our ansatz for the scattering amplitude is:

$$A_{pp}^{\overline{p}p}=P+R_f\pm R_\omega ,$$
(2)
where we introduce a tripole contribution to the DP in the simplest way (see also ):

$$P(s,t)=isg_0\left[b+\mathrm{ln}\left(i\frac{s}{s_0}\right)+c\mathrm{ln}^2\left(i\frac{s}{s_0}\right)\right]e^{\alpha ^{}bt}e^{\left(\alpha _P(t)1\right)\mathrm{ln}(\frac{is}{s_0})}.$$
(3)
The Pomeron trajectory:

$$\alpha _P(t)=\alpha _P(0)+\alpha _P^{}t+\gamma _P\left(\sqrt{t_P}\sqrt{t_Pt}\right),$$
(4)
where the two-pion threshold $`t_P=4m_\pi ^2`$.

In (2) the $`R_f,R_\omega `$ contains the subleading Reggeon contributions ($`f`$ and $`\omega `$) to the scattering amplitude:

$$Rj(s,t)=g_j\left(i\frac{s}{s_0}\right)^{\alpha _j\left(t\right)}e^{b_jt},\alpha _j\left(t\right)=1+\alpha _j^{^{}}t,j=f,\omega ;s_0=1GeV^2.$$
(5)
Recently in paper  the contribution of truncated BFKL Pomeron series to the $`\sigma _{tot}`$ of $`pp`$ and $`\overline{p}p`$ scattering was studied and it was shown that a reliable description can be obtained by using two orders in this series. As a by-product, the elastic differential cross section was obtained for the diffraction cone at low and high energies with a qualitative description of the experimental data. Contrary to , in our model we performed a simultaneous fit to the $`\sigma _{tot}`$ and $`d\sigma /dt`$ data as follows.

## 3 Comparison with experiment. Conclusion

In order to determine the parameters that control the s-dependence of $`A(s,0)`$ in a wide energy range $`4GeV\sqrt{s}1800GeV`$, we used the available data for total cross sections -. A total of 66 experimental data has been included for $`t=0`$. For the differential cross sections we selected the data at the energies $`\sqrt{s}=19;23;31;44;53;62GeV`$ (for $`pp`$scattering) and $`\sqrt{s}=31;53;62;546;1800GeV`$ (for $`\overline{p}p`$scattering). The squared 4-momentum has been limited by $`0.05GeV^2<|t|<0.5GeV^2`$, because at larger $`|t|`$ the influence of the dip region becomes visible (in particular, it can be seen quite clearly at the Collider energy $`\sqrt{s}=546GeV`$). The total of 729 experimental points have been used in the overall fit.

In the calculations we use the following normalization for the dimensionless amplitude:

$$\sigma _{tot}=\frac{4\pi }{s}ImA\left(s,t=0\right),\frac{d\sigma }{dt}=\frac{\pi }{s^2}\left|A(s,t)\right|^2.$$

The resulting fits for $`\sigma _t`$, $`\frac{d\sigma }{dt}`$ are shown in Figs. 1-3 with the values of the fitted parameters quoted in Table 1. From these figures we conclude that the multipole Pomeron model corresponding to a sum of gluon ladders up to two rungs fits the data perfectly well in a wide energy region within the diffraction cone. As result, the model gives a good behavior for the ratio $`\sigma _{el}/\sigma _{tot}`$ for $`\overline{p}p`$scattering (see Fig. 4).

The rapid decrease of the coefficients of the first three terms (approximately as $`1:\frac{1}{10}:\frac{1}{100}`$) in the Pomeron series in (3) provides the fast convergence of the series and ensures the applicability of this approximation at still much higher energies.

| $`g_0`$ | $`b`$ | $`c`$ | $`\alpha _P^{}`$ | $`\gamma _P`$ | $`\alpha _f\left(0\right)`$ | $`b_f`$ | $`g_f`$ | $`\alpha _\omega \left(0\right)`$ | $`b_\omega `$ | $`g_f`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | $`\left(GeV^2\right)`$ |  | $`\left(GeV^2\right)`$ | $`\left(GeV^1\right)`$ |  | $`\left(GeV^2\right)`$ |  |  | $`\left(GeV^2\right)`$ |  |
| $`0.253`$ | $`7.46`$ | $`0.180`$ | $`0.266`$ | $`0.137`$ | $`0.777`$ | $`3.04`$ | $`13.3`$ | $`0.524`$ | $`6.86`$ | $`7.97`$ |

Table 1: Values of the fitted parameters in the model.

In this paper we have explored only the simplest extension of the dipole Pomeron to a tripole. In fact, the scattering amplitude is much more complicated than just a simple power series in $`\mathrm{ln}s`$.

Earlier a comprehensive analysis of the $`pp`$ and $`\overline{p}p`$ diffraction cone scattering using the dipole and supercritical Pomeron models was done .

On the one hand, though we used just a simplified t-dependence in the model reasonably good results were obtained. Because the slopes of secondary Reggeons do not influence the fit sufficiently, we have fixed them at $`\alpha _f^{}=0.84`$ $`(GeV^2)`$ and $`\alpha _\omega ^{}=0.93`$ $`(GeV^2)`$, which correspond to the values of Chew-Frautschi plot. On the other hand, we included the curvature of the Pomeron trajectory that cannot be negligible. It should be also taken into account that the calculated slope of the Pomeron value is close to the generally accepted value $`\alpha _P^{}=0.25`$ $`(GeV^2)`$.

The quality of our fits $`\chi ^2/dof=2.25`$ has not reached the best fit $`\chi ^2/dof=1.3`$ obtained in . Nevertheless, we believe this difference can be significantly reduced after adding the data of Coulomb interference region $`|t|<0.05GeV^2`$ to the global fit.

Acknowledgments We thank L. Csernai, L. Jenkovszky and F. Paccanoni for discussions and BCPL, Bergen University, where this paper was completed, for hospitality and financial support. The work was supported also by the Ukrainian-Hungarian grant 2M/200-2000.


text:
# A repulsive trap for two electrons in a magnetic field

## I Introduction

Recent technological development allowed to create various types of surface superlattices for two-dimensional (2D) electron gas in semiconductor heterostructures and to investigate their transport properties in presence of magnetic field. Experiments with antidot lattices were carried out by different experimental groups (see e.g. ) and the contribution of classical periodic orbits in the resistivity peaks at certain values of magnetic field was clearly identified. In fact the size of antidots and the distance between them are relatively large, and analysis of classical trajectory dynamics can be successfully applied to understand a number of unusual transport properties in such antidot arrays . Due to nonlinearity of motion in the vicinity of antidot potential the classical dynamics can be chaotic that leads to diffusive spreading of trajectories even for perfectly periodic lattices . If the distance between antidots is large or comparable to the cyclotron radius of an electron in a magnetic field perpendicular to the lattice, then one should understand first, the properties of electron dynamics near one antidot. In the absence of electric field, the dynamics is integrable for an antidot of circular shape due to angular momentum conservation and an electron always regularly rotates around the antidot. An electric field applied in the 2D plane of the superlattice breaks the cylindrical symmetry and can lead to the electron escape from the antidot to infinity. The problem of electron dynamics in crossed electric $`\mathrm{E}`$ and magnetic $`\mathrm{B}`$ fields near a circular elastic disk (antidot) was studied in . It was shown that the dynamics can be described by a simple area-preserving map which depends only on one dimensionless parameter $`ϵ=(2\pi m/ae)(\mathrm{E}/\mathrm{B}^2)`$, where $`m,e`$ are electron mass and charge, and $`a`$ is the disk radius. For small $`ϵ<ϵ_c`$ the electron dynamics in the phase space of angular momentum $`l`$ and conjugated angle $`\varphi `$ is bounded by the invariant Kolmogorv-Arnold-Moser (KAM) curves so that electron always remains near the disk. On the contrary for $`ϵ>ϵ_c`$ the KAM curves are destroyed, global chaos sets in and the electron escapes to infinity after few collisions with the disk.

Until now the classical dynamics in antidot lattices was studied only for noninteracting electrons . In this paper for the first time we analyze the effect of Coulomb interaction between classical electrons in the vicinity of an antidot. We show that for sufficiently strong interaction between electrons their dynamics becomes chaotic. Due to that one or two electrons can escape from the antidot even in an arbitrary weak applied electric field $`E`$ that corresponds to $`ϵ_c0`$ contrary to $`ϵ_c>0`$ in absence of interaction. We determine the dependence of average escape rate $`\mathrm{\Gamma }`$ on $`ϵ`$ showing that in the limit of small electric field $`\mathrm{\Gamma }ϵ^2`$. After that we also discuss the two electron propagation in the antidot superlattice.

The paper is organized as follows. In the next Section we briefly discuss the one electron dynamics near the antidot in crossed magnetic and electric fields. In Section III the dynamics of two interacting electrons is analyzed in detail. The electron motion in the antidot superlattice is considered in Section IV. In the last Section we summarize the obtained results.

## II One electron dynamics

The dynamics of an electron in crossed electric and magnetic fields in two dimensions with one antidot is described by the Hamiltonian:

$$H_0=(𝐩e𝐀)^2/2m+U(x,y)e\mathrm{𝐄𝐫}$$
(1)
where $`𝐀`$ is the vector potential and $`U(x,y)`$ describes the antidot potential which depends only on the radius $`r=\sqrt{x^2+y^2}`$ with $`U=0`$ for $`ra`$. For convenience, following , we introduce the dimensionless variables $`\stackrel{~}{x}=x/a,\stackrel{~}{y}=y/a,\stackrel{~}{t}=t/\tau _0,\stackrel{~}{H_0}=H_0/2ϵ_F,\stackrel{~}{U}=U/2ϵ_F,\stackrel{~}{B}=B/B_0,\stackrel{~}{E}=E/E_0`$ where $`ϵ_F(v_F)`$ is the Fermi energy (velocity), $`\tau _0=(2ϵ_F/ma^2)^{1/2}=a/v_F`$, magnetic and electric fields are scaled by $`B_0=(mϵ_F)^{1/2}/ea`$ and $`E_0=2ϵ_F/ea`$ respectively. In these units a magnetic field $`B=B_0`$ gives the cyclotron radius $`R_c=a`$ for electron with energy $`ϵ_F=mv_F^2/2`$. We choose the Landau gauge $`𝐀=(By,0,0)`$. Then omitting the tildes the Hamiltonian equations of motions reads:

$$\begin{array}{c}dx/dt=v_x,dv_x/dt=Bv_ydU/dxE_x,\\ dy/dt=v_y,dv_y/dt=Bv_xdU/dyE_y\end{array}$$
(2)
where $`v_x=p_x+yB,v_y=p_y`$ . To model the antidot we chose the potential:

$$U(x,y)=U_0(1r)^6.$$
(3)
Usually we choose $`U_0`$ to be much larger than the electron energy $`H_0`$ so that this potential becomes very similar to an absolutely rigid disk with effective radius $`a_{eff}`$ about 15% smaller than $`a`$.

Far from the antidot the equations of motion are exactly solvable and give electron rotation over a circle of cyclotron radius $`R_c=v/\omega _c`$ with cyclotron frequency $`\omega _c=B`$. In addition this circle moves with the drift velocity $`v_d=E/B`$ in the direction perpendicular to electric and magnetic fields. As it was found by Berglund et al. , near the antidot, the dynamics strongly depends on the dimensionless parameter $`ϵ=v_d2\pi /\omega _c=2\pi E/B^2`$. For $`ϵ1`$ the electron scatters on the antidot and escapes to the infinity after one collision. On the contrary, the situation with not very large $`ϵ`$ is much richer . In this case electron can collide many times with the antidot and this process is described by a simple area-preserving map :

$$\begin{array}{c}\overline{\varphi }=\varphi +\pi 2\mathrm{sin}^1\beta \\ \overline{\beta }=\beta ϵ\mathrm{sin}\overline{\varphi }\end{array}$$
(4)
where bars denote the new values of variables after collision, $`\varphi `$ is the scattering angle measured in respect to the direction of drift velocity and $`\beta `$ is the scattering impact parameter divided by the antidot radius. In this way $`\beta `$ varies in the interval (-1, 1). We note that $`\beta `$ can be also considered as the orbital momentum $`l`$ of the electron divided by the maximal momentum $`l_{max}=av`$ at which electron still collides with the antidot. The real dynamics is correctly described by the map if $`R_c1`$, that corresponds to $`vB`$. For $`ϵ1`$ the variation of $`\beta `$ is bounded by the invariant KAM curves and electron is trapped near the antidot. The last KAM curve is destroyed for $`ϵ>ϵ_c0.45`$  so that the orbits with initial $`\beta 0`$ can escape from the antidot to infinity. Of course, for $`ϵ>ϵ_c`$ some islands with regular motion inside still remain, but they become very small as soon as $`ϵ`$ becomes significantly lager than $`ϵ_c`$.

To study the electron dynamics in (1) the Hamiltonian equations of motion are solved numerically by Runge-Kutta method of fourth order so that the electron energy is conserved with the relative precision better than $`10^6`$. The examples of the Poincaré cross sections constructed at $`x=0`$ and $`v_x>0`$ for trajectories trapped near the antidot is shown in Figs. 1,2. In Fig. 1 $`ϵ0.16`$ is rather small and almost all phase space is filled by integrable KAM curves. For Fig. 2, the parameter $`ϵ0.42`$ is close to $`ϵ_c`$ and KAM curves become more deformed and the chaotic component becomes visible. This case can be compared with Fig. 3 in  where the cross section for the map (4) is given for a close value of $`ϵ`$.

It is interesting to note that the map (4) can be not valid if orbits have $`|\beta |1`$ or $`|\beta |>1`$. For example, if the antidot is inside a large cyclotron circle then electron will make many many rotations before this slowly drifting circle will cross the antidot. This situation is not taken into account by the first equation in (4). An example of the electron dynamics in this case is given in Fig. 3. Here $`ϵ0.16`$ is small and the motion is still regular. We should stress that such type of trajectories separates orbits which escape to infinity and those which collide with the antidot on each cyclotron period. For the antidot superlattice with antidot spacing comparable with $`R_c`$ this type of orbits (see Fig. 3) is of special importance since these orbits can easily jump from one antidot to another leading to a global diffusion in the system. We will discuss this situation in Section III.

## III  Effects of Coulomb interaction on electron dynamics

Let us now consider how the Coulomb interaction between two electrons affects their dynamics near the antidot. In this case the Hamiltonian of the system reads

$$H=H_0(𝐩_1,𝐫_1)+H_0(𝐩_2,𝐫_2)+e^2/|𝐫_1𝐫_2|.$$
(5)
While in the free space the Coulomb interaction repels the electrons and leads to their separation the situation is more complicated in the presence of magnetic field. In the case without any antidot the total momentum of two electrons is conserved and as a result each electron rotates regularly on a cyclotron circle which in addition rotates around the center of mass of the system. Without an external electric field ($`E=0`$) the center of mass is fixed and inert whereas in the presence of the field ($`|E|>0`$) the center of mass drifts with the constant velocity $`v_d=E/B`$ but the average distance between electrons remains constant. However this electron pair can be trapped by the repulsive potential of the antidot so that the electrons will spend a long time colliding with this antidot. An example of electron dynamics in this case is shown in Fig. 4. It shows that an electron can escape from the antidot even in the situation with $`ϵ<ϵ_c`$ when without the interaction the electrons remain trapped near the antidot. In our numerical simulations we observed different cases where one or both of the electrons escape to infinity.

To investigate how the escape rate depends on the strength of an external electric field $`E`$ we studied the ensemble of 100 paths. In each path the positions and momentums of each electron are chosen randomly in the intervals $`4x,y4`$, $`2p_x,p_y2`$ in such a way that the total energy is $`H15\pm 0.5`$. We remind that the antidot with potential (3) is placed at $`(0,0)`$ and $`U_0=1000`$. The escape rate $`\mathrm{\Gamma }`$ is defined as $`\mathrm{\Gamma }=1/T`$ where $`T`$ is the time after which the distance of one of the electrons from the antidot is greater than $`R_{esc}5R_c10`$. This distance is sufficiently large and as soon as it is reached an electron escapes to infinity and never returns back to the antidot. The average value of $`\mathrm{\Gamma }`$ is obtained by averaging over $`100`$ values obtained for $`100`$ randomly chosen paths.

The dependence of the escape rate $`\mathrm{\Gamma }`$ on the strength of the applied electric field $`E`$ is presented in Fig. 5. It definitely shows that the escape takes place even at very weak electric fields with $`ϵϵ_c`$ when without Coulomb interaction electrons are forever trapped near the antidot. According to the obtained numerical data (see Fig. 5) in the limit of $`ϵ0`$ the escape rate is

$$\mathrm{\Gamma }/\omega _cϵ^2.$$
(6)
Our understanding of this dependence is based on the following argument. Due to the Coulomb interaction between electrons their dynamics in the vicinity of the antidot becomes chaotic. Therefore the phase $`\varphi `$ in the map (4) changes randomly between electron collisions with the antidot and $`\beta `$ grows diffusively with the number of collisions $`n`$ so that $`(\mathrm{\Delta }\beta )^2Dn`$ with $`D=ϵ^2/2`$. This diffusion results in the escape rate $`\mathrm{\Gamma }/\omega _cDϵ^2`$ being in agreement with the numerical data in Fig. 5.

## IV  Two electron propagation in antidot superlattice

We also studied the electron dynamics on a square antidot superlattice, when the antidot potential is given by (3) and the distance between antidots $`d>2`$. In this case our results for one electron dynamics are in qualitative agreement with the conclusions drawn in . As soon as the cyclotron radius $`R_c`$ becomes comparable with the antidot spacing $`d`$ the trajectories start to move diffusively on the whole lattice. Trapped orbits near one antidot exist only for $`2R_c<d`$ and $`ϵ<ϵ_c`$. For $`R_c>d/2`$ the cyclotron circle starts to drift in a way similar to that one shown in Fig. 3. After a time $`t_dd/v_d`$ a collision with another antidot takes place that finally originates a sequence of irregular jumps between antidots. The diffusion rate on the superlattice originated by this process can be estimated as $`D_{lat}d^2/t_ddv_d`$. This diffusion is important in the limit $`R_cd1`$. However, we note that even at $`E=0`$ at $`R_c>d/2`$ there are chaotic orbits which diffuse over the whole lattice as it was discussed in detail in  and this diffusion is dominant for $`d1`$.

It is interesting to understand how two electrons move in such a superlattice. Intuitively, one would expect that the Coulomb repulsion will separate electrons and they will not propagate together. In fact we found that it is not necessarily the case and there are regimes where two electrons propagate together. An example of such a case is shown in Fig. 6. In this case the electron pair moves with an average drift velocity $`v_dE/B`$ and the total displacement of the pair is about hundred times larger than the distance between the two electrons (Fig. 6).

The physical reason for appearance of such electron pairs is quite clear in the absence of superlattice potential. In this case, as it was discussed before at the beginning of Section III, electrons rotate one around another and propagate together, their dynamics is integrable. Then according to the KAM theorem a weak perturbation will not destroy such pairs. Indeed, in the case of Fig. 6 the antidot potential is relatively weak ($`U_0=1H/225`$) and the pair is not destroyed. We checked numerically that provided $`U_0H/2`$ the pair size starts to grow diffusively due to random scattering on a strong antidot potential, and eventually the pair is destroyed and electrons continue to propagate separately. For $`U_0H`$ the separation happens after a few collisions with antidots. On the contrary for $`U_0H`$ the life time of the classical pair becomes infinite in agreement with the KAM theorem.

## V  Conclusions

In this paper we investigated the effects of Coulomb interaction between two electrons on their classical dynamics in antidot superlattice in crossed electric and magnetic fields. We found that for weak electric field the electron pair can be trapped for a long time near an antidot even if eventually one or two electrons escape from the antidot. The escape rate $`\mathrm{\Gamma }`$ decreases proportionally to the square of electric field. This behaviour is qualitatively different from the case of noninteracting electrons which are trapped forever near the antidot in the limit of small electric field.

The study of the electron pair dynamics in the antidot superlattice showed that the Coulomb repulsion can create stable pairs propagating on a large distance. In agreement with the KAM theorem such pairs are stable when the antidot potential strength is relatively weak compared to the electron energy. On the contrary, in the opposite limit the pairs become unstable and electrons are quickly separated from one another. On the basis of this phenomenon it is possible to make a conjecture that in two dimensional heterostructures with high mobility the impurity potential is relatively weak and such electron KAM pairs will be stable and can be detected experimentally. We note that in the recent experiments  with 2D electron gas the carriers of charge $`2e`$ have been detected. It is possible that these carriers are related to the KAM pairs found in this paper.

We thank G.Casati who attracted our interest to the results found in .


text:
# Equilibrium crystal shapes in the Potts model

## I Introduction

If a binary mixture of fixed composition is brought into a coexistence state it will phase-separate into two pure phases, separated by an interface with a shape that minimizes the free energy of the system . The phase with the smaller volume will typically organize itself into a compact shape, known as the equilibrium shape. If the surface tension (excess free energy per unit interface area) is isotropic, this equilibrium shape will be a sphere. If either of the co-existing phases is crystalline, the anisotropic surface tension will lead to an aspheric equilibrium crystal shape (ECS). The ECS and the orientation dependence of the surface tension are intimately related. Once the surface tension is known for all orientations, the so-called Wulff construction  allows the generation of the ECS. The reverse is also possible, i.e., determining the surface tension as a function of orientation from the ECS; this is a procedure followed in this manuscript.

The ECS has been the subject of experimental studies . One striking feature the ECS may show is a roughening transition: the disappearance of facets (macroscopically flat surfaces) under rising temperature. In practice equilibrium shapes are rarely seen and they are hard to produce experimentally, but this does not mean they are without practical importance. For example the presence or absence of facets in the ECS influences growth properties, such as the speed of growth and the growth mode, even though growth shapes may differ strongly from the ECS. In general, the determination of the orientation dependence of the surface tension is a tough problem. Experimentally it is very difficult to measure, mainly due to equilibration problems. Numerical studies suffer from the same equilibration problems, and in the most common simulation techniques each surface orientation requires a separate simulation. Theoretical results only have been obtained for simplified models like the BCSOS model .

The prototypical model in which such properties are studied numerically, is the conserved-order-parameter Ising model with nearest-neighbor interactions, in which the total magnetization is kept constant. In the two-dimensional Ising model, the ECS is a square at zero temperature, but at any finite temperature the ECS looses all flat faces. It gradually changes into a circle when the temperature is approaching the critical temperature. The behavior is richer in the three-dimensional Ising model. At zero temperature, the ECS is a cube. At finite temperatures below the roughening temperature $`T_R`$, the ECS still has macroscopically flat faces but the corners as well as the edges are rounded. Above the roughening temperature, the ECS does no longer feature macroscopically flat faces. If the temperature is increased from the roughening temperature to the critical temperature, the ECS gradually approaches a sphere.

An extension of the Ising model is the Potts model, defined by the Hamiltonian

$$H=J\underset{i,j}{}\delta (\sigma _i,\sigma _j),$$
(1)
in which $`J`$ is the coupling constant, $`\delta `$ denotes the Kronecker delta function, and the summation runs over all pairs of nearest-neighbor sites, each having a spin with value $`\sigma =1\mathrm{}q`$. Note that the two-state Potts model is equivalent to the Ising model. The topic of this manuscript is to study how the ECS and related properties such as the roughening temperature behave with increasing number of states. Like in two dimensions the model undergoes a phase transition at a temperature we will refer to as the melting temperature $`T_m`$. Below this temperature the model has $`q`$ different phases, each of which is dominated by one of the $`q`$ possible spin values, whereas for temperatures above $`T_m`$ there is only a single phase in which on average all of the spin values are present in equal amounts. It is known that for $`q3`$ the melting transition is of first order , in contrast to the two-state (Ising) case where the transition is continuous. It is therefore not clear a priori whether there is a roughening transition when $`q3`$, nor in how far the ECS should approach a sphere when $`T`$ approaches $`T_R`$. Furthermore, one should expect the surface tension to approach a non-zero limit as the melting temperature is approached from below, whereas for continuous transitions it vanishes in this limit.

The manuscript is organized as follows. First, in section II we outline the numerical procedure to efficiently determine the ECS and show the resulting shapes for all numbers of states for which accurate numerical estimates of the melting temperature have been reported: $`q=2`$ (Ising), $`3,4`$ and $`6`$. In section III we present some theoretical considerations on the behavior of the model for large values of $`q`$. Section IV covers the measurements on the location of the roughening transition. In section V we present the orientation-dependent surface tension as extracted from the data on the ECS. We conclude with a discussion of the results and an outlook to further research.

## II Obtaining the ECS

Equilibrium shapes can be determined numerically for the $`q`$-states Potts model by forcing it into coexistence. This may be achieved by enforcing a constraint on the spin densities. The richest behavior is observed in case the maximum of $`q1`$ constraints are enforced, but mostly the resulting configurations are still determined by the orientation dependence of the surface tension, and can thus be obtained indirectly from numerical simulations with a single constraint enforced.

The constraint enforced in our simulations is the conservation of the density of state 1, while the densities of all other states are allowed to fluctuate freely. To implement this, we combined a density-conserving cluster algorithm, described in a recent article , with Glauber dynamics , constrained to never flip spins into or out of state 1. We chose the ratio of cluster updates and Glauber updates per site to be unity (i.e., a comparable amount of computational effort was spent to each). The run-times of the simulations varied between 400 000 and 800 000 Monte Carlo updates per site and we typically allowed the system to thermalize for 60 000 to 100 000 time steps. The lateral system size was 50 sites in all lattice directions, with periodic boundary conditions, and the fraction of conserved spins was typically $`0.25`$. In order to directly obtain the ECS we repeatedly took ”snapshots” of the evolving system, each time centering the cluster around the origin. Taking the autocorrelation time into account, we typically obtained 2000 independent equilibrium snapshots from one simulation. At the end we calculated the time-averaged density profile of the conserved component. Our results for the ECS were obtained from the 50% iso-density surface in this time-averaged density profile.

In figure 1, we show the equilibrium crystal shapes for several temperatures, for the case of $`q=2`$ (Ising), $`3`$, $`4`$, and $`6`$ states, respectively. These shapes are obtained from the 50% iso-density surface in the density profile, after full symmetrization, i.e. after we averaged over all 48 possible invariant mirror images of the cube.

The simulation temperatures are chosen as fractions of the estimates for the melting temperatures reported in Refs. ; these values are given in table I.

## III Behavior for large $`q`$

For large values of $`q`$ the thermal evolution of the ECS can be understood from the following arguments. First of all notice that at low temperatures the $`q`$ phases of the $`q`$-state Potts model without conservation laws will be very close to the zero temperature phases in which all spins on the lattice have equal value. In $`d`$ dimensions the free energy per site $`f`$ of these states therefore satisfies

$$f(T_{\mathrm{𝑙𝑜𝑤}})dJ.$$
(2)
The high temperature phase is completely dominated by the entropy resulting from the $`q`$ different occupation options for each site, with a resulting entropy per site of $`k_B\mathrm{ln}q`$ and free energy

$$f(T_{\mathrm{ℎ𝑖𝑔ℎ}})k_BT\mathrm{ln}q.$$
(3)
The melting temperature is obtained to a good approximation by equating these two expressions, with the result

$$\beta _mJ\frac{1}{d}\mathrm{ln}q$$
(4)
for large $`q`$, with $`\beta =1/k_BT`$. Note that for $`d=2`$ this is in agreement with the exact result  $`\beta _mJ=\mathrm{ln}(1+\sqrt{q})`$. For large $`q`$ one sees that $`Jk_BT`$ for temperatures below $`T_m`$. Therefore local excitations, typically consisting of single overturned spins, will be extremely rare in the low temperature phases, as they require an energy far exceeding the gain in free energy due to the entropy increase. Similarly the tendency to form clusters of equal spin in the high temperature phase is very weak, as the resulting gain of energy by far does not compensate the resulting loss of entropic free energy. As a result the melting transition for high $`q`$ will become very sharp.

Now we can also estimate the orientation dependence of the surface tension. In good approximation, an interface of 001-orientation (or symmetrically equivalent) is just a flat interface between two pure states, with a free energy of $`J`$ (one broken bond) per unit area. Similarly, by counting broken bonds, one finds that an interface of $`kl1`$-orientation (here $`k`$ and $`l`$ are not necessarily integers) has an energy per unit area of $`(1+k+l)/\sqrt{1+k^2+l^2}`$. For non-zero temperature there is an additional contribution from the entropy of the steps required to form such an interface. For not too large values of $`k`$ and $`l`$ this may be obtained by multiplying the step density needed to create a surface of orientation $`kl1`$ with the meandering entropy of such a step. The resulting expression is.

$$f(T)=\frac{(1+k+l)Jk_BT\mathrm{ln}\left[(k+l)^{k+l}/(k^kl^l)\right]}{\sqrt{1+k^2+l^2}}.$$
(5)
This result is independent of the value of $`q`$, especially it also holds for the Ising model. Excitations creating spins not belonging to the two phases that coexist at the interface are so rare that they may be neglected completely. For $`k/l`$ or $`l/k\mathrm{exp}(\beta J)`$ the entropy becomes dominated by thermal step fluctuations, but this concerns an extremely small range of orientations only (However, it does set the distance between facet edges). Blöte and Hilhorst  give expressions for the free energy from which the low temperature ECS may be obtained for all orientations.

The unimportance of all but the two coexisting phases implies that for large $`q`$ the 001-facets do not roughen, as $`T_m`$ is well below the roughening temperature of these facets in the Ising model (see Section IV). The equilibrium shape then remains nearly cubic for all temperatures up to $`T_m`$.

## IV Location of the roughening transition

Roughening transitions of crystal surfaces are characterized by the disappearance of macroscopically flat regions, or facets, in the ECS. Facets existing at low temperatures often disappear when the temperature increases, at a characteristic temperature known as the roughening temperature $`T_R`$ of the specific facet orientation. At this temperature the inverse radius of curvature $`(R_c)^1`$ of the crystal surface at the center of the facet jumps from zero (when a facet is present) to a non-zero value. Using renormalization-group calculations, Jayaprakash et al.  showed that the size of this jump satisfies the universal relationship

$$R_c=\frac{z_0kT_R\pi }{2\gamma _0},$$
(6)
where $`z_0`$ is the distance from the tangent plane at the facet to the center of mass of the crystal, and $`\gamma _0`$ is the surface tension of the 001-facet at $`T=T_R`$. The latter is not known exactly for any of the $`q`$-state Potts models in three dimensions. An approximation consisting of the ground state value plus the first correction term in a low temperature expansion is

$$\gamma _0=\beta J2e^{4\beta J}.$$
(7)
It will turn out that this approximation suffices for our goal: even for the smallest $`q`$-values the contribution of the correction term at $`T=T_R`$ is already smaller than the estimated systematic error in the curvature measurements.

If the center of mass of the ECS is placed at the origin, the centers of the facets that are present below the roughening temperature are located on the principal axes. To estimate the location of the roughening transition, we first measure the curvature at the six points where the ECS (obtained as in the previous section) cuts a principal axis, in the two principal directions tangential to the surface. The curvature was obtained by fitting a quadratic function $`y=y_0\frac{1}{2}y_1x^2`$; the fitted result for $`y_1/y_0`$ is the normalized inverse radius of curvature. Since we use a non-zero fitting range, the slope of the curve cannot become infinite and consequently we do not observe a jump in curvature; in fact we find quite smooth dependence of curvature on temperature, as a result of finite size, round-off and fluctuation effects. However, we made sure that our procedure to estimate $`T_R`$ is not very sensitive to this. Notably our estimate of the roughening temperature for the Ising model is close to previous estimates: Adler reported $`T_R/T_m=0.55\pm 0.02`$ , Mon et al. reported $`T_R/T_m=0.54\pm 0.02`$ , and Holzer and Wortis reported $`T_R/T_m=0.545\pm 0.004`$ .

The resulting measurements of the normalized inverse curvature are plotted in figure 2. The error bars indicate the statistical error; additional simulations for some points indicate that the systematic errors arising from the effects of thermalization, cluster shape deformations, long lived thermal excitations, etc., are larger. In the same plot, the curves described by Eq. (6) are also plotted. The intersection points of these two curves are estimates of the location of the roughening transition. The resulting values for the roughening temperature can be found in table I; the error indications are our estimates of the statistical and systematic errors combined.

In the previous section we noted that for large $`q`$ there is no roughening transition. Obviously $`q=6`$ is not large enough to observe this, but the increase of $`T_R/T_m`$ is very clear. It would be interesting to know the largest $`q`$-value for which a stable roughening transition does exist. The values of $`\beta J_R`$ quoted in table I combined with the estimate of $`T_m`$ of Eq. (4) suggest a value between 20 and 40. With the present method this would make simulations very slow, and it would also require an accurate measurement of the melting temperatures, since no literature values exist.

## V Orientation dependence of the surface tension

In general, the surface tension $`\gamma `$ can be written as $`\gamma =\gamma _0\gamma (\widehat{n})`$, where $`\gamma _0`$ is the surface tension of an interface oriented in the $`001`$, or symmetrically equivalent direction. In experiments with Pb crystals in equilibrium with their vapor, Heyraud and Métois used the inverse Wulff construction (see Ref. ) to obtain the angular part of the surface tension $`\gamma (\widehat{n})`$ as a function of orientation $`\widehat{n}`$ and the same method has been employed by Surnev et al. . Here we use it to obtain the orientation dependence of the surface tension for each ECS as obtained in section II.

The surface tension $`\gamma (\widehat{n})`$ is proportional to the distance from the center of the cluster to the tangent plane perpendicular to $`\widehat{n}`$, touching the iso-surface in the point $`\stackrel{}{X}`$. If the ECS is scaled such that the distance to the surface from the center of the shape along the lattice axes is unity, then

$$\gamma (\widehat{n})=\underset{I}{\mathrm{max}}\widehat{n}\stackrel{}{X},$$
(8)
where $`\stackrel{}{X}`$ is an element of the scaled iso-surface $`I`$. Thus, for every direction $`\widehat{n}`$, we have to find the point $`\stackrel{}{X}`$ in the iso-surface, for which the number $`\widehat{n}\stackrel{}{X}`$ is maximal. This procedure works in both two and three dimensions.

Using ECS’s obtained from fully symmetrized density profiles, we measured $`\gamma (\widehat{n})`$ along the arc in the $`11k`$-zone, which connects the $`001`$, $`111`$ and $`110`$ directions. Figure 3 shows the results.

These figures show that the angular dependence of the surface tension becomes nearly constant at increasing temperature, for all numbers of states $`q`$ considered here. Thus, for $`q6`$ the ECS approaches a sphere. For the Ising model ($`q=2`$) this was to be expected, since this model undergoes a continuous phase transition at the melting point. On the other hand, it is also clear that aspheric deviations become larger with increasing $`q`$. As we saw in Section III the ECS becomes nearly cubic for large $`q`$. Obviously $`q=6`$ still should be considered a small $`q`$-value in this context.

If the ECS is not faceted, the angular part of the surface tension will approach the point $`\mathrm{\Theta }=0`$ with zero slope. If however the ECS has $`001`$ facets, the approach to the point $`\mathrm{\Theta }=0`$ will occur with a non-zero slope, resulting in a cusp in the slope along the arc through this point . Looking for the temperature where these cusps first appear is an alternative way to measure the roughening temperature. We found this to be less accurate than the procedure described in section IV. The results were however consistent.

## VI Discussion and future research

We studied the Potts model with $`q=2`$, 3, 4, and 6 states, forced into coexistence by fixing the density of one state. The resulting ECS was studied as a function of temperature and number of states.

We found that the roughening transition, which is well known for the Ising model (equivalent to the two-states Potts model), persists for higher numbers of states, at least up to six states. The temperature $`T_R`$ at which this roughening transition occurs, measured as a fraction of the melting temperature, tends to increase with increasing number of states.

In the future, we want to study the more general and richer behavior of the ECS in case more than one quantity is conserved. For instance, in the three-state Potts model close to its melting point, with the constraint $`\rho _1\rho _2\rho _3`$, the ECS resembles the shape of two soap bubbles with a common interface; this changes under variations of temperature and the ratio $`\rho _2/\rho _3`$.

We also want to look at equilibrium shapes in constrained geometries, like fluids between parallel plates, or systems with grain boundaries.

Finally we are investigating the behavior for larger $`q`$-values with the aid of different Monte Carlo techniques.