text:
# Giant fluctuations at a granular phase separation threshold

## I Introduction

Dynamics of a system of inelastically colliding hard spheres have attracted a great deal of recent interest Kadanoff ; GG , in particular in the context of validity of kinetic theory and hydrodynamics of rapid granular flow developed in the 80-ies hydro . Hydrodynamics looks ideally suitable for a description of large-scale patterns observed in rapid granular flows: a plethora of clustering phenomena Kudrolli , vortices vort , oscillons sw1 , shocks sw2 , etc., that are difficult to understand in the language of individual particles. However, a first-principle derivation of a universally applicable continuum theory of granular gas is not a simple task, even in the dilute limit. The use of the Enskog equation, the starting point of a systematic derivation of the constitutive relations of granular hydrodynamics, is based on the Molecular Chaos hypothesis. This hypothesis is justified for not too large densities and for an ensemble of elastic hard spheres. Its use for inelastic hard spheres is not obvious, as inelasticity of the particle collisions introduces inter-particle correlations correlations . The correlations become stronger as the inelasticity of the collisions increases. On the contrary, for nearly elastic collisions, $`1r^21`$ (where $`r`$ is the coefficient of normal restitution) the correlations are small, and the Enskog equation can be safely used.

An important additional assumption, made in the process of the derivation of hydrodynamics from the Enskog equation, is scale separation. Hydrodynamics demands that the mean free path of the particles be much less than any characteristic length scale, and the mean time between two consecutive collisions be much less than any characteristic time scale described hydrodynamically. This condition should be verified, in every specific system, after the hydrodynamic problem is solved and the characteristic length and time scales determined. Again, it is safe to say that this condition can be satisfied if the particle collisions are nearly elastic Esipov ; Grossman ; Goldhirsch . Restrictive as it is, the nearly elastic limit is conceptually important just because granular hydrodynamics is expected to work here.

Another potentially important, albeit largely unexplored, limitation of the validity of granular hydrodynamics (or, rather, of any continuum approach to rapid granular flow) is due to the noise caused by the discrete nature of particles. Noise is stronger here than in classical (molecular) fluids simply because the number of particles is much smaller. In addition, noise can be amplified at thresholds of hydrodynamic instabilities as found, for example, in Rayleigh-Bénard convection of classical fluids Hohenberg .

The validity of hydrodynamic description in general, and the accuracy of constitutive relations in particular, can be conveniently checked on symmetry-breaking instabilities that are abundant in rapid granular flows. The example of a symmetry-breaking instability that we consider in this work deals with a very simple setting: a two-dimensional (2D) system of nearly elastically colliding hard spheres, confined by a rectangular box and driven by a thermal sidewall at zero gravity. The setting is described in detail in Sec. II. The basic steady state here is the “stripe state”: a stripe of enhanced density at the wall opposite to the driving wall Grossman . In the continuum language, the stripe state is uniform in the lateral direction, by which we mean the direction parallel to the driving wall. Within a certain range of parameters (delineated below), steady-state equations of granular hydrodynamics predict spontaneous symmetry breaking instability of the stripe state, when the aspect ratio of the confining box exceeds a certain threshold LMS ; Brey2 ; KM ; LMS2 . The instability leads to phase separation: the development of “droplets” (high-density domains) coexisting with “bubbles” (low-density domains). For very large aspect ratios of the box, this symmetry-breaking instability has been recently observed in event-driven molecular dynamic (EMD) simulations, and described by a phenomenological continuum model Argentina . The present work is devoted to a more detailed investigation of the phase separation instability in the range of aspect ratios comparable to the threshold value. We employ, in Sec. III, the equations of granular hydrodynamics (or rather hydrostatics) to compute the supercritical bifurcation curve for the phase separation instability. Then we report, in Section IV, on extensive EMD simulations that show that this bifurcation curve is quantitatively accurate well below and well above the threshold value of the aspect ratio. Unexpectedly, the hydrostatic theory fails in a relatively wide region of aspect ratios around the threshold value, where the system is found to exhibit giant fluctuations. In an attempt to get insight into the mechanism of this anomaly, we investigate, also in Section IV, the dependence of the magnitude of fluctuations on the total number of particles in the system. A summary and discussion of our results is presented in Section V.

## II Model system and hydrostatic equations

Let $`N`$ hard spheres of diameter $`d`$ and mass $`m=1`$ move in a 2D rectangular box $`L_x\times L_y`$. The inelasticity of particle collisions is parameterized by a constant coefficient of normal restitution $`r`$. Particle collisions with three of the walls are elastic. The fourth, thermal wall is located at $`x=L_x`$. Upon collision with it, the normal component of the particle velocity is drawn from a Maxwell distribution with temperature $`T_0`$ Grossman , while the tangential component of the particle velocity is preserved.

Working in the nearly elastic limit $`1r^21`$ and employing the Navier-Stokes hydrodynamics hydro , we introduce the number density $`n(𝐫,t)`$, granular temperature $`T(𝐫,t)`$ and mean-flow velocity $`𝐯(𝐫,t)`$. Energy input at the thermal wall can be balanced by the dissipation due to inter-particle collisions. Therefore, we assume that the system reaches a zero-mean-flow steady state $`𝐯=\mathrm{𝟎}`$, and is therefore describable by the simple momentum and energy balance equations:

$$p=\mathrm{const},(\kappa T)=I.$$
(1)
Here $`p`$ is the pressure, $`\kappa `$ is the thermal conductivity and $`I`$ is the rate of energy loss by collisions. The hydrostatic Eqs. (1) should be supplemented by constitutive relations: $`p,\kappa `$ and $`I`$ in terms of $`n`$ and $`T`$. These relations are derivable systematically only in the dilute limit hydro ; dilute . Being interested in moderate densities, we shall employ the well-known constitutive relations by Jenkins and Richman JR , that account for excluded particle volume. In the nearly-elastic limit one can neglect the inelasticity correction terms in $`p`$ and $`\kappa `$, as well as the small density gradient term, proportional to $`1r`$, in the heat flux Brey1 .

Equations (1) can be rewritten in terms of a single variable: the scaled inverse density $`z(x,y)=n_c/n(x,y)`$, where $`n_c=2/(\sqrt{3}d^2)`$ is the hexagonal close-packing density. In scaled coordinates, $`𝐫/L_x𝐫`$, the box dimensions become $`1\times \mathrm{\Delta }`$, where $`\mathrm{\Delta }=L_y/L_x`$ is the box aspect ratio. We obtain KM

$$(F(z)z)=\eta Q(z),$$
(2)
where $`F(z)=A(z)B(z)`$,

$`A(z)`$ $`=`$ $`{\displaystyle \frac{G\left[1+\frac{9\pi }{16}\left(1+\frac{2}{3G}\right)^2\right]}{z^{1/2}(1+2G)^{5/2}}},`$
$`B(z)`$ $`=`$ $`1+2G+{\displaystyle \frac{\pi }{\sqrt{3}}}{\displaystyle \frac{z(z+\frac{\pi }{16\sqrt{3}})}{(z\frac{\pi }{2\sqrt{3}})^3}},`$
$`Q(z)`$ $`=`$ $`{\displaystyle \frac{6}{\pi }}{\displaystyle \frac{z^{1/2}G}{(1+2G)^{3/2}}},`$
$`G`$ $`=`$ $`G(z)={\displaystyle \frac{\pi }{2\sqrt{3}}}{\displaystyle \frac{z\frac{7\pi }{32\sqrt{3}}}{(z\frac{\pi }{2\sqrt{3}})^2}},`$ (3)
and $`\eta =(2\pi /3)(1r)(L_x/d)^2`$ is the hydrodynamic inelasticity parameter. Introducing $`\psi (x,y)=_0^zF(z^{})𝑑z^{}`$, we arrive at the following equation:

$$^2\psi =\eta \stackrel{~}{Q}(\psi ),$$
(4)
where $`\stackrel{~}{Q}(\psi )=Q\left[z(\psi )\right]`$ (in the following the symbol $`\stackrel{~}{}`$ is omitted). The boundary conditions are

$$\frac{\psi }{x}|_{x=0}=\frac{\psi }{y}|_{x=1}=\frac{\psi }{y}|_{y=\mathrm{\Delta }/2}=\frac{\psi }{y}|_{y=\mathrm{\Delta }/2}=0.$$
(5)
Finally, the number of particles is conserved:

$$\frac{1}{\mathrm{\Delta }}\underset{\mathrm{\Delta }/2}{\overset{\mathrm{\Delta }/2}{}}\underset{0}{\overset{1}{}}\frac{dxdy}{z(\psi )}=\frac{N}{L_xL_yn_c}f.$$
(6)
The hydrostatic problem (4)-(6) is fully determined by three scaled parameters: the area fraction $`f`$, $`\eta `$, and $`\mathrm{\Delta }`$. Notice that the steady-state density distributions are independent of $`T_0`$, as the hard sphere model does not introduce any intrinsic energy scale.

## III Stripe state, symmetry-breaking instability and bifurcation curve

The trivial steady state of the system is a laterally uniform cluster of particles located at the wall $`x=0`$, opposite to the thermal wall Grossman , see Fig. 1. This state will be called the stripe state. In the language of hydrodynamics, it is described by the $`y`$-independent solution of Eqs. (4)-(6); we shall denote it by $`z=Z(x)`$, correspondingly $`\psi =\mathrm{\Psi }(x)`$.

It was predicted that, in a wide region of the parameter space $`(f,\eta ,\mathrm{\Delta })`$, the stripe state should give way, by a symmetry-breaking bifurcation (either supercritical or subcritical), to a laterally asymmetric state LMS ; Brey2 ; KM ; LMS2 . For very large aspect ratios $`\mathrm{\Delta }`$, this phase-separation instability has been observed in EMD simulations Argentina . For a laterally asymmetric steady state one can write

$$\psi (x,y)=\mathrm{\Psi }(x)+\underset{n}{}\phi _n(x)\mathrm{exp}(inky),$$
(7)
where $`\phi _n(x)=\phi _n^{}(x)`$. What happens close to the supercritical bifurcation point? Here the leading terms are those with $`n=\pm 1`$, while $`\phi _0\phi _1^2`$, $`\phi _2\phi _1^2`$, $`\phi _3\phi _1^3`$, etc. The bifurcation point itself can be found from the linear eigenvalue problem

$$\phi _{1k}^{\prime \prime }\eta Q_\mathrm{\Psi }\phi _{1k}k_c^2\phi _{1k}=0,$$
(8)
$$\phi _{1k}^{}(0)=0\text{and}\phi _{1k}(1)=0$$
(9)
that was analyzed in Refs. LMS ; Brey2 ; KM . Here

$$Q_\mathrm{\Psi }(x)=F^1dQ/dz|_{z=Z(x)}.$$

For given $`\eta `$ and $`f`$, one obtains the eigenvalue $`k=k_c(\eta ,f)`$ and corresponding eigenfunction $`\phi _{1k}(x)`$. The modes with $`k<k_c(\eta ,f)`$ are unstable. Within a spinodal interval $`f_1(\eta )<f<f_2(\eta )`$, the effective lateral compressibility of the gas is negative, and this is the mechanism of the instability KM ; Argentina . At $`\eta 1`$, there is a range of $`f`$ such that $`k_c`$ and $`\phi _{1k}(x)`$ become insensitive to the precise form of the boundary conditions at the driving wall. This is the universal ”localization regime”, when the eigenfunction $`\phi _{1k}(x)`$ is exponentially localized at the wall opposite to the driving wall LMS ; KM . The spinodal interval exists for $`\eta _c<\eta <\mathrm{}`$; it shrinks to zero at $`\eta =\eta _c344.3`$ Argentina ; KMS . It has been recently shown, for a different boundary condition at the driving wall, that the bifurcation from the stripe state to a phase-separated state is supercritical within some density interval $`f_{}(\eta )<f<f_+(\eta )`$, which is located within the spinodal interval. On each of the intervals $`f_1<f<f_{}`$ and $`f_+<f<f_2`$, the bifurcation is subcritical LMS2 .

