format: text

text:
# Sufficiently Fat Polyhedra are Not 2-CastableResearch funded in part by NSERC Canada

## 1 Introduction

Casting is a common manufacturing process where some molten substance is poured or injected into a cavity (called a *mold*), and then allowed to solidify. In many applications (see e.g. ) it is desirable to remove the cast object from the mold without destroying the mold (or, obviously the recently manufactured object). In general this requires the mold be partitioned into several parts, which are then translated away from the cast object. In the simplest case (prevalent in sand casting), a mold for polyhedron $`P`$ is partitioned into two parts using a plane. If a successful partition (i.e. both parts can be removed by translations without collisions), we say that $`P`$ is $`2`$-castable.

Guided by intuition about smooth objects, one might suspect that all convex polyhedra are $`2`$-castable. This turns out not to be the case. Bose, Bremner and van Kreveld  gave an example of a 12 vertex convex polyhedron that is not $`2`$-castable. Unfortunately the proof of non-castability relies on a computer based exhaustive search. Majhi, Gupta and Janardan  gave a simpler example with only 6 vertices; here the proof of non-castability is left as an exercise. In neither case can one draw any general conclusions (beyond the tautological) about what sort of convex polyhedra are $`2`$-castable.

In the present note we provide a general class of convex polyhedra that are not $`2`$-castable. In particular, we establish that if the polyhedron has vertices and facets in general position and is a sufficiently close approximation of a sphere, then it is not $`2`$-castable.

## 2 Background

We will actually consider a slightly more general definition of $`2`$-castability. We start with the definition of a castable polyhedron. Consider a 3-dimensional half-space $`H`$. Let $`P`$ be a polyhedron that lies in $`H`$, while one of its faces $`F`$ lies on the boundary of $`H`$. The set $`M=HP`$ is called a mold for $`P`$. We say that $`P`$ is castable with respect to a face $`F`$ (or equivalently with respect to a mold $`M`$), if we can pull $`P`$ out of $`M`$ by moving it along some vector $`d`$ without collisions (e.g. interior intersections). If there exists such a face of $`P`$, then $`P`$ is called castable. A polyhedron $`P`$ that can can be divided into $`k`$ castable parts is called $`k`$-castable. Note that for $`k=2`$ it is often required by manufacturing processes that the two halves are castable with respect to the same mutual facet. This constraint is relaxed here, however we do require that the halves are separated by a plane (i.e. have a mutual facet).

###### Definition 1.

A convex polyhedron $`P`$ is called $`(R_i,R_o)`$-fat  if there are two concentric spheres $`D_i`$ and $`D_o`$ of radii $`R_i`$ and $`R_o`$, such that $`D_iPD_o`$.

###### Definition 2.

A convex polyhedron is said to be in general position if none of its 4 points are coplanar and none of its 3 faces are parallel to a line.

Let $`P`$ be a general position $`(R_i,R_o)`$-fat polyhedron. By scaling $`P`$ we can assume $`R_i=1`$ and $`R_o=R>1`$. For the rest of the paper we use *fat* to mean $`(1,R)`$-fat. Let $`O`$ denote the center of concentric spheres from Definition 1. The following lemma gives bounds for various elements of a fat polyhedron.

###### Lemma 1.

Let $`P`$ be a fat polyhedron. Every edge of $`P`$ has length at most $`l^{}=2\sqrt{R^21}`$, every face has area at most $`S^{}=\pi (R^21)`$ and its volume is bounded $`\frac{4}{3}\pi <V(P)<\frac{4}{3}\pi R^3`$.

###### Proof.

Let $`AB`$ be an edge of $`P`$ and $`O^{}`$ be the projection of $`O`$ to $`AB`$. Since $`|OA|R`$ and $`|OO^{}|1`$, we have $`|AO^{}|\sqrt{R^21}`$, so $`|AB||AO^{}|+|O^{}B|<2\sqrt{R^21}`$.

Every face $`F`$ of $`P`$ defines a slice $`C`$ of the outer sphere $`D_o`$, indeed $`F`$ is contained within the disk $`C`$. By the previous consideration, the radius of $`C`$ is at most $`l^{}/2`$, therefore we have the bound $`S(F)\frac{\pi }{4}(l^{})^2`$ on the area of the face $`F`$.

Since $`D_iPD_o`$, we have the bounds on the volume of $`P`$ given by the lemma. ∎

The following observation is simple but important for the rest of the paper: $`l^{}`$ can be made arbitrary small by an appropriate choice of $`R`$, that is $`l^{}0`$ as $`R1`$. The following lemma gives an upper bound to the volume of a castable polyhedron.

###### Lemma 2.

Suppose that $`P`$ is castable through a face $`F`$ of area $`S`$. Let $`H`$ be the plane containing $`F`$ and $`h`$ be the maximum distance from a point $`P`$ to $`H`$. Then $`V(P)Sh`$.

###### Proof.

Let $`v`$ be the inner normal vector to $`F`$ and $`F(t)`$ be the area of $`PH+tv`$, for $`t0`$. Since $`P`$ is castable through $`F`$, the area $`F(t)`$ cannot be less then $`S`$. Thus $`V(P)=_0^hF(t)Sh`$. ∎

Let us call $`h`$ the thickness of $`P`$ with respect to $`F`$. Note that thickness is bounded by the diameter of $`P`$, thus it cannot exceed $`2R`$.

## 3 The proof of non-castability

We use the following method to prove that $`P`$ is not $`2`$-castable. We will consecutively assume that certain polyhedra are castable under some restrictions and than argue that this implies certain lower bounds on $`R`$. Assume that all the possible situations are covered and let $`R>R^{}`$ be the loosest bound on $`R`$. Then $`P`$ is not $`2`$-castable provided $`R<R^{}`$. In the following, let $`S(F)`$ denote the area of polygon $`F`$.

First assume that $`P`$ is 1-castable through some face $`F`$. Using Lemma 2 and Lemma 1, we derive the following inequality.

$$\pi (R^21)2RS(F)hV(P)>\frac{4}{3}\pi $$
(1)
This inequality implies the bound $`R>1.240011810`$.

Suppose that $`P`$ is $`2`$-castable. Let $`P`$ be sliced by a plane. Denote the larger part by $`P_1`$, the smaller by $`P_2`$ and their mutual face by $`C`$, so that $`V(P_1)V(P)/2V(P_2)`$. To simplify the presentation here, without loss of generality assume that $`O=(0,0,0)`$, $`C`$ is horizontal, $`P_1`$ lies above $`C`$ and $`P_2`$ below. Let the plane containing $`C`$ be given by the equation $`z=z_0`$. Each of $`P_1`$ and $`P_2`$ has to be $`1`$-castable.

Assume that $`P_1`$ is castable through a face $`FC`$. As before, we have

$$\pi (R^21)2RS(F)hV(P_1)\frac{V(P)}{2}>\frac{2}{3}\pi $$
(2)
Note, that the bound in (2) is looser then in (1). The numerical solution gives $`R>1.137158043`$.

Assume that $`P_1`$ is castable through $`C`$, and $`P_2`$ is castable through a face $`FC`$.

First consider the case, where $`z_00`$. Using Lemma 2 for $`P_2`$ we derive

$$\pi (R^21)2R>S(F)h>V(D_i)/2=2\pi /3$$
(3)
since $`P_2`$ contain the lower part of the inner sphere.c Solving (3) numerically we obtain $`R>1.137158043`$.

Now consider the case of $`z_0<0`$. Since $`P_1`$ contain the disk $`D_i\{z=0\}`$, we require that the diameter of the slice $`C`$ is at least $`1`$, and hence

$$\sqrt{R^2z_0^2}>1$$
(4)
Using Lemma 2 for $`P_2`$ gives

$$\pi (R^21)2\sqrt{R^2z_0^2}S(F)hV(P_2)>\pi _1^{z_0}1t^2dt=\pi (2/3+z_0z_0^3)$$
(5)
since the diameter of $`P_2`$ is at most $`2\sqrt{R^2z_0^2}`$. The numerical solution of the system of (4) and (5) gives $`R>1.07218989`$.

We summarize the arguement so far with the following lemma.

###### Lemma 3.

If $`P`$ is $`2`$-castable then one of the following conditions is true

* both $`P_1`$ and $`P_2`$ are castable through the face $`C`$.
* $`R>1.07218989`$

Assume that both $`P_1`$ and $`P_2`$ are castable through $`C`$. For each edge $`e`$ of $`C`$ consider its incident faces $`F_1`$ in $`P_1`$ and $`F_2`$ in $`P_2`$ other than $`C`$. We *mark* $`e`$, $`F_1`$ and $`F_2`$ if these faces constitute a face of $`P`$.

###### Lemma 4.

There are at most $`2`$ unmarked edges.

###### Proof.

Each unmarked edge corresponds to an edge of $`P`$ that lies in $`C`$. If we have more than two such edges, then there are at least $`4`$ vertices of $`P`$ that lie in $`C`$, which contradicts the general position assumption. ∎

Consider the set of feasible casting directions $`d=(d_x,d_y,d_z)`$ for a polytope $`P`$. Without loss of generality assume that $`d_z=1`$ for $`P_1`$ and $`d_z=1`$ for $`P_2`$.

It is known  that each face $`F`$ of $`P_1`$ ($`P_2`$) implies a linear constraint on $`d`$, namely $`(\mu ,d)0`$, where $`\mu `$ is the outward normal to $`F`$ with respect to $`P_1`$ ($`P_2`$). We restrict ourselves to the faces of $`P_1`$ ($`P_2`$) which are incident with marked edges of $`C`$. Let $`\mathrm{LP}_1`$ ($`\mathrm{LP}_2`$) be the corresponding 2-dimensional linear programs. Castability of $`P_1`$ ($`P_2`$) implies feasibility of $`\mathrm{LP}_1`$ ($`\mathrm{LP}_2`$). Note that if a face $`F_1`$ of $`P_1`$ contributes to $`\mathrm{LP}_1`$ then the incident face $`F_2`$ of $`P_2`$ contributes to $`\mathrm{LP}_2`$. The corresponding inequalities are:

$`l_1(d)`$ $`=`$ $`\mu _xd_x+\mu _yd_y\mu _z0`$ (6)
$`l_2(d)`$ $`=`$ $`\mu _xd_x+\mu _yd_y+\mu _z0`$ (7)
Let us define $`\mathrm{feas}(\mathrm{LP}_i)`$ to be the feasible region of a program $`\mathrm{LP}_i`$, i.e. the intersection of the constraints $`l_i(d)0`$ for each $`\mu `$. We study these programs in more detail. First note that feasible regions of these programs have to have inner points. For otherwise some there are either:

* Two constraints of the form $`(\mu ,d)0`$ and $`(\mu ,d)0`$ (their bounding lines coincide). Then the corresponding faces are parallel.
* Three constraints, such that their bounding lines on the plane $`d_z=\pm 1`$ intersect in a point. This means that the corresponding faces are parallel to a line.

