text:
# Brownian Dynamics Studies on DNA Gel Electrophoresis. I. Numerical Method and Quasi-Periodic Behavior of Elongation-Contraction Motions

## I Introduction

Gel electrophoresis is a most widely used technique to separate DNA according to their length $`L`$ under an electric field. Viovy (2000) The constant field method is quite efficient for DNA of small sizes. It is well known, however, that the separation does not work well for DNA of sizes longer than a certain value which depends on the field and the “pore” diameter of the gel. Norden et al. (1991) The non-constant field methods such as the pulsed-field gel electrophoresis Schwartz and Cantor (1984) have been empirically invented and they have remarkably improved the shortcoming. Nevertheless, the understandings on these gel electrophoresis methods, even on the constant field method, still remain unsatisfactory from the view point of statistical physics of polymers.

By constant field gel electrophoresis (CFGE), DNA molecules with the radius of inertia, $`R_\mathrm{I}`$, smaller than the pore diameter of gel, $`a_{\mathrm{gel}}`$, are sieved by gel. Rodbard and Chrambach (1970); Chrambach and Rodbard (1971) In this sieving regime DNA with the smaller $`R_\mathrm{I}`$ migrate the faster through the gel. In the opposite case with $`R_\mathrm{I}a_{\mathrm{gel}}`$, DNA is considered to move through a “primitive chain” Doi and Edwards (1986) which is determined by stochastic movement of its front end in pores of the gel under the field. J.Lumpkin and Zimm (1982); Lumpkin et al. (1985); Slater and Noolandi (1986) This picture, denoted as the biased reptation model, successfully describes various features of CFGE, such as the $`L`$-dependence of the migration mobility, $`\mu `$, observed experimentally Stellwagen (1985) in sufficiently weak fields. More recently, the model of the biased reptation with fluctuations has been proposed A.N. Semenov and Viovy (1995) which well describes the observed field dependence of $`\mu `$ in moderate fields. Viovy (2000) When the field is further increased, the field-dependence of $`\mu `$ begins to be saturated. In this regime, interesting behavior of an individual DNA has been observed: “hernias”-like configurations by the computer simulations Deutsch (1988); Deutsch and Madden (1989) and “quasi-periodic” behavior, i.e., “elongation-contraction” motion in a cyclic way by the real-time fluorescent microscopy experiments. Masubuchi et al. (1993); Oana et al. (1994) Indeed, the latter experimental method combined with optical or magnetic tweezers for a molecule reveals that the elasticity of DNA is originated from the entropy associated with its static conformations. Smith et al. (1992); Perkins et al. (1995); Strick et al. (1996); Austin et al. (1997) However, dynamic properties of DNA under CFGE are not described in a unified way by this picture alone.

In order to get further insights into the problem mentioned above, we have developed a Brownian dynamics (BD) method which simulates electrophoresis of DNA in a 3 dimensional (3D) space. The method has been extensively applied to the CFGE, in particular, the elongation-contraction motion of the chain under CFGE, the results of which we will discuss in the present series of papers. In the present paper I our numerical method is explained and global aspects, such as the mean mobility $`\mu `$ and time evolutions of the velocity of the center of mass, $`v_\mathrm{G}(t)`$, and of the radius of the longer principal axis, $`R_l(t)`$, are investigated. In the accompanying paper II Azuma (2001) we analyze dynamics of ‘defects’ introduced by de Gennes de Gennes (1971) in a chain under CFGE. Finally in paper III Azuma and Takayama (2001) we will examine conformational changes of the coarse-grained chain under CFGE in details, and will argue that the field-parallel component of its elongation-contraction motion is essentially understood as a deterministic dynamics of an elastic string in a 1D space with an obstacle, thereby an origin of the elasticity will be shown to be the conformational entropy of the chain.

In the Brownian dynamics (BD) simulations extensively performed in the present series of studies, DNA is modeled as a chain of spherical electrolyte beads with a constant radius $`1/2`$ in a 3D space filled with a solvent. Each bead interacts with other beads of the chain and spatially-fixed beads constituting gel by hard-core repulsive interactions. It also interacts with its two bonded beads by a non-linear elastic force which keeps the distance between the two to be less than a threshold value which we put $`\sqrt{2}`$. To solve the Langevin equation of motion under these constraints as well as under the field and random forces, we employ a BD algorithm with the Lagrangian multiplier method. It is similar to the one by Deutsch and Madden, Deutsch (1988); Deutsch and Madden (1989) but is essentially different from their method in the following respects. In their method, the time pitch to solve the equation is adjusted so as to keep a fixed distance between bonded beads under the tensile forces which are evaluated by the Lagrangian multiplier method. In our method, on the other hand, the time pitch is fixed at a certain value, which we regard as the mean period of random forces. Within each time step we evaluate paths of all beads affected by the constraints, i.e., by hard-core interactions and extremely nonlinear elastic interactions between a pair of beads. These effects are reread as the supposed constraining forces with the appropriate Lagrangian multipliers in an iterative way. The resultant constraining forces thus specified yield a new set of bead positions without any violation to the constraints. To our knowledge the present BD analysis is the first extensive simulation of DNA gel electrophoresis in a 3D space based on the Langevin equation involving purely mechanical, microscopic forces which are supposed to act on monomers of DNA.

The purposes of the present paper I are to explain our BD method and to demonstrate its results on some global properties of DNA under CFGE. The resultant mobility $`\mu `$ as a function of the field turns out to agree well with those observed experimentally once a set of three parameters, $`l^{}`$, $`\tau ^{}`$ and $`E^{}`$, are properly specified. They relate respectively the scales of length, time and field strength of the simulation and the experiments.

In the BD simulation under relatively large fields quasi-periodic time evolution is observed in $`R_l(t)`$ (see Fig. 5 below), as well as in $`v_\mathrm{G}(t)`$, as has been observed experimentally. Masubuchi et al. (1993); Oana et al. (1994) Each local maximum of $`R_l(t)`$ at $`t_{\mathrm{max}}^n`$ with $`n=1,2,\mathrm{}`$ corresponds to a conformation that a chain just gets rid of trapping due to gel. The local minimum at $`t_{\mathrm{min}:\mathrm{l}}^n`$ preceding this maximum is then considered to be the instance that, being trapped by gel, the chain starts to elongate from a certain coiled conformation, while at time $`t_{\mathrm{min}:\mathrm{r}}^n`$ of the local minimum subsequent to $`t_{\mathrm{max}}^n`$ the chain retrieves another coiled state. We establish a method to properly classify a whole time sequence of $`R_l(t)`$ into peak parts, from $`t_{\mathrm{min}:\mathrm{l}}^n`$ to $`t_{\mathrm{min}:\mathrm{r}}^n`$ and the rest. We call the former parts deterministic ranges and the latter stochastic ones. We then investigate statistics of periods of the deterministic ranges, $`D_nt_{\mathrm{min}:\mathrm{r}}^nt_{\mathrm{min}:\mathrm{l}}^n`$, and those of the stochastic ranges, $`(1/\lambda )_nt_{\mathrm{min}:\mathrm{l}}^{n+1}t_{\mathrm{min}:\mathrm{r}}^n`$. It is found that {$`D_n`$} obey a Gaussian distribution with the mean $`D\overline{D_n}M`$ and the width $`\mathrm{\Delta }D`$, while {$`(1/\lambda )_n`$} do a Poisson distribution with the mean $`1/\lambda \overline{(1/\lambda )_n}M^0`$ for large $`M`$, where $`M`$ is the number of beads of the chain and the overline stands for the averages over $`n`$, i.e., peak structures. These results, combined with the observation that in each elongation-contraction motion the chain moves to the field direction by the distance proportional to $`M`$, yield $`\mu M^0`$ and explain why CFGE cannot separate DNA according to their sizes $`L(M)`$ for large $`L`$.

The organization of the present paper is as follows. In the next section we explain our BD method. The mobility $`\mu `$ of the chain obtained by the BD method is presented and compared with the experimental results in Sec. III. The method of analysis of peak structures in $`R_l(t)`$ is explained, and statistical nature of the deterministic and stochastic ranges is discussed in Sec. IV. The last section is devoted to further discussions of our results.

## II The Brownian Dynamics Method

A chain of spherical beads in a continuous 3D space, shown schematically in Fig. 1, is investigated by a new type of Brownian dynamics (BD) calculation. The hard core interaction is assumed in order for centers of any pair of beads not to come closer than 1 (the excluded volume effect), where the radius of beads is set $`1/2`$. Furthermore distances between neighboring beads, $`l`$, are restricted to $`l<\sqrt{2}`$. We suppose that this is realized by the extremely nonlinear elastic interaction between the neighboring beads: it is infinite for $`l>\sqrt{2}`$ and zero otherwise. These interactions assure the self-avoiding walk of the chain.

The BD method we adopt in the present work is to solve the following Langevin equation of motion for the chain,

$$\zeta \dot{𝒙}_i=q𝑬_\mathrm{b}+\underset{\alpha }{}\underset{j}{}𝑺_{ij}^\alpha +\underset{\alpha }{}\left[𝑭_i^\alpha 𝑭_{i1}^\alpha \right]+𝒇_i,$$
(1)
where $`𝒙_i`$ is the position of the center of the $`i`$-th bead, and $`\zeta `$, $`q`$, and $`𝑬_\mathrm{b}`$ are the viscosity, the charge of a bead and the (bare) external field, respectively. The second and the third terms in r.h.s. are the constraining forces required for the beads to satisfy the conditions,

$`𝒙_{ij}=`$ $`𝒙_j𝒙_i`$ $`>1`$ (2)
$`𝒍_i=`$ $`𝒙_{i+1}𝒙_i`$ $`<\sqrt{2}.`$ (3)
The vectors $`𝑺_{ij}^\alpha `$ and $`𝑭_i^\alpha `$ are respectively written as $`\theta _{ij}^\alpha [𝒙_i𝒙_j]`$ and $`\varphi _i^\alpha [𝒙_i𝒙_{i+1}]`$ by introducing the Lagrangian multipliers $`\theta _{ij}^\alpha `$ and $`\varphi _i^\alpha `$. The last term in eqn. 1 stands for the random force from a solvent. Its distribution is characterized by the following averages

$`𝒇_i`$ $`=`$ $`\mathrm{𝟎}`$ (4)
$`𝒇_i(t)𝒇_j(t^{})`$ $`=`$ $`2\zeta k_\mathrm{B}T\mathrm{\Delta }t\delta _{ij}\delta (tt^{}).`$ (5)
Here, $`k_\mathrm{B}`$, $`T`$, and $`\mathrm{\Delta }t`$ are the Boltzmann constant, the absolute temperature, and the average period of the random force, respectively.