As we have already noted, the present work focuses on the phase separation via a supercritical bifurcation. To obtain the asymptotics of the supercritical bifurcation curve close to onset, one should go to the second order of the perturbation theory and take into account, in Eq. (7), the terms $`n=0,\pm 1`$ and $`\pm 2`$. In this way one obtains three linear ordinary differential equations, presented in Ref. LMS2 , where the same problem was solved for a different boundary condition at the driving wall. The solvability condition for these equations Iooss  yields the bifurcation curve: $`A`$ versus $`k_c^2k^2`$. The amplitude $`A`$ can be uniquely defined by the relation

$$\phi (x)=A\mathrm{\Phi }_0(x)+A|A|^2\delta \phi (x),$$

where $`\mathrm{\Phi }_0(x)`$ is the solution of Eqs. (8) and (9) such that $`\mathrm{\Phi }_0(0)=1`$, while $`\delta \phi (x)=𝒪(1)`$. This yields

$$A\left(k_c^2k^2\right)=CA|A|^2,$$

where $`C=`$const. The trivial solution $`A=0`$ describes the stripe state, while the nontrivial one, $`k_c^2k^2=C|A|^2`$ describes the bifurcated state. The constant $`C`$ can be computed numerically. $`C>0(<0)`$ corresponds to supercritical (subcritical) bifurcation. We present here the resulting bifurcation curve for $`Y_c`$, the (normalized) $`y`$-coordinate of the center of mass of the granulate

$$Y_c=\frac{_0^1𝑑x_{\mathrm{\Delta }/2}^{\mathrm{\Delta }/2}yn(x,y)𝑑y}{\mathrm{\Delta }_0^1𝑑x_{\mathrm{\Delta }/2}^{\mathrm{\Delta }/2}n(x,y)𝑑y},$$
(10)
Let us fix $`\eta `$ and $`f`$ and treat $`\mathrm{\Delta }`$ as the control parameter. When $`\mathrm{\Delta }`$ is slightly larger than $`\mathrm{\Delta }_c=\pi /k_c(f)`$, only the fundamental mode $`k=\pi /\mathrm{\Delta }`$ is unstable, and the bifurcation curve has the form

$$\left|Y_c\right|=\mathrm{{\rm Y}}\left(\mathrm{\Delta }\mathrm{\Delta }_c\right)^{1/2}.$$
(11)
Here

$$\mathrm{{\rm Y}}=\frac{2^{3/2}f_0}{C^{1/2}\mathrm{\Delta }_cf},f_0=2_0^1𝑑x\frac{\mathrm{\Phi }_{01}}{Z^2F},$$

and $`\mathrm{\Phi }_{01}(x)`$ is the solution of initial-value problem for Eq. (8) with the initial conditions $`Y(0)=1`$ and $`Y^{}(0)=0`$. Equation (11) assumes $`C>0`$: a supercritical bifurcation. We have found that, at fixed $`\eta `$, $`C>0`$ on an interval $`f_{}(\eta )<f<f_+(\eta )`$ that lies within the spinodal interval $`(f_1,f_2)`$. On the intervals $`f_1<f<f_{}`$ and $`f_+<f<f_2`$ the coefficient $`C`$ becomes negative which indicates a subcritical bifurcation. The solid line in Fig. 6 shows the supercritical bifurcation curve (11) for $`\eta =11,050`$ and $`f=0.025`$. Here $`\mathrm{\Delta }_c0.514`$ and $`\mathrm{{\rm Y}}0.142`$.

When $`\mathrm{\Delta }`$ is well above $`\mathrm{\Delta }_c`$, the weakly nonlinear theory is invalid, and a numerical solution of the fully nonlinear hydrostatic problem (4)-(6) is needed for the determination of $`\left|Y_c\right|`$. An alternative approach is a hydrodynamic simulation, that is a numerical solution of the hydrodynamic equations. Numerical simulations of this type were done in Ref. LMS2  for a different version of constitutive relations Grossman  and a different boundary condition at the driving wall. It was observed that the phase-separation instability produces multiple clusters whose further dynamics proceed as gas-mediated competition and coarsening. Direct merging of clusters can also occur. The final symmetry-broken state, as observed in the hydrodynamic simulations, is always a single, almost densely packed stationary 2D cluster coexisting with gas (or dilute bubble coexisting with denser fluid). The cluster is located in one of the system’s corners (unless periodic boundary conditions are used). This scenario was confirmed in a hydrodynamic simulation of the present system (for $`\eta =11,050`$, $`f=0.025`$ and $`\mathrm{\Delta }=3`$) done by E. Livne Livne . A density map of the hydrodynamic final state in this case is shown in Fig. 2D. The steady-state value $`\left|Y_c\right|0.265`$, obtained in this simulation, is shown by the circle in Fig. 6.

## IV EMD simulations

### IV.1 Simulation method, parameters and diagnostics

We put the predictions of the granular hydrostatics into test by doing extensive EMD simulations of this system. Most of the simulations were done with $`N=210^4`$ particles: hard disks of diameter $`d=1`$ and mass $`m=1`$. The thermal wall temperature is $`T_0=1`$, so the scaled time unit is $`d(m/T_0)^{1/2}=1`$. A standard event-driven algorithm Alder  was used. Two of the hydrodynamic parameters, $`\eta =11,050`$ and $`f=0.025`$, were fixed in all simulations, while $`\mathrm{\Delta }`$ was varied in the range of $`0.1<\mathrm{\Delta }<3`$. This was achieved by varying $`L_x`$, $`L_y`$ and $`r`$. Indeed, for a fixed $`\eta `$, $`f`$, $`\mathrm{\Delta }`$ and $`N`$ the coefficient of normal restitution,

$$r=1\frac{\sqrt{3}\eta f\mathrm{\Delta }}{\pi N},$$
(12)
and the system’s dimensions,

$$L_x=\left(\frac{\sqrt{3}N}{2f\mathrm{\Delta }}\right)^{1/2}\text{and}L_y=\left(\frac{\sqrt{3}N\mathrm{\Delta }}{2f}\right)^{1/2},$$
(13)
are uniquely determined. For the values of the parameters that we used, $`r`$ was always in the range of nearly elastic collisions: $`r0.977`$. The initial spatial distribution of the particles was (statistically) uniform, while the initial velocity distribution was Maxwell’s with the wall temperature $`T_0=1`$. The center-of-mass coordinate $`Y_c(t)`$ was used as a quantitative probe of the lateral asymmetry of the system. Before taking the steady-state measurements we waited until transients died out. This was monitored by the time-dependence of the average kinetic energy of the particles (that first decayed and then approached an almost constant value) and by the time-dependence of the center of mass itself, see below. Selected movies of these simulations can be downloaded from http://bioinf.charite.de/kies/giantfluctuations/.

### IV.2 Final states at different $`\mathrm{\Delta }`$

The EMD simulations showed that, at aspect ratios well below the threshold value of $`\mathrm{\Delta }=\mathrm{\Delta }_c0.512`$, the final state is a (weakly fluctuating) stripe state. The number density profile versus $`x`$, found in the simulations, compares very well with the hydrostatic solution (see Fig. 1), while $`Y_c(t)`$ stays close to zero. Notice that the Jenkins-Richman constitutive relations JR , that we used in this comparison, do not include any fitting parameters. Therefore, well below the instability threshold in $`\mathrm{\Delta }`$, the hydrostatic solution yields a quantitatively accurate leading-order description of the system.

At aspect ratios well above the instability threshold we always observed several clusters nucleating at the wall opposite to the driving wall. The cluster dynamics (Fig. 2 A to C) proceeds as gas-mediated competition and coarsening (sometimes as direct mergers) of clusters, in accord with hydrodynamic simulations LMS2 . As time increases, the number of clusters goes down, and only one dense cluster, fluctuating around its average position in one of the two corners, opposite to the thermal wall, finally survives. Fig. 2C shows a snapshot of the final state for $`\mathrm{\Delta }=3`$. For comparison, Fig. 2D shows a density map of the final steady state obtained by E. Livne in a hydrodynamic simulation for the same hydrodynamic parameters. The center-of-mass position $`Y_c`$ of the steady state agrees well with the average-in-time center-of-mass position, measured in the EMD simulations, as shown by the circle in Fig. 6. This indicates that, well above the instability threshold, the hydrostatic theory describes the steady states of the system well. We can also refer the reader to the recent EMD simulation results for very large aspect ratios Argentina . As no appreciable fluctuations around a broken-symmetry steady state were reported, one can safely assume that the broken-symmetry steady states observed in Ref. Argentina  should be also describable by the hydrostatic theory.

The system behavior changes dramatically, however, as the aspect ratio $`\mathrm{\Delta }`$ approaches $`\mathrm{\Delta }_c`$. We found that, in a wide region of $`\mathrm{\Delta }`$ around $`\mathrm{\Delta }_c`$, the final state of the system exhibits large-amplitude irregular oscillations, as dense clusters at the wall opposite to the driving wall nucleate, move in the lateral direction, dissolve and reappear. Figure 3 shows a typical sequence of snapshots from an EMD simulation for $`\mathrm{\Delta }=1`$.

Figure 4 shows the time history of the center-of-mass coordinate $`Y_c`$ for six different values of $`\mathrm{\Delta }`$. One can see that, in a wide region of intermediate $`\mathrm{\Delta }`$, the center-of-mass coordinate $`Y_c(t)`$ shows large-amplitude irregular oscillations. Noticeable are multiple zero crossings of $`Y_c(t)`$ at aspect ratios above the hydrodynamic bifurcation point $`\mathrm{\Delta }_c0.512`$ (Fig. 4 c-e). Smaller but still significant irregular oscillations are also observed below $`\mathrm{\Delta }_c`$, as if the system persistently tends to break the lateral symmetry there. The hydrostatic picture is recovered when one moves farther away, in any direction, from the region of $`\mathrm{\Delta }\mathrm{\Delta }_c`$. Indeed, Fig. 4e shows that zero crossings of $`Y_c(t)`$ occur less often for $`\mathrm{\Delta }=1.3`$, than for $`\mathrm{\Delta }=0.7`$ or $`1`$. At still larger $`\mathrm{\Delta }`$ (Fig. 4f) no zero crossings are observed for any reasonable simulation time, and $`Y_c`$ fluctuates around a constant value that is very close to that predicted by the hydrostatic theory (and shown by the circle in Fig. 6).