Let $`d`$ and $`e`$ be inner points in $`\mathrm{feas}(\mathrm{LP}_1)`$ and $`\mathrm{feas}(\mathrm{LP}_2)`$ respectively. Then for every $`l_1\mathrm{LP}_1`$, we have $`l_1(d)<0`$. So $`l_2(d)=l_1(d)>0`$ and $`d`$ satisfies no constraints of $`\mathrm{LP}_2`$. Consider the ray $`r=(d)+\lambda (e(d))`$ starting at $`d`$ towards $`e`$. The segment between $`d`$ and $`e`$ of this ray has to intersect every bounding line in $`\mathrm{LP}_2`$, since one of its endpoints satisfies all the constraints while the other satisfies none of them. Therefore the remainder of $`r`$ (beyond $`e`$, $`\lambda 1`$) can not intersect any any of the bounding lines of $`\mathrm{LP}_2`$. We conclude that $`\mathrm{feas}(\mathrm{LP}_2)`$ is unbounded with respect to $`r`$. Similarly we can prove, that $`\mathrm{feas}(\mathrm{LP}_1)`$ is also unbounded. Suppose that $`\mathrm{feas}(\mathrm{LP}_1)`$ is unbounded along a ray $`r^{}=p+\lambda v`$, $`\lambda 0`$.

Consider the plane containing the convex polygon $`C`$ illustrated in Figure 1. Suppose there exists a marked edge $`e`$ with the outward normal $`n`$, such that $`(n,v)>0`$. Let $`\mu `$ be the outward normal of the corresponding face of $`P`$, whose projection to $`\{z=0\}`$ is $`n`$. Then

$$\mu _x(p_x+\lambda v_x)+\mu _y(p_y+\lambda v_y)+\mu _z=\lambda (n,v)+\mu _xp_x+\mu _yp_y+\mu _z0$$

for every $`\lambda 0`$, which is a contradiction. So all marked edges have outward normals $`n`$ such that $`(n,v)0`$. Consider the edges of $`C`$ with the outward normal $`n`$, such that $`(n,v)>0`$. All such edges are unmarked and form a chain since $`C`$ is convex. Define the segments $`QS`$ and $`RT`$ to be the segments parallel to the vector $`v`$ and touching the interior circle at the points $`Q`$ and $`R`$ and the points $`S`$ and $`T`$ lie on the exterior circle (see Figure 1). The first and last edge of $`C`$, that intersect the interior of $`QRTS`$ have to belong to this chain. So the chain connects $`QS`$ and $`RT`$, hence its length is at least $`2`$. But we know that there are at most two unmarked edges, thus we one of them has to be longer then $`1`$. This means that $`l^{}=2\sqrt{R^21}>1`$, an thus $`R>\sqrt{5/4}>1.118033989`$. Bringing all of the bounds on $`R`$ together, we conclude with

###### Theorem 1.

Let $`P`$ be a $`(R_i,R_o)`$-fat polyhedron in general position (i.e. no four vertices of which lie on a plane and no three faces are parallel to a line). Then $`P`$ is not $`2`$-castable if $`R_o/R_i<1.07218989`$.


text:
# Mott Domains of Bosons Confined on Optical Lattices

A considerable amount of work has been done in the last decade to determine the ground state phase diagram of correlated bosons on a lattice described by the “boson-Hubbard” Hamiltonian. On-site repulsion can produce a Mott insulating phase at commensurate fillings, with a quantum phase transition to a superfluid as the density is shifted or the interaction strength weakened. Longer range interactions can cause charge density wave, stripe, or even supersolid order. Extensions to disordered systems have allowed the detailed study of the interplay of randomness and interactions in quantum systems.

Recently, the trapping of atoms on optical lattices has given another experimental realization of these bosonic phases. However, the quadratic confining potential present in addition to the regular “lattice” potential, leads to a number of fundamentally new, and open, issues: (i) Does the confining potential preclude the formation of Mott regions by providing a continuous, unbounded, distribution of local site energies? (ii) If an insulating phase still exists, how is it characterized? (iii) What are the quantitative values of the trap curvature and interaction strength that support Mott phases? These questions are largely unaddressed in the literature.

In this paper we report the first quantum Monte Carlo simulation of the one-dimensional boson-Hubbard model in a confining quadratic potential and provide a quantitative map of the state diagram. To our knowledge, the only non-mean-field work on this problem is on very small systems (five particles). We find that the trap changes the physics fundamentally from that found in earlier simulations and subsequent analytic studies. For example, the vanishing of the global compressibility, discussed at length in  and reported in recent mean-field studies of the boson-Hubbard model in the context of optical lattices , but which ignores the confining potential, is absent. Other recent papers  on bosons in optical traps have likewise emphasized similarities to the physics in the absence of a confining potential, and have used values of the unconfined lattice critical coupling to compare with experimental data.

We find the inhomogeneous potential resulting from the confining trap is a crucial feature in discussing Mott regions at incommensurate fillings , and the physics of bosons on optical lattices generally. The Mott insulating regions exist above a threshold interaction strength, even without the commensurate filling required in the non-confined case. This is a consequence of the inhomogeneous distribution of boson densities which allows extended Mott domains with, for example, one particle per site to coexist with regions of other density and is a unique feature which distinguishes the behavior of the confined model. Similarly, the global compressibility is non-zero for all densities, including commensurate ones. The different regimes of the state diagram can be characterized by the topology of the density and local compressibility profiles. Therefore, in this context, the (global) density, $`\rho `$, defined as the total number of particles divided by the system size, loses its meaning since the particles are not uniformly distributed and the “system size” itself is ill-defined. Adding particles can push bosons deeper into the confining potential at the edges of the system thus changing the “size”.

We review briefly the properties of the ground state of the $`d`$-dimensional non-confined boson Hubbard model. This model has two phases: a Mott insulator at commensurate fillings and sufficiently strong interactions, and a superfluid elsewhere . The critical behavior is of two types: mean field for transitions induced by tuning the density, and of the $`(d+1)`$ dimensional XY universality class when the interaction strength is swept at fixed commensurate filling. One of the key new results of this paper will be that this special status of commensurate filling is lost in the case of a confining potential since commensuration itself well defined only locally..

We study the Hamiltonian

$`H`$ $`=`$ $`t{\displaystyle \underset{i}{}}(a_i^{}a_{i+1}+a_{i+1}^{}a_i)`$ (1)
$`+`$ $`V_0{\displaystyle \underset{i}{}}n_i(n_i1)+V_c{\displaystyle \underset{i}{}}(iL/2)^2n_i,`$ (2)
at zero temperature. Here $`t`$ measures the boson kinetic energy, $`V_0`$ the on-site repulsion, $`V_c`$ the curvature of the quadratic confining potential and $`L`$ the number of sites. In the presence of the trap, the value of $`L`$ should be chosen such that for the given trap curvature, the bosons do not see the edge of the system and therefore do not leak out. Our simulations were done with the world-line quantum Monte Carlo algorithm in the canonical ensemble . The chemical potential $`\mu =E/N`$ is obtained by differentiating numerically the energy with respect to the particle number . In the presence of a confining potential, it is important to measure the local density of bosons, $`n_i=a_i^{}a_i`$ as well as the local compressibility, $`\kappa _i=n_i/\mu _i=\beta [n_i^2n_i^2]`$.

Figure 1 shows the evolution of the local boson density with increasing total occupancy of the lattice. At low fillings the density profile is smooth, with an inverted parabolic shape reflecting the confining potential. Above a critical filling of about 30 bosons (for this choice of $`V_0`$ and $`V_c`$) a plateau with a local filling of one boson per site develops in the density profile which is analogous to the Mott structure of $`N`$ vs $`\mu `$ in the unconfined model. This plateau indicates the presence of an incompressible, insulating region where $`\kappa _i`$, as defined above, drops to a small but finite value (Fig. 2) which vanishes for $`V_0\mathrm{}`$. Here, this behavior of $`\kappa _i`$ will be taken as the signal for a Mott region. As the density is increased further, the plateau widens spatially. But when the energy cost of extending the plateau to increasingly large values of the confining potential becomes prohibitive, the occupancy begins to exceed one at the center of the lattice, indicating a breakdown of Mott behavior there, but not everywhere. Increasing the filling further, for example $`N_b=116`$, eventually produces a second Mott region in the center of the system with a local filling of two bosons per site without destroying totally the first Mott region. Four slices from Fig. 1 are shown in Fig. 2 along with the local compressibilities. It is clear that at higher boson numbers, richer structures where the local compressibility vanishes at several locally commensurate densities can occur. Mean-field work in two dimensions shows a similar coexistence of Mott and superfluid regions.

A central feature of the Mott phase transition of the unconfined boson-Hubbard model is global incompressibility: A charge gap opens up, i.e. the density gets “stuck” at $`\rho =\mathrm{integer}`$ for a range of chemical potentials $`\mu `$. One might, then, crudely interpret the spatial dependence of the local density in the confined case, Fig. 2, as rather analogous to the chemical potential dependence in the unconfined case . This assumption, and the $`\rho `$ versus $`\mu `$ curve in the unconfined case, allow us to calculate the site at which a Mott domain is entered or exited. These are shown as arrows in Fig. 2. However it is vital to emphasize that while the confined system has locally incompressible regions, the global compressibility is never zero which is seen clearly in Fig. 3. The main figure should be contrasted with the nonconfined case (inset).

An important difference in the behavior of the local $`\kappa _i`$ is especially evident in one dimension where, in the unconfined case, the global compressibility diverges as the Mott lobe is approached, $`\kappa |\rho 1|^1`$ for the first lobe. Here, instead, we find $`\kappa _i(n_i1)`$, as shown in Fig. 4. The origin of these differences is of course that the global compressibility, $`\kappa =\beta (_{ij}n_in_jn^2)`$, probes density correlations at all length scales. In the unconfined case, contrary to the confined system, the establishement of the Mott phase is a true quantum phase transition: It happens collectively throughout the system and the correlation length diverges. There is, however, an interesting “universality” in the trapped system as a Mott region is approached. That is, the values of $`\kappa _i`$ are the same even as the total filling and the on-site repulsion are varied (see Fig. 4). The same behavior is observed for the $`n=2`$ locally incompressible phase.

Sets of runs such as those shown in Fig. 1 allow us to determine the state diagram as a function of boson filling and interaction strength for a given trap curvature. This is shown in Fig. 5. Because of the absence of true phase transitions, we have referred to Fig. 5 as a state diagram rather than phase diagram. As the filling is increased at fixed interaction strength, one crosses from a smooth density profile to one which has locally incompressible “Mott” domains, if $`V_0`$ is large enough. Further increase in the filling ultimately leads to the formation of regions at the well center where $`n_i>1`$. In Fig. 5, region A admits only locally incompressible regions with $`n_i=1`$, as in Fig. 2(b) and (c). Region B has $`n_i>1`$ surrounded by incompressible regions, Fig. 2(d). Region C is where the central part of the system has an $`n_i=2`$ incompressible region which, when the boundary of the system is approached, falls off to a shoulder of $`n_i=1`$ “Mott” region before reaching zero density. The $`n_i=1`$ and $`n_i=2`$ incompressible regions are separated by compressible regions. Region D is where the center of the system is compressible $`n_i>2`$, bounded by $`n_i=2`$ which in turn is bounded by $`n_i=1`$ incompressible regions. Region E has no incompressible regions.

Note that the values of $`V_0`$ at which the A and C regions in Fig. 5 are entered are of the same order as those of the first two true Mott lobes in the nonconfined case . This is consistent with the experimental results on three-dimensional optical lattices which appear to be in agreement with the expected value in the nonconfined case. Furthermore, the narrowness of region C could help understand why the experiments have not shown signs of the $`n=2`$ Mott region, even though $`n_i=2.5`$ in the core of the system.