To solve the one step evolution of eqn. 1, thereby the terms $`\{𝑺_{ij}^\alpha ,𝑭_i^\alpha \}`$ with the Lagrangian multipliers are determined, we examine the full Langevin equation with the inertia term and the forces arising from the hard core interaction and the nonlinear elastic interaction mentioned above. We solve the equation in a time interval of $`\mathrm{\Delta }t>t>0`$, and take the overdamped limit, $`m/\zeta 0`$, where $`m`$ is the supposed mass of the bead. The constraints of eqns. 2 and 3 are supposed to be satisfied by positions of the beads at $`t=0`$, denoted by $`\{𝒙_i^0\}`$. At this instance $`𝒇_i`$ and $`q𝑬_\mathrm{b}\mathrm{\Delta }t`$ are added as impulsive forces with the integrated strength $`\stackrel{~}{𝒇}_i`$ where $`\stackrel{~}{𝒇}_i=𝒇_i+q𝑬_\mathrm{b}\mathrm{\Delta }t`$. When the constraints are discarded, the solution in the overdamped limit is described by the relation between the displacement $`\mathrm{\Delta }𝒙_i`$ in this time interval and the force as

$$\mathrm{\Delta }𝒙_i=\frac{1}{\zeta }\stackrel{~}{𝒇}_i,$$
(6)
which is just what is obtained by simply integrating eqn. 1 without the second and the third terms.

In the case that positions of a pair of beads $`𝒙_i^{01}𝒙_i^0+\mathrm{\Delta }𝒙_i`$ and $`𝒙_j^{01}𝒙_j^0+\mathrm{\Delta }𝒙_j`$ violate the constraint of eqn. 2, we solve the full Langevin equation for this pair with the hard core scattering, which can be either elastic or inelastic. We then obtain, again in the overdamped limit, positions $`𝒙_i^1`$ and $`𝒙_j^1`$ at $`t=\mathrm{\Delta }t`$ as shown in Fig. 2a in case of the inelastic scattering. The result is reread as the consequence of the constraining force, $`𝑺_{ij}^1\theta _{ij}^1[𝒙_i^1𝒙_j^1]`$, which is determined by the conditions $`𝒙_i^1𝒙_i^0=(\stackrel{~}{𝒇}_i+𝑺_{ij}^1)/\zeta `$ and $`𝒙_j^1𝒙_j^0=(\stackrel{~}{𝒇}_j𝑺_{ij}^1)/\zeta `$ (note that $`𝑺_{ji}^1=𝑺_{ij}^1`$). This constraining force $`𝑺_{ij}^1`$ just corresponds to the one in the second term of eqn. 1. In the case $`𝒙_i^{01}`$ and $`𝒙_{i+1}^{01}`$ violate eqn. 3, the nonlinear elastic force is supposed to reduce the distance $`𝒍_i`$ to $`\sqrt{2}`$. This effect is represented by the constraining force $`𝑭_i^1(=𝑭_{i+1}^1)`$, which is added to the third term of eqn. 1. Its Lagrangian multiplier $`\varphi _i^1`$ is determined in such a way that $`𝒙_i^1𝒙_{i+1}^1=\sqrt{2}`$ with $`𝒙_i^1=𝒙_i^0+𝑭_i^1\mathrm{\Delta }t`$ and $`𝒙_{i+1}^1=𝒙_{i+1}^0+𝑭_{i+1}^1\mathrm{\Delta }t`$ (see Fig. 2b). If $`𝒙_i^{01}`$ involves violation to more than one constraints of eqns. 2 and/or 3, the above-mentioned procedures are applied to evaluate the constraining forces for each pair of beads with the common $`𝒙_i^{01}`$. If, on the other hand, $`𝒙_i^{01}`$ does not involve any violation to eqns. 2 and 3, we put $`𝒙_i^1=𝒙_i^{01}`$.

In the configuration $`\{𝒙_i^1\}`$, it is generally expected that new pairs of beads may violate the constraints since we have tried to remedy the violations in $`\{𝒙_i^{01}\}`$ in a pairwise way. Then we repeat the above-mentioned procedure, in which $`\stackrel{~}{𝒇}_i`$ is replaced by $`\stackrel{~}{𝒇}_i+𝑪_i^1`$, where $`𝑪_i^1=_j𝑺_{ij}^1+𝑭_i^1𝑭_{i1}^1`$ is the constraining force determined by the first procedure. If a configuration without any violation to the constraints is obtained by the $`\alpha _{\mathrm{con}}`$ times repetition of this procedure, we regard $`\{𝒙_i^{\alpha _{\mathrm{con}}}\}`$ as the solution of eqn. 1 evolved from $`\{𝒙_i^0\}`$ by one unit of time $`\mathrm{\Delta }t`$. In this manner the multi-scattering processes in the interval of $`\mathrm{\Delta }t`$, including those due to the nonlinear elastic interaction, can be taken into account, and the corresponding forces $`\{𝑺_{ij}^\alpha ,𝑭_i^\alpha \}`$ in eqn. 1 are determined.

Although our BD method is based on the overdamped Langevin equation with the Lagrangian multipliers, it distinctly differs from the one due to Deutsch and Madden. Deutsch (1988); Deutsch and Madden (1989) In the latter the equidistant condition on all the pairs of beads is imposed, and the constraining tensile forces associated with it are evaluated at the instance when random forces are applied. Then the time interval $`\mathrm{\Delta }t`$ is adjusted so that the configuration at $`t+\mathrm{\Delta }t`$ in fact obeys the equidistant condition within accuracy of $`0.1\%`$. In contrast, our BD method takes into account a similar, but looser, constraint on the bond length as well as the excluded volume effect of the beads by explicitly examining the supposed multiple scattering within the interval $`\mathrm{\Delta }t`$ which is fixed. Therefore there is no finer time unit other than $`\mathrm{\Delta }t`$ in solutions of our BD method. It is the mean period of random forces $`\{𝒇_i\}`$.

In the present work we simulate dynamics of a chain in a continuous 3D space confined into a cubic box with a volume $`360^3`$. The periodic boundary condition is imposed to each direction of the box. The gel is represented by a network of immobile bars arranged so that they form a simple cubic lattice. Each bar consists of tightly connected beads with a radius $`1/2`$, and its interaction with a bead of the mobile chain is computed according to the constraint of eqn. 2. We employed three lattices of gel with different lattice distances $`a_{\mathrm{gel}}=10`$, $`17`$, and $`20`$, which are considered to correspond to the ’pore’ diameter of agarose gel. Norden et al. (1991) The lengths of chains, $`M`$, used in the BD calculation are $`M=30`$, $`40`$, $`60`$, $`80`$, $`160`$, $`240`$, and $`320`$. We have performed typically $`64`$ runs whose initial configurations are given independently. The direction of the field is fixed at $`(1,1,1)`$. The time is measured in unit of $`\mathrm{\Delta }t`$, i.e., 1 BD step.

Here we want to emphasize that the present BD method correctly reproduces the expected static and dynamic properties of a real (not phantom) polymer in a 3D space under the vanishing field: $`R_IM^{\nu _\mathrm{F}}`$ with $`\nu _\mathrm{F}0.59\pm 0.01`$ and the Einstein relation $`D_G\mu /M`$ are confirmed (see Fig. 3 in Sec. III). Here, $`\nu _\mathrm{F}`$, $`R_I`$, and $`D_G`$ are the Flory exponent, the radius of inertia, and the diffusion constant of the center of mass, respectively.

The number of the multiple scattering, $`\alpha _{\mathrm{con}}`$ above mentioned, is expected to be at most of the order of $`M^2`$. It has turned out that $`\alpha _{\mathrm{con}}`$ becomes relatively larger when the chain is in an extended configuration under gel electrophoresis. We quote here typical figures for a chain with $`M=80`$: $`\alpha _{\mathrm{con}}4\times 10^1`$, $`5\times 10^2`$ for $`E=0`$, $`0.032`$, respectively.

## III Mobility

The most basic quantity that has been commonly measured in CFGE of DNA is the long time average of velocity divided by the electric field, or the mobility $`\mu `$. In the present work, we have implemented long time BD simulations with a highly efficient algorithm using a massive parallel computer and have eventually obtained $`\mu `$, hereafter denoted by $`\mu _{\mathrm{BD}}`$, with enough statistical precision for all the set of parameters investigated. Here we only mention typical figures of computation to get $`\mu _{\mathrm{BD}}`$ of a $`M=240`$ chain for each run: more than 150 hours are consumed on a single processor of SGI2800 system to pass $`10^6`$ BD steps. In the present work the following values of the parameters are chosen: $`\mathrm{\Delta }t=1`$, $`\zeta =10`$ and $`k_\mathrm{B}T=2`$. The results are shown in Fig. 3 with the experimental data due to Heller et alC. Heller and Viovy (1994) Here and hereafter the value of field $`E`$ stands for $`qE_\mathrm{b}`$ of eqn. 1.

The simulated $`\mu _{\mathrm{BD}}`$ at $`E=0`$ in Fig. 3 is evaluated through the Einstein relation

$$\frac{\mu _{\mathrm{BD}}^0}{M}=D_\mathrm{G}\underset{t\mathrm{}}{lim}\frac{1}{6t}(𝑹_\mathrm{G}(t)𝑹_\mathrm{G}(0))^2|_{E=0},$$
(7)
where $`𝑹_\mathrm{G}(t)`$ is the position of the center of mass of the chain. The diffusion constant $`D_\mathrm{G}`$ in the r.h.s. is evaluated by the BD simulation in a vanishing field in the time range $`t=(510)\times 10^6`$. The resultant $`\mu _{\mathrm{BD}}^0`$ turns out to be compatible with the value of $`lim_{E0}\mu _{\mathrm{BD}}`$ as seen in the figure. This confirms that the Einstein relation in fact holds in the present BD simulation for the chain.

In order to compare the BD results with the experimentally observed data, we have to specify proper conversion of values of the parameters in the simulation and the experiment. Let us first consider the length scale of a chain. In the accompanying paper II we examine in detail the embodiment of ‘defects’ due to de Gennes de Gennes (1971) in our BD chain and obtain the mean distance $`n`$ of adjacent ’defects’ along the chain under $`E=0`$ as $`n5.7`$. This value, in unit of the mean distance between neighboring beads, $`l`$, is almost independent of the set of parameters $`a_{\mathrm{gel}}`$ and $`M`$ we have examined. Regarding $`n`$ as an estimate of the ‘persistent length’ in the BD simulation, $`𝒫_{\mathrm{BD}}`$, we put $`𝒫_{\mathrm{BD}}=cn`$, where $`c`$ is an adjustable parameter whose value is nearly equal to unity (note that $`\sqrt{2}l1`$). The conversion factor of the length-scale $`l^{}`$ is then determined by