To better characterize the fluctuation-dominated region, we computed the probability distribution function $`P(|Y_c|)`$ of different values of $`|Y_c|`$ in a statistical steady state, that is, after transients die out. The stationarity of the remaining data was tested by dividing the respective time interval into three sub-intervals and checking that the differences in $`P(|Y_c|)`$ for the sub-intervals are small and not systematic. The probability distribution $`P(|Y_c|)`$ is shown, at different $`\mathrm{\Delta }`$, in Fig. 5. At $`\mathrm{\Delta }\mathrm{\Delta }_c`$ the maximum of $`P(|Y_c|)`$ is at $`|Y_c|=0`$, and it is relatively narrow. Correspondingly, there is no symmetry breaking there, the fluctuations are relatively small, and the hydrostatic theory yields an accurate leading-order description. At $`\mathrm{\Delta }\mathrm{\Delta }_c`$, the maximum of $`P(|Y_c|)`$ is at a non-zero $`|Y_c|`$. This is a clear manifestation of symmetry-breaking: a dense cluster develops in one of the corners away from the driving wall. The probability distribution $`P(|Y_c|)`$ is also quite narrow here, the fluctuations are relatively small, and there is a good agreement between the hydrostatic theory and EMD-simulations. On the contrary, in a wide region of $`\mathrm{\Delta }`$ around $`\mathrm{\Delta }_c`$, the probability distribution $`P(|Y_c|)`$ is very broad, and the hydrostatic theory breaks down. By following the position of the maximum of $`P(|Y_c|)`$ at different $`\mathrm{\Delta }`$ (see Fig. 6), one can see that the symmetry-breaking transition occurs somewhere in the region of $`0.3<\mathrm{\Delta }<1.0`$. Because of the extreme flatness and broadness of the probability distribution $`P(|Y_c|)`$ in this region, a more accurate estimate of the position of the maximum of $`P(|Y_c|)`$ requires a much better statistics (that is, a much longer simulation time) than we could afford in this series of simulations simulationtime .

Noticeable in Fig. 6 is a systematic discrepancy, within the wide fluctuation-dominated region, between the positions of the maxima of $`P(|Y_c|)`$ and the hydrostatic bifurcation curve computed in Sec. III. We even cannot exclude a change in the character of bifurcation caused by the fluctuations (apparently without shifting the bifurcation point). Indeed, the maxima of $`P(|Y_c|)`$ at $`\mathrm{\Delta }=1.0`$, $`1.3`$ and $`2.0`$ appear to lie on a straight line passing through the theoretical transition point $`\mathrm{\Delta }_c0.5`$. As $`\mathrm{\Delta }`$ increases further, the discrepancy between the positions of the maxima of $`P(|Y_c|)`$ and the theoretical bifurcation curve goes down accuracy . Importantly, the fluctuation-dominated region $`0.3<\mathrm{\Delta }<1.0`$ does include the hydrostatic transition point $`\mathrm{\Delta }_c0.5`$.

We should stress that the failure of hydrostatics is observed at intermediate values of the aspect ratio $`\mathrm{\Delta }`$, when the hydrodynamic parameters $`\eta `$ and $`f`$, and the number of particles $`N`$, are fixed. In view of Eq. (12), while increasing $`\mathrm{\Delta }`$, one increases the inelasticity of particle collisions $`1r`$. That the hydrostatic theory fails at intermediate values of the inelasticity, and improves at small enough, or large enough inelasticities, excludes the inelasticity itself as the reason for the failure.

### IV.3 Simulations with different $`N`$

We did a series of simulations with different number of particles $`N`$ in order to verify the hydrodynamic scaling and investigate the $`N`$-dependence of the (relatively weak) fluctuations well below and well above $`\mathrm{\Delta }_c`$. These additional simulations were done for $`\mathrm{\Delta }=0.1`$ and three values of $`N`$: $`510^3`$, $`10^4`$ and $`1.510^4`$, and for $`\mathrm{\Delta }=3.0`$ and $`N=410^4`$.

When varying $`N`$ at fixed $`\mathrm{\Delta }`$, we kept the hydrodynamic parameters $`\eta =11,050`$ and $`f=0.025`$ constant. Therefore, if the hydrostatic equations provide a correct leading-order theory of the steady states far below and far above $`\mathrm{\Delta }_c`$, the time-averaged steady state values of $`Y_c`$ should become $`N`$-independent for large enough $`N`$. Figure 7 shows $`Y_c`$ versus time for $`\mathrm{\Delta }=0.1`$ at the four different values of $`N`$. One can see that, in all these cases, the average value of $`Y_c`$ is close to zero as expected, while fluctuations are relatively small. Figure 8 shows the dynamics of $`Y_c(t)`$ for $`\mathrm{\Delta }=3`$ and two different values of $`N`$: $`210^4`$ and $`410^4`$. Here the symmetry-breaking is evident, as a dense cluster develops in a corner. With a moderate accuracy determined by the relatively high level of fluctuations of $`Y_c`$, the average values of $`Y_c`$ at late times are close to each other. Therefore, well below and well above $`\mathrm{\Delta }_c`$ the hydrodynamic scaling is obeyed.

Simulations with fixed scaled parameters $`\eta `$, $`f`$ and $`\mathrm{\Delta }`$, but different $`N`$ can also help in identifying the mechanism of breakdown of the hydrostatic theory at aspect ratios around $`\mathrm{\Delta }_c`$. Indeed, it is natural to interpret the giant oscillations, shown in Fig. 4c-e, in terms of a strong coupling between the two bifurcated states predicted by the hydrostatic theory. One possible scenario of this coupling (which we call Scenario I) relies on the discrete-particle noise, unaccounted for by granular hydrodynamics. Below $`\mathrm{\Delta }_c`$, the discrete-particle noise is expected to cause fluctuations, that is to broaden the distribution of $`Y_c`$ as indeed observed in Fig. 5. If Scenario I is correct, the standard deviation $`\sigma `$ of $`Y_c(t)`$ from its average value should vanish as $`N`$ goes to infinity, at fixed hydrodynamic parameters $`\eta `$, $`f`$ and $`\mathrm{\Delta }`$.

Another possibility (Scenario II) is that the fluctuations persist in the limit of $`N\mathrm{}`$. If this is the case, the dominating mechanism of fluctuations has a purely hydrodynamic nature and should be explainable by a full hydrodynamic analysis (as opposed to our hydrostatic analysis, and to the simplified hydrodynamic simulations that used a model Stokes friction instead of the full viscosity). Here the coupling between the two symmetry-broken states may be due to either an unstable hydrodynamic mode (Scenario IIa), or a weakly damped mode (Scenario IIb). In Scenario IIb, $`\sigma `$ should vanish, as $`N\mathrm{}`$, if one waits for a sufficiently long time. Therefore, to distinguish between the two sub-scenarios, one should, in addition to the limit of $`N\mathrm{}`$, take the limit of $`t\mathrm{}`$.

Obviously, one is unable to take any of these two limits in actual EMD simulations, where the maximum achievable values of $`N`$ and $`t`$ are limited by the available computer resources. So what was observed in our EMD simulations with different $`N`$? Figures 7 and 9 show what happens well below $`\mathrm{\Delta }_c`$, when $`N`$ increases from $`\mathrm{5\hspace{0.17em}000}`$ to $`\mathrm{20\hspace{0.17em}000}`$. One can see from Fig. 7 that, as $`N`$ grows, the high-frequency components of the fluctuations do decrease, but the low frequency component does not show any pronounced decrease. Overall, the fluctuation spectrum moves towards the lower frequencies. As the result, a good resolution of the low-frequency part of the power spectrum requires longer and longer simulations (which rapidly become prohibitively long). This introduces an additional, non-trivial constraint on simulations with a large number of particles. A similar situation occurs well above $`\mathrm{\Delta }_c`$. Figure 8 does indicate that $`\sigma `$ goes down as $`N`$ goes up from $`\mathrm{20\hspace{0.17em}000}`$ to $`\mathrm{40\hspace{0.17em}000}`$. However, one also observes that, as $`N`$ grows, the role of the low-frequency components of the fluctuations increases.

Hydrodynamics provides a hint for the mechanism of the “red shift” of the power spectrum with an increase of $`N`$. There are four hydrodynamic modes in the system: two acoustic modes, the entropy mode and the shear mode. The frequencies of the acoustic modes are the highest, as they are determined by the “ideal” (non-dissipative) terms in the hydrodynamic equations, and they scale like the inverse system size. The frequencies of the entropy and shear modes are much lower, as they are determined by the transport coefficients: the heat conduction, viscosity and inelastic loss rate, and they scale like the inverse square of the system size. In the units of $`d=m=T_0=1`$, and at fixed hydrodynamic parameters $`\eta `$, $`f`$ and $`\mathrm{\Delta }`$, a larger $`N`$ implies a larger system, see Eqs. (13). Correspondingly, as $`N`$ increases, the characteristic frequencies of the entropy/shear modes go down much faster than those of the acoustic modes. Therefore, it seems likely that one of these modes is responsible for the low-frequency components of the fluctuations. A related issue is that, in contrast to the hydrostatic problem (1), the full time-dependent hydrodynamic problem has an additional scaled parameter: $`d/L_x`$. This parameter describes the role of the dissipative terms compared to the “ideal” terms in the hydrodynamic equations. As it is clear from Eq. (13), when increasing $`N`$ at constant $`\eta `$ and $`f`$, one reduces this additional parameter. Therefore, as $`N`$ increases, the low-frequency shear/entropy modes should become more and more persistent. As these modes are not necessarily broad-band, $`\sigma `$ might cease to provide a good characterization of the system at large $`N`$.

Still, if one continues following $`\sigma `$ as $`N`$ increases, one observes (see Fig. 9) that $`\sigma `$ decreases much slower than the classic dependence $`N^{1/2}`$ characteristic of equilibrium systems. If one attempts to interpret the decrease of $`\sigma `$ with an increase of $`N`$ in terms of an empiric power law, one obtains an exponent $`0.23`$, instead of the classical value of $`1/2`$ for equilibrium systems. Importantly, we did reproduce the classical $`N^{1/2}`$ scaling of $`\sigma `$ in a control series of simulations with the same $`f`$ and $`\mathrm{\Delta }`$, but with $`\eta =0`$ (elastic collisions). Moreover, a good quantitative agreement was obtained with a theoretical result for $`\sigma `$ that directly follows from the classic expression for the density correlation function in equilibrium LL . We also found that, for the same total number of particles $`N`$, the fluctuation levels in the elastic case are significantly lower than in the inelastic case. That is, well below $`\mathrm{\Delta }_c`$, the fluctuations, though much smaller than those observed for $`\mathrm{\Delta }\mathrm{\Delta }_c`$, are still large compared to the elastic case.

To summarize this subsection, our simulations with different $`N`$ strongly indicate that the hydrostatic equations provide a correct leading order theory of this system well below and well above $`\mathrm{\Delta }_c`$. On the other hand, the simulations proved to be insufficient for determining the mechanism of giant fluctuations that we observed in this system at $`\mathrm{\Delta }\mathrm{\Delta }_c`$. We cannot even be sure at this point whether the fluctuations (or, more precisely, their low-frequency components) persist or not as $`N\mathrm{}`$.

## V Summary and Discussion

The main results of this work can be summarized in the following way. Granular hydrostatics, in combination with simplified hydrodynamic simulations, correctly predict the phase separation instability in this prototypical driven granular system. Well above and well below the critical value of the aspect ratio $`\mathrm{\Delta }_c`$, the hydrostatic theory describes the steady state of the system well. However, in a wide region of aspect ratios around $`\mathrm{\Delta }_c`$ the system is dominated by fluctuations, and the hydrostatic theory fails. The fluctuation levels are anomalously high even relatively far from the hydrostatic bifurcation point, and they certainly do not exhibit the classic $`N^{1/2}`$ scaling with the number of particles $`N`$.