One of the interesting experimental results is how rapidly coherence is re-established when $`V_0`$ is suddenly reduced from a value large enough to have produced large incompressible regions. It was argued that the characteristic time is of the order of the tunneling (hopping) time between sites . This is entirely consistent with the picture we present here. If, for example, the density profile is as in Fig. 2(c) when $`V_0`$ is suddenly reduced to a very small value, the system will evolve to a profile like in 2(a) albeit with higher local density in the center. This is accomplished by particles near the edges hopping towards the center, and is greatly accelerated by the fact that the trap is much lower near the center than near edges. In addition, since there is only one center of nucleation (the geometric center of the system) there is no slowing down due to competition at domain walls where different nucleation zones meet.

In summary, we have discussed the nature of locally incompressible “Mott” insulating behavior in a one-dimensional system of interacting bosons in a confining potential at $`T=0`$. We conclude that because of the destruction of translation invariance, great care should be taken in drawing on the analogy with the unconfined case at a fundamental level. For one thing, it can support “Mott” behavior off commensurate fillings. While “Mott” regions still exist, the critical properties are completely altered: Incompressible regions are established in a very localized way and not at all critically in the usual sense. These localized regions grow or shrink with $`V_0`$, but are always in co-existence with other regions, some compressible and some which might be incompressible but at higher local integer filling. In that sense, the formation of these “Mott” regions is not a true quantum critical phenomenon as it is in the unconfined case.

While we have focussed on the new qualitative physics which results from the confining potential, it is important to emphasize that experiments on one and two dimensional trapped systems are currently underway . For these, our paper should provide specific quantitative predictions for the critical ratios of interaction strength to kinetic energy and trap curvature, as a function of density. We are currently undertaking these comparisons. In addition, in the absence of traps, the phase diagram is qualitatively the same in one and two dimensions . We expect this to be true in the confined case too. In fact, initial simulation results in two dimensions show this to be true.

Note Added: A QMC study of the three-dimensional system appeared recently in which the authors also conclude, like we do, that one cannot characterize globally the transitions discussed here. They discuss a signal that can be used experimentally to study the transition.

We gratefully acknowledge financial support from NSF DMR 9985978, a CNRS-NSF co-operative grant (France-USA), a PROCOPE grant (France-Germany), the Stichting FOM and the Swiss National Science Foundation. We thank R. Kaiser, C. Miniatura, T. Pfau and A.B.N. Sue for useful conversations.


text:
# Elastic behavior in Contact Dynamics of rigid particles

## I Introduction

In computational physics one can distinguish two different validation tasks, which have to be solved in order to make simulations a useful research tool. First one must prove the validity of a simulation model by comparing its results to laboratory experiments, and second, equally important, one must assess the systematic errors due to the practical implementation in order to tell, how precisely the simulation results reflect the theoretical properties of the simulation model. In this paper we address the second type of validation problem for the simulation technique of contact dynamics, which was developed about 10 years ago Jean and Moreau (1992); Moreau (1994); Jean (1999) with the aim to investigate granular media Herrmann et al. (1998) in the limit of high rigidity of the particles at a high packing density. This method has been successfully applied and reproduces experiments (see e.g. Daudon et al. (1997)). However, its systematic errors and the computational effort to keep them tolerably small have not been investigated in detail before.

This is in marked contrast to other discrete elements methods for granular media (cf. e.g. Wolf (1996)), in particular the soft particle *molecular dynamics* simulation model, which has been widely used since more than 20 years Cundall and Strack (1979). During this time possible pitfalls such as the *detachment effect*Luding et al. (1994) and the *brake failure effect*Schäfer and Wolf (1995) could be discovered, analyzed and hence avoided.

Contact dynamics simulations have been applied to study a large variety of questions in dense granular systems, where excluded volume interactions and static friction, so called unilateral constraints, are believed to be essential Radjai et al. (1996a, b, 1998, 1999); Nouguier et al. (2000), but it should be mentioned that such constraints arise also in other areas like virtual reality, engineering, especially in robotics, and operations research, where the numerical treatments are similar Lynch (1994); Stewart and Trinkle (1996).

In a system of perfectly rigid particles the sound velocity would be infinite. This is in principle borne out by the contact dynamics simulation model. However, its practical implementation will normally give rise to sound-waves in the granular material, as we are going to show in the following, even if each single collision is modeled as being perfectly inelastic. Our aim is to elucidate this artifact after introducing briefly the principles of this method.

## II The Contact Dynamics Method

First, let us point out the basic difference between molecular dynamics (MD) on one side and contact dynamics (CD) on the other. Both have in common the integration of Newton’s equation of motion where the occurring forces are due to external fields (gravity) or are – more important – *contact forces*, i.e. caused by contacts between particles or their contacts with confining walls.

The spirit of MD is to calculate the contact forces according to their *cause*, i.e. the (usually microscopic) deformation of the contact region and the involved velocities. Since the full treatment of every particle as a deformable body would render the simulation of a large number of particles exceedingly time consuming, a lot of models exist, how to replace the deformation by the local overlap Schäfer et al. (1996), the latter being a virtual quantity obtained from the undeformed shapes.

The principles of CD are different: Here the contact forces are calculated by virtue of their *effect*, which is to fulfill certain constraints. Typically, such a constraint is the volume exclusion of the particles or the absence of sliding due to static friction. As can be seen immediately, this problem cannot be solved locally: In a cluster of particles where many contacts are simultaneously present, the force at one contact depends on adjacent contacts and so on. In that case the aim is to find a global force-network, which is consistent with the constraints at all contacts. The method to carry out this calculation is often called the *solver* in this context, which is commonly an iterative scheme, as the one we describe in the following section. In order to make our points very clear, we perform an analytical investigation for a one dimensional example but can prove the existence of the discovered effect in two dimensions as well.

### II.1 The dynamical equations

As CD is designed to obey excluded volume constraints exactly, particle collisions lead to discontinuous velocity changes (“shocks”), i.e. to nonsmooth mechanics Moreau and Panagiotopoulos (1988). Therefore, higher order terms than employed in the Euler integration scheme are of no use. For the $`i`$th particle’s positions $`x_i`$, this reads

$$x_i(t+\mathrm{\Delta }t)=x_i(t)+v_i(t+\mathrm{\Delta }t)\mathrm{\Delta }t,$$
(1)
and correspondingly for its velocity $`v_i`$ we have

$$v_i(t+\mathrm{\Delta }t)=v_i(t)+\frac{F_i(t+\mathrm{\Delta }t)}{m}\mathrm{\Delta }t,$$
(2)
where $`F_i`$ is the total force acting on the particle, $`m`$ its mass and $`\mathrm{\Delta }t`$ is the time step.

A remark about the Euler scheme (2) being *implicit* is in order: Whereas in conventional MD the choice of evaluating the force for the previous or for the new configuration (i.e. at $`t`$ or $`t+\mathrm{\Delta }t`$ respectively) is merely a matter of stability of the integration scheme (cf. e.g. Press et al. (1993)), the constraints in CD can only be imposed on the *new* configuration. In this sense, the integration scheme of CD is inevitably of implicit type. (Taking into account the configurational change during a time step consistently (*fully implicit* integration Jean (1999)) leads to difficulties which are analogous to implicite schemes in MD, when the forces have to be evaluated for the yet unknown new configuration. In one dimension, though, these difficulties do not arise.)

#### II.1.1 One contact

We now turn to the force on the $`i`$’th particle, $`F_i`$, occurring in Eq. (2): In order to determine it one has to know the contact forces between the grains. In CD they are calculated from the condition that the constraints must not be violated. In one dimension and if one disregards rotations, this is simply the excluded volume contraint. To give a specific example, let us consider two particles with equal masses $`m`$ subjected to constant external forces (Fig. 1). A contact force $`R`$ (which is the reaction force due to the constraint) is active only, if interpenetration needs to be prevented. Otherwise, i.e. if the gap $`g`$ would remain non-negative (no overlap) anyway, it takes on its minimal value, $`R=0`$ (this is expressed in the CD literature as *Signorini’s condition*).

Without the repulsion $`R`$ between the particles the gap at the end of a time step would be given by

$$g^{}=g_t+(v_{2,t}v_{1,t})\mathrm{\Delta }t\frac{F_1^{\text{ext}}+F_2^{\text{ext}}}{m}\mathrm{\Delta }t^2$$
(3)
according to (1) and (2). However, the excluded volume constraint requires that $`g_{t+\mathrm{\Delta }t}=\mathrm{max}\{g^{},0\}`$. In order that this results from (1) and (2), the contact force $`R_{t+\mathrm{\Delta }t}=\mathrm{max}\{R^{},0\}`$ must be taken into account, where

$`R^{}`$ $`=`$ $`{\displaystyle \frac{mg^{}}{2\mathrm{\Delta }t^2}}`$ (4)
$`=`$ $`{\displaystyle \frac{m}{2\mathrm{\Delta }t}}\left({\displaystyle \frac{g_t}{\mathrm{\Delta }t}}+v_{1,t}v_{2,t}\right)+{\displaystyle \frac{F_1^{\text{ext}}+F_2^{\text{ext}}}{2}}`$ (5)
Note that this scheme corresponds to a completely inelastic collision (i.e. with the so-called *restitution coefficient* being zero), which is accomplished in two time steps: At first the gap closes, in the next step also the relative velocity vanishes. (Finite restitution coefficients can also be incorporated into this algorithm Radjai (1998).)

The above determination of the contact force has a drawback, though: If $`g_t<0`$ occurs due to a previous inaccuracy, then the *elimination* of this overlap is accompanied by a surplus of kinetic energy. Therefore, mostly the *quasi-inelastic shock* formula Jean (1999)

$$R^{}=\frac{m}{2\mathrm{\Delta }t}\left(\frac{g_t^{\text{pos}}}{\mathrm{\Delta }t}+v_{1,t}v_{2,t}\right)+\frac{F_1^{\text{ext}}+F_2^{\text{ext}}}{2}$$
(6)
is used instead of (5), where $`g_t^{\text{pos}}=\mathrm{max}\{0,g_t\}`$. That means, negative gaps are treated differently from Eq. (5) in such a way that an already existing overlap is not eliminated but only its further growth is inhibited. Hence the inelastic shock law (6) is in a way even “more inelastic” than the original law (5), because it avoids overlap correcting impulses which could destroy stable equilibrium states.

#### II.1.2 Many contacts

We now address the question how to solve the problem of the constraint forces if we consider not only one, but many contacts at the same time. Fig. 2 shows such a system, where $`R_i`$ denotes the contact force between the particles $`i`$ and $`i+1`$, and for the sake of simplicity we have no external force acting on them. Furthermore, we will concentrate on the situation where the neighboring pairs are permanently in contact (i.e. $`g_{i,t}^{\text{pos}}=0`$) which corresponds to compressed dense packings.

For the $`i`$th contact in this setup, $`R_{i1}`$ and $`R_{i+1}`$ play the role of the external forces in Eq. (6), though they are not constant, their values are not known for the next time step:

$$R_{i,t+\mathrm{\Delta }t}=\frac{m}{2\mathrm{\Delta }t}\left(v_{i,t}v_{i+1,t}\right)+\frac{R_{i1,t+\mathrm{\Delta }t}+R_{i+1,t+\mathrm{\Delta }t}}{2}.$$
(7)
Obviously the contact force is coupled to the neighboring contacts and through those to further contacts. In a similar way in higher dimensions, large numbers of contacts are coupled within clusters defined by the contact network. Hence, the determination of the proper reaction forces becomes a global problem.

A standard way in CD is to solve this problem iteratively: In one iteration step we calculate the forces according to the constraint conditions pretending that the corresponding neighboring contacts already exhibit the right forces. In that way the process traverses the list of contacts many times until satisfactory convergence is reached (for the question of the convergence cf. the works Jourdan et al. (1998); Jean and Pratt (1985)).