$$𝒫_{\mathrm{exp}}=l^{}𝒫_{\mathrm{BD}}=cl^{}n,$$
(8)
where $`𝒫_{\mathrm{exp}}`$ is the persistence length experimentally measured. With $`𝒫_{\mathrm{exp}}`$ 60nm at $`T_{\mathrm{exp}}`$ 300K taken from Ref. Volkmuth et al., 1994, we obtain $`cl^{}`$ 10nm.

Another length scale of importance is the mean pore size of gel, which is rather difficult to be accurately specified even in the experiments. Viovy (2000) Here we simply suppose that DNA in $`1\%`$ agarose corresponds to a BD chain in the gel with $`a_{\mathrm{gel}}=17`$. With $`𝒫_{\mathrm{BD}}`$ above fixed, this value of $`a_{\mathrm{gel}}`$ is about three times larger than that of $`𝒫_{\mathrm{BD}}`$.

Next let us relate $`E`$ in the BD simulation to the experimental field, $`E_{\mathrm{exp}}`$, in V/cm. For this purpose we impose that the ratios $`\mathrm{\Delta }\epsilon /k_\mathrm{B}T`$ in the BD calculation and the experiment be identical, where $`\mathrm{\Delta }\epsilon `$ is the energy that a part of chain of a length equal to the persistence length $`𝒫`$ gains when it moves a distance of $`𝒫`$ under the electric field:

$$\frac{\sigma E_{\mathrm{exp}}𝒫_{\mathrm{exp}}^2}{k_\mathrm{B}T_{\mathrm{exp}}}=\frac{E𝒫_{\mathrm{BD}}^2}{k_\mathrm{B}T_{\mathrm{BD}}},$$
(9)
or

$$E_{\mathrm{exp}}=E^{}E;E^{}=\frac{k_\mathrm{B}T_{\mathrm{exp}}𝒫_{\mathrm{BD}}^2}{2\sigma 𝒫_{\mathrm{exp}}^2}5\times 10^2c^2\mathrm{V}/\mathrm{cm}.$$
(10)
where $`\sigma `$ is the density of charge of DNA, and $`k_\mathrm{B}T_{\mathrm{BD}}=2`$ has been used. Quoting further the experimental value of $`\sigma `$ 0.6e<sup>-</sup>/Å, Volkmuth et al. (1994) as well as $`cl^{}`$ 10nm at $`T_{\mathrm{exp}}`$ 300K in eqn. 8, we obtain the last expression for $`E^{}`$ in the above equation. As for the ordinates of Fig. 3 we may relate $`\mu _{\mathrm{BD}}`$ to $`\mu _{\mathrm{exp}}[\mathrm{cm}^2/\mathrm{V}\mathrm{s}]`$ by

$$\mu _{\mathrm{exp}}=\mu ^{}\mu _{\mathrm{BD}};\mu ^{}=\frac{l^{}}{\tau ^{}E^{}}2\times 10^9\frac{c}{\tau ^{}}[\mathrm{cm}^2/\mathrm{V}\mathrm{s}],$$
(11)
where $`\tau ^{}`$\[s\] is the real time which corresponds to one discretized time step $`\mathrm{\Delta }t`$ in the BD calculation with $`\zeta =10`$ (or $`\mathrm{\Delta }t/\zeta =0.1`$).

The two sets of data of $`\mu `$ in Fig. 3 are drawn as follows. The data $`\mu _{\mathrm{exp}}`$ are plotted directly using the units indicated on the right ordinate and the upper abscissa, whereas the units for $`\mu _{\mathrm{BD}}`$ are indicated on the left ordinate and the lower abscissa with $`c^2=1.60`$ and $`1.26`$ respectively for $`M=240`$ and $`320`$. The scales of the two ordinates are related by means of $`\tau ^{}=1.4\times 10^6`$. We see from the figure that the field dependence of $`\mu _{\mathrm{BD}}`$, particularly that of $`M=240`$, well coincide with that of $`\mu _{\mathrm{exp}}`$ of DNA with the corresponding length we have specified, i.e., $`4.3`$ and $`6.5`$ kbp DNA for $`M=240`$ and $`320`$, respectively. This is quite satisfactory if we take into account the fact that only three adjustable parameters are used to draw the two sets of data in Fig. 3. We therefore consider that the present BD model catches up even quantitatively essences of DNA dynamics under CFGE.

When we compare $`\mu _{\mathrm{BD}}`$ with $`\mu _{\mathrm{exp}}`$ further in detail, however, there are some unsatisfactory features in the BD results. As compared with $`\mu _{\mathrm{exp}}`$ of $`6.5`$ kbp DNA, an $`E`$-independent branch of $`\mu _{\mathrm{BD}}`$ has not been ascertained in the weak field limit for the $`M=320`$ chain. One possible reason for this may be the insufficient time of simulation (up to $`t=10^6`$ in finite fields) for this length of chain in weak fields. Thus, within the limited CPU time, it is a rather hard task for the present BD simulation to judge the field dependence of $`\mu _{\mathrm{BD}}`$, proportional to $`E`$ or $`E^2`$, in weak fields. Hervet and Bean (1987); Holmes and Stellwagen (1990); T. Duke and Semenov (1994) In the strong field regime, on the other hand, $`\mu _{\mathrm{BD}}`$ tend to saturate, while, $`\mu _{\mathrm{exp}}`$ are still significantly increasing with $`E`$. These different tendencies are seen more clearly if we compare $`\mu _{\mathrm{exp}}`$ measured in $`10`$V/cm with the corresponding $`\mu _{\mathrm{BD}}`$ (not seen). The reason of this discrepancy is not clear at the moment.

We have also examined $`\mu _{\mathrm{BD}}`$ of the chains in gels with $`a_{\mathrm{gel}}=20`$ and $`10`$. If only the volume ratio of gel is considered, they correspond roughly to $`0.75\%`$ and $`3\%`$ agaroses, respectively. Qualitatively, $`\mu _{\mathrm{BD}}`$ becomes larger and its $`E`$-dependence becomes weaker for the larger $`a_{\mathrm{gel}}`$ as expected, but these tendencies are quantitatively weaker than those observed experimentally. The origins of the discrepancy may be attributed to our model for gel, i.e., a perfectly rigid, regular jungle gym.

## IV Elongation-Contraction Motion

### IV.1 Motion of an individual chain

In fields larger than a certain value depending on $`M`$ and $`a_{\mathrm{gel}}`$, chains in the present 3D BD simulation exhibit quasi-periodic behavior, i.e., they exhibit elongated and contracted shapes alternatively, as observed in the experiments Masubuchi et al. (1993); Oana et al. (1994); Larsson and Åkerman (1995) and in the previous 2D BD models. Deutsch (1988); Matsumoto and Doi (1994); Masubuchi et al. (1995) It is here emphasized that such quasi-periodic behavior has been observed within the field range where we have discussed the $`E`$-dependence of $`\mu `$ in Fig. 3. We show a typical conformational change of a chain observed in one MD run of CFGE with the parameters $`M=240`$, $`E=0.032`$, and $`a_{\mathrm{gel}}=17`$ in Fig. 4. Within a time window of the figure the elongation-contraction motion occurs three times at around $`t=`$2, 6 and 11 $`\times 10^5`$. At ranges between them centered around $`t=4`$ and 8 $`\times 10^5`$ a chain is in a rather compact form. In this range it can happen, though not so frequently, that the front and rear ends exchange their roles as seen around $`t=9\times 10^5`$ in the figure. It is also pointed out here that the front end of a chain moves in the field direction with almost a constant rate.

### IV.2 Quasi-periodic behavior of $`R_l(t)`$ and $`v_\mathrm{G}(t)`$

In Fig. 5 we show the time evolution of the radius of longer principal axis, $`R_l(t)`$, and the velocity of the center of mass, $`v_\mathrm{G}(t)`$, in the MD run whose chain conformational change is shown in Fig. 4. The fluctuation in $`R_l(t)`$ is notable, showing successive $`\mathrm{\Lambda }`$-shaped peaks. We also note that just before the maxima of $`R_l(t)`$, minima are observed in $`v_\mathrm{G}(t)`$. These features are qualitatively in good agreement with those observed in the experiment. Oana et al. (1994)

The time evolution of $`R_l(t)`$ and $`v_\mathrm{G}(t)`$ in Fig. 5 looks quasi-periodic. Actually it was found by the experiment that the autocorrelation function of $`R_l(t)`$, $`C_{\mathrm{RR}}(t)`$, exhibits a damped oscillation. Oana et al. (1994) This is also the case for our present BD data. In Fig. 6 $`C_{\mathrm{RR}}(t)`$ evaluated for various $`M`$ are presented. For $`M80`$, $`C_{\mathrm{RR}}(t)`$ exhibit the first undershoot below the line of $`C_{\mathrm{RR}}(t)=0`$, and even the first overshoot is clearly seen for $`M160`$. From these results we may conclude that, in chains with $`M160`$ under $`E=0.032`$, elongation-contraction motions occur quasi-periodically. In the inset of Fig. 6 we show nearly linear dependence of the period of oscillation $`\tau =4t_0`$ on $`M`$, where $`t_0`$ is defined as the time of the first intercept of $`C_{\mathrm{RR}}(t)`$ with the abscissa. Oana et al. (1994)

### IV.3 Analyses on peak structures of $`R_l(t)`$

Oana et al. Oana et al. (1994) argued that the steady-state time evolution of DNA under CFGE is classified into two types of behavior. One is the elongation-contraction motion which looks apparently deterministic. We call the time branches in which such a motion is undergoing the deterministic ones (ranges centered at around $`t=`$2, 6 and 11 $`\times 10^5`$ in Fig. 4). The other is what is observed between two deterministic branches. We call them stochastic branches. Oana et al. approximately described the quasi-periodic behavior of $`R_l(t)`$ in terms of a simplified function: peaks in $`R_l(t)`$ are approximated by $`\mathrm{\Lambda }`$-shape branches whose widths are assumed to obey a Gaussian distribution, and others by constants whose duration times are assumed to obey a Poisson distribution. They calculated $`C_{\mathrm{RR}}(t)`$ by means of this model function for various cases and remarked that when $`\lambda D`$ becomes larger than unity a damped oscillation in $`C_{\mathrm{RR}}(t)`$ becomes prominent. Here $`1/\lambda `$ and $`D`$ stand for the mean period of the stochastic branches and the mean width of the deterministic ones, respectively.