Though we are unable to pinpoint the mechanism of excitation of the giant fluctuations, we can suggest two different scenarios for their origin. In Scenario I the fluctuations are driven by discrete particle noise. Indeed, it is well known that discrete particle noise can drive relatively large fluctuations in the vicinity of thresholds of hydrodynamic instabilities Hohenberg  and non-equilibrium phase transitions Haken . Fluctuations of this type should vanish as one increases indefinitely the number of particles in the system, keeping the hydrodynamic parameters constant. Unfortunately, our simulations with different $`N`$, but fixed $`\eta `$, $`f`$ and $`\mathrm{\Delta }`$, have been insufficient to prove or disprove this scenario.

A difficulty with Scenario I is that the fluctuations are so big in so wide a region of aspect ratios. No anomaly of this type has been observed in any other symmetry-breaking instability of granular flow, even with much smaller numbers of particles. As an example, let us consider for a moment the same system, but introduce gravity in the $`x`$ direction. Now the granular gas is heated from below, and the system exhibits another symmetry-breaking instability: thermal convection, similar to the Rayleigh-Bénard convection of classical fluids. The transition to convection occurs via a supercritical bifurcation Ramirez ; HMD ; KMconvection . Though EMD simulations of thermal granular convection Ramirez  involved only $`N=2,300`$ particles (which is much less than $`N=210^4`$ used in the present work), a sharp supercritical bifurcation was observed, in agreement with a hydrodynamic analysis HMD ; KMconvection . By comparison, the giant fluctuations, observed in a wide region of $`\mathrm{\Delta }`$ in the present work, are an anomaly, as one needs some (hydrodynamic?) mechanism of strong amplification of the discrete-particle noise.

If Scenario I proves to be correct, the corresponding theory can be developed in the framework of Fluctuating Hydrodynamics LL , generalized to granular gases in the limit of nearly elastic collsions. Fluctuating Hydrodynamics is a Langevin-type theory that takes into account the discrete character of particles by adding delta-correlated noise terms in the momentum and energy equations LL . Fluctuating Hydrodynamics is by now well established for classical fluids in 3D, including non-equilibrium states Hohenberg ; Mansour . We should mention here that the 2D case has an additional difficulty. The coupling of fluctuations here is anomalously strong, even in the elastic case: the transport coefficients diverge in the thermodynamic limit, except for a sufficiently dilute gas 2D . Therefore, one can hope to generalize the Fluctuating Hydrodynamics to the 2D gas of inelastic hard spheres in the dilute limit unfortunately . Close to the phase separation threshold, the dilute limit holds with a reasonable accuracy. It would be interesting to investigate the phase separation problem in 3D, where important differences in the fluctuation behavior may occur.

Alternatively, in Scenario II the low-frequency component of the giant fluctuations has a purely hydrodynamic origin and is driven either by a presently unknown hydrodynamic instability (Scenario IIa), or by a long-lived transient mode (Scenario IIb). Effects of these type are obviously missed by a hydrostatic analysis. They may have also been missed by the time-dependent hydrodynamic simulation Livne  that employed a model Stokes friction, rather than the hard-sphere viscosity, to accelerate the convergence to a steady state. If Scenario II is correct, the low-frequency component of the fluctuations should be observable in hydrodynamic simulations with the true hard-sphere viscosity. These simulations, therefore, should be an important next step in the analysis of this fascinating problem.

We are very grateful to E. Livne for doing the hydrodynamic simulation. We thank A. Barrat, I. Goldhirsch and E. Khain for useful discussions. This research was supported by Deutsche Forschungsgemeinschaft (Grant PO 472/5-2), by the Israel Science Foundation (grant No. 180/02), by the Russian Foundation for Basic Research (grant No. 02-01-00734) and by Deutscher Akademischer Austauschdienst.


format: text

text:
# Escaping free-energy minima

## I Methodology

We shall assume here that there exist a finite number of relevant collective coordinates $`s_i,i=1,n`$ where $`n`$ is a small number, and we consider the dependence of the free energy $`(s)`$ on these parameters. Practical examples of appropriate choices of these variables will be given below. The exploration of the FES is guided by the forces $`F_i^t=/s_i^t`$. In order to estimate these forces efficiently, we introduce an ensemble of $`P`$ replicas of the system, each obeying the constraint that the collective coordinates have a pre-assigned value $`s_i=s_i^t`$, and each evolved independently at the same temperature $`T`$. Since the $`P`$ replicas are statistically independent, the estimate of thermodynamic observables (e.g. the forces on the constraints) is improved with respect to an evaluation on a single replica, and it can be parallelized in a straightforward manner. The constraints are imposed on each replica via the standard methods of constrained molecular dynamicsallen  by adding to the Lagrangean a term $`_{i=1,n}\lambda _i(s_is_i^t)`$, where $`\lambda _i`$ are Lagrange multipliers. Averaging over the time and over the replicas we can evaluate the derivative of the free energy relative to the $`s_i^t`$’s as $`F_i^t=\lambda _i_t`$Constraints1 ; Constraints2 . For the sake of simplicity we neglect here the small kinematic correction term discussed in Ref. Constraints2 , although this can easily be added. Taking a cue from the work of Gear and Kevrikidiskevrikidis , we use these forces, determined by a microscopic dynamics, to define a coarse-grained dynamics in the space of the $`s_i`$’s. In our case the definition of this dynamics is rather arbitrary and designed only to explore the FES efficiently. The dynamics is defined from the discretized evolution equation:

$$\sigma _i^{t+1}=\sigma _i^t+\delta \sigma \frac{\varphi _i^t}{\left|\varphi ^t\right|}$$
(1)
In eqn. (1) we have introduced the scaled variables $`\sigma _i^t=s_i^t/\mathrm{\Delta }s_i`$ and the scaled forces $`\varphi _i^t=F_i^t\mathrm{\Delta }s_i`$, where $`\mathrm{\Delta }s_i`$ is the estimated size of the FES well in the direction $`s_i`$, $`\left|\varphi ^t\right|`$ is the modulus of the n-th dimensional vector $`\varphi _i^t`$ and $`\delta \sigma `$ is a dimensionless stepping parameter. It should be stressed that eqn. (1) has the form of a steepest descent equation in the direction given by the forces $`\varphi _i`$, and does not imply any real dynamical evolution. The index $`t`$ is only used to label the states. After the collective coordinates are updated using eqn. (1), a new ensemble of replicas of the system with values $`\sigma _i^{t+1}`$ is prepared, and new forces $`F_i^{t+1}`$ are calculated for the next iteration. At the same time the driving forces are evaluated from the microscopic Hamiltonian in short standard microscopic molecular dynamics runs. Since there is no dynamical continuity between the $`P`$ replicas at different iterations, one can use large values of $`\delta \sigma `$ and move very efficiently in the space of the collective coordinates.

Clearly eqn. (1) alone cannot guarantee an efficient exploration of the FES, nor is it useful to determine the FES. This task can be achieved if we replace the forces in eqn. (1) with a history-dependent termHuber ; wang :

$$\varphi _i\varphi _i\frac{}{\sigma _i}W\underset{t^{}t}{}\underset{i}{}e^{\frac{|\sigma _i\sigma _i^t^{}|^2}{2\delta \sigma ^2}}$$
(2)
where the height and the width of the Gaussian W and $`\delta \sigma `$ are chosen in order to provide a reasonable compromise between accuracy and efficiency in exploring the FES, as we will show below. The component of the forces coming from the Gaussian will discourage the system from revisiting the same spot and encourage an efficient exploration of the FES. As the system diffuses though the FES, the Gaussian potentials accumulate and fill the FES well, allowing the system to migrate from well to well. After a while the sum of the Gaussian terms will almost exactly compensate the underlying FES well.

A typical example of this behaviour is shown in Figure 1, in which a dynamics in eqn. (1) is used to explore a one-dimensional potential energy surface (PES) $`V(\sigma )`$ with three wells. The dynamics starts from a local minimum that is filled by the Gaussians in $``$ 20 steps. Then the dynamics escapes from the well from the lowest energy saddle point, filling the second well in $``$ 80 steps. The second highest saddle point is reached in $``$ 160 steps, and the full PES is filled in a total of $``$ 320 steps. Hence, in the case of this example, since the form of the potential is known, it can be verified that, for large t and if the width of the Gaussians is sufficiently small with respect to the length of a typical variation of V,

$$\underset{t^{}t}{}We^{\frac{|\sigma \sigma t^{}|^2}{2\delta \sigma ^2}}V(s)$$
(3)
modulus an additive constant. The same property can be verified for any multidimensional potential if the variables $`\sigma _i`$ are allowed to vary in a finite region. Hence, we assume that eqn. (3) holds also if the method is used to explore the FES, and that the free energy can be estimated from eqn. (3) for large t. This procedure resembles the algorithm recently proposed by Wang and Landauwang , in which a time-dependent bias is introduced in order to modify the density of states to produce locally flat hystograms.

We observe an interesting isomorphism between our dynamics and the ordered deposition of the beads of a polymer chain on the FES. In fact we can regard $`\sigma _i^t`$ as the position of a monomer in the n-dimensional parameter space. Each monomer is held at a distance $`\delta \sigma `$ from the previous one and monomers repell each other with the Gaussian term of eqn. (2). Neglecting terms of order $`\delta \sigma ^2`$ each monomer is placed by eqn. (1) in the position of free enegy minimum compatible with the requirement of minimum overlap with previous beads. Hence, the polymer chain generated in this manner fills the wells in the FES. This isomorphism helps to clarify how our approach works and why it can be made more precise by reducing $`\delta \sigma `$.

The efficiency of the method in filling a well in the PES (or in the FES) can be estimated by the number of Gaussians that are needed to fill the well. This number is proportional to $`\left(1/\delta \sigma \right)^n`$, where n is the dimensionality of the problem. Hence, the efficiency of the method scales exponentially with the number of dimensions involved. If n is large, the only way to obtain reasonable efficiencies is to use Gaussians with a size comparable to that of the well. On the other hand, a sum of Gaussians can only reproduce features of the FES on a scale larger than $`\delta \sigma `$. A judicious choice of $`\mathrm{\Delta }s_i`$, W and $`\delta \sigma `$ will ensure the right compromise between accuracy and sampling efficiency, and the optimal height and width of the Gaussians are determined by the typical variations of the FES.

If prior information on the nature of the free-energy well is not available, the scaling parameters $`\mathrm{\Delta }s_i`$ are chosen by performing short coarse-grained dynamic runs without bias potential (see eqn. (1)). In such a case the system moves around the minimum. The scaling parameters are chosen so that the elongations have approximately the same value in all directions. This amounts to an empirical form of preconditioning which makes the FES minimum nearly spherical in n dimensions and easy to fill with n-dimensional Gaussians.

The metastep length $`\delta \sigma `$ determines the efficiency of the method in a manner that is exponential in $`n`$ and it should be made as large as possible. However, the dynamics in eqn. (1) is able to reconstruct details of the FES only on the length scale of $`\delta \sigma `$. Moreover, an over-large value of $`\delta \sigma `$ would cause the system to jump irregularly from one basin to the other. The value of $`\delta \sigma `$ is chosen requiring the metatrajectory to remain localized in the FES minimum if the bias potential is not applied. With this choice of $`\delta \sigma `$, a single step of the dynamics in eqn. (1) cannot lead a system from a FES minimum to another, and the initial state at each new iteration can be generated from the last MD step of the previous iteration without creating major overlaps. Since the shake algorithm is used to impose a new set of $`\sigma _i^{t+1}`$ values, large forces on the constraints are generated at the beginning of each microscopic dynamics, and the initial part of the trajectory has to be discarded for the calculation of $`F_i^{t+1}`$.