For our one-dimensional chain of particles, this means that Eq. (7) simply gets the meaning of an *assignment* of the right hand side to the left hand side, and one iteration step consists of applying this assignment sequentially once to each contact. The order of this sequence is preferably random (with the pattern changing for every sweep), in order not to create any bias in the information spreading111We call this procedure “random sweep” as opposed to “random sequential” where the choice of the contact to be updated would be random *and independent* of the previous choices.. With each sweep a globally consistent solution is approached, until finally a chosen convergence criterion is fulfilled.

In the next section we shall show, what kind of consequences the local update scheme presented here has on large time and length scales.

## III The large scale description

In order to analyze the coarse grained behavior of the microscopic equations derived in the last section we can also regard them as the discretized form of a continuum description, making a treatment in terms of partial differential equations (PDE) possible. In order to obtain the corresponding PDEs, we consider the particle index $`i`$ as space variable $`x`$ and replace the differences of consecutive quantities by derivatives, the error term for the first and second order derivatives being of first and second order, respectively. E.g. $`v_{t+\mathrm{\Delta }t}v_t\mathrm{\Delta }t_tv`$ and $`R_{i+1}+R_{i1}2R_id^2_x^2R`$, where $`d`$ is the particle diameter.

### III.1 The relaxation of the contact forces

While the continuum versions of the updates (1) and (2) can be obtained straight forwardly, the force change (7) lacks a time variable, for during the force-iteration, being just a calculation, no physical time passes. Hence, to be able to describe this force-development as well, let us introduce a fictitious time $`t^{}`$ with time interval $`\mathrm{\Delta }t^{}`$ for one iteration-step. With this, the continuum version of Eq. (7) reads:

$`_t^{}R`$ $`=D_x^2R\beta _xv`$ (8)
$`\text{with}D`$ $`=q{\displaystyle \frac{d^2}{\mathrm{\Delta }t^{}}}`$ (9)
$`\beta `$ $`=q{\displaystyle \frac{md}{\mathrm{\Delta }t\mathrm{\Delta }t^{}}}`$ (10)
$`\text{and}q`$ $`={\displaystyle \frac{1}{2}}`$ (11)
This analytic form clearly reveals the nature of the iteration loop: The reaction forces relax towards the solution in a diffusive way. (Note that the $`_xv`$ term is constant in $`t^{}`$, it only depends on $`x`$.)

The introduction of the constant $`q`$ reflects a subtlety regarding the sequential character of the update discussed in the appendix. In fact, since the PDE (8) describes the change of the whole field $`R(x)`$ *at once*, given its actual value at time $`t^{}`$, it corresponds to a *parallel update* (in the sense that the right hand side of Eq. (7) always employs the values $`R_i`$ from the *beginning* of the iteration sweep, not the freshly updated ones). In appendix A we shall show, though, that a random sweep update instead of a parallel one only renormalizes the value of $`q`$ to about $`0.8`$ while leaving the form of the PDE untouched.

### III.2 Sound waves

To connect the velocity update, whose continuum version follows immediately from Eq. (2) as

$$_tv=\frac{d}{m}_xR,$$
(12)
to the force update, we must relate the “iteration time” $`t^{}`$ to the physical time $`t`$. Although, depending on the convergence criterion, there can be in principle a varying number of iterations during one physical time step $`\mathrm{\Delta }t`$, we assume for simplicity this number $`N_I`$ being fixed. (Actually, in practice this crude “criterion” is sometimes applied.)

Hence, with $`\mathrm{\Delta }t=N_I\mathrm{\Delta }t^{}`$, we can express everything in terms of the physical time:

$`_tR`$ $`=D_x^2R\beta _xv,`$ (13)
$`\text{and}D`$ $`=qN_I{\displaystyle \frac{d^2}{\mathrm{\Delta }t}}`$ (14)
$`\beta `$ $`=qN_I{\displaystyle \frac{md}{\mathrm{\Delta }t^2}}`$ (15)
With the equations (12) and (13) we obtained two coupled PDEs. We can combine them to arrive at a wave equation with an additional damping term:

$$_t^2R=c^2_x^2R+_t\left(D_x^2R\right)$$
(16)
The sound velocity appearing is of finite value:

$$c=\sqrt{qN_I}\frac{d}{\mathrm{\Delta }t}$$
(17)
This equation indicates that the CD simulation of the particle chain, as presented in the previous section, can lead to sound propagation like in an elastic medium, which however contradicts the conception of perfect rigidity. The constraint conditions applied at the contacts should in principle prohibit overlaps, i.e. prohibit elastic deformation of the grains. It can be seen that this deviation from the perfect rigidity enters at the force relaxation: A finite number of iterations means a finite range for the information spreading and thus yields systematic errors in the calculated reaction forces. As a consequence, the finite $`N_I`$ involves soft particles and a finite sound velocity $`c\sqrt{N_I}`$. Note that in the limit of an infinite $`N_I`$, the exact value of the forces is reached, which corresponds to the case $`c\mathrm{}`$, as it should be for rigid particles.

#### III.2.1 Dispersion

Performing a Fourier transformation on Eq. (16), one obtains the properties of the different wave modes. The oscillation frequency $`\omega `$ of the wave number $`k`$ is

$$\omega \left(k\right)=k\sqrt{c^2\frac{D^2k^2}{4}}.$$
(18)
That means, $`\omega \left(k\right)`$ becomes zero at a critical wave number

$$k_c=\frac{2c}{D}\frac{1}{\sqrt{N_I}},$$
(19)
and waves with $`k`$ larger than $`k_c`$ (short wave lengths) are over-damped. The damping time $`\tau (k)`$ for the oscillating modes is given by:

$$\tau (k)=\frac{2}{Dk^2}$$
(20)
We derived the dispersion relation (18) in the continuum limit which is a good approximation for small wave numbers, but not near to the border of the Brillouin zone ($`k_{\text{Br}}=2\pi /d`$), where the effect of the spatial discreteness can be strong. However, increasing the number of the iterations sufficiently, $`k_c`$ becomes small compared to $`k_{\text{Br}}`$. Actually, for $`N_I10`$ the formula (18) works well not only for small wave numbers but for all the oscillating modes, as could be verified numerically.

#### III.2.2 Numerical confirmation

In order to confirm the results of this section, we performed the following numerical experiment: The starting configuration of the simulation consists of an array of $`50`$ discs and an immobile wall, the geometry can be seen in Fig. 3. Initially the gap between the wall and the leftmost particle is one disc diameter ($`d`$), the gap between the particles is zero and the array has zero velocity. Starting from $`t=0`$ a constant external force ($`F^{\text{ext}}`$) is acting on the rightmost particle which accelerates the array towards the wall (only horizontal motion takes place). As simulation parameters we chose $`N_I=40`$ and $`F^{\text{ext}}=0.05dm\mathrm{\Delta }t^2`$.

The collision with the wall induces a relative motion of the grains and generates sound waves in the array. After a transient period the grains remain permanently in contact (the whole array is pressed against the wall by $`F^{\text{ext}}`$). Since the different wave modes have different relaxation time, after a while only the largest wave length mode survives. This wave length is four times the system size because the wall represents a fixed boundary while the right side is free. Since the wave length is given, the oscillating frequency and the damping time can be calculated from Eq. (18) and Eq. (20), respectively. For comparison with the simulation we measured the motion of the rightmost particle. The expected motion is a damped oscillation

$$x(t)=x_0+A\mathrm{exp}\left(t/\tau \right)\mathrm{sin}\left(\omega t+\varphi \right),$$
(21)
where the offset $`x_0`$, the amplitude $`A`$ and the phase shift $`\varphi `$ have to be fitted (in contrast to $`\omega `$ and $`\tau `$) for a comparison. In Fig. 4 the measured data (dots) and the fitted curve can be seen. It shows that the simulation is in good agreement with our continuum description.

### III.3 Global Elasticity

It is instructive to compare our test-system to its simplest MD counterpart where the contact forces depend linearly on the local kinematic variables, i.e. the so called *linear spring/dashpot model*

$$R_i=\kappa \left(x_{i+1}x_id\right)\gamma \left(v_{i+1}v_i\right)$$
(22)
with the spring stiffness $`\kappa `$ and the damping coefficient $`\gamma `$. Employing again the updates (1) and (2) for the positions and velocities, respectively, the continuum limit yields the same type of PDE as Eq. (13) with its coefficients being inherited from Eq. (22):

$$_tR=\frac{\gamma d^2}{m}_x^2R\kappa d_xv$$
(23)
This allows us to relate the physical MD model parameters to the “technical” CD parameters:

$`\kappa `$ $`=qm{\displaystyle \frac{N_I}{\mathrm{\Delta }t^2}}`$ (24)
and
$`\gamma `$ $`=qm{\displaystyle \frac{N_I}{\mathrm{\Delta }t}}`$ (25)
This equivalence shows that on large scales the CD chain should behave identical to its MD counterpart, e.g. it will exhibit a global shrinkage proportional to an external compressive load. Note that a real congruence can be expected only for the collective behavior but not on the level of the contacts. In the CD method, as explained above, the contact forces are not related to the overlaps, which must merely be regarded as due to the incompleteness of the force-calculation and in fact are stochastic quantities because of our random update procedure. Only on scales larger than the grain size, where the fluctuations of these local “deformations” are averaged out, the behavior can be smooth like in an elastic medium, as is shown in Fig. 4.

In sections III.2 and III.3, our calculation was based on the assumption of a constant number of iterations for every time step, and due to this premise the analytical treatment became simple and directly comparable to the corresponding simulation. We should keep in mind, though, that the application of a convergence criterion involves a fluctuating $`N_I`$ (i.e. it can vary from time step to time step), and therefore steps with a different “stiffness” are mixed during the integration of motion. Consequently, the behavior of the CD method is more complex in detail, but qualitatively the results for the constant $`N_I`$ remain relevant also here. (E.g. the mechanism resulting in soft particles or the way how shock-waves can arise with finite velocity.)

## IV 2D simulation

After the analysis of the regular 1D system, the important question arises whether the behavior is similar in higher dimensions and less regular systems. Hence, we performed CD simulations with two-dimensional random packings of discs and observed the same “elastic” waves (even transversal modes were found).

The simulation presented here consists of $`1000`$ discs with radii distributed uniformly between $`r_{\text{min}}`$ and $`r_{\text{max}}=2r_{\text{min}}`$, the mass of each disc being proportional to its area. Fig. 5 shows the geometry: The base and the two side-walls are fixed while the upper piston is mobile.

Starting from a loose state, we compressed the system and waited until the packing reached an equilibrium state (the compression force $`F`$ applied on the piston was kept constant). The simulation was carried out without gravity and with a Coulomb friction coefficient of $`0.05`$ for all the disc-disc and disc-wall contacts (cf. Kadau et al. (2002)).

After the packing was relaxed completely, we generated sound waves by increasing the compression force abruptly to $`F+\mathrm{\Delta }F`$. After a transient period only one standing wave mode survives (both the wavenumber vector and the collective motion are vertical), where the piston, representing a free boundary, oscillates with a relatively large amplitude. We measured the vertical position of the piston versus time and found that the data can be fitted by an exponentially damped sine function (Fig. 6).