Following the classification due to Oana et al., we analyze our BD results by making use of an algorithm which appropriately picks up peaks of $`R_l(t)`$. For this purpose we first pick up time sets $`\{t_{\mathrm{min}:\mathrm{l}}^i,t_{\mathrm{max}}^i,t_{\mathrm{min}:\mathrm{r}}^i\}`$, i.e., the times of preceding local minimum of the $`i`$-th maximum, the $`i`$-th maximum itself, and the subsequent minimum. Then we select out peaks whose height relative to their adjacent minima is larger than a certain threshold value. The latter is chosen in such a way that the sum of periods of peak branches selected surmounts more than $`50\%`$ of the total duration of the observation. This procedure divides a whole time sequence of $`R_l(t)`$ into the two types of branches whose numbers are equal to each other. It picks up ensemble of various peaks but not a limited numbers of peaks of special shapes. The number of peaks thus selected is more than a hundred for each set of parameters.

### IV.4 Statistics of deterministic and stochastic branches

For the set of times $`\{t_{\mathrm{min}:\mathrm{l}}^n,t_{\mathrm{max}}^n,t_{\mathrm{min}:\mathrm{r}}^n\}`$ of peaks picked up by the procedure described above we examine the distribution of peak widths (or periods of deterministic branch) $`D_nt_{\mathrm{min}:\mathrm{r}}^nt_{\mathrm{min}:\mathrm{l}}^n`$ and that of duration times between neighboring peaks (periods of stochastic branch) $`1/\lambda _nt_{\mathrm{min}:\mathrm{r}}^nt_{\mathrm{min}:\mathrm{l}}^{n1}`$. The results for $`M=240,E=0.032`$ and $`a_{\mathrm{gel}}=17`$ are shown in Fig. 7. They numerically confirm the basic assumptions of Oana et al., i.e., $`\{D_n\}`$ nearly obey a Gaussian distribution, while $`\{1/\lambda _n\}`$ do a Poisson distribution. The mean and the variance of the former are given by $`D=3.3\times 10^5`$ and $`\mathrm{\Delta }0.2D`$, respectively, while the mean width of the latter is given by $`\lambda D1.7`$. The sum $`D+1/\lambda `$ should be equal to $`\tau =4t_0`$ evaluated from $`C_{\mathrm{RR}}(t)`$ in Sec. IV.2, which is verified within our numerical accuracy. The ratio $`\lambda D1.6`$ fulfills the condition for a damped oscillation to be observed in $`C_{\mathrm{RR}}(t)`$.

The ratio $`\lambda D`$ is found to decrease below unity for $`M`$ which is smaller than a certain value. In Fig. 8, $`D`$ and $`1/\lambda `$ are plotted against $`M`$ for $`E=0.032`$ and $`a_{\mathrm{gel}}=17`$. The results strongly suggest that $`DM^1`$ and $`1/\lambda M^0`$ for large $`M`$. In Fig. 9, the $`M`$-dependence of ratios $`\lambda D`$ and $`\mathrm{\Delta }/D`$ are presented. As expected from the results shown in Fig. 8, $`\lambda D`$ is an approximately linearly increasing function for large $`M`$. On the other hand, $`\mathrm{\Delta }/D`$ is almost constant ($`0.3`$), indicating that the Gaussian distribution of $`\{D_n\}`$ is scaled by the mean $`D`$ alone. It should be noted that, for $`M>100`$ where $`\lambda D>1`$ holds, $`C_{\mathrm{RR}}(t)`$ clearly exhibits a damped oscillation as has been observed in Fig. 6. Also the result $`4t_0M`$ for $`M>100`$ shown in the inset of Fig. 6 is in accordance with the fact that $`\lambda D`$ is larger than unity for such $`M`$. For sufficiently large $`M`$ the deterministic branches dominate a whole sequence of the time evolution, and so $`4t_0DM`$ holds. These results consistently indicate that there exists a crossover between the regimes with and without elongation-contraction motions at around $`M100`$ in CFGE with the parameters studied, i.e., $`E=0.032`$ and $`a_{\mathrm{gel}}=17`$.

We have also examined the field dependence of $`D`$ and $`1/\lambda `$ for a fixed $`M(=240)`$, the results of which are shown in Fig. 10. In strong fields $`E>0.06`$, $`DE`$ becomes to be saturated, while $`E/\lambda `$ decreases with increasing $`E`$. As a consequence, $`D\lambda `$ is an increasing function of $`E`$, indicating that the deterministic branches, or in other words, the elongation-contraction motions survive and even dominate a whole evolution of the chain in sufficiently strong fields.

Lastly, in Fig. 11, the averaged time evolution of individual peaks of $`R_l(t)`$ are plotted for various $`M`$ and $`E`$. The abscissa and ordinate are normalized by the period $`\tau 4t_0`$ and the averaged peak height, respectively. It is seen that with these normalizations, the data with different $`M`$ and $`E`$ almost lie on a universal curve during the period from $`0.6\tau `$ to $`0.4\tau `$. The results obtained above such as $`DM,E^1`$, $`\mathrm{\Delta }D(M)`$ can be derived from the scaling behavior of $`R_l(t)`$ in Fig. 11. Such scaling behavior of $`R_l(t)`$ including the fact that the normalized peak exhibits an anti-symmetric shape is just what was observed experimentally. Oana et al. (1994) It was mentioned in Ref. Oana et al., 1994 that the ratio of slope of the $`\mathrm{\Lambda }`$-shape immediately after the maximum to that before it is approximately constant ($`3`$) for all conditions they employed. The corresponding ratio of the BD result is $`2`$. The origin of this quantitative difference is not clear at the moment.

## V Discussions

We have performed extensive simulations on the DNA constant-field gel electrophoresis (CFGE) by means of a new Brownian dynamics (BD) method which we have recently developed. Our BD method is based on a generalized bead-spring model for a polymer in a 3D space with rigid, regularly arrayed obstacles. By the word ‘generalized’ we mean that the ‘spring force’ is not linear but extremely nonlinear: it is infinite when the distance $`l`$ between two beads connected by the ‘spring’ is larger than $`\sqrt{2}`$, while it is zero when $`l<\sqrt{2}`$. Also in the model we strictly take into account of the excluded volume effects between any pair of beads whose diameter is put unity. These forces are supposed to be of a microscopic origin, i.e., molecular forces acting on monomers of a polymer. It should be emphasized that there involves no characteristic scale in these interactions except for the maximum distance between the connected beads ($`\sqrt{2}`$) relative to the bead diameter (unity).

In a vanishing field the viscosity $`\zeta `$ (or $`\mathrm{\Delta }t/\zeta `$ with $`\mathrm{\Delta }t`$ being the time unit of the BD dynamics) and the temperature $`T`$ in eqn. 1 give rise to a characteristic scale of the model. The two parameters specify the viscous force and the random force from a solvent and are related to each other by eqn. 5. As is discussed in the accompanying paper II, the characteristic scale is the ‘persistence length’ of the chain. Actually, with $`\zeta =10,\mathrm{\Delta }t=1`$ and $`k_\mathrm{B}T_{\mathrm{BD}}=2`$ the mean distance between two connected beads ($`cl^{}`$ in eqn. 8) turns out to be about one six-th of the persistence length. Equating the latter to the persistence length of real DNA, the mean bead distance is roughly estimated as 10nm. Thus a bead in our BD model represents a potion of DNA consisting of a few tens of base-pairs. We may call our BD model with these values of parameters a bead-spring model of a semi-microscopic level. We will argue in paper III that the resultant chain dynamics, when it is coarse-grained in time and space, can be interpreted as dynamics of a charged, elastic string, and that the elasticity is due to the conformational entropy of the original, semi-microscopic chain.

A value of the electric field, $`E=qE_b`$ in eqn. 1, is reasonably converted to that of the experimental CFGE on DNA as discussed in Sec. III. By the procedure adopted there the conversion of the field scale is essentially governed by that of the length scale as expressed in eqn. 10. The parameter $`c`$ in eqn. 10 has been adjusted for the chains with different sizes, but within the range that the results are compatible with the condition $`\sqrt{2}l1`$ where $`l`$ is the mean distance of connected beads. Then the $`E`$-dependence of mobility $`\mu _{\mathrm{BD}}`$ has been shown to agree semi-quantitatively with that of $`\mu _{\mathrm{exp}}`$ once absolute magnitudes of the simulational and experimental $`\mu `$’s are adjusted at the vanishing field limit. This procedure fixes the conversion of the time scales of the simulation and the experiment. To our knowledge, such a semi-quantitative coincidence between the simulation and the experiment has not been so far achieved.

We have also demonstrated that our BD simulation on CFGE in fields stronger than a certain crossover value reproduces the quasi-periodic evolution of a chain whose characteristic aspects are quite similar to those experimentally observed. Actually, as proposed by Oana et al. Oana et al. (1994) based on their experimental observations, a whole sequence of time evolution of the chain can be divided into two branches according to behavior of the longer radius of the chain $`R_l(t)`$. One is the deterministic branch where $`R_l(t)`$ exhibits a distinct peak structure which corresponds to the elongation-contraction motion of the chain, and the other is the stochastic branch between two neighboring deterministic branches where the chain exhibits a rather compact form.

We have studied statistics of time durations of the two branches and have found that, as assumed by Oana et al.Oana et al. (1994) those of the deterministic branch obey a Gaussian distribution with mean $`D`$ and variance $`\mathrm{\Delta }`$, while those of the stochastic branch do a Poisson distribution with the mean $`1/\lambda `$. Furthermore it is found that, in a fixed field, $`DM^1,\mathrm{\Delta }D`$ and $`1/\lambda M^0`$ for large $`M`$. This means that the deterministic branches become more dominant in a whole time sequence, and that the quasi-periodicity in the chain dynamics becomes more distinct (see Fig. 6) for larger $`M`$. These results, combined with the observation that the chain moves to the field direction by the distance proportional to $`M`$ in each deterministic branch, yield $`\mu M^0`$. This explains why CFGE cannot separate DNA according to their size $`L(M)`$ for large $`L`$. We have also checked that, for a fixed $`M`$ which is sufficiently large, the deterministic branches become dominant as $`E`$ is increased at least within the field range of our BD simulation (see Fig. 10).

By the present BD simulation a crossover between the regimes with and without elongation-contraction motions is expected to occur for a chain with $`M100`$ in $`E=0.032`$ and $`a_{\mathrm{gel}}=17`$. According to the length and field conversions described in Sec. III, this roughly corresponds to DNA of 2kbp in $`1\%`$ agarose in $`E_{\mathrm{exp}}4\mathrm{V}/\mathrm{cm}`$. Experimently, on the other hand, the direct measurement of the elongation-contraction motion by the fluorescent microscopy has been rather limited; the shortest is for T7 DNA (38kbp), Larsson and Åkerman (1995) but mostly for T4 DNA (166kbp) Oana et al. (1994) and T2 DNA (164kbp). Larsson and Åkerman (1995) But, there exist another experiment which strongly indicates the occurrence of the elongation-contraction motion, i.e., the antiresonance of mobility in the filed inversion gel electrophoresis. In such a experiment the antiresonance has been clearly observed for DNA with 9.42kbp. T. Kobayashi and Ogawa (1990) Therefore, our criterion mentioned above is considered to be a reasonable estimate for the elongation-contraction motion to occur.