The height of the Gaussians $`W`$ can be estimated as follows. We notice that we want to reach a situation in which the modified FES is flat. In such a case the forces coming from the FES and that coming from the Gaussian approximately balance each other. If we require the maximum value of the Gaussian forces to be smaller than the typical FES force, we arrive at the relation $`W/\delta \sigma e^{\frac{1}{2}}=\alpha F_i^2^{\frac{1}{2}}`$ with $`\alpha <1`$. In all cases so far studied a value of $`\alpha `$ close to $`.5`$ allows a fast escape from the local minima in the FES without significant loss of the underlying structure.

In the future adaptive ways of determining all these parameters should be considered in order to adapt the procedure in an optimal manner to the local shape of the FES.

## II Results

### II.1 Dissociation of NaCl in water

We first discuss an application of the method to the dissociation of a NaCl molecule in water. The most stable state is the dissociated one, while the undissociated contact ion pair corresponds to a metastable local minimum. Here we shall study how the system escapes from this metastable state. We model the system using the AMBER95 force fieldamber , and solvate the NaCl complex in 215 TIP3P water molecules with periodic boundary conditions. The electrostatic interaction between classical atoms is taken into account by the P3M methodphilippe . As collective coordinates, we use the distance $`r_{NaCl}`$ between Na and Cl, and, in order to take into account the dynamics of the solvation shell during the dissociation, we also use the electric fields $`V_{Na}`$ and $`V_{Cl}`$ on the Na and on the Cl due to the water molecules within $`6.5`$ Å$`,`$ of the ions. $`V_{Na}`$ and $`V_{Cl}`$ depend significantly on the hydration pattern around the Na and on the Cl ion. If, for example, a hydrogen bond to one of the two ions is formed or broken, the field on the ion changes by several kcal/mol. A dynamics of the form (1) was performed on the system, with $`\delta \sigma =0.25`$ and $`W=0.3`$ kcal/mol. The scaling parameters $`\mathrm{\Delta }s_i`$ were 0.53 Å, 32 and 32 kcal/mol for $`r_{NaCl}`$, $`V_{Na}`$ and $`V_{Cl}`$ respectively. The forces in eqn. (1) were evaluated by short MD runs with a time step of 0.7 fs performed on 6 replicas of the system. The replicas were then equilibrated for 200 MD steps at 300 K, and the forces evaluated by averaging over the following 500 MD steps. Starting from a distance of 2.6 Å, corresponding to a contact pair, the ions dissociate after 120 coarse-grained iterations. The dynamics was then continued for another 350 steps imposing a maximum separation of the ion pair of 5 Å, in order to study the structure of the free energy surface around the transition state (see Figure 2). Seven recrossings for the coarse-grained dynamics were observed. The transition state is located $`r_{NaCl}4.02`$ Å, $`V_{Na}120`$ kcal/mol and $`V_{Cl}81.5`$ kcal/mol (see Figure 2). The overall topology of the FES confirms the importance of the solvent degrees of freedom for describing the reaction, since the transition region is in transverse orientation relative to the axis of the coordinatesgeissler , and hence neither the distance nor the fields alone can provide an exhaustive description of the dissociation. The transition barrier (S in Figure 2), estimated with thermodynamic integrationConstraints2  along $`r_{NaCl}`$, is of 3.3 kcal/mol. The one-dimensional free energy as a function of $`r_{NaCl}`$ can also be obtained by integrating the three-dimensional FES ($`F_3`$) with respect to $`V_{Na}`$ and $`V_{Cl}`$ as $`F_1=1/\beta \mathrm{log}\left[𝑑V_{Na}𝑑V_{Cl}\mathrm{exp}(\beta F_3)\right]`$. This leads to a barrier of 3.4 kcal/mol. Hence, the method is able to reproduce the one-dimensional free energy profile estimated with thermodynamic integration.

We have verified that these collective coordinates identify the transition state and the reaction coordinate also in the dynamical sensechandler . To this effect we prepared a set of replicas of the system with the critical values of the collective coordinates. We then followed the procedure described in Ref.geissler , assigning to the particles values of the velocity chosen from a Maxwellian distribution at 300 K, and following the trajectory backwards and forwards in time. By this procedure, we found that, within a time scale of 200 fs, the overwhelming majority of the replicas fell into one of the two basins of attraction without any recrossing. Moreover, the time-reversed dynamics leads systematically to the opposite basin of attractionchandler . Finally, if we prepare the systems with the SHAKE algorithmallen  for a given value of the collective coordinates, and if the particle velocities are not reassigned before letting the system evolve freely, the time needed to go in either of the attraction basins is of the order of a few ps, and several recrossings are observed. This is due to the fact that the SHAKE algorithm imposes zero velocities in the direction of the constraint and it is only after thermalization, which takes place in a time scale of picoseconds, that the collective coordinates acquire sufficient velocity to evolve towards the FES minima. This is a strong confirmation that we have identified the correct transition state.

In order to characterize the solvent structure during the dissociation, we have computed the Na-Water pair correlation function at the transition state S. The coordination number, obtained by integrating the pair correlation function up to a distance of 3 Å, is 5.28. This value is consistent with the picture depicted in Ref. geissler , in which the sodium at the transition state is shown to be typically only five-fold coordinated.

### II.2 Isomerization of alanine dipeptide in water

As a second application, we used our method to explore the free energy surface of an alanine dipeptide as a function of the backbone dihedral angles $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ bolhuis . The system is described by the all-atom AMBER95 force field amber  and solvated in 287 TIP3P water molecules. The parameters of the coarse-grained dynamics of eqn. (1) are $`\delta \sigma =0.25`$ and $`W=0.25`$ kcal/mol, while the scaling parameters $`\mathrm{\Delta }s_i`$ are $`60^o`$ for both $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$. The other settings are the same as for the simulation of the Na-Cl complex. The FES of alanine dipeptide in water has been studied in detail by other authorsbolhuis ; apostolakis , who identified several minima. Our dynamics (1) is able to capture the overall features of the FES extremely quickly.

Starting from a configuration close to the state $`C_{7eq}`$bolhuis ; apostolakis , the dynamics fills this attraction basin in approximately 45 steps. The minimum is localized at $`\mathrm{\Phi }80^o`$ and $`\mathrm{\Psi }150^o`$. The estimated free energy at the saddle point is 2 kcal/mol higher than the free energy in $`C_{7eq}`$. This compares to a value of 2.1 kcal/mol obtained in Ref. bolhuis  by thermodynamic integration. Our estimate is obtained with 45 steps of coarse-grained dynamics on 6 replicas, corresponding to a total of only $`120`$ ps of mycroscopic dynamics. After another 150 steps, the dynamics fills the attraction basin of the state $`\alpha _r`$, which corresponds to the equilibrium configuration of the system, and recrosses to $`C_{7eq}`$. Finally, the dynamics visits the $`\mathrm{\Phi }>0`$ region of the configuration space after $`300`$ steps. In Figure 3 we plot the free energy as a function of $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ after 40, 100 and 200 steps of coarse-grained dynamics. Also small details of the free energy landscape, like the location of the saddle points and the presence of secondary minima (e.g. the state $`C_5`$ apostolakis ), can be detected graphically after a few steps of coarse-grained dynamics. The one-dimensional free energy profile as a function of $`\mathrm{\Psi }`$ obtained by analytic integration of the FES reproduces the profile obtained in Ref. bolhuis  within 0.2 kcal/mol for all the values of $`\mathrm{\Psi }`$.

Concerning the choice of the collective coordinates for this system, it should be remarked that the dihedral angles $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ do not provide a complete description of the dialanine isomerization reactionbolhuis , and the real reaction coordinate should take into account also the solvent degrees of freedom. Still, our results show that the method is able to reproduce the overall feature of the FES even if relevant degrees of freedom are not explicitly included in the collective coordinate space. This shows that the neglected degrees of freedom, although relevant for determining the reaction coordinate, are associated with relatively small barriers and are sampled efficiently during the microscopic dynamics for each value of $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$, in spite of the relatively short runs. This is due to the use of the replicas and to the ability of the dynamics to retrace the same region of parameter space during the dynamics.

## III Discussion

The method is based on the combination of the ideas of coarse-grained dynamicskevrikidis  in the space defined by a few collective coordinates, and on the introduction of a history-dependent biasHuber ; wang . Constructing a dynamics on a free energy surface that depends on a few collective coordinates, we simplify enormously the complexity of the problem, which depends exponentially on the number of degrees of freedom. The FES is much smoother than the underlying potential energy surface and it is topologically simpler, with a greatly reduced number of local minima. The history-dependent bias prevents the system from vistiting regions it has already explored. This term is crucial for the efficiency of the method: for example, a Langevin dynamics in the space of the collective coordinates would provide scant information on the underlying FES and it would be slower in locating global minima. This has been checked by performing a Langevin dynamics for the dissociation of a NaCl molecule in water, with the same collective variables introduced in Section 2 and adding to eqn. (1) a Gaussian noise at a temperature $`T`$ that is gradually reducedannealing . The number of evaluations of the forces required to converge to the minimum of the FES (i.e. the dissociated state) is well above 1000, if the temperature is decresed with a logarithic schedule, as required in order to ensure quasi-ergodicity in the collective coordinate spaceannealing1 . A history dependent bias potential as defined in eqn. (2) but applied in a regular MD simulation without coarse graining the dynamics would be efficient in finding escapes from the local minima, at it has been shown elsewhereHuber , but it would not provide quantitative informations about the FES. In particular, equation (3) holds only because the dynamic is performed in the collective coordinates space using the derivatives of the free energy. Hence, both the coarse-grained dynamics on the FES and the history-dependent bias are essential ingredients of the method, and only their combination allows an efficient and accurate determination of the FES.

Another advantage of our search algorithm is that it can be easily parallelized over the $`P`$ replicas, thus optimally exploiting present-day computer architectures and, moreover, does not require a very accurate determination of the forces. In fact, in model calculations on an analytic PES, we have found that the algorithm is effective even for noise levels as large as 20-30 % of the forces, since the height and the width of the gaussians is such that the coarse-grained dynamics explores a region of space of the size of a Gaussian several times. and the system is only discouraged (and not prevented) from revisiting the same spot of configuration space. Hence, the errors in the forces tend to even out during the evolution. The dynamics (1) is ”self-healing”, i.e. it is in principle capable to compensate the effect of complitely wrongly located gaussians or of a wrong additional bias potential. We can therefore use quite short runs and a small number of replicas to evaluate the forces. This is at variance with thermodynamic integrationallen  and the blue moon ensemble methodConstraints1 ; Constraints2 , where a very precise determination of the forces is needed and where treating multidimensional reaction coordinates is a daunting task. Other techniques, like the umbrella sampling or bias potential methods allen ; Voter1 ; Voter2 ; VandeVondele1 , allow an accurate estimate of the free energy surface in many dimensions. The drawback of these methods is that they are efficient only if a good approximation for the FES is known a priori, while this is not required in our method.

A last advantage of the method is its ability to provide qualitative information on the free energy of a system in a very short time. For example, the overall topology of a free energy surface can be determined by very few coarse-grained dynamics steps using ”large” Gaussians. Subsequently, our ”qualitative” knowledge of the FES can be improved using smaller Gaussians, eventually reducing the dimensionality of the problem by exploiting the topological information obtained with the large Gaussians.