Here, in contrast to the 1D case, also $`\omega `$ and $`\tau `$ are fit parameters, since, due to the different geometry, the values (18) and (20) cannot be adopted, but, because of the random structure of the system, a more complex treatment is required for a quantitative description. However, we checked the most important relation, namely that the scaling properties of $`\omega `$ and $`\tau `$ remain valid also for the 2D random system, that is $`\omega \sqrt{N_I}`$ and $`\tau N_I^1`$, which means that the artificial visco-elasticity of the particles depends on the number of the iterations in the same manner as we showed for the 1D chain.

## V Discussion

The artificial elasticity found in CD simulations was analyzed. We showed that the systematic errors of the force calculation can lead to a collective elastic behavior, even though single contacts are assumed perfectly rigid and perfectly inelastic. For the 1D chain of particles, we could, starting from the “microscopic” laws, reproduce the numerical results analytically, including the dependence of the effective stiffness and viscous dissipation of the contacts on the computational parameters ($`N_I,\mathrm{\Delta }t`$).

Besides elucidating the origin of the elastic behaviour, the coarse grained description reveals important characteristics of the CD method, which were less obvious on the discrete level. We saw that using the iterative solver, the proper contact forces are approached in a diffusion like manner which is a crucial information concerning the computational time. The conception of perfectly rigid particles requires that the calculated forces are consistent even for contacts far from each other. Therefore, the “diffusion length” $`\sqrt{D\mathrm{\Delta }t}`$ must be larger than the linear system size $`L`$, which defines a lower boundary for the number of iterations of the order $`(L/d)^2`$. The same condition is obtained if we want to avoid other consequences of the effective softness. E.g. if we want the sound to travel a larger distance than $`L`$ during one time step (i.e. $`c>L/\mathrm{\Delta }t`$) or if we would like all possible wave modes to be overdamped (i.e. with $`\lambda _c>L`$ no oscillations). All these cases are equivalent, one is forced to apply a relatively large number of iterations: $`N_IL^2`$. Going further on this line, the scaling with respect to the number of particles $`n`$ can be determined: One step of the iteration consists of as many force-updates as there are contacts, which is proportional to $`n`$. Therefore, the computational effort of one time step scales with the particle number like $`nN_I`$, which is $`n^2`$ in 2D or $`n^{5/3}`$ in 3D. Therefore a large CD-simulation is computationally more costly than MD, where the computational effort scales like $`n`$. This is the price for simulating rigid particles without getting elasticity artifacts, which cannot be done with MD.

To avoid this super-linear scaling when dealing with large systems, we can also accept the finite stiffness by keeping $`N_I`$ constant independently of $`n`$. Then, besides gaining a running time of order $`n`$, of course, elastic deformations and sound waves can arise with an increasing number of the particles, and consequently they have to be monitored. We want to mention the idea, though, that in certain situations advantage can be taken of the artifact. E.g. when being applied deliberately, Coulombian friction can be combined with global elasticity easily; this way considerable computational time could be saved and even better performance than MD could be achieved.

## VI Acknowledgments

This work was supported by DAAD and Federal Mogul Technology GmbH. Partial support by OTKA T029985, T035028 is acknowledged.

## Appendix A $`q`$ for the random sweep update

To find the proper value of the constant $`q`$ appearing in sec. III, two more things have to be taken into account: Firstly the sequential type of the force update, and secondly that the order is random. The latter results in a stochastic force relaxation, i.e. the change of $`R_i`$ in one iteration sweep is a stochastic variable. In order to obtain a similar equation as (8) we shall determine the average value $`\mathrm{\Delta }R_i`$, where the average for site $`i`$ is meant to be taken over all possible update sequences.

Before going any further, let us introduce a few notations:

* For $`N`$ being the total number of contacts, the mapping $`u:\{1,\mathrm{},N\}\{1,\mathrm{},N\}`$ denotes the order of the update sequence; i.e. if the contact labeled $`i`$ is updated before $`j`$, then $`u_i<u_j`$.
* Throughout this appendix, the notation $`R_i`$ means the value from the beginning of the iteration sweep, the value at the end is $`R_i+\mathrm{\Delta }R_i`$.
* We define $`\delta R_i`$ as the change according to a parallel update (cf. Eq. (7)), i.e.

$`\delta R_i`$ $`={\displaystyle \frac{m}{2\mathrm{\Delta }t}}\left(v_iv_{i+1}\right)`$
$`+{\displaystyle \frac{R_{i1}+R_{i+1}2R_i}{2}},`$ (26)
as opposed to the *total* change $`\mathrm{\Delta }R_i`$.

Given that, it can easily be seen how $`\mathrm{\Delta }R_i`$ depends on the update order. If e.g. the site $`i`$ is updated earlier than its neighbors (i.e. $`u_{i1}>u_i<u_{i+1}`$), then $`\mathrm{\Delta }R_i=\delta R_i`$, but in the case $`u_{i1}>u_i>u_{i+1}<u_{i+2}`$, we get $`\mathrm{\Delta }R_i=\delta R_i+\delta R_{i+1}/2`$ (because for contact $`i`$, Eq. (7) employs the already updated force at contact $`i+1`$). Similarly (but with less probability) even very far contacts can contribute to $`\mathrm{\Delta }R_i`$, which can be summarized in the following way:

$$\mathrm{\Delta }R_i=\delta R_i+\underset{r}{}\frac{1}{2^r}\left(p_r\delta R_{ir}+p_r\delta R_{i+r}\right)$$
(27)
Here, $`p_r`$ is the probability that $`\mathrm{\Delta }R_i`$ contains information from the update of a contact at distance $`r`$, that is $`p_r=P(u_i>u_{i+1}>\mathrm{}>u_{i+r})`$. (This definition holds true for contributions from contacts with labels higher than $`i`$, but due to left-right symmetry the same value is inevitably obtained for the corresponding lower ones.)

The value of $`p_r`$ can be obtained from the following combinatorial consideration: Given an index $`i`$ and a distance $`r`$, we can classify the set of all update orders into groups such that the sequences in one group differ only in the permutations of the elements $`u_j`$, $`iji+r`$. Such a group contains $`(r+1)!`$ sequences, but only one of them satisfies the condition $`u_i>u_{i+1}>u_{i+2}>\mathrm{}>u_{i+r}`$. Since all update sequences are equally probable, the value of $`p_r`$ is equal to $`1/(r+1)!`$.

The factor $`(2^r(r+1)!)^1`$, relating the contacts $`i`$ and $`i+r`$, decays faster than exponentially; already for $`r=8`$, it drops below $`10^6`$. Therefore, the sum in Eq. (27) reaches only the immediate vicinity of contact $`i`$, such that, for our large wavelength considerations, the approximation $`\delta R_{i+r}\delta R_i`$ can be applied. This allows us to calculate the average change of the contact force:

$`\mathrm{\Delta }R_i`$ $`=\delta R_i\left(1+2{\displaystyle \underset{r=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{2^r(r+1)!}}\right)`$
$`=\delta R_i\left(4\sqrt{e}5\right).`$ (28)
Thus, it is shown that the random sweep results in a larger change of $`R_i`$ than the parallel update. Eq. (28) provides also the sought value of the parameter $`q`$ as

$$q=\frac{4\sqrt{e}5}{2}0.797,$$
(29)
which completes the continuum description given in sec. III.


text:
# REMARKS ON CP ASYMMETRIES IN 𝐷⁰/(𝐷⁰)̄→𝐾_𝑆𝜋⁺𝜋⁻

## Abstract

We consider the interference of resonant amplitudes leading to the final state $`K_S\pi ^+\pi ^{}`$ in $`D^0/\overline{D^0}`$ decays. Each of these amplitudes consists of both Cabibbo allowed and doubly Cabibbo suppressed transitions. The role of strong phase arising out of Breit-Wigner resonant propagators is emphasised. Invoking the $`\mathrm{\Delta }S=\mathrm{\Delta }Q`$ rule, $`K_S`$ in the final state is identified as the mass eigenstate of superposed weak eigenstates $`K^0`$ and $`\overline{K}^0`$. A nonzero CP asymmetry appears to be possible in three body decays.

preprint: BARI-TH/432-2002

As immense activities are on at charm/B factories searching for CP violating effects, the study of weak dynamics of charm sector acquires renewed interest. The CP asymmetry occurs if there exist nonvanishing weak and strong phases provided by a pair of amplitudes which are distinct by both phases with respect to one another, through interference . For instance, in the two body $`D^0`$ decays into $`K^\pm \pi ^{}`$, the strong phase arises due to strong interaction effects such as rescattering and final state interaction. But, such a strong phase is small, for example it is about $`13^{}`$ in the model of Buccella et al , thus resulting in a CP asymmetry of $`O(10^3)`$.

Alternatively, one can look at three body modes. The three body final state is reached by more than one way, namely, nonresonant and resonant modes. There are many resonant modes: That is, $`D^0M_1RM_1[M_2M_3]_R`$, where $`R`$ stands for resonance and $`M`$’s for mesons. Given two resonant modes that lead to the same final state, there arise distinct strong phase between them due to the Breit-Wigner (BW) resonant propagator and the angular momentum quantum numbers of the resonance and of its decay products . Thus, a strong phase is nontrivial.

Besides, the correlation between the net strangeness (produced)In the strange decays, this is change in strangeness. and the net charge (involved), namely, the $`\mathrm{\Delta }S=\mathrm{\Delta }Q`$ rule. Unlike in strange decays where violation strangeness-charge symmetry means the mixing of $`K_S`$ and $`K_L`$ , in charm decays into strange final state, this rule qualifies $`K_S`$ is the mixed state of the weak eigen states $`K^0`$ and $`\overline{K}^0`$. This rule implies the coherent superposition of CA and DCS transitions.

In this letter, we look at the significance of strong phase that arises due to BW propagators of intermediate resonances which lead to a common final state in $`D^0`$ three body decays, as an extension of the Atwood and Soni’s proposal for $`B`$ decays to charm sector and the importance of invoking the $`\mathrm{\Delta }S=\mathrm{\Delta }Q`$ rule thereof. For $`D^0K_S\pi ^+\pi ^{}`$, we consider the resonant modes of $`K^{}(890)^\pm ,\rho ^0(770)`$ and $`f_0(990)`$:

$`D^0`$ $``$ $`K^{}\pi ^+\overline{K}^0\pi ^+\pi ^{}`$ (1)
$`D^0`$ $``$ $`\overline{K}^0\rho ^0\overline{K}^0\pi ^+\pi ^{}`$ (2)
$`D^0`$ $``$ $`\overline{K}^0f_0\overline{K}^0\pi ^+\pi ^{}`$ (3)
These resonances show up in the Dalitz plot of BABAR data .

Let us write the amplitude for the resonant decay $`D^0K_S\pi ^+\pi ^{}`$ mediated by a resonance $`i`$ as

$$M_i=A_i\mathrm{\Pi }_iB_i$$
(4)
where $`A_i`$ is the weak part of the amplitude containing the CKM phase, $`\mathrm{\Pi }_i=[sm_i^2+i\mathrm{\Gamma }_im_i]^1`$ the BW propagator that provides the strong phase and the strong coupling $`B_i=(16\pi m_i^3\mathrm{\Gamma }_i/\lambda ^{1/2}(m_i^2,m_1^2,m_2^2))^{1/2}`$. In order to extract the strong phase, we rewrite $`\mathrm{\Pi }_i`$ as

$$\mathrm{\Pi }_i=\stackrel{~}{\mathrm{\Pi }}_ie^{i\delta _i}$$
(5)
The strong phase $`\delta `$ is given by the width and mass of the resonance.