Lastly two further comments are in order. In the whole parameter ranges we have examined, including the ones where the deterministic branches sufficiently dominate the stochastic ones, so-called ’hernias’-like configurations of a chain has been scarcely observed in the present BD simulation. This implies that they do not play a significant role in determining the saturation of $`\mu `$ at least for chains with moderate sizes we have studied. Another interesting observation has been already pointed out in Sec. IV.2: the front end of a chain moves in the field direction with almost a constant rate (Fig. 4). In relation to this, we note that the average of $`v_\mathrm{G}(t)`$ over the deterministic branches and that over the stochastic branches are found to almost coincide with each other in a whole range of field we have examined. We shall interpret these results in paper III where the details of the chain dynamics in the deterministic branches is discussed.

###### Acknowledgements.

The authors wish to thank M. Doi and Y. Masubuchi for useful discussions on their experimental and simulational results. The computation in this work has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo, and those of the Computer Center of University of Tokyo.

## References

FIG. 1

FIG. 2

FIG. 3

FIG. 4

FIG. 5

FIG. 6

FIG. 7

FIG. 8

FIG. 9

FIG. 10

FIG. 11


format: text

text:
# Theory for high spin systems with orbital degeneracy

\[

## Abstract

High-spin systems with orbital degeneracy are studied in the large spin limit. In the absence of Hund’s coupling, the classical spin model is mapped onto disconnected orbital systems with spins up and down, respectively. The ground state of the isotropic model is an orbital valence bond state where each bond is an orbital singlet with parallel spins, and neighbouring bonds interact antiferromagnetically. The possible relevance to the transition metal oxides are discussed.

\]

In many transition metal oxides, there is an orbital degeneracy in the electron occupation energy, resulting in rich and novel magnetic phenomena.  The orbital ordering and orbital density wave have been observed experimentally in a family of manganites. In these systems, the spin coupling depends on the electron’s orbital occupations. The simplest model to describe spin-$`\frac{1}{2}`$ systems with two-fold orbital degeneracy is the SU(4) model where the Hund’s rule coupling and the orbital anisotropy are neglected. Some transition metal oxides, such as manganites and lanthanum vanadium oxides, have higher spins, which has attracted much theoretical interests.  Mathematical models for high spin systems are generally complicated. The Hund’s rule coupling favors ferromagnetic (FM) spins, and the anisotropy of the electron hopping integrals breaks orbital SU(2) symmetry. While the phase diagrams of these complex spin systems are rich, there is lack of a tractable method to systematically study these systems. In this paper, we study a class of spin systems with two-fold orbital degeneracy in the limits of large spin. Our model can be derived from a lattice of ions with high spins ferromagnetically coupled to the strongly correlated electrons with double orbital degeneracy. In the classical limit of spin, such a system can be mapped onto two decoupled black and white sub-systems, representing classical spin up and down, respectively. This provides an efficient method to study the ground states. We apply this method to study the case where the electron hopping integrals are isotropic and use the linearized spin wave theory to examine the effects of the quantum spin fluctuation. The ground state in 1D chain is a spin disordered valence bond state, in which each bond is an orbital valence bond (orbital singlet, parallel spins) and neighboring bonds interact antiferromagnetically. The ground states in a 2D square lattice and in a 3D cubic lattice are spin antiferromagnetically ordered orbital valence bond (OVB) state. We discuss possible relevance of the theory to some transition metal oxides.

We start with an Hamiltonian describing correlated electrons coupled to the spins of local ions,

$`H`$ $`=`$ $`t{\displaystyle \underset{ij,\alpha ,\sigma }{}}\left(c_{i\alpha \sigma }^{}c_{j\alpha \sigma }+c_{j\alpha \sigma }^{}c_{i\alpha \sigma }\right)`$ (2)
$`+\text{ }U{\displaystyle \underset{i,\alpha \sigma \alpha ^{^{}}\sigma ^{}}{}}n_{i\alpha \sigma }n_{i\alpha ^{}\sigma ^{}}{\displaystyle \underset{i,\alpha }{}}J_H𝐒_i^{(e)}𝐒_i^{(ion)}`$
In the above equations, $`\alpha =1,2`$ are the orbital indices and $`\sigma `$ is the electron spin. The sum over $`ij`$ runs all the nearest neighbor (n.n.) pairs. The electronic part of the Hamiltonian is a generalized Hubbard Hamiltonian with two-fold orbital degeneracy, where the Hund’s rule coupling and anisotropy of the electron hopping integrals are neglected. $`𝐒_i^{(e)}`$ and $`𝐒_i^{(ion)}`$ are the spin operators of electron and ion at site i, respectively. We consider $`J_H0`$, consistent with the Hund’s rule. In the strong coupling limit, $`Ut,sJ_H`$, with $`s`$ the quantum number of the ion spin, the projection perturbation theory may be applied to study low energy physics of the system. At the filling of one electron per site, and up to the order of $`t^2/U`$, this leads to an effective Hamiltonian,