MP would like to thank Yannis Kevrikidis for a very illuminating discussion on the coarse-grained method, and David Chandler for having patiently shared his understanding of the statistical mechanics of transition paths. AL thanks Joost VandeVondele for many precious suggestions.


text:
# Classical information and Mermin’s non-technical proof of the theorem of Bell

(1 Beckman Institute, Department of Electrical Engineering and Department of Physics, University of Illinois, Urbana, Il 61801
2 Beckman Institute, Department of Statistics and Department of Mathematics, University of Illinois,Urbana, Il 61801
)

## Abstract

We show that Mermin’s reasoning against our refutation of his non-technical proof for Bell-type inequalities is of limited significance or contains mathematical inconsistencies that, when taken into account, do not permit his proof to go forward. Our refutation therefore stands.

We discuss two recent notes of Mermin ,  that deal with our refutation  of Mermin’s non-technical proof  of Bell’s inequalities and other publications related to this discussion -. We first show, that the latest note  does not add any substance to our original refutation , . Secondly, we give a more detailed explanation of our original refutation.

Mermin discusses in his latest note  the admissibility of certain classical information that can be exchanged between the stations $`S_1`$ and $`S_2`$ in Einstein- Podolsky-Rosen (EPR) type experiments. All classical information is permitted except “ …any information whatever about the setting it has randomly been given in that run.” Mermin writes: “Let us turn Hess and Philipp upside down and explore the extent to which Bell’s theorem survives, not only if, following Hess and Philipp, we take advantage of properties of the detectors correlated by the time on local synchronized clock’s, but even if we allow further correlation of the detectors through direct straightforward ongoing classical communication between them.” What Mermin does not appreciate here is, that our introduction of time as an independent variable in Bell’s functions $`A,B`$ adds not only time as a variable but also adds the set of all functions of time and settings. Indeed, we have stated repeatedly and clearly that our extension of Bell’s parameter space is the addition of functions of time and settings, namely the addition of $`\mathrm{𝐭𝐢𝐦𝐞}`$ and $`\mathrm{𝐬𝐞𝐭𝐭𝐢𝐧𝐠}`$ dependent parameter random variables $`\lambda _{𝐚,t}^{}`$, $`\lambda _{𝐛,t}^{}`$, $`\lambda _{𝐜,t}^{}`$ for station $`S_1`$ and $`\lambda _{𝐚,t}^{}`$, $`\lambda _{𝐛,t}^{}`$, $`\lambda _{𝐜,t}^{}`$ for station $`S_2`$. The structure of these functions may be constituted such as to carry information about some of the history of how and when the settings were and are actually chosen and thus may contain information on the setting of the given run. Therefore these functions (parameter random variables) cannot be communicated during the course of each run to the other station. Mermin’s addition of classical information is then of very limited significance to our enlargement of the parameter space and we are back to the previous discussion whether our addition is sufficient to refute Mermin’s model in the first place. We add here further explanations to this discussion which should help to better understand our initial refutation of Mermin’s proof , .

To facilitate the discussion, we attempt to provide a one to one correspondence of Mermin’s notation and ours. We cannot complete our reasoning with Mermin’s all too abbreviated non-technical notation and therefore need to take this step. We also proceed now in two stages. We first consider only the parameters  of the original publication of Mermin which are all parameters “that in each run both particles carry to their detectors” and only later discuss the enlarged set of parameters. We therefore consider random variables $`A=\pm 1`$ in station $`S_1`$ and $`B=\pm 1`$ in station $`S_2`$ that describe the potential outcome of spin measurements and are indexed by instrument settings that are characterized by three-dimensional unit vectors $`𝐚,𝐛,𝐜`$ in both stations. Mermin introduces no precise counterpart for the values that these random variables $`A,B`$ may assume as distinguished from the actual outcomes in form of green and red detector flashes (the data). These green and red flashes correspond then also to the values that $`A,B`$ may assume: $`+1`$ which we can identify with green and $`1`$ which we can identify with red.

The key assumption of Bell  and also of Mermin in his original non-technical proof  is that the random variables $`A,B`$ depend only on the setting in the respective station and on another random variable $`\mathrm{\Lambda }`$ that carries the information with the particles that are emitted from a common source. Because of Einstein locality and the particular way EPR experiments are performed , the parameter random variables $`\mathrm{\Lambda }`$ are independent of the settings. Mermin uses instead of $`\mathrm{\Lambda }`$ instruction sets e.g. GGR meaning flash green for settings $`𝐚,𝐛`$ (which Mermin actually labels $`1,2`$) and flash red for setting $`𝐜`$ (labelled $`3`$ by Mermin). In our and Bell’s notation this means that for a particular value $`\mathrm{\Lambda }^{}`$ that the variable $`\mathrm{\Lambda }`$ may assume and that corresponds to the specific instruction set GGR we have $`A(𝐚,\mathrm{\Lambda }^{})=A(𝐛,\mathrm{\Lambda }^{})=+1`$ and $`A(𝐜,\mathrm{\Lambda }^{})=1`$. The following Table 1 summarizes the eight possible instruction sets and the nine possible different $`AB`$ products which are used by Mermin as a model for EPR-experiments.

According to Mermin’s point (i) , the columns of $`A_𝐚B_𝐚`$, $`A_𝐛B_𝐛`$ and $`A_𝐜B_𝐜`$ all have entries $`+1`$. It is clear, and this is the main point of Mermin’s argument, that each of the 8 rows of the Table 1 contains at least five entries $`+1`$ and at most 4 entries $`1`$. Each pair of settings occurs with probability $`1/9`$. Denote by $`p_1,\mathrm{},p_8`$ the probabilities that the instruction set $`\mathrm{\Lambda }`$ of rows $`1,\mathrm{},8`$ is carried by both particles in a given run. Then, no matter of whether or not we consider actual outcomes or potential outcomes, the average over all possibilities obeys

$$\frac{1}{9}[9p_1+9p_8+(54)(p_2+p_3+p_4+p_5+p_6+p_7)]\frac{1}{9}$$
(1)
as $`_{l=1}^8p_l=1`$, thus contradicting point (ii) of references , . This represents Mermin’s non-technical argument in a more formal way. As can easily be seen from Table 1, it is crucial that $`A_𝐚`$ in Table 1 is the same in each of the three entries where it appears and that $`A_𝐚=B_𝐚`$ (and similarly for $`A_𝐛`$ etc.).

We now turn to the expanded instruction set including time and setting dependent instrument parameters. We emphasize that the time dependence is crucial and all our reasoning rests on it. If it could be shown that EPR experiments contain no time dependencies and all elements of a corresponding model are also independent of time, then we have no argument. Let us therefore assume that the experiments contain time dependencies and that the theory needs to include time and setting dependent parameter random variables. It becomes then a question of which mathematical objects correspond to Mermin’s instruction sets. Following Mermin we identify now the instruction sets with instructions or operations that turn $`A`$ into $`\pm 1`$ and $`B`$ into $`\pm 1`$, i.e. in Mermin’s newest definition the instruction sets are the random variables $`A`$ and $`B`$ . We can now again consider Table 1 with the same notation except that the instruction sets refer to the expanded instruction sets. However, we have lost now the clear distinction that Einstein locality necessitates for source and instrument parameters, respectively. More importantly, it can no longer be guaranteed that $`A_𝐚`$ in Table 1 is the same in each of the three entries (and similar for $`A_𝐛`$ etc.). Furthermore, we have $`A_𝐚=B_𝐚`$ only for the same time. This was explained by us in Eqs.(1) and (2) of reference . However, this point contains some subtle reasoning and we give therefore the following elaboration. For a theory to provide valid conclusions, special care must be taken when combining, adding and/or counting different elements (actual or imagined) that are mutually exclusive. In the case that we consider there are two further subtleties. First, different elements of the theory are mutually exclusive only at the same time and second, the measurement time itself appears as a random variable. This puts restrictions on the possible measurement times which now must be different for different settings. If they are chosen to be equal, Table 1 will contain mutually exclusive alternatives that cannot simultaneously be used in mathematical operations such as counting, adding or averaging. We will return to this point below.

The instruction sets become now also time and setting dependent and the joint frequency of occurrence of the setting dependent instrument parameters may now be different in different columns of Table 1. We see that Mermin’s proof, that is based on the equality of all $`A_𝐚`$, $`A_𝐛`$ and $`A_𝐜`$ (and similarly for the $`B`$’s) in Table 1, comes to a halt. His argument  that the expanded instruction sets for a given single experiment in one time interval must be one of his eight instruction sets is true but meaningless. It has no consequences for the statistics and the associated probability measure and therefore does not guarantee Mermin’s way of counting $`+1`$ or $`1`$.

One way one still could attempt to proceed with the proof is to assume that EPR experiments are equivalent to other experiments that are made all at the same time. This does not only involve counterfactual reasoning but contradicts the assumed fact of possible time dependencies. The only other way to proceed with Mermin’s proof is to add to the possible outcomes of the Table 1 the outcomes for all the nine settings as they would have been obtained when taken at the same time and to include them into the counting of positive and negative $`AB`$ products. This procedure, however, again does not only contain counterfactual reasoning, but contradicts the facts because then the theory $`\mathrm{𝐜𝐨𝐮𝐧𝐭𝐬}`$ also impossible outcomes by adding mutually exclusive alternatives. This, in turn, has the consequence that the set of elements so counted is nine times larger than the set of experimental results. Any reasonable procedure to establish a theoretical model for an experiment will require a one to one correspondence of what is added/counted/averaged in theory to what is added/counted/averaged in the experiment. We conclude that Mermin’s proof, as well as those of other Bell type inequalities, apply only to a stationary situation that may have nothing to do with past and present EPR experiments .

We finish this discussion with an example of what can happen if imagined alternatives are added (or counted) that involve, simultaneously (time!), mutually exclusive alternatives. Mutually exclusive events are discussed in every probability text well before the concept of countable additivity of probability measures is introduced. When probabilities of mutually exclusive events are added, it is understood that the simultaneous occurrence of mutually exclusive events is impossible. Our example below has nothing to do with quantum mechanics. It just illustrates the rules of standard probability theory that must be obeyed. Consider a coin with a little magnet inside and a hidden bigger magnet under ground. The experimenter can set the underground magnet N or S and correspondingly the coin will show with higher probability head or tail, respectively, in any given experiment. If one performs a large number of coin tosses, the result will be biased according to the choices that are made for the underground magnet (N, S). Now in analogy to the above addition of nine terms for each single run of the experiment, construct a theory in which one simply adds or counts, simultaneously, the mutually exclusive alternatives of head or tail. Count 1 if a head shows and count 0 if a tail shows in each of the potential outcomes. Then for each toss the count of potential heads plus the count of potential tails equals 1. Because we have twice as many potential outcomes than actual tosses, the likelihood of heads equals exactly $`1/2`$. Hence such a theory necessarily concludes that this is a fair game. This theory, however, has combined simultaneously mutually exclusive alternatives. It necessarily counts twice as many elements than any given set of experimental results contains and is therefore not admissible. This is, of course, immediately obvious to anyone.

As mentioned above, the reasoning becomes more subtle when time plays a role in the random variables and the elements become mutually exclusive at the same time, hence impossible. Then care must be exercised in the possible choices of the random variables related to measurement time or the same mistake as outlined in the above example could be made. One way to proceed safely is to label the measurement times by the actual setting that is chosen e.g. $`t_{\mathrm{𝐚𝐛}}`$, $`t_{\mathrm{𝐚𝐜}}`$ etc. and regard measurement times with different label as different random variables. Such labelling clearly shows the special role of measurement time as random variable.