We take into account four amplitudes due to them: both the CA and DCS transitions of each. The amplitude for $`D^0K_S\pi ^+\pi ^{}`$ is then expressed after factoring out the strong and weak phases as

$`M`$ $`=`$ $`e^{i\gamma _1}\left[e^{i\delta _i}M_i^C+e^{i\delta _j}M_j^C\right]+e^{i\gamma _2}\left[e^{i\delta _i}M_i^S+e^{i\delta _j}M_j^S\right]`$ (6)
$`\overline{M}`$ $`=`$ $`e^{i\gamma _1}\left[e^{i\delta _i}\overline{M}_i^C+e^{i\delta _j}\overline{M}_j^C\right]+e^{i\gamma _2}\left[e^{i\delta _i}\overline{M}_i^S+e^{i\delta _j}\overline{M}_j^S\right]`$ (7)
where $`\gamma _{1,2}`$ stands for the CA and DCS weak phase, $`\delta _{i,j}`$ the strong phase and the superscript $`C`$ and $`S`$ for CA and DCS transitions. The amplitude $`|M|`$ is resultant of coherent superposition of the four amplitudes. As notation, $`K^{}`$ is identified with $`i`$ and $`\rho ^0`$ and $`f_0`$ with $`j`$.

Then we have the asymmetry as

$$a_{CP}=\frac{2W\mathrm{sin}\mathrm{\Delta }\mathrm{sin}\varphi }{1+X+4Y\mathrm{cos}\mathrm{\Delta }+2Z\mathrm{cos}\varphi +2W\mathrm{cos}\mathrm{\Delta }\mathrm{cos}\varphi }$$
(8)
where $`\varphi =|\gamma _1\gamma _2|`$ is the weak phase and $`\mathrm{\Delta }=|\delta _i\delta _j|`$ the strong phase and

$$W=R_3+R_1R_2,X=R_1^2+R_2^2+R_3^2,Y=R_1+R_2R_3,Z=R_2+R_1R_3$$
(9)
with $`R_1,R_2`$ and $`R_3`$ respectively the ratio of $`M_j^C,M_i^S`$ and $`M_j^S`$ with respect to $`M_i^C`$:

$`R_1`$ $`=`$ $`{\displaystyle \frac{f_\rho \overline{K}^0|(\overline{s}c)_{VA}|D^0}{f_\pi K^{}|(\overline{s}c)_{VA}|D^0}}{\displaystyle \frac{\widehat{\mathrm{\Pi }}_\rho B_\rho }{\widehat{\mathrm{\Pi }}_K^{}B_K^{}}},R_2={\displaystyle \frac{f_K^{}\pi |(\overline{d}c)_{VA}|D^0}{f_\pi K^{}|(\overline{s}c)_{VA}|D^0}},`$ (10)
$`R_3`$ $`=`$ $`{\displaystyle \frac{f_{K^0}\rho ^0|(\overline{d}c)_{VA}|D^0}{f_\pi K^{}|(\overline{s}c)_{VA}|D^0}}{\displaystyle \frac{\widehat{\mathrm{\Pi }}_\rho B_\rho }{\widehat{\mathrm{\Pi }}_K^{}B_K^{}}}`$ (11)
where in (11) factorisation approximation is applied for the weak amplitude. Similarly, $`R`$’s are for $`f_0`$. We note caution that not much is known learly about $`f_0`$. However, we have treated this on par with $`\rho ^0`$. The width of these resonances are 50, 150 and 50-100 MeV respectively of $`K^{}(890)`$, $`\rho ^0(770)`$ and $`f_0(980)`$. In this note, for calculational purpose, the width of $`f_0`$ is chosen as 75 MeV.

The strong phase $`\mathrm{\Delta }`$ is determined as a function of $`s`$. As an order of magnitude, the weak phase $`\varphi `$ is chosen as $`0.4\times 10^3`$ as $`\varphi arg(V_{cd}V_{us}^{}/V_{cs}V_{ud}^{})`$. The ratios of the weak matrix elements in (11) are assumed to be $`O(1)`$ as they may turn out to be so in the $`SU(3)`$ limit. The CP asymmetry is shown in Fig. (1) as a function of $`s`$. The strong phase $`\mathrm{\Delta }`$ is about $`57^{}`$ for $`K^{}`$ and $`\rho ^0`$ and about $`64^{}`$ for $`K^{}`$ and $`f_0`$ at $`s=m_K^{}^2`$. At $`s=m_K^{}^2`$, the amplitudes are in phase with respect to one another, yielding constructive interference of each pair leading to the same final state. The minima of $`a_{CP}`$ for $`K^{}\rho `$ and $`K^{}f_0`$ occurs respectively at $`s=0.75`$ and $`0.85`$ GeV2 respectively, whereas the maxima are respectively at $`m_\rho ^2`$ and $`m_{f_0}^2`$.

$`\mathrm{\Delta }S=\mathrm{\Delta }Q`$ rule: At quark level, $`D^0K_S\pi ^+\pi ^{}`$ proceeds via $`cs\overline{d}u(q\overline{q})`$, internal $`W`$-emission, and $`cd\overline{s}u(q\overline{q})`$, external $`W`$-emission, corresponding to the CA and DCS transitions respectively. While the net charge $`\mathrm{\Delta }Q`$ is zero in both, the strangeness production $`\mathrm{\Delta }S`$ is $`1`$ and $`+1`$ respectively. Thus, for $`D^0K_S\pi ^+\pi ^{}`$, the $`\mathrm{\Delta }S=\mathrm{\Delta }Q`$ is satisfied by (6) if and only if the mass eigenstate $`K_S`$ in the final state is the superposition of the weak eigenstates $`K^0`$ and $`\overline{K}^0`$. In other words, in order to reach this final state, the coherent superposition of the CA and DCS transitions is the basic requirement.

In the case of $`D^\pm K_S\pi ^\pm \pi ^0`$, the intemediate resonances play similar role in giving rise to a strong phase. As noted in , the mixing of final state neutral $`K`$ mesons brings in an additional phase factor of $`3.3\times 10^3`$. The role of this extra contribution to the CP asymmetry would eventually be significant.

To conclude, we considered the effects arising out of the BW propagators. These effects are large and thus lead to a possible large CP asymmetry. Application of the strangeness-charge symmetry rule reveals the structure of $`K_S`$ in the final state and is expected to shed more light on the possible new physics which would be responsible for a reasonable CP violating effects. A precise determination of CP aymmetry depends on the contribution of angular momenta of the resonces and of their decay products.

###### Acknowledgements.

The author is greatly indebted to P. Colangelo and A. Palano for useful discussions, S. Pakvasa for fruitful comments and P. Majewski for helps.


text:
# Reentrant Melting of Soliton Lattice Phase in Bilayer Quantum Hall System

## Abstract

At large parallel magnetic field $`B_{}`$, the ground state of bilayer quantum Hall system forms uniform soliton lattice phase. The soliton lattice will melt due to the proliferation of unbound dislocations at certain finite temperature leading to the Kosterlitz-Thouless (KT) melting. We calculate the KT phase boundary by numerically solving the newly developed set of Bethe ansatz equations, which fully take into account the thermal fluctuations of soliton walls. We predict that within certain ranges of $`B_{}`$, the soliton lattice will melt at $`T_{\mathrm{KT}}`$. Interestingly enough, as temperature decreases, it melts at certain temperature lower than $`T_{\mathrm{KT}}`$ exhibiting the reentrant behaviour of the soliton liquid phase.

When the interlayer spacing $`d`$ is comparable to the mean particle distance, the bilayer quantum Hall system at total filling factor $`\nu _{\mathrm{tot}}=1`$ can exhibit the quantum Hall effect due to the strong interlayer correlations. In the presence of large parallel magnetic field $`B_{}`$, the ground state of the system is known to form a uniform soliton lattice (SL) phase made of the periodic array of the phase solitons of the field variable $`\theta (𝐫)`$, which represents the relative phase difference between electrons in two layers. Using the isospin language , the effective energy functional for the bilayer system can be written by

$$E(\theta )=d^2r\left\{\frac{1}{2}\rho _s|\theta |^2\frac{t}{2\pi \mathrm{}^2}\mathrm{cos}[\theta (𝐫)+Qx]\right\},$$
(1)
where $`\rho _s`$ is the isospin stiffness, $`t=t_0e^{Q^2\mathrm{}^2/4}`$ the interlayer tunneling amplitude, $`\mathrm{}=(\mathrm{}c/|e|B_{})^{1/2}`$, and $`Q=2\pi /L_{}`$ with $`L_{}=\mathrm{\Phi }_0/B_{}d`$, which defines the length associated with one flux quantum $`\mathrm{\Phi }_0`$ enclosed between two layers. We set the magnetic length $`\mathrm{}=1`$. This continuum model is valid, when the field $`\theta (𝐫)`$ varies smoothly over the lengths $`\mathrm{},L_{}`$. With the increase of $`B_{}`$, the system exhibits a quantum phase transition from commensurate (C) to incommensurate (I) phase at $`Q_c=(2/\pi )(2t/\pi \rho _s)^{1/2}`$ at $`T=0`$. The incommensurate phase at $`Q>Q_c`$ describes the uniform SL phase at zero temperature.

In the paper, we study the thermodynamics of the SL phase through a mapping of the 2D statistical mechanical model in Eq. (1) to the 1D quantum sine-Gordon (QSG) Hamiltonian. When one neglects the compactness of the field variable $`\theta `$, which will be reasonably valid for $`T<(\pi /2)\rho _s`$, the 1D QSG model is exactly soluble by the Bethe ansatz method. The Bethe ansatz solution provides the thermodynamic CI phase boundary $`Q_c(T)`$. Depending on the values of $`C=t/(32\pi \rho _s)`$, the critical value $`Q_c(T)`$ will increase (decrease) with temperature for $`C>1(C<1)`$. Within the incommensurate phase, the ground state of the system is the uniform SL phase at $`T=0`$. As $`T`$ increases to about $`(\pi /2)\rho _s`$, the compactness of the angle variable $`\theta `$ will become important, which introduces dislocations to the system leading to the KT melting of the SL phase. Recently Hanna et al. have studied the melting of the SL phase by analyzing the elastic moduli based on the rigidity of the zero temperature ground state. However, the renormalizaion of the elastic moduli due to thermal fluctuations of soliton walls has not been taken into account. Furthermore, as one approaches to the CI phase boundary, the soliton density will vanish. In order to understand the melting of the SL phase near the CI phase boundary, it is crucial to include the thermal fluctuations and the soliton density variations as well.

In the paper, we develop a new set of Bethe ansatz equations, which fully take into account both the thermal fluctuations of the soliton walls and the density variations. By numerically solving the Bethe ansatz equations, we have calculated both the compression modulus $`\kappa _{xx}`$ and the shear modulus $`\kappa _{yy}`$ of the SL phase. By making an asymptotic expansion near the CI phase boundary, we have shown that $`\kappa _{xx}`$ goes like $`(2\tau /M)^{1/2}\left(QQ_c(T)\right)^{1/2}\rho _s`$ and $`\kappa _{yy}(8M\tau ^3)^{1/2}\left(QQ_c(T)\right)^{1/2}\rho _s`$, where $`M`$ is the soliton mass which varies with $`\tau `$ and $`\tau =T/(8\pi \rho _s)`$. This asymptotic behaviour explicitly confirms the predictions made by Coppersmith et al.. Based on the elastic moduli thus obtained, $`T_{\mathrm{KT}}`$ is calculated and compared with the zero temperature estimates. We predict the following reentrant behaviour of the soliton liquid phase. At certain values of $`Q`$ slightly below $`Q_c(0)`$, the system initially stays at the C phase. As $`T`$ increases, it makes a transition to the soliton liquid phase. With the further increase of $`T`$, the soliton liquid solidifies to the SL phase due to the rapid increase of the soliton density, which subsequently melts reentering the soliton liquid phase.