$`H_p`$ $`=`$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{i,j}{}}\left(2𝐒_i^{(e)}𝐒_j^{(e)}+1/2\right)\left(2𝐓_i𝐓_j+1/2\right)`$ (4)
$`{\displaystyle \underset{i}{}}J_H𝐒_i^{(e)}𝐒_i^{(ion)}`$
where $`𝐓_i`$ is the electron orbital operator, $`J=4t^2/U`$. Furthermore, in the case $`sJ_H>t^2/U`$, the electron and ion spins strongly bind ferromagnetically to form a state with total spin $`s_t=s+1/2`$. The effective Hamiltonian within this Hilbert space of total spin $`s_t`$ can be obtained by applying the projection operator,

$$P=\underset{i}{}\frac{2𝐒_i^{(e)}𝐒_i^{(ion)}+(s+1)}{2s+1}.$$
(5)
The new coupled spin-orbital model is found to be ,

$`H_{eff}`$ $``$ $`PH_pP`$ (6)
$`=`$ $`{\displaystyle \frac{J}{2s_t^2}}{\displaystyle \underset{ij}{}}\left(𝐒_i𝐒_j+s_t^2\right)\left(𝐓_i𝐓_j+1/4\right)`$ (7)
with $`𝐒_i=𝐒_i^{(e)}+𝐒_i^{(ion)}`$ the total spin operator with maximal eigenvalue $`s_t`$ at site i.

The model has SU(2)$`\times `$SU(2) symmetry, representing rotational invariance in both spin and orbital spaces. A special case is at $`s_t=1/2`$, corresponding to $`𝐒_i^{(ion)}=0`$. In that case, the Hamiltonian possesses a higher SU(4) symmetry and the model has been studied extensively. In this paper, we focus on another limit where $`s_t1`$. We start with the classical spin to replace $`𝐒_i`$ by a classical vector of length $`s_t`$ described by two angles $`\theta _i`$ and $`\varphi _i`$. The Hamiltonian then reads,

$$H_c=J\underset{ij}{}\mathrm{cos}^2\frac{\mathrm{\Theta }_{ij}}{2}\left(𝐓_i𝐓_j+1/4\right)$$
(8)
where $`\mathrm{\Theta }_{ij}`$ is the angle between the two spin vectors. To investigate the ground state, we use the variational principle to find the equations for $`\theta _i`$ and $`\varphi _i`$: $`\delta H_c/\delta \theta _i=0`$ and $`\delta H_c/\delta \varphi _i=0.`$ We see that $`\varphi _i=\varphi _j`$ and $`\theta _i=\theta _0`$ or/and $`\theta _0+\pi `$ at all the sites are solutions of the equations. Below we will consider these solutions. In 1D it can be shown that these solutions give the lowest energy . In 2D or 3D we speculate these solutions contain the lowest energy states, although a rigorous proof is absent. We set $`\theta _0=0`$ below, and label all the lattice sites with $`\theta _i=0`$ by blacks, and all the sites with $`\theta _i=\pi `$ by whites. The bond Hamiltonian is then reduced to

$$H_c(ij)=\{\begin{array}{cc}J(𝐓_i𝐓_j+1/4),\hfill & \text{if }\theta _i=\theta _j;\hfill \\ 0,\hfill & \text{if }\theta _i=\theta _j\pm \pi .\hfill \end{array}$$
(9)
Therefore, the coupled spin-orbital system is decomposed into disconnected black and white sites or blocks (collection of the connected same colored sites). In the same colored block, all the spins are parallel but the orbitals interact antiferromagnetically; and the interaction vanishes between different blocks. The total energy of the system is then the simple sum of the these blocks. This greatly simplifies the calculations.

Let us first consider a two-site problem. The ground state of Eq.(8) is an orbital singlet with total spin $`2s_t`$. We shall call this 2-site state as an orbital valence bond (OVB), whose energy is $`0.5J`$. There are two competing terms in Eq.(8). One is the spin-orbital coupled term, $`\mathrm{cos}^2(\mathrm{\Theta }_{ij}/2)𝐓_i𝐓_j`$, which favors FM spins and antiferromagnetic (AFM) orbitals. The other is the spin interaction term $`\mathrm{cos}^2(\mathrm{\Theta }_{ij}/2)/4`$, which favors AFM spins. Because of this competition, the ground state of more than 2-sites is generally not a uniform FM spin state as we will see explicitly below. In the 1D chain, we divide the chain into blacks and whites and calculate the lowest energies of these spin configurations. In particular, we consider the configurations with alternating black and white segments of n-sites (see Fig. 1a). The total energy of the chain is the sum of these independent segments. The ground state energies of these segments are calculated by using exact numerical diagonalization method from $`n=2`$ to $`16`$. The results are plotted in Fig. 2. The energies of the even n segments are lower than those of the odd n ones. The ground state is an OVB solid with alternating spins up and down with an energy $`0.25J`$ per bond, much lower than that of the uniform FM spin state corresponding to $`n\mathrm{}`$. The energy of the latter state can be deduced from the result of the Bethe ansatz solution, and it is $`(12\mathrm{ln}2)J/2=0.193J`$ per bond. Our numerical results approach to the exact result rapidly as n increases. Including the spin quantum fluctuation, we expect that the 1D AFM spin long range order be destroyed based on the dimensionality consideration and the calculation of the linearized spin wave theory. The system is then described by a spin-$`2s_t`$ chain of the OVBs with AFM Heisenberg coupling $`J/(32s_t^2)`$. Since $`2s_t`$ is an integer, the ground state is a Haldane’s gap state. The quantum ground state has two-fold degeneracy corresponding to the translational invariance by one lattice constant in the OVBs.

In a 2D square lattice, we consider various spin configurations including alternating colored 2-site bond state (Fig.1b), alternating 4-site plaquette state (Fig.1c), stripes and the uniform FM states. The ground states are found to be highly degenerate. The lowest energy state in Fig.1b is an OVB state, and the lowest energy state in Fig. 1c is a plaquette orbital singlet state whose orbital is given by  t$`(12)`$t$`(34)`$t$`(14)`$t$`(23)`$, with t(12) representing an orbital singlet of sites 1 and 2. The OVB, the plaquette orbital, and their mixed states are degenerate, and have energy $`0.125J`$ per bond. The uniform FM spin state has a much high energy $`<T_iT_j+1/4>_{2D}=0.085J`$ per bond as estimated from the known result of 2D spin-1/2 model . The OVB state is also found to have the lowest energy in 3D cubic lattice. The large degeneracy of the ground states found in the classical spin limit in higher dimension is removed when the spin quantum fluctuation is included. We have used the spin wave theory to calculate the energy correction to the classical spin states in 2D square lattice, and found it to be $`0.01157J/s_t`$ for the OVB state, and $`0.01012J/s_t`$ for the plaquette orbital state. Therefore the quantum spin fluctuation favors the OVB state. In 2D and 3D, the AFM spin long range order is expected to survive from quantum fluctuation. We thus conclude that the ground state of Eq.(7) in 2D and 3D is the OVB state with AFM ordered spins between the n.n. bonds.

We now turn to discuss the possible relevance of our theory to some transition metal oxides. We first examine the ground state of cubic vanadate LaVO3. The observed magnetic order of the vanadate is of C-type AFM phase (FM along the c-axis and AFM in the a-b plane), arising from a practically undistorted structure above the Néel temperatures. To explain the unusual magnetic ordering, Khaliullin et al. recently proposed a spin-orbital Hamiltonian for the oxides. In the limit $`\eta =J_H/U0`$ (the notations are the same as in their paper) and along a given cubic axis, their model is reduced to Eq.(7) in the present paper with $`s_t=1`$. The corresponding two equivalent orbitals for the coupling along c-axis are, dzx and dyz, for example. Khaliullin et al. compared the ground state energies of the C-type and G-type (AFM in all three directions) AFM phases. They concluded that the C-type AFM is more stable at $`\eta =0`$, and is further stabilized at $`\eta >0.`$ In that work, the comparison between the C-type AFM and the OVB state was not included. As we discussed earlier, the ground state is an OVB state at $`\eta =0`$ in the large spin limit. We believe the experimentally observed C-type AFM in that compound is stabilized by the Hund’s coupling. Indeed we have compared the energies of C-type AFM and OVB state in their model (Eq.(1) of Ref.), and found that the C-type AFM has lower energy for $`\eta >\eta _c`$ with $`\eta _c0.06`$ It will be interesting to examine the higher order spin fluctuation or to use numerical techniques to verify if the large spin limit applies to spin-1 system.

In the second example, we consider the spin-orbital model for LaMnO3 assuming the lattice is not distorted. This type of model has been studied by many authors and the ground state is A-type AFM (AFM along c-axis and FM within a-b planes). Here we wish to point out that the A-type ordering may be obtained by treating the Hund’s rule coupling as a small perturbations. The effective Hamiltonian describing spin $`S=2`$ and the two degenerate eg orbitals of the Mn-ion is given in the limit $`J_H0`$,

$$H=\frac{J}{2s_t^2}\underset{ij}{}\left(𝐒_i𝐒_j+s_t^2\right)\left(\tau _i^\gamma +\frac{1}{2}\right)\left(\tau _j^\gamma +\frac{1}{2}\right)$$
(10)
where $`\gamma `$ is along $`ji`$ and the orbital operators $`\tau _i^\gamma =\mathrm{cos}(2m_\gamma \pi /3)𝐓_i^z\mathrm{sin}(2m_\gamma \pi /3)𝐓_i^x`$ ($`m_\gamma =1,2,3`$) with eigenvalues $`\pm \frac{1}{2}`$. The ground state of this Hamiltonian in the large spin limit is highly degenerate. As the expectation value of $`\left(\tau _i^\gamma +\frac{1}{2}\right)\left(\tau _j^\gamma +\frac{1}{2}\right)`$ is always not less than zero, the low bound for the ground state energy is zero. Thus the G-, C-, and A-type AFM are all degenerate ground states with zero energy. When a small Hund’s rule coupling is introduced, the high degeneracy in the ground state is removed. To simplify our problem, we keep the perturbation term due to the Hund’s rule coupling up to the order $`JJ_H/U`$, $`\mathrm{\Delta }H=(J/2s_t^2)\left(𝐒_i𝐒_j+s_t^2\right)h`$, where $`h`$ is of order $`sJ_H/U1`$. Including $`\mathrm{\Delta }H`$, the A-type AFM spin state is most favorable in energy. In that state, the orbitals are C-type antiferro-orbitally ordered: d$`_{x^2z^2}`$/d$`_{y^2z^2}`$ alternative in the a-b plane, and ferro-orbital along the c-axis. This state is favored because of the energy gain of $`\mathrm{\Delta }H`$ from the bonds within the a-b planes. Our calculation shows that the superexchange interaction including the Hund’s coupling leads to a spin A-type and orbital C-type structure. This state might be realized in KCuF$`_3.`$ Note that the orbital ordering pattern is different from that experimentally observed in LaMnO3 with a d$`_{3x^2r^2}`$/d$`_{3y^2r^2}`$ alternating ordering. The observed orbital ordering is due to strong Jahn-Teller distortion , but not due to the superexchange .

Finally, we examine an extended model of Eq.(7) in a cubic lattice to include the anisotropic coupling in different axis and the Hund’s rule coupling,

$$H=\underset{i,a}{}J_\alpha \left(𝐒_i𝐒_{i+a}+s_t^2\right)(𝐓_i𝐓_j+\frac{B}{4})$$
(11)
where $`J_\alpha `$ is the coupling strength along the direction $`\alpha =a,b,c`$. At the symmetric point $`B=1`$ and $`J_a=J_b=J_c,`$ the model is reduced to Eq.(7) . In the presence of $`J_H,`$ $`B<1`$ (see in $`\mathrm{\Delta }H`$ in the previous example). We compare the ground state energies of various states including OVB state, the uniform FM spin state, and the A-, C-, and G-type AFM. In Fig. 3, we plot the phase diagram for the lowest energy states in the parameter space B and $`\eta =J_{ab}/J_c.`$ The OVB state is found to be stable ground state in a finite parameter space around the symmetry point. It will be interesting to find an experimental realization of the OVB state. Since the OVB state break the translational symmetry, we expect the electron-lattice interaction may accompany a crystal structure phase transition to result in shorter and longer bonds. At a larger $`\eta `$ and smaller $`B`$, C-type AFM spin state is realized, in which the orbitals are disordered liquid. At a smaller $`\eta `$ and smaller $`B`$, A-type AFM becomes stable, in which the orbitals has a long range AFM order in the x-y plane, and FM along the z-axis. The phase diagram is thus very rich, reminiscent of some of the features in the transition metal oxides.

In conclusion, we have examined the ground state of a coupled high-spin orbital model from the large spin limit point of view. A novel orbital valence bond state is found to be stable in a certain parameter region, and may be realized in future experiments. The anisotropic model contains both C-type and A-type AFM spin states among others. Our method should be useful in study of more realistic models when Hund’ rule coupling is weak. The large spin approach emphasizes orbital quantum fluctuation over spins, and special care is needed to apply the method to spin-1/2 case. In the SU(4) limit, there is a permutation symmetry between spin and orbital, and we do not expect the method will work. The quantum fluctuation may drive the two degenerated spin and orbital valence bond states to a quantum state with more intricate correlations. In fact, from the Bethe ansatz solution, the 1D chain has a gapless SU(4) liquid ground state, which is not dimerized.

We would like to thank Jian Wang and Michael Ma for useful discussions. This work was in part supported by the RGC Grant No. HKU 7088/01P from Hong Kong, and by the US DOE Grant No. DE/FG03-01ER45687, and by the Chinese Academy of Sciences.


text:
# Turning the Cosmological Constant into Black Holes

Permanent address: \]Physique Théorique et Mathématique,Université Libre de Bruxelles,Campus Plaine C.P. 231, B–1050 Bruxelles, Belgium.

(November 4, 2001)

## Abstract

It is known that there is a quantum mechanical tunneling process which, through nucleation of a membrane, induces a transition between two de Sitter spaces, lowering the cosmological constant. It is shown in this paper that a different, new, membrane nucleation process exists, which, in addition to lowering the cosmological constant, leaves a black hole behind. Once a black hole is present, the relaxation of the cosmological constant may proceed via an analog of the old process, which decreases the black hole horizon area, or via the new process, which increases it.

Microscopic physics would appear to predict an enormously large value for the cosmological constant $`\mathrm{\Lambda }`$ Weinberg . However, the observed value appears to be very small and positive obs . This makes it very appealing to consider $`\mathrm{\Lambda }`$ as a variable that can evolve dynamically, and relax through quantum mechanical tunneling.

The simplest way of implementing this view is, in four spacetime dimensions, through the coupling to the gravitational field of a 3–form gauge potential Aurilia ; Duff ; Freund ; Hawking ; HennCT . In addition to the gauge potential it is natural to bring in also its source, which is a two dimensional membrane that sweeps a three dimensional history. There exists then a quantum mechanical process – analogous to pair creation by an electric field in two spacetime dimensions – which nucleates a membrane and induces a transition between two de Sitter spacetimes, the final one having a value of the cosmological constant lower than the initial one BrownCT ; BP ; Feng .

In the case of pair production, the energy of the pair is provided by the decrease in the energy of the electric field in the space between the pairs. Much in the same way, the energy of the membrane which is nucleated may be accounted for, in physical terms, by considering the reduction of the cosmological constant inside the bubble as a lowering of the vacuum energy.

Now, with the dynamical cosmological constant, one has a vast energy sea at one’s disposal. It is then natural to ask whether one could tap this reservoir to get from it more than just an expanding membrane, which escapes to infinity, and whose only effect is to lower the cosmological constant . Could one perhaps transmute the cosmological background energy into “localized matter” ?. And what localized matter can be more natural in the context of gravitation theory than a black hole ?. The purpose of this note is to answer this question in the affirmative: One may also nucleate a membrane which contracts to form a black hole, leaving outside of the hole a lower cosmological constant.

To make the point in the most economical manner we will consider membrane nucleation around a pre–existent black hole of mass $`M_+`$ and will then study, in particular, the case $`M_+0`$. This is useful because $`M_+0`$ is the generic situation and the limit $`M_+0`$ is somewhat singular, so that some important features become blurred. The case $`M_+0`$ was considered in CTplett  in the context of black hole thermodynamics. However, the new process presented in this paper, as well as some important features of the analog of the old process treated in BrownCT , were not discussed there.

We consider an (Euclidean) spacetime element of the form

$$ds^2=f^2(r)dt^2+f^2dr^2+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right).$$
(1)
The antisymmetric field strength tensor takes the form

$$F_{\mu \nu \lambda \rho }=(dA)_{\mu \nu \lambda \rho }=E\sqrt{g}ϵ_{\mu \nu \lambda \rho }.$$
(2)
The history of the membrane will divide the spacetime in two regions, one which will be called the interior, labeled by the suffix “-” and the other the exterior, labeled by the suffix “+”. The boundary in question may be described by the parametric equations

$$r=R(\tau )t_\pm =T_\pm (\tau ),$$
(3)
where $`\tau `$ is the proper time, so that its line element reads

$$ds^2=d\tau ^2+R^2(\tau )\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right),$$
(4)
with

$$1=f_\pm ^2\left(R(\tau )\right)\dot{T}_\pm ^2+f_\pm ^2\left(R(\tau )\right)\dot{R}^2.$$
(5)
In the “+” and “-” regions the solution of the field equations read

$`f_\pm ^2`$ $`=`$ $`1{\displaystyle \frac{2M_\pm }{r}}{\displaystyle \frac{r^2}{l_\pm ^2}},`$ (6)
$`E_\pm ^2`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}\left({\displaystyle \frac{3}{l_\pm ^2}}\lambda \right).`$ (7)
The actual cosmological constant $`\mathrm{\Lambda }=3/l^2`$ is thus obtained by adding $`\lambda `$ (normally taken to be negative), coming from “the rest of physics” and not subject to change, and the contribution $`4\pi E^2`$, which is subject to dynamical equations. The discontinuities in the functions $`f^2`$ and $`E`$ across the membrane are given by