We present now our argument in a still more technical way to bring the problems of Mermin’s reasoning into focus. Consider the random variables $`A_𝐢=A_𝐢(\mathrm{\Lambda })`$ and $`B_𝐣=B_𝐣(\mathrm{\Lambda })`$ with $`𝐢,𝐣=𝐚,𝐛,𝐜`$ and which, just as above, assume the values $`\pm 1`$. Let $`\mu `$ be the probability measure that assigns weight $`p_l`$ to the $`\mathrm{\Lambda }`$ in the $`l`$’th row. Then according to Mermin’s (i) and (ii) we have

* $$A_𝐣(\mathrm{\Lambda })=B_𝐣(\mathrm{\Lambda })\text{ , }𝐣=𝐚,𝐛,𝐜$$

and
* the average over all products $`A_𝐢B_𝐣`$ is about 0.

The average for actual as well as potential outcomes equals about:

$$\frac{1}{9}(A_𝐚(\mathrm{\Lambda })+A_𝐛(\mathrm{\Lambda })+A_𝐜(\mathrm{\Lambda }))(B_𝐚(\mathrm{\Lambda })+B_𝐛(\mathrm{\Lambda })+B_𝐜(\mathrm{\Lambda }))𝑑\mu $$
(2)
If we now consider the original instruction sets that are independent of the setting then we have by Mermin’s condition (i)

$$A_𝐚(\mathrm{\Lambda })+A_𝐛(\mathrm{\Lambda })+A_𝐜(\mathrm{\Lambda })=B_𝐚(\mathrm{\Lambda })+B_𝐛(\mathrm{\Lambda })+B_𝐜(\mathrm{\Lambda })$$
(3)
for any given value of $`\mathrm{\Lambda }`$ (which could be e.g. $`\mathrm{\Lambda }^{}`$). This gives for Eq.(2)

$$\frac{1}{9}(A_𝐚(\mathrm{\Lambda })+A_𝐛(\mathrm{\Lambda })+A_𝐜(\mathrm{\Lambda }))^2𝑑\mu \frac{1}{9}$$
(4)
contradicting Mermin’s condition (ii).

However, if we use now the enlarged instruction set which includes both time and setting dependent parameters, then the $`AB`$ products that involve different pairs of settings involve different times and the parameter random variables $`\mathrm{\Lambda }`$ in Eq.(3) may all be different. In fact, just as above, it can no longer be guaranteed that $`A_𝐚`$ is the same in each of the three products $`A_𝐚B_𝐚`$, $`A_𝐚B_𝐛`$ and $`A_𝐚B_𝐜`$ and that $`A_𝐚=B_𝐚`$, independent of time. Again, Mermin’s proof comes to halt.

One final comment. It is obvious that the above argument still works for an arbitrary probability measure $`\mu `$, that is for any probability distribution governing $`\mathrm{\Lambda }`$ of the original instruction set . However, if the enlarged instruction set is considered, Mermin’s proof cannot be completed unless it can be guaranteed that $`A_𝐚`$ (or at least its probability distribution) is the same in each of the three products $`A_𝐚B_𝐚`$, $`A_𝐚B_𝐛`$ and $`A_𝐚B_𝐜`$. More specifically one needs to take into account that there are $`\mathrm{𝐧𝐢𝐧𝐞}`$ joint distributions in operation namely $`\mu _{\mathrm{𝐚𝐚}}`$, $`\mu _{\mathrm{𝐚𝐛}}`$, $`\mu _{\mathrm{𝐚𝐜}}`$,…, $`\mu _{\mathrm{𝐜𝐜}}`$. Thus, unless it can be guaranteed that the first marginal distribution of the three distributions $`\mu _{\mathrm{𝐚𝐚}}`$, $`\mu _{\mathrm{𝐚𝐛}}`$ and $`\mu _{\mathrm{𝐚𝐜}}`$ is the same, the proof cannot go forward.

We conclude that our refutation of Mermin’s non-technical proof for the theorem of Bell stands.

Acknowledgement: the work was supported by the Office of Naval Research N00014-98-1-0604, the MURI supported by ONR and by NSF Grant CCR-0121616. We thank Salvador Barraza-Lopez for providing sources regarding language, valuable discussions and help with the manuscript and Istvan Berkes for valuable discussions.


text:
# Polarization of superfluid turbulence

## Abstract

We show that normal fluid eddies in turbulent helium II polarize the tangle of quantized vortex lines present in the flow, thus inducing superfluid vorticity patterns similar to the driving normal fluid eddies. We also show that the polarization is effective over the entire inertial range. The results help explain the surprising analogies between classical and superfluid turbulence which have been observed recently.

PACS 67.40Vs, 47.37.+q

Our concern is the experimental evidence that turbulent helium II appears similar to classical hydrodynamic turbulence. For example, Smith et al found that the temporal decay of helium II turbulence behind a towed grid is the same as that expected in an ordinary fluid. Maurer & Tabeling observed the classical Kolmogorov $`k^{5/3}`$ dependence of the energy spectrum on the wavenumber $`k`$ in helium II continuously agitated by rotating propellers. In other experiments it was found that if helium II is forced at high velocity along pipes and channels the same pressure drops are detected which are observed in a classical liquid. Furthermore, when a sphere moves at high velocity in helium II, the same drag crisis is measured that occurs if the fluid is air.

These apparently classical results are surprising because helium II is a quantum fluid. According to Landau’s two-fluid theory, it consists of the intimate mixture of an inviscid superfluid component and a viscous normal fluid component. The latter is similar to a classical Navier - Stokes fluid, so, when made turbulent, it consists of eddies of various sizes and strengths. On the contrary quantum mechanics constrains the rotational flow of the superfluid to discrete quantized vortex lines, each vortex characterized by the same quantum of circulation $`\mathrm{\Gamma }=9.97\times 10^4cm^2/sec`$. Unlike what happens in a classical Euler fluid, superfluid vortex lines can reconnect with each other. They also interact with the normal fluid via a linear mutual friction force, and, when helium II is made turbulent, they form a disordered, apparently random tangle. Clearly the experimental results described above call for an explanation in terms of the basic physical ingredients of the problem (normal fluid eddies and superfluid vortices) and their interaction. What is remarkable is that these classical aspects of helium II turbulence are observed to be independent of temperature, whereas the relative proportion of superfluid and normal fluid is a strong function of temperature.

A possible solution of the puzzle is that both fluids independently obey the classical Kolmogorov scaling. A few years ago Nore et al solved the Gross-Pitaevskii equation for a tangled system of vortices and noticed that the energy spectrum follows a $`k^{5/3}`$ dependence. Unfortunately the effect only took place over a small range in $`k`$ space and further calculations performed using the vortex filament model did not reproduce it. Another difficulty with this argument is that one has also to explain why the mutual friction coupling between the two fluids does not change the shape of the spectra which the two fluids would have independently of each other.

A second explanation of the experiments has been suggested by Vinen. He argued that at spatial scales larger than the average separation $`\delta `$ between the quantized vortex lines, the superfluid and normal fluid are coupled by a small degree of polarization of the almost random tangle of superfluid vortex lines. If that is the case, on these scales helium II behaves as a single fluid of density $`\rho =\rho _s+\rho _n`$, consistently with the experiments ($`\rho _s`$ and $`\rho _n`$ are the superfluid and normal fluid density respectively).

The aim of this letter is to support Vinen’s argument with quantitative evidence of polarization. First we shall introduce some simple models which, although very idealized, capture the essential physics of polarization. Secondly we shall look for evidence of polarization by direct numerical simulation.

Our first model is concerned with the reaction of superfluid vortices to a normal fluid shear. Consider a row of point vortices of alternating circulation $`\pm \mathrm{\Gamma }`$ initially set along the $`x`$ axis at distance $`\delta `$ from each other. We assume that the normal fluid is $`𝐯_n=V_ncos(ky)\widehat{𝐱}`$. The governing equation of motion of a vortex point is

$$\frac{dy}{dt}=\pm \alpha V_ncos(ky),$$
(1)
where $`\alpha `$ is a known temperature dependent friction coefficient. The solution of Eq. (1) corresponding to the initial condition $`y(0)=0`$ is

$$y(t)=\frac{2}{k}\left(\frac{\pi }{4}+tan^1(e^{\pm \alpha kV_nt})\right).$$
(2)
Given enough time ($`t\mathrm{}`$), positive and negative vortices will reach stable locations $`y_{\mathrm{}}=\pm \pi /2k`$ respectively. In a turbulent normal flow, however, the shear does not last longer than few times the turnover time $`\tau 1/\omega _n1/kV_n`$. Since $`\alpha `$ is small (it ranges from $`0.037`$ at $`T=1.3\mathrm{K}`$ to $`0.35`$ at $`T=2.16\mathrm{K}`$) we have $`y(\tau )\pm b`$ where $`b=\alpha /k`$. Within the lifetime of the shear we have thus created a separation $`2b`$ between positive and negative vortices, that is to say we have polarized the initial configuration. The velocity of this Karman - vortex street is approximately $`V_s\pi \mathrm{\Gamma }/2\delta `$ for $`k\delta <<1`$ in the direction along the $`x`$ axis where the normal fluid (which induced the polarization in the first place) is stronger. The result suggests that it is not necessary to create extra vortex lines to generate a superfluid pattern that mimics the normal fluid one - rearranging existing vortices is enough. We also notice that the induced polarization is proportional to $`\alpha `$.

Our second model is concerned with the expansion of favourably oriented superfluid vorticity. We think of the superfluid vortex tangle as a collection of vortex rings of radius approximately equal to the average separation of vortices in the tangle, $`R_0\delta `$. We assume for simplicity that the rings are on the $`x,y`$ plane with equal numbers of rings oriented in the $`\pm z`$ directions. Depending on whether they have positive or negative orientation, the rings move along $`\pm z`$ with self induced speed given by

$$V_{R_0}=\frac{\mathrm{\Gamma }}{4\pi R_0},$$
(3)
where $`=\mathrm{ln}(8R_0/a_0)1/2`$ is a slowly varying term and $`a_010^8\mathrm{cm}`$ is the vortex core radius. Now we apply a normal fluid velocity $`V_n`$ in the $`z`$ direction. The radius $`R`$ of a ring is determined by

$$\frac{dR}{dt}=\frac{\gamma }{\rho _s\mathrm{\Gamma }}(V_nV_R),$$
(4)
where $`\gamma `$ is a known friction coefficient and $`\gamma /\rho _s\mathrm{\Gamma }\alpha `$ at almost all temperatures of interest. Eq. (4) shows that $`V_n`$ selectively changes radius and velocity of vortex rings moving in opposite directions. A ring which grows (shrinks) by an amount $`\delta R=\alpha \delta t(V_nV_{R_0})`$ in time $`\delta t`$ slows down (speeds up) by an amount $`\delta V_R=V_{R_0}\delta R/R_0`$. In this way a superflow is induced in the same direction of the normal fluid which induced the polarization in the first place. A simple estimate of the spatial averaged magnitude of this superflow yields $`V_s3V_{R_0}\delta R/R_0`$.