Following the effective energy functional of Eq. (1), the low temperature thermodynamics of the system can be defined by the statistical partition function $`𝒵=[𝒟\theta ]e^{E(\theta )/T}`$. One can map the statistical ensemble summation into quantum transfer matrix by identifying $`𝒵=lim_R\mathrm{}\mathrm{Tr}\left[e^{R\widehat{H}}\right]`$,

$$\widehat{H}=𝑑x\left[\frac{T}{2\rho _s}\mathrm{\Pi }^2+\frac{\rho _s}{2T}(_x\theta )^2\frac{t}{2\pi T}\mathrm{cos}\left(\theta +Qx\right)\right].$$
(2)
If one neglects the compactness of the angle variable $`\theta (𝐫)`$, the canonical conjugate momentum $`\mathrm{\Pi }(𝐫)`$ can be easily integrated out leading to the corresponding classical partition function. We make the following change of variables $`\varphi (𝐫)=\sqrt{\rho _s/T}\left[\theta (𝐫)+Qx\right]`$, which leads to the quantum Hamiltonian of the 1D sine-Gordon model with external $`U(1)`$ coupling,

$$\widehat{H}=𝑑x\left[\frac{1}{2}\mathrm{\Pi }^2+\frac{1}{2}(_x\varphi )^2A\frac{\beta }{2\pi }_x\varphi +2\mu \mathrm{cos}\left(\beta \varphi \right)+\frac{1}{2T}\rho _sQ^2\right].$$
(3)
Here $`\mu =t/(4\pi T),\beta =\sqrt{T/\rho _s},\mathrm{and}A=2\pi Q\rho _s/T`$. The compactness of $`\theta (𝐫)`$ would restrict $`\mathrm{\Pi }(𝐫)`$ to the integer values. Hence the ground state energy $``$ of the 1D QSG model corresponds to $`/T`$, where $``$ is the free energy of the 2D statistical mechanical system.

The ground state of the 1D QSG model is exactly soluble by the Bethe ansatz method. There is a competition between the finite soliton mass $`M`$ which prefers the commensurate phase and the external field $`A`$ coupled to the topological charge which favors the incommensurate phase. This is a ‘quantum’ version of the CI transition, where the soliton mass includes full ‘quantum’ fluctuation effects. With external $`U(1)`$ field, a soliton can be created with energy cost $`\mathrm{\Delta }=M\mathrm{cosh}\mathrm{\Theta }A`$. The soliton mass $`M`$ is given by

$$M=\frac{2\mathrm{\Gamma }\left(\frac{p}{2}\right)}{\sqrt{\pi }\mathrm{\Gamma }\left(\frac{p+1}{2}\right)}\left[\pi \mu \frac{\mathrm{\Gamma }\left(\frac{1}{p+1}\right)}{\mathrm{\Gamma }\left(\frac{p}{p+1}\right)}\right]^{\frac{1}{2}(p+1)},$$
(4)
where $`p=T/(8\pi \rho _sT)`$ and $`\mathrm{\Theta }`$ represents the rapidity of the soliton.

For $`A<M`$, the ground state will be the vacuum: the commensurate phase. However for $`A>M`$, a nontrivial vacuum will arise, since the energy cost to create a soliton can be negative. In this case, the ground state will be described by the soliton condensations: The SL phase. The CI transition will occur at $`A=M`$. By equating the two quantities, one can obtain the exact CI phase boundary $`Q_c(T)`$

$$Q_c(T)=\frac{8}{\sqrt{\pi }}\frac{\tau \mathrm{\Gamma }\left(\frac{\tau }{2(1\tau )}\right)}{\mathrm{\Gamma }\left(\frac{1}{2(1\tau )}\right)}\left[C\frac{\mathrm{\Gamma }(1\tau )}{\mathrm{\Gamma }(1+\tau )}\right]^{\frac{1}{2(1\tau )}},$$
(5)
where we have neglected the $`Q`$-dependence of $`t`$, since the continuum model is only valid for $`Q1`$. The parameter $`C`$ is given by $`t/(4T_0)`$. The reduced temperature variable $`\tau `$ is between $`0`$ and $`1`$. As $`\tau 0`$, one can reproduce the classical limit: $`Q_c(0)=16\sqrt{C}/\pi `$.

Now we want to study the melting of the SL phase by numerically solving the Bethe ansatz equations. Using the fact that the exact scattering matrix is known for the 1D QSG model, one can calculate the ground state energy and the soliton density $`\overline{n}_s(Q,T)`$ as a function of $`Q`$ and $`T`$. The Bethe ansatz equation is given by

$$2\pi \rho (\mathrm{\Theta })=M\mathrm{cosh}\mathrm{\Theta }+_B^B𝑑\mathrm{\Theta }^{}\phi (\mathrm{\Theta }\mathrm{\Theta }^{})\rho (\mathrm{\Theta }^{}),$$
(6)
where $`\rho (\mathrm{\Theta })`$ is the density of solitons between $`\mathrm{\Theta }`$ and $`\mathrm{\Theta }+d\mathrm{\Theta }`$ per unit length and the integral kernel $`\phi (\mathrm{\Theta })`$ is given by

$$\phi (\mathrm{\Theta })=_{\mathrm{}}^{\mathrm{}}𝑑\omega e^{i\omega \mathrm{\Theta }}\frac{\mathrm{sinh}\left(\frac{\pi (p1)\omega }{2}\right)}{2\mathrm{cosh}\left(\frac{\pi \omega }{2}\right)\mathrm{sinh}\left(\frac{\pi p\omega }{2}\right)}.$$
(7)
In terms of the density function, one can calculate the ground state energy as follows

$$\mathrm{\Delta }=_B^B𝑑\mathrm{\Theta }(M\mathrm{cosh}\mathrm{\Theta }A)\rho (\mathrm{\Theta }),$$
(8)
where the ground state is determined by imposing the condition: $`\mathrm{\Delta }/B=0`$. By making an asymptotic expansion of Eq. (6) near the CI phase boundary, one can confirm that the soliton density exhibits a power-law behaviour: $`\overline{n}_s(2M/\tau )^{1/2}\left(QQ_c(T)\right)^{1/2}/(2\pi )`$.

As $`T`$ increases to about $`(\pi /2)\rho _s`$, the compactness of the variable $`\theta (𝐫)`$ becomes important, which introduces topological defects (dislocations) into the system leading to the KT melting of the SL phase at temperatures much below the CI transition temperature. Although the 1D QSG Hamiltonian does not involve dislocation excitations, one can get a reasonable estimate for $`T_{\mathrm{KT}}`$ by analyzing the elastic moduli of the SL at finite temperature. The elastic moduli $`\kappa _{ij}`$ are defined as follows

$$\kappa _{ij}=\frac{^2}{Q_s^iQ_s^j},$$
(9)
where $`Q_s^x`$ is given by $`2\pi /L_s`$, where $`L_s`$ is the lattice constant of the SL phase along the $`\widehat{x}`$ direction, $`Q_s^y=\mathrm{tan}\varphi Q_s^x`$ with $`\varphi `$ the tilt angle of the soliton walls, and $`i,j=x,y`$. In order to obtain the elastic moduli of the system, one needs to calculate the energy of the system by varying the soliton density from the ground state value and/or tilting the solitons from the vertical orientations. We first calculate the compression modulus $`\kappa _{xx}`$. By shifting the fermi momentum from $`B`$ to $`B+ϵ`$, one can increase (decrease) the soliton density for positive (negative) values of $`ϵ`$. The soliton density $`\rho (\mathrm{\Theta })`$ is a function of both $`ϵ`$ and $`\mathrm{\Theta }`$. For small $`ϵ`$, one can expand the soliton density up to quadratic order in $`ϵ`$

$$\rho (\mathrm{\Theta })\overline{\rho }(\mathrm{\Theta })+ϵ\rho _b^{}(\mathrm{\Theta })+\frac{ϵ^2}{2}\rho _b^{\prime \prime }(\mathrm{\Theta }),$$
(10)
where $`\overline{\rho }(\mathrm{\Theta })`$ is the soliton density profile of the ground state, $`\rho _b^{}=\rho /ϵ|_{ϵ=0}`$, and $`\rho _b^{\prime \prime }=^2\rho /ϵ^2|_{ϵ=0}`$. At the lowest order of $`ϵ`$, we obtain the usual Bethe ansatz equation for the ground state as shown in Eq.(6) with the condition $`\mathrm{\Delta }/ϵ=0`$. Up to the order of $`ϵ^2`$, we obtain the following two equations

$$2\pi \rho _b^{}(\mathrm{\Theta })=[\phi (\mathrm{\Theta }B)+\phi (\mathrm{\Theta }+B)]\overline{\rho }(B)+_B^B𝑑\mathrm{\Theta }^{}\phi (\mathrm{\Theta }\mathrm{\Theta }^{})\rho _b^{}(\mathrm{\Theta }^{}),$$
(11)
$`2\pi \rho _b^{\prime \prime }(\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \frac{}{\mathrm{\Theta }}}\left[\phi (\mathrm{\Theta }+B)\phi (\mathrm{\Theta }B)\right]\overline{\rho }(B)+[\phi (\mathrm{\Theta }B)+\phi (\mathrm{\Theta }+B)]\left[{\displaystyle \frac{\overline{\rho }(B)}{B}}+2\rho _b^{}(B)\right]`$ (12)
$`+`$ $`{\displaystyle _B^B}𝑑\mathrm{\Theta }^{}\phi (\mathrm{\Theta }\mathrm{\Theta }^{})\rho _b^{\prime \prime }(\mathrm{\Theta }^{}).`$ (13)
The energy of the system is given as follows: $`\mathrm{\Delta }=\mathrm{\Delta }\overline{}+\frac{ϵ^2}{2}\chi _b`$, where $`\mathrm{\Delta }\overline{}`$ is the ground state energy and $`\chi _b`$ is given by

$$\chi _b=2M\mathrm{sinh}B\overline{\rho }(B)+2(M\mathrm{cosh}BA)\left[\frac{\overline{\rho }(B)}{B}+2\rho _b^{}(B)\right]+_B^B𝑑\mathrm{\Theta }(M\mathrm{cosh}\mathrm{\Theta }A)\rho _b^{\prime \prime }(\mathrm{\Theta }).$$
(14)
The $`Q_s^x`$ can be calculated as follows: $`Q_s^x=2\pi _{Bϵ}^{B+ϵ}𝑑\mathrm{\Theta }\rho (\mathrm{\Theta })\overline{Q}_s^x+\alpha _bϵ+𝒪(ϵ^2)`$, where $`\overline{Q}_s^x=2\pi _B^B𝑑\mathrm{\Theta }\overline{\rho }(\mathrm{\Theta })`$ is the ground state value and $`\alpha _b=4\pi \overline{\rho }(B)+2\pi _B^B𝑑\mathrm{\Theta }\rho _b^{}(\mathrm{\Theta })`$. The compression modulus $`\kappa _{xx}`$ is given as follows: $`^2/Q_{s}^{x}{}_{}{}^{2}=T\chi _b/\alpha _b^2`$.