$`f_{}^2\dot{T}_{}f_+^2\dot{T}_+`$ $`=`$ $`\mu R.`$ (8)
$`E_+E_{}`$ $`=`$ $`q.`$ (9)
Here $`\mu `$ and $`q`$ are the tension and charge on the membrane respectively (the action is given in (16) below). Eq. (9) follows from integrating Gauss’ law for the antisymmetric tensor across the membrane, whereas Eq. (8) represents the discontinuity in the extrinsic curvature of the membrane when it is regarded as embedded in either the “-” or the “+” spaces Israel . In writing these equations the following conventions have been adopted, and will be maintained from here on: $`(i)`$ The coordinate $`t`$ increases anticlockwise around the cosmological horizon, $`(ii)`$ the variable $`\tau `$ increases when somebody traveling along the curve leaves the interior on his right side.

Equation (8) may also be interpreted as the first integral of the equation of motion for the membrane, which is thus obtained by differentiating it with respect to $`\tau `$ (“equations of motion from field equations”). Hence, satisfying (6-9) amounts to solving all the equations of motion and, therefore, finding an extremum of the action.

We are interested in a solution of the equations of motion of the membrane, which is a closed orbit, that will be interpreted as an “instanton”, leading to a reduction of the cosmological constant accompanied by the production of a black hole.

To this end we examine the “radii of formation” of the membrane, namely, the values of $`R`$ for which $`\dot{R}=0`$. This is most efficiently analyzed if one uses the proper time condition (5) to rearrange (8) to read

$`\mathrm{\Delta }M`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\alpha ^2\mu ^2\right)R^3\mu R^2f_+^2\dot{T}_+,`$ (10a)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(\alpha ^2+\mu ^2\right)R^3\mu R^2f_{}^2\dot{T}_{},`$ (10b)
where

$$\alpha ^2=\frac{1}{l_+^2}\frac{1}{l_{}^2}\mathrm{\Delta }M=(M_{}M_+).$$
(11)
The graph of $`\mathrm{\Delta }M`$ as a function of $`R`$ when $`\dot{R}=0`$ is given in Fig. 1. There are two branches which merge smoothly, corresponding to taking both signs of the square root of $`\dot{T}_\pm ^2`$ in (5) (much as $`x=\pm \sqrt{1y^2}`$ gives the smooth circle $`x^2+y^2=1`$).

We see that for each $`\mathrm{\Delta }M`$ within the range shown in fig. 1, there are two roots, which correspond to the simultaneous creation of two membranes, which have opposite polarities, i.e., one has charge $`q`$, the other $`q`$. In order for the solution to represent an extremum of the action associated with tunneling, these two roots must be the two turning points of a single, time symmetric, Euclidean orbit. Now, it may be verified by studying the orbits which solve (10), that this does not happen for the generic $`\mathrm{\Delta }M`$. However, there are two special “mass gaps” for which the orbit indeed closes. They correspond to the limiting cases when $`(i)`$ the smaller membrane sits on the horizon of the initial black hole, $`(ii)`$ the larger membrane sits on the final cosmological horizon. In the limiting case, $`M_+=0`$, discussed in BrownCT  the smaller membrane has zero radius, disappearing from the problem. However, in the generic case, $`M_+0`$, two membranes are formed simultaneously. The closed orbits and their properties are discussed in Fig. 2.

The joining of the Euclidean and Lorentzian sections is not devoid of subtleties, and will be discussed in detail in GHT . The resulting processes may be described as follows. At a certain moment two membranes are created and the cosmological constant is lowered in the space in between them. For the old process, the first membrane materializes at the radius $`R_0`$ between the black hole and cosmological horizons, and proceeds to expand, “collapsing outwards”. At the same time, the second membrane materializes at a radius $`r_+`$, the value of $`r`$ corresponding to the initial black hole horizon. With the creation of the membranes, the radius of the black hole horizon decreases ($`r_{}<r_+`$), while that of the cosmological horizon increases ($`r_{}>r_{++}`$). The second membrane collapses inwards beyond the new black hole horizon. Note that this second membrane can exist only if there is initially a black hole. The inside of the first membrane is the “interior” in this case.

For the new process the roles of the two horizons are interchanged. The first membrane materializes at the corresponding $`R_0`$ and proceeds to contract, collapsing inwards. At the same time, the second membrane materializes right on the final cosmological horizon. This time the black hole radius increases and the cosmological horizon radius decreases. The second membrane expands beyond the cosmological horizon. In this case the first membrane collapses onto the black hole, or forms one if initially there is none. Now the outside of the first membrane is the “interior”.

Note that with the existence of the new process the symmetry between cosmological and black hole horizons, which had been lost with the presence of the old process alone, is restored. The only remaining asymmetry is that one can consider the absence of an initial black hole, but one cannot allow, in this context, for the absence of an initial cosmological horizon.

The tunneling processes change both parameters $`M`$ and $`l`$ of the Schwarzschild–de Sitter solution. For a given change in the cosmological constant, the change in the mass is obtained by setting, in Eq. (10), $`R`$ equal to either the initial black hole horizon, $`r_+`$ (old process), or to the final cosmological horizon, $`r_{}`$ (new process). This yields

$`\mathrm{\Delta }M`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\alpha ^2\mu ^2)r_+^3\text{(old process)},`$ (12)
$`\mathrm{\Delta }M`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\alpha ^2+\mu ^2)r_{}^3\text{(new process)}.`$ (13)
In particular, when $`M_+=0`$, Eq. (12) gives $`M_{}=0`$, in accordance with BrownCT . However, for the new process, one obtains for $`M_{}`$ the value,

$$M_0=\frac{1}{2}\frac{\left(\frac{\mathrm{\Lambda }_+}{3}\frac{\mathrm{\Lambda }_{}}{3}+\mu ^2\right)}{\left(\frac{\mathrm{\Lambda }_+}{3}+\mu ^2\right)^{3/2}}$$
(14)
for the mass of the black hole spontaneously formed out of de Sitter space, when the cosmological constant decreases from $`\mathrm{\Lambda }_+`$ to $`\mathrm{\Lambda }_{}`$, through nucleation of a pair of membranes of tension $`\mu `$.

Having established the existence of the instantons, we pass to discuss the probability for each of the two processes. The general rule is that, in the semiclassical approximation, the probability is the exponential of the Euclidean action. However, since there are thermal effects at play due to the presence of event horizons, one must distinguish those from the underlying quantum mechanical probability for the process. This is equivalent to subtracting the action of the background.

When there is only one horizon at play, as it happens for the old process in the absence of an initial black hole, the identification of the thermal effects can be done neatly since there exists a complete action principle for the coupled system of the gravitational field, the antisymmetric field, and the membrane. The action, appropriate for keeping fixed the initial cosmological constant, $`\mathrm{\Lambda }_+`$, is then equal to the sum of four contributions. They are the “bulk” gravitational and antisymmetric tensor hamiltonian actions, the membrane action (which includes couplings with the gravitational and antisymmetric tensor fields), and the surface term one fourth of the cosmological horizon area. The surface term, is the horizon entropy, and it stems from the thermal nature of the problem. The rest of the action may be thought of as yielding the quantum mechanical probability for the process itself. Furthermore, the on–shell value of the bulk Hamiltonian action turns out to be zero, and thus the action for the background is just one fourth of the area. Therefore, in the semi–classical approximation, the quantum mechanical probability for the process is given by the exponential of the membrane action alone.

$$P=\mathrm{exp}\left[I_{\text{membrane}}\right],$$
(15)
with

$$I_{\text{membrane}}=\frac{\mu }{4\pi }V_3+\frac{3}{8\pi }\alpha ^2V_4.$$
(16)
Here $`V_3`$ is the $`3`$–volume of the history of the membrane and $`V_4`$ is the $`4`$–volume of its interior. The $`4`$–volume term may be understood as arising from the standard minimal coupling through application of Stoke’s theorem. The coefficient $`3\alpha ^2/8\pi `$ comes from writing

$`{\displaystyle \frac{3}{8\pi }}\alpha ^2`$ $`=`$ $`{\displaystyle \frac{3}{8\pi }}\left({\displaystyle \frac{1}{l_+^2}}{\displaystyle \frac{1}{l_{}^2}}\right)={\displaystyle \frac{1}{2}}\left(E_+^2E_{}^2\right)`$
$`=`$ $`\left(E_+E_{}\right){\displaystyle \frac{1}{2}}\left(E_++E_{}\right)=qE_{\text{av}},`$
where $`E_{\text{av}}=(E_++E_{})/2`$ is the field on the membrane.