Our third model is concerned with the rotation of existing superfluid vorticity, We represent a superfluid vortex line as a straight segment pointing away from the origin and study how its orientation is changed by a normal fluid rotation about the $`z`$ axis. Using spherical coordinates $`(r,\theta ,\varphi )`$, we assume that the vortex is initially in the plane $`\theta =\pi /2`$. The normal fluid’s velocity is $`𝐯_n=(0,0,\mathrm{\Omega }r\mathrm{sin}\theta )`$ and the motion of the vortex segment is determined by

$$\frac{d\theta }{dt}=\alpha \mathrm{\Omega }sin(\theta ),$$
(5)
and $`dr/dt=0`$ and $`d\varphi /dt=0`$. The solution to Eq.(5) is

$$\theta (t)=2tan^1(e^{\alpha \mathrm{\Omega }t}),$$
(6)
with $`r`$ and $`\varphi `$ constant. Given enough time, the vortex segment will align along the direction of the normal fluid rotation ($`\theta 0`$ for $`t\mathrm{}`$), but the lifetime $`\tau `$ of the eddy is only of the order of $`\tau 1/\mathrm{\Omega }`$, so the vortex can only turn to the angle

$$\theta (\tau )\pi /2\alpha .$$
(7)
Despite the smallness of the angle, the effect is sufficient to create a net polarization of the tangle in the direction of the normal fluid’s rotation, provided that there are enough vortices. The following argument shows how this is possible. The normal fluid is like a classical viscous Navier - Stokes fluid, and, if left to itself, its spectrum $`E_k`$ would obey Kolmogorov’s law

$$E_k=Cϵ^{2/3}k^{5/3},$$
(8)
Eq. (8) is valid in the inertial range $`1/\mathrm{}_0<k<1/\eta `$ in which big eddies break up into smaller eddies, transferring energy to higher and higher wavenumbers without viscosity playing a role. Here $`k`$ is the magnitude of the three dimensional wavevector, $`ϵ`$ is the rate of energy dissipation per unit mass, $`C`$ is a constant of order unity, $`\mathrm{}_0`$ is the integral length scale (the scale at which energy is fed into the energy cascade) and $`\eta `$ is the Kolmogorov scale at which kinetic energy is dissipated by the action of viscosity. In reality the normal fluid is not alone but is forced by the quantized vortex filaments. We know little of the effects of this forcing (it has been studied only for very simple geometries) so, for lack of further information, we assume that the classical relation Eq. (8) is valid for the normal fluid.

The quantity which is often used to describe the intensity of the superfluid vortex tangle is the vortex line density $`L`$ (length $`\mathrm{\Lambda }`$ of vortex line per volume $`𝒱`$) because it is easily measured by detecting the attenuation of second sound. From $`L`$ one infers the average separation between vortices, $`\delta L^{1/2}`$. The quantity $`\mathrm{\Gamma }L`$ can be interpreted as the total rms vorticity of the superfluid. Note that the net amount of superfluid vorticity in a particular direction can be much less than $`\mathrm{\Gamma }L`$ (even zero, if the tangle is randomly oriented).

The key question is whether, as a result of mutual friction, sufficient quantized vortex lines can re-orient themselves within a normal fluid eddy of wavenumber $`k`$ so that the resulting net superfluid vorticity matches the vorticity $`\omega _k`$ of that eddy. The process must take place in a time scale shorter than the typical lifetime of the eddy, which is of the order of few times the turnover time $`1/\omega _k`$. At this point we use the result of the third model, for which the initial condition $`\theta (0)=\pi /2`$ represents the average case. In doing so we remark that, if the initial orientation of the vortex is toward the origin rather than away from it, then the vortex segment turns to $`\theta (\tau )\pi /2+\alpha `$ rather than $`\pi /2\alpha `$ but still contributes to positive vorticity in the $`z`$ direction. Therefore in the time $`1/\omega _k`$, re-ordering of existing vortex lines creates a net superfluid vorticity $`\omega _s`$ of the order of $`\alpha L\mathrm{\Gamma }/3`$ in the direction of the vorticity $`\omega _k`$ of the driving normal fluid eddy of wavenumber $`k`$. Since $`\omega _k`$ is approximately $`\omega _k=\sqrt{k^3E_k}`$ we have $`\omega _k=C^{1/2}ϵ^{1/3}k^{2/3}`$. Matching $`\omega _s`$ and $`\omega _k`$ would then require

$$\frac{1}{3}\alpha \mathrm{\Gamma }\mathrm{L}\mathrm{C}^{1/2}ϵ^{1/3}\mathrm{k}^{2/3}.$$
(9)
The normal fluid vorticity increases with $`k`$ and is concentrated at the smallest scale ($`k1/\eta `$), so a vortex tangle with a given value of $`L`$ may satisfy the above equation only up to a certain critical wavenumber $`k_c`$. Substituting $`ϵ=\nu _n^3/\eta ^4`$ where $`\nu _n`$ is the normal fluid’s kinematic viscosity (the viscosity of helium II divided by $`\rho _n`$), we obtain $`(\delta /\eta )=C^{1/4}(\alpha /3)^{1/2}(\mathrm{\Gamma }/\nu _n)^{1/2}(\eta k_c)^{1/3}`$. If $`k_c1/\eta `$ then

$$\frac{\delta }{\eta }=C^{1/4}\left(\frac{\alpha }{3}\right)^{1/2}\left(\frac{\mathrm{\Gamma }}{\nu _n}\right)^{1/2}$$
(10)
In the temperature range of experimental interest $`\mathrm{\Gamma }/\nu _n`$ ranges from $`0.43`$ at $`T=1.3\mathrm{K}`$ to $`5.86`$ at $`T=2.15\mathrm{K}`$, so $`\delta /\eta =𝒪(1)`$ and we conclude that matching of the two vorticities ($`k_c1/\eta `$) is possible throughout the inertial range.

Due to the computational cost, numerical simulations of superfluid turbulence do not produce vortex tangles dense enough to cover the range $`k<1/\mathrm{}`$ and determine unambiguously the dependence of the energy spectrum on $`k`$ in this range. To make progress in the problem and confirm the above arguments we study the reaction of the superfluid vortex tangle to a single scale ABC normal flow given by $`𝐯_n=(A\mathrm{sin}(kz)+C\mathrm{cos}(ky),B\mathrm{sin}(kx)+A\mathrm{cos}(kz),C\mathrm{sin}(ky)+B\mathrm{cos}(kx))`$ where $`k=2\pi /\lambda `$ is the wavenumber, $`\lambda `$ is the wavelength and $`A`$, $`B`$ and $`C`$ are parameters. ABC flows are solutions of the Euler equation and the forced Navier - Stokes equation and have been used as idealized model of eddies in fluid dynamics, magneto-hydrodinamics and superfluid hydrodynamics. For the sake of simplicity we take $`A=B=C`$ and $`\lambda =1`$.

Following Schwarz, we represent a superfluid vortex filament as a space curve $`𝐬=𝐬(\xi ,t)`$ where $`\xi `$ is arclength and $`t`$ is time. Neglecting a small transverse friction coefficient, the curve moves with velocity

$$\frac{d𝐬}{dt}=𝐯_{si}+\alpha 𝐬^{}\times (𝐯_n𝐯_{si}).$$
(11)
where $`𝐬^{}=d𝐬/d\xi `$ and the self induced velocity $`𝐯_{si}`$ is given by the Biot - Savart integral

$$𝐯_{si}=\frac{\mathrm{\Gamma }}{4\pi }\frac{(𝐫𝐬)\times \mathrm{𝐝𝐫}}{|𝐫𝐬|^3}.$$
(12)
The calculation is performed in a cubic box of volume $`𝒱=1cm^3`$ with periodic boundary conditions. The numerical technique is standard and the details of our algorithm, including how to perform vortex reconnections, have been published elsewhere.

We start the calculation with $`N=50`$ superfluid vortex rings set at random positions and orientation and integrate in time at a variety of temperatures ($`\alpha =0.1`$, $`0.5`$ and $`1.0`$) and normal fluid’s velocities ($`A=0.01`$, $`0.1`$, $`1.0`$ and $`10.0cm/sec`$). The vortex length ($`\mathrm{\Lambda }=76.8cm`$ at $`t=0`$) increases or decreases depending on whether the ABC flow is strong enough to feed energy into the normal fluid via instabilities of vortex waves (for example, for $`\alpha =1.0`$, the final length $`\mathrm{\Lambda }`$ is as high as $`781.3cm`$ at $`A=10.0cm/sec`$, and as low as $`56.67cm`$ at $`A=0.01cm/sec`$) The rings interact with each other and with the normal fluid, get distorted, reconnect, and soon an apparently random tangle is formed (see figure 1).

The quantity $`<\mathrm{cos}(\theta )>=<𝐬^{}\widehat{𝝎}_n>`$, which we monitor during the evolution, gives us the tangle - averaged projection of the local tangent to a vortex in the direction of the local normal fluid vorticity, $`\widehat{𝝎}_n=(1/\omega _n)𝝎_n`$, where $`𝝎_n=\times 𝐯_n`$. At $`t=0`$ $`<\mathrm{cos}(\theta )>=0`$ due to the random nature of the initial state, and it is apparent from figure 2 that $`<\mathrm{cos}(\theta )>`$ increases with time, no matter whether $`\mathrm{\Lambda }`$ decreases or increases.

The results are analyzed in figure 3. From the simple models descrived above we expect that the polarization induced by the normal fluid vorticity is proportional to $`\alpha `$. We also know from the discussion above that we should restrict the analysis to times $`t<\tau `$ where $`\tau =1/\omega _n`$ with $`\omega _n=\sqrt{3}Ak`$ is the lifetime of the normal fluid eddy, which we assume to be the same as the eddy’s turnover time. It is apparent from the figure that, no matter whether the tangle grows or decays, approximately the same polarization takes place for $`t/\tau <1`$.

In conclusion we have put the theory of superfluid turbulence on firmer ground. Using simple models which capture the essential physical mechanisms of polarizaton and then a numerical simulation, we have shown that, within the lifetime of a normal fluid eddy of wavenumber $`k`$, superfluid vortex lines can rearrange themselves so that the superfluid vorticity and the normal fluid vorticity are aligned. Provided that enough vortex lines are present, vorticity matching should take place over the entire inertial range, up to wavenumbers $`k`$ of the order of $`1/\mathrm{}`$.

Our result has theoretical and experimental implications. Numerical simulations of vortex lines driven by normal fluid turbulence show a $`k^1`$ superfluid energy spectrum in the accessible region $`k1/\delta `$. More intense (hence computationally expensive) vortex tangles should be investigated to explore the region $`k<<1/\delta `$ where, as a consequence of our result, we predict the classical $`k^{5/3}`$ dependence. Clearly an important issue which must be investigated is the nonlinear saturation of the polarization process. On the experimental side, our result supports the use of helium II to study classical turbulence. This has been done recently by Skrbek et al who exploited the physical properties of liquid helium to study the decay of vorticity on an unprecedented wide range of scales.

ACKNOWLEDGMENTS

We are indebted to W.F. Vinen who suggested to study some of the above described models of the interaction between normal fluid and quantized vortex lines.

FIGURE CAPTIONS

Figure 1. Vortex configuration at $`t=0.123sec`$ for $`\alpha =0.5`$ and $`A=1.0cm/sec`$.

Figure 2. Average polarization $`<\mathrm{cos}(\theta )>`$ versus time $`t`$ computed for different values of $`A`$ and $`\alpha `$.

Figure 3. $`<\mathrm{cos}(\theta )>/\alpha `$ versus scaled time $`t/\tau `$ where $`\tau `$ is the eddy’s lifetime (same symbols as in Figure 1).