Now we calculate the shear modulus $`\kappa _{yy}`$. To obtain the shear modulus, one needs to fix $`Q_s^x`$ and vary $`Q_s^y`$ alone. Global shift of the rapidity will not only tilt the soliton walls but also change the soliton density along the $`\widehat{x}`$-direction. In order to vary $`Q_s^y`$ alone, one needs to change the fermi momentum $`B`$ to $`B+ϵ`$ and rapidity by $`\eta `$ simultaneously along the following trajectory: $`ϵ=\alpha _\mathrm{T}\eta ^2`$ for some constant $`\alpha _\mathrm{T}`$ . By making a perturbative expansion in $`\eta `$, we obtain the following two Bethe ansatz equations up to $`\eta ^2`$,

$$2\pi \rho _s^{}(\mathrm{\Theta })=[\phi (\mathrm{\Theta }B)\phi (\mathrm{\Theta }+B)]\overline{\rho }(B)+_B^B𝑑\mathrm{\Theta }^{}\phi (\mathrm{\Theta }\mathrm{\Theta }^{})\rho _s^{}(\mathrm{\Theta }^{}),$$
(15)
$`2\pi \rho _s^{\prime \prime }(\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \frac{}{\mathrm{\Theta }}}\left[\phi (\mathrm{\Theta }+B)\phi (\mathrm{\Theta }B)\right]\overline{\rho }(B)+\left(\phi (\mathrm{\Theta }B)+\phi (\mathrm{\Theta }+B)\right)\left[{\displaystyle \frac{\overline{\rho }(B)}{B}}+2\rho _s^{}(B)\right]`$ (16)
$`+`$ $`{\displaystyle _B^B}𝑑\mathrm{\Theta }^{}\phi (\mathrm{\Theta }\mathrm{\Theta }^{})\rho _s^{\prime \prime }(\mathrm{\Theta }^{}),`$ (17)
where one can easily notice that $`\rho _s^{}(\mathrm{\Theta })`$ is an odd function of $`\mathrm{\Theta }`$ and $`\rho _s^{\prime \prime }(\mathrm{\Theta })`$ an even function. The energy of the system along the trajectory is given by

$$\mathrm{\Delta }=\mathrm{\Delta }\overline{}+\frac{\eta ^2}{2}\chi _s,$$
(18)
where $`\chi _s`$ is given by

$$\chi _s=2M\mathrm{sinh}B\overline{\rho }(B)+2(M\mathrm{cosh}BA)\left[\frac{\overline{\rho }(B)}{B}+2\rho _s^{}(B)\right]+_B^B𝑑\mathrm{\Theta }(M\mathrm{cosh}\mathrm{\Theta }A)\rho _s^{\prime \prime }(\mathrm{\Theta }).$$
(19)
Since $`Q_s^y=\mathrm{tan}\varphi \overline{Q}_s^x=\eta \overline{Q}_s^x`$, the shear modulus $`\kappa _{yy}`$ is given by the following relation: $`^2/Q_{s}^{y}{}_{}{}^{2}=T\chi _s/(\overline{Q}_s^x)^2`$. The off-diagonal elements of $`\kappa _{ij}`$ are zero. Using the asymptotic expansion, one can show that near the CI phase boundary, $`\kappa _{xx}`$ goes like $`(2\tau /M)^{1/2}\left(QQ_c(T)\right)^{1/2}\rho _s`$, and $`\kappa _{yy}(8M\tau ^3)^{1/2}\left(QQ_c(T)\right)^{1/2}\rho _s`$.

In Murphy et al.’s experiment, for equally populated layers, the experimental values of the various parameters are given as follows: $`\rho _s0.35K,t_01.2K,\mathrm{}126\AA `$, and $`d200\AA `$. For the above sample, the parameter $`C`$ can be estimated to be about $`0.033`$. In the inset of Fig. (1), the open circles represent $`\kappa _{xx}`$ and the open squares $`\kappa _{yy}`$ at fixed value of $`Q0.919`$, which obey the correct asymptotic behaviour near the CI phase boundary. We have confirmed that both $`\kappa _{xx}`$ and $`\kappa _{yy}`$ approach to $`\rho _s`$ at large $`Q`$. The KT transition temperature $`T_{\mathrm{KT}}`$ can be estimated by solving the following equation: $`k_BT/(\pi /2)\rho _s=[\kappa _{xx}(T)\kappa _{yy}(T)]^{1/2}/\rho _s`$. In Fig. (1), the closed circles represent $`(\kappa _{xx}(T)\kappa _{yy}(T))^{1/2}/(16\rho _s)`$ as a function of $`\tau `$. The intersections with the solid line with slope $`1`$ locate the positions of $`T_{\mathrm{KT}}`$. It can be shown that for $`Q_{\mathrm{min}}<Q<Q_c(0)`$, the equations have two solutions, where $`Q_{\mathrm{min}}0.9,Q_c(0)0.927`$. In Fig. (2), we plot the KT melting temperature as a function of $`Q`$. The closed circles represent the KT phase boundary from our Bethe ansatz calculation, the open squares from Hanna et al. based on the zero temperature value of the elastic moduli, and the solid curve the CI phase boundary. For $`Q_{\mathrm{min}}<Q<Q_c(0)`$, the system initially stays at the C phase. As $`T`$ increases, it makes a transition to the soliton liquid phase. With the further increase of $`T`$, the soliton liquid becomes the SL phase due to the rapid increase of the soliton density upon entering the CI phase boundary. Subsequently it melts exhibiting the reentrant behaviour of the soliton liquid phase. We have confirmed that at large $`Q`$, $`T_{\mathrm{KT}}`$ approaches to $`(\pi /2)\rho _s`$ . Read argues that at $`T=0`$, quantum fluctuations are not important, since the domain ‘sheets’ are marginally rough. However at finite $`T`$, thermal length scale $`\xi _\beta =(\mathrm{}v/k_BT)`$ becomes finite. Hence near the CI phase boundary and/or KT phase boundary, the large distance thermal fluctuations will be still important.

To summarize, we have studied the melting of the SL phase by numerically solving the newly developed set of Bethe ansatz equations, which fully take into account both the thermal fluctuations of the soliton walls and the density variations. Based on the elastic moduli thus obtained, $`T_{\mathrm{KT}}`$ is calculated. We predict that the system will exhibit the reentrant behaviour of the soliton liquid phase within certain ranges of $`Q`$. Our proposed KT phase boundary can be experimentally observed by measuring the longitudinal resistance of the high purity bilayer sample. Towards the completion of the paper, we became aware that the preprint cond-mat/0201343 also studied the bilayer quantum Hall system based on the exact solution of the 1D sine-Gordon model. In our paper, we mainly focus on the role of the compactness of the angle variables and explicitly calculate the KT phase boundary.

K.M. acknowledges C. Hanna for useful discussions and allowing his data for the publication in our paper. This work was supported by grant No. R01-1999-00018 from the interdisciplinary Research program of the KOSEF. K.M. acknowledges the financial support by Yonsei University Research Fund of 2000.


text:
# Cost of QCD simulations with 𝑛_𝑓=2 dynamical Wilson fermions

## 1 INTRODUCTION

Realistic QCD simulations with dynamical fermions require to operate beyond the $`\rho `$ decay threshold and closer to the continuum limit. High energy physics can profit from lattice QCD as soon as systematical and statistical errors are comparable to or smaller as experimental ones .

In order to make best usage of the next generation QCD machines concerning these physical goals, we should employ the most appropriate lattice actions and simulation algorithms. Hence, cost predictions are needed for different lattice discretizations and algorithms beyond $`\frac{m_{\text{PS}}}{m_\text{V}}=0.5`$.

As suggested by R. Kenway (Edinburgh)111The Wuppertal and Edinburgh groups are members of the working group “Algorithms” supported by the EU network HPRN-CT-2000-00145 Hadrons Lattice QCD., chairman of the Lattice 2001 panel discussion, we determine the costs of the SESAM and T$`\chi `$L experiments to estimate the costs of future hybrid Monte Carlo simulations with two degenerate flavors of standard Wilson fermions.

## 2 SESAM/T$`\chi `$L SIMULATIONS

SESAM/T$`\chi `$L has generated 10 ensembles of full QCD vacuum configurations with $`O(5000)`$ HMC trajectories each, at $`\beta =5.6`$ and $`5.5`$ in the region $`0.57<\frac{m_{\text{PS}}}{m_\text{V}}<0.85`$ with two flavors of Wilson fermions. The lattice sizes are $`16^3\times 32`$ (SESAM) and $`24^3\times 40`$ (T$`\chi `$L), corresponding to physical sizes (from $`\rho `$ mass) of $`1.372(36)`$ fm (SESAM) and of $`1.902(34)`$ fm (T$`\chi `$L) after chiral extrapolation. We used the standard hybrid Monte Carlo algorithm , boosted by BiCGStab as solver , ll -SSOR preconditioning  and the educated guess procedure (chronological inversion method) . Running primarily on APE100 systems at DESY/Zeuthen, DFG/Bielefeld and INFN/Rome, the total costs of the simulations amount to about 0.06 Tflops-yrs.

In table 1 we give important quantities from SESAM/T$`\chi `$L. We have exploited both o/e and SSOR preconditioned fermion actions . The latter series are used for the cost analysis, as SSOR shows better scaling behavior.

## 3 SCALING FITS

Figure 1a presents fits to the average number of iterations, $`N_{\text{it}}`$, of the ll -SSOR preconditioned BiCGStab solver, as function of $`1/m_{\text{PS}}a`$. Note that the $`\beta =5.6`$ result scales much better than $`\beta =5.5`$. Presumably, smoother gauge fields allow for better ll-SSOR preconditioning.

Figure 1b shows fits to the integrated autocorrelation times of the series of $`1/N_{\text{it}}`$, the inverse number of iterations. $`1/N_{\text{it}}`$ is related to the minimal eigenvalue of $`M^{}M`$. Its autocorrelation time is comparable to that of the topological charge. Again, we observe a strong $`\beta `$ dependence.

## 4 COST FORMULA

Comparing the sustained CPU time for one HMC trajectory on APE100 with the costs for the iterative solver allows to determine the normalization of the cost function. Multiplication with the autocorrelation time gives the effort to generate one statistically independent configuration at given $`\beta `$ as function of $`1/m_{\text{PS}}a`$. We assume the volume to scale like $`(L/a)^5`$ (unlike $`(L/a)^{4.55}`$ of Ref. ) and the temporal lattice extent to be twice as large as the spatial extent $`L`$.

$`\beta =5.6:N_{\text{flops}}=2.3(7)10^7\left({\displaystyle \frac{L}{a}}\right)^5\left({\displaystyle \frac{1}{am_{\text{PS}}}}\right)^{2.8(2)},`$
$`\beta =5.5:N_{\text{flops}}=1.6(4)10^7\left({\displaystyle \frac{L}{a}}\right)^5\left({\displaystyle \frac{1}{am_{\text{PS}}}}\right)^{4.3(2)}.`$ (1)
## 5 EXTRAPOLATIONS

The CP-PACS quenched simulations  achieved finite $`a`$ results for light hadrons with errors $`<1\%`$ and continuum result with errors between $`1`$ and $`3\%`$. Linearly extrapolating (1) in $`a`$ and $`\frac{m_{\text{PS}}}{m_\text{V}}`$, we find the upper bounds to the CPU time (see table 2) needed to carry out an analogous simulation with $`n_f=2`$ Wilson fermions.