When two horizons are present in the instanton, as it happens for the old process with an initial black hole or for the new process, always, there is no complete action principle. This is because one cannot include appropriate surface terms in the action which ensure that, on–shell, there is no conical singularity at either horizon (the time periods do not match). This problem is not a peculiarity of the membrane, or of the antisymmetric tensor. It is already present for the pure gravitational field. Thus, already the action for the background is not well defined. In physical terms this is a reflection of the fact that the black hole and cosmological horizons are not in thermal equilibrium GibbHawk . However, in view of the previous discussion, it would seem justified to assume that Eq. (15) remains valid as the quantum mechanical probability for membrane nucleation for each of the two processes.

Since the main purpose of this communication is to point out that the new process exists, we will not dwell here on the actual evaluation of the semiclassical tunneling probability. We just point out that it will be necessary to appeal to approximations to make the problem tractable, because in all cases except for the tunneling between two de Sitter spaces discussed in BrownCT , the geometry has a low degree of symmetry. For the old process with $`M_+=M_{}=0`$, the instanton consists of two $`4`$–spheres joined on a $`3`$–sphere. However, in all other cases, for the new and the old processes, one has two $`S_2\times S_2`$’s joined on a $`S_1\times S_2`$. Furthermore, the $`S_2`$’s are just topological, not metrical spheres.

A fuller account of the results presented in this paper and their possible implications is now in preparation, and will be published elsewhere GHT . However, even at this early stage, it is hard to refrain from speculating that the low value of the cosmological constant might point to the presence of mass in the universe, in the form of black holes created by its relaxation.

###### Acknowledgements.

We thank Valery Rubakov for an enlightening discussion. This work was funded by an institutional grant to CECS of the Millennium Science Initiative, Chile, and also benefits from the generous support to CECS by Empresas CMPC. AG gratefully acknowledges support from FONDECYT grant 1010449 and from Fundación Andes. AG and CT acknowledge partial support under FONDECYT grant 1010446. The work of MH is partially supported by the “Actions de Recherche Concertées” of the “Direction de la Recherche Scientifique - Communauté Française de Belgique”, by IISN - Belgium (convention 4.4505.86) and by the European Commission RTN programme HPRN-CT-00131, in which he is associated to K. U. Leuven.


text:
# Gravitational clustering in N-body simulations

## Introduction

The study of the features of gravitational clustering is one of the most challenging problems of classical physics, and an interesting example of interplay between statistical mechanics and astrophysics.

From the point of view of statistical mechanics, it is very hard to study the properties of an infinite system of self-gravitating particles. This is mainly due to the long range nature of gravitational potential, which is not shielded by the balance of far away charges, as e.g. in a plasma. Therefore all scales contribute to the evaluation of the potential energy of a particle. As a consequence of such difficulties, a satisfying thermodinamic equilibrium treatment of such systems is still lacking.

Most studies on self gravitating systems have been focusing on clustering properties of well-defined physical systems: from interstellar cold molecular clouds to globular clusters, to clusters and superclusters of galaxies. The range of scales on which such processes take place is impressive (from $`10^1pc`$ to $`10^8pc`$). This also implies interactions with physical processes of different types, depending on the scale (from e.g. turbulence in cold molecular clouds to cosmological expansion above galaxy scale).

In spite of the large amount of analytical and numerical studies on gravitational clustering, it is very difficult to retrieve a clear picture of the statistical properties of self-gravitating systems, exactly because of the richness of physical processes which they take into account.

What has especially raised our interest in the subject was the results of a series of analyses of space correlations in galaxy catalogues by Pietronero et al. Pietro . They revealed that galaxies form a fractal set with well defined properties at least up to a certain scale. The discovery of the fractal characteristics of large scale structures in the universe raises a number of new and extremely appealing questions, which have neverthless been almost completely neglected. For example we can ask whether gravity is a self organizing process, or if the formation of fractal properties depend on the choice of a particular set of initial conditions, or even one has to invoke some other physical process.

In a more general perspective, we intend to study in detail the characteristics of the evolution of self gravitating systems, from the point of view of statistical mechanics. In particular we focused on the evolution of spatial correlations in self-gravitating systems, and on the possibility that they may develop some kind of self similar spatial or temporal feature.

N-body simulations represent in fact a valuable experimental setting for a careful investigation of the dynamics of clustering. Current astrophysical simulations have reached a high level of refinement, both in resolution and in the number of different physical processes which they take into account. Such characteristics allow them to study in great detail the single physical problem for which they are developed. On the other hand they don’t allow a clarification the common rôle and the peculiarities of gravitational interaction.

Our approach has been completely different: we start from the simplest possible model and study in detail what happens during the clustering. Even such a simple model gives raise to a rich phenomenology, which has not been studied in a systematic way up to now.

## Theoretical issues on gravitational clustering

The peculiar form of the gravitational potential makes the statistical properties of self-gravitating systems a very difficult subject to study. Two classes of problems arise: those due to the short range (i.e. $`r0`$) divergence and those due to the long range (i.e. $`r+\mathrm{}`$) behaviour.

The former is not uncommon, since it is the same problem which arises in electromagnetism. The divergence would cause, e.g., the Boltzmann factor to diverge in the limit $`r0`$. A typical prescription is to put a small distance cut-off in the potential. The physical nature of this cut-off may be due to many effects, e.g. the dynamical emergence of angular momentum barriers.

The long range behaviour is of much more concern and is, in fact, the problem. It is an easy exercise to verify that the energy of a particle in an infinite self-gravitating system diverges. This causes the energy to be non-extensive. As a consequence, a thermodynamical limit is not achieved, since as the number of particles goes to infinity while keeping the density constant, the energy per particle diverges. Strangely enough, such a problem has not been fully appreciated by many physicists in the field (see e.g. Padma ), as they try to avoid the long range divergence by putting the system in a box “as it is usually done with ordinary gas”. In fact, the difference is that in ordinary gas, when confining the system in a box, the energy per particle is equal to a constant plus a surface term that goes to zero in the thermodynamical limit. In self-gravitating systems, due to non extensivity, the energy per particle is neither a constant, nor the surface term goes to zero (in fact, it is of the same order of magnitude as the potential energy due to particles belonging to the system).

Another very interesting consequence, which is often not appreciated, is that the thermodynamical definition of temperature, as the parameter which controls the equilibrium of the system, doesn’t hold for a self gravitating system, since one cannot divide a system into smaller subsystems with the same thermodynamic properties of the larger system.

So far we have discussed the problems one encounters when one tries to build a theory for thermodynamic equilibrium. However we are much more interested in what happens out of equilibrium, during the evolution of a system. For example it is highly probable that any fractal state could not survive when the system would reach the equilibrium. Problems due to long range interactions are present all the same in most approaches, but they are more subtle to be put in evidence.

## N-body simulations

N-body simulation is one of the main tools for the study of gravitational clustering. They provide in fact the possibility to test theoretical models and to study the evolution of selfgravitating systems, a feature that is very hard to obtain from observations. Several different algorithms are used (see nbdrev  for a review).

A brute force method to evaluate forces acting on particles in a simulation requires $`O(N^2)`$ computations, a feature which forbids them to be useful if large number of particles are required. The tree algorithm (see barnes , and alltree  for a detailed explanation) is an approximation scheme which involves $`O(NlogN)`$ computations, and whose level of accuracy can be tuned to match the problem requirements.

It is based on the decomposition of the simulation volume into a tree-like structure. Particles are grouped with well-defined rules in a hierarchy of clusters. The contribution of distant clusters to the force on a particle is taken into account by a multipole expansion of the potential generated by them. This approximation greatly reduces the computational expense of a simulation.

The algorithm has to be elaborated for our purposes by adding periodic boundary conditions, which allow a statistically equivalent environment for all the particles, and also include the long range interactions.

In our simulation we started from the simplest possible model:

i) random (white noise) initial positions of particles;

ii) no cosmological expansion;

iii) zero initial velocities;

iv) equal mass particles.

The main tool we use to analyse spatial correlation of the system is to measure the conditional average density $`\mathrm{\Gamma }^{}(r)`$. It is defined as $`\mathrm{\Gamma }^{}(r)=<N(r)>/V(r)`$ where $`<N(r)>`$ is the average number of points (not considering the center) in a sphere of radius $`r`$ and volume $`V(r)`$ centered on each point of the system.

## Results

All the simulations we have performed have the same average number density, $`\overline{n}=1`$ which sets the physical scale of the system.

Figure 1 shows the evolution of $`\mathrm{\Gamma }^{}(r)`$ with time for a simulation of $`32^3`$ particles. This picture shows some interesting features, and allow us to focus on the different aspects of gravitational clustering. We can identify the region left of the dash-crossed line as a regime of mainly two particles interactions. In fact for such small scales, we find on average less than one particle in each sphere centered on each particle of the system. This region shows a very nice power law behaviour which seems stable for quite a long time, Such behaviour is not reported in the literature, as far as we know. The growing front of the correlations shows beautiful features of “temporal self-similar behaviour”, in the sense that if we rescale length scales with a function of time the curves collapse. The scaling factor is exponential. It is interesting to notice that while the exponential time dependence can be naively retrieved (with the right coefficients) using the solutions for the linear theory of the growth of density perturbation in a self gravitating fluid, the behaviour of the correlations is not well described by such solutions. In fact they predict a selfsimilar (in time) growth of the initial correlations of density fluctuations, until $`\delta n/n<1`$. On the other hand, an inspection of the fig.1 shows that such growth takes place at large number density contrasts (i.e. $`\delta n/n>1`$), and does not correspond to a pure amplification of the initial conditions. Anyway, such theory is considered very well verified in cosmological simulations. There may be several explanations. One possibility is that some of the approximations used to derive such solutions break earlier without expansion. Another interesting remark could be that in our simulations we don’t suppress the granular features of particles with large smoothing length, as it is customary in many cosmological simulations (see e.g. Theis ). Such an observation paves the way to a very interesting issue which we are currently investigating: whether the granularity is an essential feature to the understanding of gravitational clustering, or which is the limit of validity of the fluid description of a self gravitating medium. Another region of interest would be the one between the two we have just described. With present data, we cannot resolve the behaviour of correlations in such regime, due to small numbers of particles. Since simulations are independent of finite size effects until $`\mathrm{\Gamma }^{}(r)`$ starts differing from $`\overline{n}`$ on scales of the order of the system size, the more particles we have, the larger the system size, the farther in time we can observe the evolution of correlations. For this reason we are now working on a parallelisation of our code, which would allow us to simulate systems with more than a million of particles.