text:
# The surface brightness of dark matter: unique signatures of neutralino annihilation in the Galactic halo.

## I Introduction

Determining the nature of dark matter is of fundamental importance to both Astronomy and Particle Physics. Both theory and observational data currently favour a universe with a matter density that is dominated by non-baryonic particles. Many candidates have been proposed: some are known to exist, others are more speculative (e.g., Ref.  and references therein). Structure formation in a universe dominated by cold dark matter (CDM) has been extensively tested against observations and the model has proven highly successful at reproducing the large scale properties and distributions of galaxies . On the non-linear scales of galactic halos it remains to be confirmed wether the model can successfully reproduce the observational data .

Direct detection in the laboratory is the ultimate technique for verifying the existence of dark matter particles (see, Ref. ). However, even the most popular candidate for dark matter, the neutralino, has a cross section that spans many orders of magnitude and the current laboratory searches are only just becoming sensitive to the cosmologically interesting parameter range. Presently, Astronomical observations provide the best insights into the nature of the dark matter, furthermore direct detection relies on the existence of a smooth component of dark matter.

Within the next few years indirect detection of neutralinos will provide interesting constraints on their possible cross-section and masses. Neutralino-neutralino annihilation produces observable photons (as well as a host of other particles) that may be observed as a diffuse gamma-ray background from the halo surrounding the Milky Way as discussed in Refs. , and more recently, in Refs. .

Renewed interest in these predictions has recently arisen because of an unexplained component of diffuse high energy photons in the Egret data (e.g. Ref. ), and also the possibility of an excess from the center of the Galaxy itself , unexpected clumpy emission and the unresolved “discrete sources” . Progress in this area will result from several new and sensitive gamma-ray surveys such as GLAST  and VERITAS .

The efficiency of the annihilation process is strongly dependent on the local density and the cross-section of the neutralino. Many authors have calculated the expected flux from the Galactic halo using simple models for the expected mass distribution of neutralinos within the Galaxy (see, e.g., Refs. ) or from its satellites .

Advances in computational Cosmology have lead to several recent breakthroughs that have direct relevance to the detection of dark matter. In particular, the numerical resolution that can be achieved using parallel computational techniques is now sufficient to study the internal structure of dark matter halos that form within a cosmological context. The results of these simulations have important implications for indirect (and direct) detection of dark matter candidates. Most significantly for particle-particle annihilation, we are now confident that the central density profile of CDM halos follows a singular power law down to small scales . Thus we may expect a point like source of mono-chromatic gamma-rays emanating from the center of the Milky Way, where the annihilation rate will be very high.

A second fundamental prediction of the CDM model is that previous generations of the merging hierarchy survive within halos . Halos that accrete into larger systems may be tidally stripped of most of their mass, however their dense central cores survive and continue to orbit within the parent halos. This may present some problems for the CDM model since the predicted number of satellites within the Milky Way’s halo is 50–100 times as many as observed . If the CDM model is correct, then only a fraction of these satellites must have formed stars and most of the substructure remains as dark objects within the Galactic halo.

The possibility of an enhanced gamma-ray background from dark matter substructure was explored by Bergström et al. , who made simple assumptions as to the mean density and abundance of such clumps. We can now use the high resolution N-body simulations to directly measure these quantities. The simulations also allow us to study the influence of the halo shapes on the diffuse gamma-ray background and the intensity of the central halo emission that arises from the singular dark matter density profiles. This paper is organized as follows. In Section II we explore the gamma-ray background that results from the smoothly distributed component of dark matter using both analytical and simulated halos. In Section III we focus on the substructure within halos. Our conclusions are summarised in Section IV.

## II The sky distribution of the gamma-ray background

In what follows we will consider a flux of photons (or other particles) that are a by-product of the annihilation of dark matter particles within the smooth component of dark matter that surrounds the Galaxy. It is not our intention to discuss the details of neutralino interactions, a complete overview on these processes (and super-symmetric matter in general) can be found in Refs. .

### A Model neutralino halos

We calculate the gamma-ray flux along a given line of sight through a spherically symmetric galactic halo using:

$$\varphi (\psi )=\frac{K}{4\pi }_{Lineofsight}\rho ^2(l)𝑑l(\psi )$$
(1)
where $`\psi `$ is the angle between the direction of galactic center and observation; $`\rho `$, the density of dark matter at distance $`l`$ from the observer. We have summed up the dependence of the flux on neutralino mass and interaction cross section in the constant $`K`$. This is enough scope for the present discussion - it is straightforward to take our results and input a neutralino cross-section, $`<\sigma v>`$, and mass, $`M_\chi `$ to determine the absolute gamma-ray flux (where $`K`$ is defined to be $`<\sigma v>/M_\chi ^2`$). Our results can also be used to infer the sky distribution of other products of the annihilation, such as neutrinos or positrons.

The line of sight distance, $`l`$, is related to the radial distance from the halo center, $`r`$, via

$`r^2=l^2+R_o^22lR_o\mathrm{cos}(\psi )`$

where $`R_o`$ is our galacto-centric distance, taken here to have the IAU standard value of $`R_o=8.5`$kpc , and $`\psi `$ is related to galactic coordinates ($`\mathrm{}`$, $`b`$) through

$`\mathrm{cos}(\psi )=\mathrm{cos}(\mathrm{})\mathrm{cos}(b)`$.

For the halo density profile, $`\rho (r)`$, we take the latest results from the highest resolution numerical simulations of galactic halos carried out to-date . These authors simulated 6 different galactic mass halos with force resolution of $`0.5`$ kpc and mass resolution of $`10^6M_{}`$. (Throughout the paper we will use the Hubble constant value of $`H_o=100h\mathrm{kms}^1\mathrm{Mpc}^1`$ and $`h=0.5`$; as adopted for the simulations.) The best fitting density profile to these data is (subscript moore):

$$\rho _{moore}(r)=\frac{\rho _{moore}^{}}{(r/a)^{1.5}(1+(r/a)^{1.5})}.$$
(2)
Where $`r`$ is the distance from the halo center, $`a=r_{200}/c_{moore}`$ the scale radius for halos of mass $`1\times 10^{12}M_{}`$. The virial radius of our fiducial Galactic halo, $`r_{200}300`$ kpc, is defined as the radius of a sphere at which the mean overdensity is 200 times the cosmological mean density. (A central density profile of slope -1.5 on galactic scales was also found by Jing & Suto  and confirmed as an asymptotic slope by Ghigna et al. .)

We also compare this profile with that determined by Navarro et al.  using a sequence of lower resolution studies (subscript nfw) (the main difference being that the central dark matter density profile has a slope of -1):

$$\rho _{nfw}(r)=\frac{\rho _{nfw}^{}}{(r/a)(1+r/a)^2}$$
(3)
and the modified isothermal profile with a constant density core (subscript is):

$$\rho _{is}(r)=\frac{\rho _{is}^{}}{[1+(r/a)^2]^{3/2}}.$$
(4)
The scale radius, $`a`$, is determined directly from the numerical simulations, except for the modified isothermal model which we normalise to match the observational rotation curve data (as in Ref. ); $`a_{is}=24.3`$ kpc, $`a_{nfw}=27.7`$ kpc and $`a_{moore}=33.2`$ kpc (this radius is directly related to the concentration parameter, $`c=\mathrm{r}_{200}/a`$). We normalise each density profile such that the peak circular velocity, $`v_{peak}=200`$ km/s (the maximum of the $`v_c=\sqrt{GM/r}`$ curve), which gives: $`\rho _{is}^{}=4.96\times 10^6\mathrm{M}_{}`$ kpc-3, $`\rho _{nfw}^{}=5.11\times 10^6\mathrm{M}_{}`$ kpc-3 and $`\rho _{moore}^{}=1.64\times 10^6\mathrm{M}_{}`$ kpc-3. We plot the effective circular velocity profiles and density profiles of these model halos in Fig. 1(a) and Fig. 1(b) respectively.

In Fig. 2 we plot the flux, $`\varphi `$, along the line of sight through a spherical Milky Way halo using the above density profiles as the observer looks towards the Galactic center at $`\psi =0^{}`$, to the Galactic anticenter at $`\psi =180^{}`$. As expected, the central annihilation flux depends strongly on the form of the inner density profile. At an angle of five degrees from the Galactic center, the ratio of fluxes from the three different profiles, moore:nfw:is are 1000:100:1.

The peak central value depends upon the distance from the Galactic center that we are willing to consider integrating from – the flux slowly diverges for the density profile in Eq. (2). However, within a given radius, most of the neutralinos would have self annihilated leaving a tiny constant density core. We can estimate the size of this core using $`(n\sigma v)^1=t_h`$, where $`t_h10`$ Gyrs is the Hubble time. Taking a typical cross section, $`\sigma v=10^{30}\mathrm{cm}^3`$ s-1, and adopting the Moore et al. density profile we find that the annihilation radius within the Milky Way is approximately $`4\times 10^7`$ parsecs $`10^{12}r_{200}`$.

The total flux that arises within 5 degrees of the Galactic center using the Moore et al. density profile is a factor of 20 larger than that found using the NFW profile (both integrated down to the annihilation radius calculated above).

### B Comparison with high resolution CDM simulations

We can use the numerical simulations to compare directly with the above predictions that were obtained assuming spherical symmetry. We refer the reader to Moore et al. (Ref. ) for details of the numerical simulations.111Images, data and movies of these dark matter simulations can be downloaded from http://www.nbody.net. To construct the expected gamma-ray sky maps we choose a simulated dark matter halo at a redshift z=0 that has a peak circular velocity of $`200\mathrm{k}\mathrm{m}\mathrm{s}^1`$ and a total mass, within the virial radius, $`r_{200}=300\mathrm{k}\mathrm{p}\mathrm{c}`$, of $`1\times 10^{12}\mathrm{M}_{}`$. This simulated halo is from the “Local Group” simulation and is close to our fiducial Milky Way cold dark matter halo that we adopted in the previous section.

N-body simulations attempt to simulate a collisionless fluid of dark matter using discrete massive particles. We calculate the local density at the position of each particle by averaging over its nearest 64 neighbours. The observer is placed 8.5 kpc from the halo center (defined using the most bound particle in the simulation) and we sum up the flux of annihilation products along each line of sight using the discrete equivalent to Eq. ( 1):

$$\mathrm{\Phi }(\mathrm{},b)=\frac{K}{\mathrm{\Omega }}\underset{LOS}{}\rho _i^2(\mathrm{},b)\mathrm{\Delta }r_i(\mathrm{},b)$$
(5)
where $`\mathrm{}`$, $`b`$ are galactic longitude and latitude respectively. The flux is binned in angular windows of size $`\mathrm{\Omega }=1^{}\times 1^{}`$ and in the radial direction in fixed increments $`\mathrm{\Delta }r_i=`$ 1kpc.

The simulated dark matter halos are typically flattened oblate or prolate systems . We do not know a-priori in which axis the stellar disk would be located, therefore we show two all-sky maps using the same dark matter halo but viewed using two different locations for the “observer”: Figure 3(a) and Fig. 3(b) has the observer located on the short and long axes respectively. Both of these plots show the enhanced brightening towards the halo center, as well as some clumpy substructure in the halo itself. Note that both the central halo and the centers of the substructure halos are artificially “dimmed” in these plots due to the numerical resolution $`0.5`$ kpc, which sets a maximum density that can be resolved. The non-spherical shape of the halo is also clearly evident by inspecting the plots with different observer positions.

Recent estimates for the shape of the Milky Way’s halo (see, e.g.,  and references therein), suggest that it may be flattened with a short/long axis ratio of 0.5. An independent estimate from the orbit of the Sagittarius debris stars yields a nearly spherical dark matter halo . The simulated halo that we have chosen to analyse represents a typical prolate CDM halo with a short to long axis ratio of 0.4, and intermediate to long axis ratio of 0.5.

It is straightforward to estimate the effects of flattened dark matter halos by modifying Eq. ( 1) to accommodate triaxial shaped bodies. The simplest way to achieve this is to change from spherical coordinate $`r`$ to

$`\xi ^2=\frac{x^2+y^2}{b^2}+\frac{z^2}{c^2}`$

where $`b>c`$ for the oblate case, and $`b<c`$ for prolate and we leave $`z`$ as the axis of symmetry. A 2d visualization of these 3d shapes is illustrated in Fig. 4.

In Fig. 5 we plot spherical, oblate (2:1) and prolate (2:1) versions of the integral in Eq. ( 1) using the Moore et al. (1999) density profile. The observer is located on a plane parallel to the axis of symmetry, again at a distance $`R_o=8.5`$kpc from the center of the halo. The halo shape leads to little difference towards the Galactic center, but at the anti-center prolate halos can be 100 times brighter than oblate halos.

We can also compare the predicted angular flux with that measured directly from the N-body simulation. The annihilation flux is averaged in ten degree bins from the simulated dark matter halo, along a great circle from from the galactic center to its anti-center. This direct measurement of the flux is also plotted (as points) in Fig. 5. These data are particularly noisy due to the large numbers of substructure clumps in the simulation - the spike at $`\psi =125^{}`$ is due to a massive dark clump that happens to lie exactly along this chosen line of sight.

## III Substructure

### A Enhancement of global flux due to substructure

Cold dark matter substructure clumps have singular density profiles that will be a significant source of annihilation products. The velocities and spatial distribution of dark matter substructure is unbiased with respect to the smooth dark matter background . Therefore, to first order, substructure increases the global sky brightness in any given direction. However, the details depend on how much substructure survives within the solar radius and also on how far down the mass function substructure halos form and survive.

First we will estimate the annihilation flux from clumps of dark matter that are known to exist in the Galactic halo i.e. the dark matter halos that surround the Magellanic Clouds and dwarf spheroidal galaxies. In fact, high-Energy gamma-ray emission from the Large Magellanic Cloud (LMC) was detected with EGRET by Sreekumar et al.  in 1992 (although the origin of this emission was reported to be the interaction of cosmic rays with interstellar matter).

We estimate the average flux, $`\mathrm{\Phi }_{AV}`$, from the dark matter halos that surround some of the principal structures in the Local Group: The “Andromeda Galaxy”; M31 ($`v_{peak}=200\mathrm{k}\mathrm{m}\mathrm{s}^1`$ at a distance of 700 kpc), The Large and Small Magellanic Clouds ($`v_{peak}=70\mathrm{k}\mathrm{m}\mathrm{s}^1`$ and $`v_{peak}=40\mathrm{k}\mathrm{m}\mathrm{s}^1`$, respectively; both at a distance of 50 kpc), Draco ($`v_{peak}=10\mathrm{k}\mathrm{m}\mathrm{s}^1`$ at a distance of 50 kpc) and a small dark matter clump ($`v_{peak}=2\mathrm{k}\mathrm{m}\mathrm{s}^1`$ at a distance of 10 kpc). A sketch of the geometry is given in Fig. 6.

The total flux from a substructure halo at distance $`R_c`$ from the observer is

$$\mathrm{\Phi }_{TOT}(R_c)=\frac{K}{R_c^2}\rho ^2(r)r^2𝑑r.$$
(6)
by considering the central $`\mathrm{\Delta }\mathrm{\Omega }=1^{}\times 1^{}`$ patch over each clump, we define the maximum integration limit in 6 above and the average flux is then

$$\mathrm{\Phi }_{AV}=\frac{\mathrm{\Phi }_{TOT}}{\mathrm{\Delta }\mathrm{\Omega }}$$
(7)
(we set $`\mathrm{\Delta }\mathrm{\Omega }`$ in steradians, so we may compare directly with the smooth flux of Section II).

For the dark matter distribution within the substructure clumps we use the Moore et al. profile, which provides a good fit to the smallest, well-resolved substructure halos. The concentration of CDM halos is a function of mass  and for the density profile in Eq. 2 this can be written:

$$c_{moore}102(\frac{M_{\mathrm{𝑣𝑖𝑟}}}{1h^1M_{}})^{0.084}.$$
(8)
This defines the scale radius of each substructure clump: $`a_{\mathrm{M31}}=33.3\mathrm{kpc},a_{\mathrm{LMC}}=6.7\mathrm{kpc},a_{\mathrm{SMC}}=3.1\mathrm{kpc},a_{\mathrm{Draco}}=0.5\mathrm{kpc},a_{\mathrm{Tiny}}=0.05\mathrm{kpc}`$.

The integral in Eq. 6 diverges as $`r0`$ for the density profile that we are using, however, even the smallest substructure halos will have a maximum density set by the radius within which most of the neutralinos would have self annihilated. We therefore present results for the average flux from these clumps as a function of the minimum integration radius $`\mathrm{R}_{\mathrm{min}}/a`$ in Fig. 7, where $`a`$ is the scale radius as defined above.

For comparison, we plot the range of background emission at the Galactic anti-center as the shaded line in Fig. 7. The Tiny clump is only marginally visible above the background flux (depending on whether or not the Galactic halo is prolate or oblate) whereas most of the subhalos are easily visible. Also for comparison we have plotted the flux from the inner region of the galaxy which is the brightest of these sources.

Although the Galactic halo is expected to contain just a few clumps more massive that the Magellanic Clouds, there are many thousands of smaller mass objects. The mass function of substructure is a power law close to $`dn(m)/dmm^{1.9}`$ or in terms of circular velocity $`dn(v_c)/dv_cv_c^{3.8}`$). Above a circular velocity $`v_{peak}=10`$ km s-1 and $`1`$ km s-1 we expect the galactic halo to host roughly 1000 and $`5\times 10^5`$ substructure halos respectively. Future simulations should be able to measure how far down the mass function substructure halos can survive as well as to determine their central density profiles. (We note that the highest resolution simulation to date resolved the substructure within a dark matter mini-halo of mass $`10^7M_{}`$. The force resolution was 10 parsecs and the mass resolution was 10M allowing substructure with peak circular velocities as low as a few hundred meters per second to be resolved. The survival of substructure continues even down to this scale, where the slope of the power spectrum is close to -3.)

We calculate the total flux from substructure using Monte-Carlo techniques. First we generate a list of peak circular velocities and positions of $`5\times 10^5`$ substructure halos in the range of 1–70km s-1. (Distances are randomly selected using the Moore et al. density profile and peak circular velocities are randomly assigned from a power law distribution scaling as $`v^{3.8}`$). For each lump, we estimate its total flux as in the previous cases, integrating Eq. 6 with a density profile scaled according to Eq. 8 for the concentration.

In the absence of further constraints on the possible value for $`\mathrm{R}_{\mathrm{min}}`$, we use the same criteria as before and choose it to be a fixed fraction of the virial radius, $`\mathrm{R}_{\mathrm{min}}10^{12}r_{200}`$. This corresponds to a mean density of $`10^{22}M_{}`$ kpc-3 for the galactic halo. The results are not too sensitive to the value of the minimum integration radius as is apparent from inspection of Fig. 7. The total flux is then averaged over the entire sky and we repeat this process in order to estimate the variance. The cumulative distribution of flux above a given substructure peak circular velocity ($`\mathrm{\Sigma }\mathrm{\Phi }_{TOT}`$) is plotted for ten of these random halo realisations in Fig. 8.

It is evident from this plot that the effects of including the entire mass spectrum of substructure is quite dramatic and boosts the expected flux from the smooth halo by several orders of magnitude. However, most of the flux arises from the subhalos with circular velocities larger than 10 km s-1. Extrapolating to very small halos would not change the total flux by a large factor.

To quantify the brightening of the background due to substructure, we have calculate the average flux due to all clumps with $`v_{peak}>1`$ km s-1 within a spherical halo. The point plotted in Fig. 7 represents this contribution to the flux, where the error-bar is the 1$`\sigma `$ variation among the different Monte-Carlo models. From this plot we see that the flux due to substructure is over two orders of magnitude brighter than the smooth background from a spherical halo. We note that one needs to observe a fairly large fraction of the sky ($`>100`$ square degrees) to ensure a significant number of clumps lie in the field of view. (Also note that the variance at high peak circular velocities is due to the proximity of the largest few dark matter substructures, however, the mean total flux converges to similar values for each Monte-Carlo model.)

### B The flux due to substructure in prolate and oblate halos

Not only is the mean flux at a given position on the sky dominated by substructure halos, the spatial distribution of flux across the sky will be determined by the substructure. The convergence study by Ghigna et al. (Ref. ) shows that substructure halos trace the global mass distribution of the halo. Therefore, we can use the N-body simulations to generate Monte-Carlo distributions of substructure halos and construct all-sky maps of the expected gamma-ray flux. We take a random particle from the simulation and assign a circular velocity from a distribution $`dn(v_c)/dv_cv_c^{3.8}`$. For each sub-halo we calculate its total annihilation flux and then repeat the process until we have 500,000 halos above a circular velocity of 1$`\mathrm{kms}^1`$.

Figure 9(a) and Fig. 9(b) show the resulting sky distribution of flux from sub-halos binned in one degree bins where the observer has been placed in the short and long axis of the simulation respectively. Large substructure halos, such as the Magellanic Clouds in our own halo, will contain its own gravitationally bound sub-halos which leads to clustering of gamma-ray emission in the all sky maps.

Future observations may only be able to make deep strip maps therefore in Fig. 10 we have binned the flux along lines of constant galactic $`l`$ and $`b`$ with the observer placed in the short and long axis of the global density distribution. From these plots we can see that the emission from substructure peaks at the galactic center, as one would expect, this effect is not that different for spherical halos than for prolate or oblate halos.

Since the substructure traces the global mass distribution, a prolate halo would also have a prolate distribution of satellites. Therefore we can study the variation of flux within smooth prolate or oblate halos to examine how the background flux from substructure can be used to quantify the halo shape. We calculate the observed flux as a function $`l`$ and $`b`$ for spherical, prolate and oblate flattened 2:1 and 3:1 geometries. In each case, the density profile is taken from Eq. 2 and again the observer is placed in either the short (Figure 11(a)) or long (Figure 11(b)) axis of symmetry.

These plots show how the distribution of flux on the sky can vary significantly depending on the shape of the density distribution and on where the observer is situated within the halo.

### C The distribution of point sources

Individual substructures may be observed and quantified if the resolution of the telescope is sufficient. However, all of the past and present observations would only detect substructure as unresolved point sources. The distribution of their fluxes (and spatial distribution on the sky) may be used to rule out alternative origins, such as extra-galactic sources. In Fig. 12 we plot the cumulative distribution of point sources above a given flux within one degree square bins. The two curves consider substructure with peak circular velocities larger than 10 km s-1 and 1 km s-1. The number density of the brightest sources in the sky scales as $`NF^{0.7}`$.

Higher resolution simulations are vital to quantify how much substructure survives within the galactic halo, how it is spatially distributed and to quantify the internal structure of surviving substructure. However, Fig. 12 gives an idea of what to expect if an all sky survey is carried out that is capable of detecting the brightest substructure halos.

## IV Conclusions

Numerical simulations that follow the growth of structure within a universe dominated by neutralinos (cold dark matter) have achieved a resolution that allows their global structure and internal structure to be quantified. The density profiles, shapes of dark matter halos, abundance and properties of dark matter substructure, all play an important role in determining the absolute surface brightness of observable products from dark matter annihilation.

We have used the results from the highest resolution simulations ever performed of CDM halos to examine the expected all-sky distribution of gamma-rays from neutralino annihilation. Substructure can boost the expected flux significantly over that originating from a smooth dark matter halo. Thus, gamma-ray observations, such as EGRET data, may already have the potential of constraining a large part of the parameter range of the neutralino cross-sections. The distinguishing shapes of CDM halos and the unique spatial and flux distribution of point sources from substructure within the Galactic halo should allow a unique identification of observational data with dark matter.

## Acknowledgments

The authors would like to thank Prof. Arnold Wolfendale for numerous discussions and suggestions that have improved the quality of this work. Carlos Calcáneo–Roldán continues his research thanks to the generous support from the People of México, through a grant by CONACyT. Ben Moore is a Royal Society research fellow. Computations were carried out as part of the Virgo consortium.

Figure 1

Figure 2

Figure 3(a)

Figure 3(b)

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9(a)

Figure 9(b)

Figure 10

Figure 11

Figure 12


format: text

text:
# Phase Transitions of Hard Disks in External Periodic Potentials: A Monte Carlo Study

## I Introduction

The liquid -solid transition in systems of particles under the influence of external modulating potentials has recently attracted a fair amount of attention from experiments , theory  and computer simulations . This is partly due to the fact that well controlled, clean experiments can be performed using colloidal particles  confined between glass plates (producing essentially a two -dimensional system) and subjected to a spatially periodic electromagnetic field generated by interfering two, or more, crossed laser beams. One of the more surprising results of these studies, where a commensurate, one dimensional, modulating potential is imposed, is the fact that there exist regions in the phase diagram over which one observes re -entrant  freezing/melting behaviour. As a function of the laser field intensity the system first freezes from a modulated liquid to a two dimensional triangular solid – a further increase of the intensity confines the particles strongly within the troughs of the external potential, making the system quasi -one -dimensional which increases fluctuations and leads to re -melting.

Our present understanding of this curious phenomenon has come from early mean -field density functional  and more recent dislocation unbinding  calculations. The mean field theories neglect fluctuations and therefore cannot explain re -entrant behaviour. The order of the transition is predicted to be first order for small laser field intensities, though for certain combinations of external potentials (which includes the specific geometry studied in the experiments and in this paper) the transition may become second order after going through a tricritical point. In general, though mean field theories are applicable in any dimension, the results are expected to be accurate only for higher dimensions and long ranged potentials. The validity of the predictions of such theories for the system under consideration is, therefore, in doubt.

A more recent theory  extends the dislocation unbinding mechanism for two -dimensional melting  for systems under external potentials. For a two -dimensional triangular solid subjected to an external one -dimensional modulating potential, the only dislocations involved are those which have their Burger’s vectors parallel to the troughs of the potential. The system, therefore, maps onto an anisotropic, scalar Coulomb gas (or XY model)  in contrast to a vector Coulomb gas  for the pure $`2d`$ melting problem. Once bound dislocation pairs are integrated out, the melting temperature is obtained as a function of the renormalized or “effective” elastic constants which depend on external parameters like the strength of the potential, temperature and/or density. Though explicit calculations are possible only near the two extreme limits of zero and infinite field intensities one can argue effectively that a re -entrant melting transition is expected on general grounds quite independent of the detailed nature of the interaction potential for any two -dimensional system subject to such external potentials. The actual extent of this region could, of course, vary from system to system. In addition, these authors predict that the auto -correlation function of the Fourier components of the density (the Debye -Waller correlation function) decays algebraically in the solid phase with a universal exponent which depends only on the geometry and the magnitude of the reciprocal lattice vector.

Computer simulation results in this field have so far been inconclusive. Early simulations  involving colloidal particles interacting via the Derjaguin, Landau, Verwey and Overbeek (DLVO) potential  found a large re -entrant region in apparent agreement with later experiments. On closer scrutiny, though, quantitative agreement between simulation and experiments on the same system (but with slightly different parameters) appears to be poor . Subsequent simulations  have questioned the findings of the earlier computation and the calculated phase diagram does not show a significant re -entrant liquid phase.

Motivated, in part, by this controversy, we have investigated the freezing/melting behaviour of an unrelated system subjected to similar modulating external potentials. In this paper we have computed the phase behaviour of a two dimensional hard disk system in an external potential. The pure hard disk system is rather well studied  by now and the nature of the melting transition in the absence of external potentials reasonably well explored. Also, there exist colloidal systems with hard interactions  so that, at least in principle, actual experiments using this system are possible. Finally, a hard disk simulation is relatively cheap to implement and one can make detailed studies of large systems without straining computational resources. The main outcome of our calculations, the phase diagram, is shown in Fig. 1. We have shown results from our simulation of a system of $`N=1024`$ hard disks (diameter $`\sigma `$) of density $`0.86<\rho ^{}(=\rho \sigma ^2)<0.91`$ and the amplitude of the external potential $`0<V_0^{}(=\beta V_0)<1000`$. Within our range of densities, one has a clear signature of a re -entrant liquid phase showing that this phenomenon is indeed a general one as indicated in Ref. .

The rest of our paper is organized as follows. In Section II we specify the model and the simulation method including details of the finite size analysis used. In Section III we present our results for the order parameter and its cumulants with a discussion on finite size and hysteresis effects. We also present results for the specific heat, order parameter susceptiblity and correlation functions which further illustrates the nature of the phase transitions in this system. In Section IV we discuss our work in relation to the existing literature on this subject, summarize and conclude.

## II Model and method

### A The Model

We study a system of N hard disks of diameter $`\sigma `$ in a two dimensional box of size $`S_x\times S_y`$ ($`S_x/S_y=\sqrt{3}/2`$) interacting with the pair potential $`\varphi (r_{ij})`$ between particles $`i`$ and $`j`$ with distance $`r_{ij}`$,

$$\varphi (r_{ij})=\{\begin{array}{cc}\mathrm{}\hfill & r_{ij}\sigma \hfill \\ 0\hfill & r_{ij}>\sigma \hfill \end{array}$$
(1)
In addition a particle with coordinates $`(x,y)`$ is exposed to an external periodic potential of the form:

$$V(x,y)=V_0\mathrm{sin}(2\pi x/d_0)$$
(2)
The constant $`d_0`$ in Eq.(2) is chosen such that, for a density $`\rho =N/S_xS_y`$, the modulation is commensurate to a triangular lattice of hard disks with nearest neighbor distance $`a_s`$: $`d_0=a_s\sqrt{3}/2`$ (see Fig. 2). The only parameters which define our system are the reduced density $`\rho \sigma ^2=\rho ^{}`$ and the reduced potential strength $`V_0/k_BT=V_0^{}`$, where $`k_B`$ is the Boltzmann constant and $`T`$ is the temperature.

### B The Method

#### 1 Numerical Details

We perform NVT - Monte Carlo (MC) simulations  for the system with interactions given by Eqs. (1) and (2) for various values of $`V_0^{}`$ and $`\rho ^{}`$. Averages $`<>`$ of observables have been obtained with the canonical measure. In order to obtain thermodynamic quantities for a range of system sizes, we have analyzed various quantities within subsystems as shown in Fig. 3. We have used $`<>_L`$ to denote averages in subsystems. The subsystems are of size $`L_x\times L_y`$ where $`L_x`$ and $`L_y`$ are chosen as $`L_y=La_s`$ and $`L_x=L_y\sqrt{3}/2`$ consistent with the geometry of the triangular lattice.

Most of the simulations described below have been done for $`N=1024`$ particles unless otherwise indicated. Phase transitions have been studied in most cases by starting in the ordered solid and reducing $`\rho ^{}`$ for fixed $`V_0^{}`$. Runs where the density $`\rho ^{}`$ is increased were also performed in a few cases.

A typical simulation run with 4 $`\times `$ $`10^7`$ Monte Carlo steps (MCS) (including 1.5 $`\times 10^7`$ MCS for relaxation) took about 50 CPU hours on a PII/500 MHz PC. At high values of $`V_0^{}`$ in addition to ordinary (local) MC moves we also used ‘through-moves’, by which particle placements in neighboring troughs are tried. Besides producing faster equilibration, including such moves ensures that the formation of dislocations for large $`V_0^{}`$ and $`\rho ^{}>\sqrt{3}/2`$ ($`d_0<\sigma `$) are not artificially hindered since particles can bypass each other — this is impossible with purely local MC moves.

To guarantee good equilibration and averaging, we simulated only systems up to N=1600 particles in the region of the phase boundary. Systems of N=4096 and N=16384 were used only once in Fig. 11, where the interesting region is clearly in the liquid phase and equilibration is much easier.

#### 2 Order parameter

The nature of the fluid-solid phase transition in two dimensions has been a topic of controversy throughout the last forty years . It is well known that true long range positional order is absent in the infinitely large system due to low energy long wavelength excitations so that translational correlations decay algebraically. According to the dislocation unbinding mechanism  the two dimensional solid (with quasi long ranged positional order) first melts into a “hexatic” phase with no positional order but with quasi long ranged orientational order signified by a non -zero bond orientational order parameter $`\psi _6=\mathrm{exp}(i6\theta )`$ where $`\theta `$ is the angle of a bond and the sum is over all distinct bonds. A liquid, with no bond orientational order either ($`\psi _6=0`$) is produced by a second Kosterlitz -Thouless (KT)  transition from the hexatic.

In an external periodic field given by Eq.(2), however, the bond orientational order parameter is non -zero even in the fluid phase . This is because for $`V_0^{}0`$ we have now a “modulated” liquid, in which local hexagons consisting of the six nearest neighbors of a particle are automatically oriented by the external field. Thus $`<\psi _6>`$ is non -zero both in the (modulated) liquid and the crystalline phase and it cannot be used to study phase transitions in this system. The order parameters corresponding to a solid phase are the Fourier components of the (non -uniform) density $`\rho (\stackrel{}{r})`$ calculated at the reciprocal lattice points $`\{\stackrel{}{G}\}`$. This (infinite) set of numbers are all zero (for $`\stackrel{}{G}0`$ ) in an uniform liquid phase and non -zero in a solid. We restrict ourselves to the star consisting of the six smallest reciprocal lattice vectors of the two -dimensional triangular lattice. In the modulated liquid phase that is relevant to our system, the Fourier components corresponding to two out of these six vectors, viz. those in the direction perpendicular to the troughs of the external potential, are non -zero . The other four components of this set consisting of those in the direction $`\stackrel{}{G}_1`$ (as defined for the ideal crystal in Fig. 2), and those equivalent to it by symmetry, are zero in the (modulated) liquid and non -zero in the solid (if there is true long range order). We therefore use the following order parameter:

$$\psi _{G_1}=\left|\underset{j=1}{\overset{N}{}}\mathrm{exp}(i\stackrel{}{G}_1\stackrel{}{r}_j)\right|$$

where $`\stackrel{}{r}_j`$ is the position vector of the $`j^{th}`$ particle. Note that though the order parameter $`<\psi _{G_1}>`$ decays to zero with increasing system size in the 2 -d solid — quasi long ranged order — this decay, being weak, does not hinder us from distinguishing, in a finite system, a modulated liquid from the solid phase with positional order in the $`\stackrel{}{G}_1`$ direction.

#### 3 Cumulants

We have determined phase transition points by the order parameter cumulant intersection method . The fourth order cumulant $`U_L`$ of the order parameter distribution is given by:

$$U_L(V_0^{},\rho ^{})=1\frac{<\psi _{G_1}^4>_L}{3<\psi _{G_1}^2>_L^2}$$
(3)
In case of a continuous transition close to the transition point the cumulant is only a function of the ratio of the system size $`La_s`$ and the correlation length $`\xi `$: $`U_L(La_s/\xi )`$. Since $`\xi `$ diverges at the critical point the cumulants for different system sizes intersect in one point: $`U_{L_1}(0)=U_{L_2}(0)=U^{}`$. Even for first order transitions these cumulants intersect  though the value $`U^{}`$ of $`U_L`$ at the intersection is not universal any more. The intersection point can, therefore, be taken as the phase boundary regardless of the order of the transition. This is useful since the order of the melting transition in $`2d`$ either in the absence  or with  external potentials is not unequivocally settled.

In order to map the phase diagram we systematically vary the system parameters $`V_0^{}`$ and $`\rho ^{}`$ to detect order parameter cumulant intersection points which are then identified with the phase boundary. It should be noted that though the order parameter (defined for long range positional order) vanishes  with increasing system size in the crystalline phase, its cumulants are well defined and can be used to determine phase boundaries. For large $`L`$ the cumulants approach the value 2/3 in the solid phase and 1/3 in the liquid  so that they are guaranteed to intersect ! For very large $`V_0^{}`$ we do not find an unique point of intersection for $`U_L`$, instead the cumulants for various values of $`L`$ collapse onto a single curve. In this case the onset of the collapse is taken as the “intersection” density. It is curious to note that this behaviour is, in fact, typical of the anisotropic XY model . In this case although the order parameter cumulants have an intersection point, the value of the cumulant at the intersection differs for various anisotropies and drifts towards a limiting value at zero anisotropy. The intersection “point” therefore changes to a “line” of intersections for different system sizes and for small anisotropies. In our system, for the large $`V_0^{}`$ we see similar behaviour.

## III Results and Discussion

### A Order parameter and cumulants

In Fig. 4 we present data for the average order parameter $`<\psi _{G_1}>`$ and its cumulants as functions of the density for $`V_0^{}=0.05`$ and $`V_0^{}=0.5`$ calculated within various subsystems. In both cases $`<\psi _{G_1}>_L`$ and $`U_L`$ increase with $`\rho ^{}`$ with a sharpening of the structures for increasing $`L`$. As discussed above, we observe that for any density increasing subsystem size $`L`$ depresses $`<\psi _{G_1}>_L`$. The cumulant $`U_L`$, on the other hand, approaches limiting values (2/3 for solid and 1/3 for liquid). The values of the cumulants are higher for larger $`L`$ in the ordered (solid) phase and vice versa in the disordered (modulated liquid) phase thus resulting in an intersection point — the transition density.

In Fig. 5 we compare the density dependence of the average order parameter $`<\psi _{G_1}>`$ (calculated over the entire system) for different $`V_0^{}`$ and for the same system size ($`N=1024`$). With increasing $`V_0^{}`$ the turning point in $`<\psi _{G_1}>(\rho ^{})`$ is shifted to lower densities and then for even larger $`V_0^{}`$ values to higher densities. This indicates, already, that the system prefers having smaller transition densities for intermediate values of $`V_0^{}`$ compared to smaller and higher $`V_0^{}`$ values — i.e. we have a re -entrant transition.

In Fig. 6 we show a systematic study of $`<\psi _{G_1}>_L`$ and $`U_L`$ as a function of $`V_0^{}`$ at the density $`\rho ^{}=0.89`$ for different $`L`$\- values. Maxima in the $`<\psi _{G_1}>_L`$\- and $`U_L`$\- curves are found near $`V_0^{}=2`$. Again we note that the $`<\psi _{G_1}>_L`$\- values decrease with increasing $`L`$ (see Fig. 6(a)). The cumulants $`U_L`$, on the other hand increases with $`L`$ for intermediate values of $`V_0^{}`$ (the ordered, solid phase) and decreases with $`L`$ for either large or small $`V_0^{}`$ (the disordered, liquid phase) resulting in intersection points indicating two consecutive phase transitions (see Fig. 6(b)).

If $`V_0^{}`$ is increased to $`20`$ the value of the cumulant at intersection $`U^{}`$ is shifted upwards, see Fig. 7(a). For very high $`V_0^{}`$\- values the cumulant curves for different $`L`$ merge on the high density side, see Fig. 7(b) (see discussion in Section II B 3). In Fig. 8 the cumulant intersection values are shown as a function of $`V_0^{}`$, where for large $`V_0^{}`$-values the value at the onset of the merging is shown. We observe that $`U^{}`$ is not an universal number but, nevertheless, goes to a limiting value for large $`V_0^{}`$ .

In Fig. 9 we show $`<\psi _{G_1}>`$ as a function of the density with $`V_0^{}=0.5`$ for different $`N`$-values. The general features of $`<\psi _{G_1}>`$ as discussed above is retained though there is a shift of the turning point to slightly higher densities with increasing $`N`$. The effect on the phase diagram is discussed in Section III C.

### B Susceptibility, specific heat, finite size effects and hysteresis

In addition to $`<\psi _{G_1}>_L`$ and $`U_L`$ we have computed the order parameter susceptibility $`\chi _{G_1}`$ and the specific heat for different system- and subsystem-sizes.

The order parameter susceptibility $`\chi _{G_1}`$ is defined as :

$$k_BT\chi _{G_1}=L^2\left[\left(\psi _{G_1}\right)^2\psi _{G_1}^2\right]$$
(4)
In Fig. 10(a) we show $`\chi _{G_1}`$ as a function of $`\rho ^{}`$ at $`V_0^{}=0.05`$ for different $`L`$\- values. The increase of $`\chi _{G_1}`$ with increasing $`L`$ signals the presence of a phase transition in the density range where the transition has been found by cumulant intersection techniques ($`\rho _t^{}0.896`$). In Fig. 10(b) $`\chi _{G_1}`$ is shown for the same system size ($`N=1024`$) and various $`V_0^{}`$-values. We note that the density of the $`\chi _{G_1}`$\- maxima are smallest for the intermediate value of $`V_0^{}`$ which again show that for these $`V_0^{}`$\- values the transition density is lowest. Compared to the cumulant intersection values, $`\chi _{G_1}`$ maxima are located at slightly smaller densities (see also Section III C) which may be due to finite size effects, which often show the feature that phase transition points in finite systems are shifted to slightly different values depending on the observable under investigation. In particular one expects (and we get) a shift towards parameter values in the disordered region (here a liquid, i.e. low densities) for the order parameter and the susceptibility as compared to the cumulant intersection parameters.

We have also calculated the specific heat $`C_V`$ (Fig. 11) as a function of the density for $`V_0^{}=0.2`$ with $`N=4096`$ and $`N=16384`$. For a second order transition, the maximum of the specific heat scales with the system size as $`C_V^{max}L^{\alpha /\nu }`$ where $`\alpha `$ and $`\nu `$ are critical exponents. For a first order transition, on the other hand $`C_V^{max}L^d`$ where $`d`$ is the dimensionality ($`=2`$ in our case). We, however, do not see any of this behaviour. In contrast, the specific heat is relatively featureless. Although it shows a peak, surprisingly, the height of this peak is almost insensitive to system size. This is a strong indication that the phase transition we observe is unconventional and is KT like . Further, as expected for such transitions, the maximum does not lie at the density where the cumulants intersect and it would be incorrect to identify specific heat maxima with the phase boundary (see discussion in Section IV).

In order to study the effect of the path taken through the parameter space on the location of the cumulant intersection densities, we compared $`U_L`$ as a function of the density for $`V_0^{}=0.5`$ as obtained from two runs. In Fig. 4(d) we have already presented the data for a run where the density was decreased systematically. In Fig. 12 we present data for a second run where the density was increased instead. We find negligibly small hysteresis effects on the cumulants as well as on the value of intersection density. This shows that the transition points are not affected much by the path through the parameter space.

### C The Phase Diagram

We have obtained the phase diagram of the system for $`.86<\rho ^{}<.91`$ and $`0<V_0^{}<1000`$. For each density and $`V_0^{}`$\- value we computed cumulants $`U_L`$ for a range of subsystem sizes $`L`$ and located intersection points which we identify with the phase boundary. The resulting locus of the transition points is shown in the phase diagram, see Fig. 1. At very small values of $`V_0^{}`$ we find good agreement of our transition densities with the melting densities ($`\rho _m0.91`$) known from literature  on the pure hard disk solid ($`V_0^{}=0`$). The values of the transition density initially drop and subsequently rise as $`V_0^{}`$ increases. The minimum transition density is found for $`V_0^{}12`$. These transition points separate a high density solid from a low density modulated liquid. Thus, at a properly chosen density, we observe an initial freezing transition followed by a re -entrant melting at a higher $`V_0^{}`$\- value. Such an effect had been found earlier in experiments on colloidal systems in an external laser field .

In order to quantify finite size effects on the phase diagram, we have computed the transition points for different system sizes. The resulting phase diagrams are shown in Fig. 13. We note that due to residual finite size effects with increasing system size all transition points are slightly shifted to higher densities, the structure of the phase diagram with a pronounced minimum at intermediate values of $`V_0^{}`$ is not affected by this shift.

### D Correlation functions

The Debye-Waller- correlation function is defined as follows:

$$C_{\stackrel{}{G}_1}(\stackrel{}{R})=<e^{i\stackrel{}{G}_1(\stackrel{}{u}(\stackrel{}{R})\stackrel{}{u}(0))}>$$

where $`\stackrel{}{R}`$ points to the elementary cell of the ideal lattice, and $`\stackrel{}{u}(\stackrel{}{R})`$ is the deviation of the actual particle position from the ideal lattice: $`\stackrel{}{r}=\stackrel{}{R}+\stackrel{}{u}(\stackrel{}{R})`$. In this case we have chosen the direction of $`\stackrel{}{R}`$ to lie along the $`y`$ axis (i.e. along the troughs of the potential).

We have also computed the spatial correlation function $`g(y)`$ which is the pair correlation function in the $`y`$-direction. We compute it in the following way: for a particle i, $`g(y)dy`$ number of particles j for which: $`|y_iy_j|[y,y+dy]`$ and $`|x_ix_j|<d_0/2`$, normalized so that $`g(y)1`$ as $`y\mathrm{}`$.

These correlation functions are plotted in Fig. 14 as functions of $`y`$. The Debye- Waller correlation function $`C_{G_1}(y)`$ and the correlation function $`g(y)`$ along the potential valley are compared in Fig. 14(a) at a density just below the transition. We see that the decay of the maxima of $`g(y)`$ as function of $`y`$ is similar to the decay in $`C_{G_1}(y)`$. The decay of $`C_{G_1}(y)`$ is analyzed in more detail in Fig. 14(b) for parameter values in the liquid and in the solid phase. In the liquid phase the decay is exponential while in the solid region it is algebraic: $`C_{G_1}(y)y^{\eta _{G_1}}`$. Taking the data-points in the crystal which are closest to the phase boundary for each $`V_0^{}`$, we get $`\eta _{G_1}`$ in the range of $`0.20\mathrm{}0.27`$. The exponent $`\eta _{G_1}`$ is predicted  to be universal and equal to $`1/4`$ for our geometry, so this value is consistent with our numerical results.

### E Scaling behavior

We next try to determine the order of the phase transitions encountered in this system for various values of $`V_0^{}`$. In order to investigate this issue we studied the scaling behavior of the order parameter, susceptibility and the order parameter cumulant near the phase boundary for a small ($`0.5`$) and a large ($`1000`$) $`V_0`$. From finite size scaling theory (for an overview see Ref. ) we expect these quantities to scale as  :

$$<\psi _{G_1}>_LL^bf(L/\xi )$$
(5)
$$\chi _Lk_BTL^cg(L/\xi )$$
(6)
$$U_Lh(L/\xi )$$
(7)
Here $`b=\beta /\nu `$, $`c=\gamma /\nu `$ (for critical scaling) and $`f`$, $`g`$, $`h`$ are scaling functions. The correlation length $`\xi `$ diverges as $`\xi (1\rho /\rho _c)^\nu `$ for an ordinary critical point, while for a KT- transition we have an essential singularity and $`\xi \mathrm{exp}(a(1\rho /\rho _c)^{\stackrel{~}{\nu }})`$.

According to general arguments given in Ref., we expect that for a finite lattice, the identification of the properties of our system with those of the anisotropic XY model should improve with increasing $`V_0^{}`$. Indeed, for large $`V_0^{}`$, scaling according to the KT- theory seems to be supported by our data. In Fig. (15) we have plotted the left hand sides of Eqs.(7) and (6) versus $`L/\xi `$ for $`V_0^{}=1000`$, where data points for $`0.86\rho ^{}0.898`$ have been considered and $`\rho _c^{}=0.902`$, obtained by cumulant intersection. In order not to introduce an unwarranted bias, we have separately considered (a) ordinary critical scaling and (b) a KT scaling forms and adjusted the values of the parameters $`b`$, $`c`$ and $`\nu `$ till we obtained collapse of our data onto a single curve determined by a least square estimator. Though good collapse of our data is observed both in (a) and (b), the numerical values for $`\stackrel{~}{\nu }`$, $`\eta =2b`$ and $`c=2\eta `$ for KT scaling ($`b0.138`$, $`c1.70`$, $`\stackrel{~}{\nu }0.44`$) are close to the predicted values  ($`b=\eta /2=1/8`$, $`c=1.75`$, $`\stackrel{~}{\nu }=0.5`$).

The situation is less straight -forward for $`V_0^{}=0.5`$. The critical parameters were obtained in this case for densities $`0.85\rho ^{}0.876`$, with $`\rho _c^{}=0.878`$. In Fig. 16 (a) the data collapse looks slightly better than in (b), such that relying on this data alone one may conclude that KT- scaling in this region of the phase diagram seems less likely. It must be kept in mind though that for small values of $`V_0^{}`$ in a finite system, the analysis of the data would be complicated by crossover effects. Strictly for $`V_0=0`$ we do not have a correspondence with the anisotropic XY model but rather with a vector Coulomb gas with a different set of exponents. Our results for the numerical values of the parameters are summarized in Table I.

In summary, from the scaling analysis in Figs. (16) and (15) a KT- scenario at least for large $`V_0^{}`$\- values seems likely. This is supported by the behavior of the cumulants as well (see Sect. III A). A more precise classification of the phase transitions with the present data and system sizes is not easy. This topic is left for future work, in particular we plan to compute the elastic properties of the system by a method recently developed for the hard- disk system  and to test the KT- predictions .

## IV Summary and conclusion

In summary, we have calculated the phase diagram of a two dimensional system of hard disks in an external sinusoidal potential. We found freezing followed by re -entrant melting transitions over a significant region of the phase diagram in tune with previous experiments on colloids  and with the expectations of a recent dislocation unbinding theory . One of the main features of our calculation is the method used to locate phase boundaries. In contrast to earlier simulations  which used either the jump of the order parameter, or specific heat maxima to locate the phase transition, we have used the more reliable cumulant intersection method. It must be noted that the specific heat in this system does not show a strong peak at the phase transition density so that its use may lead to confusing results. This, in our opinion, may be the reason for part of the controversy in this field. It is possible that earlier simulations which used smaller systems and no systematic finite size analysis may have overlooked this feature of $`C_V`$ which becomes apparent only in computations involving large system sizes. We have shown that finite size scaling of the order parameter cumulants as obtained from sub-system or sub-block analysis, on the other hand, yields an accurate phase diagram.

What is the order of the phase transition? We know that  for the pure hard disk system in two -dimensions this question is quite difficult to answer and our present understanding  is that hard disks in two dimensions show a KTHNY transition. This transition, however, lies very close (in an appropriately expanded parameter space) to a first order boundary so that crossover effects may be significant. The present system has an imposed external periodic potential which stabilizes the hexatic phase  and an (anisotropic) KT transition  is expected. Our results show several features which suggest that this is, perhaps, what we have. Though we have discussed these observations in the rest of the paper, we list them below for clarity:

* The behaviour of the value of the cumulants at intersection $`U^{}`$ is similar to an earlier work  on the anisotropic XY system which shows a KT transition.
* The specific heat is relatively featureless and does not scale with system size in a fashion expected of a true first order or conventional continuous transition.
* The decay of the correlation functions is similar to what is predicted  for an anisotropic scalar Coulomb gas.
* For large $`V_0^{}`$\- values the scaling of the order parameter, the susceptibility and the cumulant may be described by the KT- theory.

Of course, in order to resolve this issue unambiguously yet larger simulations are required. Also, we need to compute elastic properties  of this system in order to compare directly with the results of Ref. . Work along these lines is in progress.

Before we end we would like to point out that after completion of this work and prior to the submission of this manuscript we received a preprint  where the same system as ours has been studied using simulations. The phase diagram obtained by these authors is similar to ours (thus confirming and corroborating our results), though there exist some quantitative differences. These differences may be attributed to the absence of systematic finite size scaling in the latter work.

## V Acknowledgment

We are grateful for many illuminating discussions with C. Bechinger and K. Binder. One of us (S.S.) thanks the Alexander von Humboldt Foundation for a Fellowship. Support by the SFB 513 and granting of computer time from the NIC and the HLRS is gratefully acknowledged.

| $`V_0^{}`$ | $`b`$ | $`c`$ | $`\nu `$ | $`b`$ | $`c`$ | $`\stackrel{~}{\nu }`$ | $`a`$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 1000 | 0.13(1) | 1.68(5) | 2.25(25) | 0.138(8) | 1.70(5) | 0.44(3) | 1.05(25) |
| 0.5 | 0.152(5) | 1.65(6) | 1.06(13) | 0.170(12) | 1.83(4) | 0.38(10) | 1.0(2) |

Table I Parameters in the scaling plots (Figs. (16) and (15)) for $`V_0^{}=0.5`$ and $`V_0^{}=1000`$. The first three parameter columns are for critical scaling, the last four for KT- scaling.


text:
# Total 4He Photoabsorption Cross Section Revisited: Correlated HH versus Effective Interaction HH

## Abstract

Two conceptually different hyperspherical harmonics expansions are used for the calculation of the total 4He photoabsorption cross section. Besides the well known method of CHH the recently introduced effective interaction approach for the hyperspherical formalism is applied. Semi-realistic NN potentials are employed and final state interaction is fully taken into account via the Lorentz integral transform method. The results show that the effective interaction leads to a very good convergence, while the correlation method exhibits a less rapid convergence in the giant dipole resonance region. The rather strong discrepancy with the experimental photodisintegration cross sections is confirmed by the present calculations.

PACS numbers: 21.45.+v, 25.20.Dc, 31.15.Ja, 21.30.Fe

The 4He photodisintegration in the giant dipole resonance region is a particularly interesting reaction. An understanding of this process in terms of a microscopic calculation is not only a challenge in few-body physics, but could also lead to a deeper insight in the physics of more complex nuclei. Because of the obvious difficulties in describing correctly the four-nucleon dynamics, calculations with realistic interactions are still lacking. In the most advanced calculation  semi-realistic NN potential models were employed and the complicated four-nucleon final state interaction was treated exactly applying the Lorentz integral transform (LIT) method . The semi-realistic potential models lead to rather realistic results for the total three-nucleon photoabsorption cross section  and the 4He inverse energy weighted sum rule (see below). Hence one might expect a rather realistic description of the 4He giant dipole resonance cross section, which is situated close to the breakup threshold. However, the obtained cross sections show a considerably more pronounced giant dipole resonance than seen in the experimental results published in the last two decades (see ). On the other hand the experimental situation is not yet completely settled. Older photoabsorption data (see discussion in Ref. ) and a more recent determination of the photoabsorption cross section via photon scattering  show a stronger giant dipole peak. A new round of experiments presently carried out at Lund will hopefully help to clarify this unsatisfying situation.

Besides clarification on the experimental side it is also necessary to check the obtained theoretical result. In a recent calculation of the 4He photoabsorption with the effective interaction hyperspherical harmonics (EIHH) method  small deviations from the above mentioned calculation of Ref.  in a correlated hyperspherical harmonics (CHH) approach were found. The slightly different results are most probably due to not fully convergent HH expansions. Therefore it is the aim of the present work to study the convergence in both cases in greater detail extending the calculations to higher order terms. On the one hand it will allow establishing the correct 4He total photoabsorption cross section with semi-realistic NN potential models. On the other hand it will show which is the most efficient HH approach. This is very important in view of calculations of the 4He photoabsorption with realistic forces.

The rate of convergence of an HH expansion is generally rather slow in nuclear physics problems. In particular the short range repulsion of the NN interaction leads to high momentum components in the nuclear wave function which can be parametrized only by including a very large number of HH basis functions $`H_n`$. The convergence can be improved introducing proper short range two-body correlation functions $`f(r_{ij})`$,

$$H_n(\rho ,\mathrm{\Omega })\underset{i,j}{}f(r_{ij})H_n(\rho ,\mathrm{\Omega }),$$
(1)
where $`\rho `$ and $`\mathrm{\Omega }`$ denote hyperradius and hyperangle, while $`r_{ij}`$ is the relative distance of particles $`i`$ and $`j`$. The function $`f(r)`$ can be obtained from the solution of the two-body Schrödinger equation, since at short distances the role of other particles is rather unimportant. Though such correlations lead to a considerable improvement  one still needs in general a rather large number of HH terms in order to reach convergence. The convergence can be improved considerably if one introduces longer range correlations . However, these two-body long range correlations are less under control, because correlations among more particles become more important. In addition, different from short range correlations, they change considerably from ground to continuum states. In this respect long range ground state correlations would not be appropriate for the description of LIT states $`\stackrel{~}{\mathrm{\Psi }}`$ in electromagnetic disintegrations of nuclei, since these states contain information about the continuum states.

A quite different approach is the HH effective interaction method. The two-body Hilbert space is divided in two subspaces with projection operators $`P`$ and $`Q`$ ($`P+Q=1`$, $`PQ=0,`$ $`dim_P=N_p`$). In an HH calculation $`P`$ and $`Q`$ spaces are realized via the hyperangular quantum number $`K`$ ($`P`$ space: all HH states with $`KK_{max}`$). Applying the Lee-Suzuki similarity transformation  one decouples $`P`$ and $`Q`$ space interactions in such a way that one obtains an effective interaction in the $`P`$ space,

$$V_P=P\left[\underset{i<j}{}V_{i,j}\right]_{eff}P.$$
(2)
The two-body Hamiltonian with this potential has exactly the same eigenvalues as the lower $`N_p`$ eigenvalues of the bare interaction acting on the full $`P+Q`$ space. Due to the effective two-body interaction one yields an enormous improvement of the convergence for the ground state energies of nuclei with A=3-6 . The introduction of such two-body effective interactions in few-body calculations was first made for the harmonic oscillator basis . In comparison the HH basis offers further advantages, e.g., the presence of collective HH coordinates allows to construct a medium affected two-body effective interaction which leads to a better convergence.

Since one can interpret the effective interaction as a kind of momentum expansion one cannot expect that high momentum components are included in a correct way in the wave function if one works with a small number of HH functions. On the other hand one should find a much improved convergence for observables which contain little information on high momentum components. Thus it is no surprise that also for the nuclear radii an extremely good convergence was observed in Ref. . Similar good convergence results are expected to hold for observables that are governed by not too high momentum components.

Both methods, the correlation (CHH) and the effective interaction approach (EIHH), will be used in the following study of the $`\gamma +^4`$He $`X`$ reaction.

We calculate the nuclear 4He total photoabsorption cross section in the dipole approximation

$$\sigma _{tot}(E_\gamma )=4\pi ^2(e^2/\mathrm{}c)E_\gamma R(E_\gamma )$$
(3)
with

$$R(E_\gamma )=𝑑f|f|D_z|\mathrm{\Psi }_0|^2\delta (E_fE_0E_\gamma ),$$
(4)
where $`D_z`$ is the third component of the nuclear dipole operator $`\stackrel{}{D}`$, while $`E_0`$ and $`E_f`$ denote the eigenvalues of the nuclear Hamiltonian for ground and final states, respectively, and $`E_\gamma `$ is the photon energy. The dipole approximation is well established in low-energy photonuclear reactions. Deviations from the exact result are expected to be very small for the total cross section, since also the important meson exchange current contribution is implicitly taken into account within the dipole approximation (Siegert’s theorem).

With the LIT method the cross section is calculated indirectly. In a first step the LIT of the response function $`R(E_\gamma )`$, i.e.

$$L(\sigma _R,\sigma _I)=𝑑E_\gamma \frac{R(E_\gamma )}{(E_\gamma \sigma _R)^2+\sigma _I^2}=\stackrel{~}{\mathrm{\Psi }}|\stackrel{~}{\mathrm{\Psi }},$$
(5)
is determined via the asymptotically vanishing LIT state $`\stackrel{~}{\mathrm{\Psi }}`$, which fulfills the following boundstate like differential equation

$$(HE_0\sigma _R+i\sigma _I)|\stackrel{~}{\mathrm{\Psi }}=D_z|\mathrm{\Psi }_0.$$
(6)
The second step of the method consists in the inversion of $`L`$ in order to obtain $`R(E_\gamma )`$ (see, e.g. ).

Formally, one can rewrite equation (5) in the following form

$$L(\sigma _R,\sigma _I)=\underset{\nu }{}\frac{|\mathrm{\Psi }_0|D_z|\stackrel{~}{\mathrm{\Psi }}_\nu |^2}{(\stackrel{~}{E}_\nu E_0\sigma _R)^2+\sigma _I^2},$$
(7)
where $`|\stackrel{~}{\mathrm{\Psi }}_\nu `$ ($`\stackrel{~}{E}_\nu `$) are eigenfunctions (eigenvalues) of $`H`$ in a truncated space. They are obtained with the same boundary conditions as $`|\mathrm{\Psi }_0`$. Generally speaking, these states are either bound states or pseudo-resonance states. It is clear that the low energy (small $`\sigma _R`$) behaviour of $`L`$ is dominated by the positions of the lowest eigenvalues $`\stackrel{~}{E}_\nu `$.

In the following we compare results obtained with CHH and EIHH methods for the 4He photodisintegration. Different from the above mentioned calculation of Ref.  correlations are also introduced for the ground state wave function. Our CHH and EIHH calculations agree very well for ground state energy (-30.69 MeV (CHH), -30.71 MeV (EIHH)) and rms matter radius (1.421 fm (CHH), 1.422 fm (EIHH)) being a bit different from those of Ref.  (29.24 MeV, 1.43 fm). However, there is no significant change of the photoabsorption cross section due to the more precise bound state. In fact one finds a small reduction of the response function $`R(E_\gamma )`$, but due to the higher binding energy the decrease is compensated by the increased value for $`E_\gamma `$ (see Eq.(3)).

In Fig.1 we show the convergence patterns of the transform $`L`$ using the MT-I/III potential  as NN interaction. For the CHH case, depicted in Fig. 1a, one observes a very rapid convergence for $`K_{max}`$ values from 1 to 7 ($`K_{max}`$=7 was adopted in Ref. ), but for higher $`K_{max}`$ the convergence is slowed down considerably and $`K_{max}=11`$ does not yet lead to a completely convergent result. The EIHH results of Fig.1b show a much nicer convergence behavior, particularly at low $`\sigma `$. This can be understood in view of the fact that the effective interaction method is, as mentioned, a kind of momentum expansion and thus brings an enormous acceleration to the convergence of the lowest eigenvalues which dominate the low-energy cross section (see Eq.(7)). Besides the CHH results in Fig.1a we illustrate the EIHH transform for $`K_{max}=11`$. One sees that the results are very similar and studying the convergence patterns one might expect that the converged CHH result will come quite close to the EIHH result. In the comparison one should also not forget that small differences might remain even for the converged results, since both calculations are carried out in completely different ways. In particular the CHH calculation is numerically less accurate, since it includes a nine-dimensional Monte Carlo integration.

In Fig.2 we illustrate the results for the total 4He photoabsorption cross section obtained from the inversion of the LITs of Fig.1 with $`K_{max}=7,\mathrm{\hspace{0.17em}9},\mathrm{\hspace{0.17em}11}`$. Again one has a very nice convergence for the EIHH case, while the CHH results are not yet completely convergent. One sees that for increasing $`K_{max}`$ the peak of the CHH cross section is shifted to lower energies in direction of the EIHH peak.

In Fig.3 we show the EIHH cross section results with MT-I/III and TN  potentials ($`K_{max}=11`$) in comparison to experimental data. Though the theoretical results are somewhat different from those of Ref.  the comparison with experiment is not improved. One still observes a considerably stronger giant dipole peak than that found in photoabsorption experiments. At lower energies the situation has changed a bit, since the new results exhibit a stronger deviation from experiment. On the other hand it is evident that the experimental result of Ref. , where the total photoabsorption cross section is extracted from Compton scattering via dispersion relations, agrees much better with the theoretical results similarly to the already mentioned older 4He photoabsorption data.

It will be very interesting to get a further clarification of the experimental cross section from the new experiments at Lund. Since the semi-realistic potential models lead to a rather good result for the 4He rms radius, which is the dominant ingredient in the sum rule for the inverse energy weighted cross section, one may think that there is not much space to change the theoretical cross section. In case of the photodisintegration of the three- nucleon systems one finds only a 10 % reduction of the peak cross section due to more realistic NN interactions and three-nucleon forces . On the other hand it would be extremely interesting if the difference between realistic and semi-realistic interactions for the photonuclear cross section is much higher in the four-nucleon than in the three-nucleon case.


text:
# A sandpile model with tokamak-like enhanced confinement phenomenology

## Abstract

Confinement phenomenology characteristic of magnetically confined plasmas emerges naturally from a simple sandpile algorithm when the parameter controlling redistribution scalelength is varied. Close analogues are found for enhanced confinement, edge pedestals, and edge localised modes (ELMs), and for the qualitative correlations between them. These results suggest that tokamak observations of avalanching transport are deeply linked to the existence of enhanced confinement and ELMs.

preprint: HEP/123-qed

The introduction of the sandpile paradigm- into magnetized plasma physics (fusion-, magnetospheric-, and accretion disk-: for recent reviews see Ref.) has opened new conceptual avenues. It provides a framework within which observations of rapid nondiffusive nonlocal transport phenomena can be studied; recent examples include analyses of auroral energy deposition derived from global imaging, and of electron temperature fluctuations in the DIII-D tokamak, both of which involve avalanching. Insofar as such phenomena resemble those in experimental sandpiles or mathematically idealized models thereof, they suggest that the confinement physics of macroscopic systems (plasma and other) may reflect unifying underlying principles.

In this paper we present results suggesting that this unity may extend to some of the most distinctive features of toroidal magnetic plasma confinement: enhanced confinement regimes (“H-modes”), edge localised modes (“ELMs”), steep edge gradients (“edge pedestal”), and their observed phenomenological and statistical correlations – for recent quantitative studies, see for example Refs. and references therein. An important question is whether the L to H transition necessarily reflects a catastrophic bifurcation of confinement properties, or can be associated with a monotonic change in the character of the turbulence.Transitional behaviour resembling that observed in fusion plasmas has been found in other sandpile algorithms: for example in Ref  a local sheared flow region is imposed which reduces the characteristic avalanche length thereby affecting confinement; and in Ref , changes in the redistribution rule lead to changes in profile stiffness. We show that key elements of the observed phenomenology emerge naturally from a simple one-dimensional sandpile model, that of Chapman , which incorporates other established models as limiting cases. This centrally fueled (at cell $`n=1`$) model’s distinctive algorithmic feature relates to the local redistribution of sand at a cell (say at $`n=k`$) where the critical gradient $`z_c`$ is exceeded: the sandpile is instantaneously flattened behind the unstable cell over a length $`L_f`$, embracing the cells $`n=k(L_f1),k(L_f2),\mathrm{},k`$; and this sand is conservatively relocated to the cell at $`n=k+1`$. $`L_f`$ defines the length scale over which the most rapid redistribution occurs; in a plasma context this could be considered a proxy for turbulent correlation length or eddy size. In Ref. the sandpile is explored for all regimes $`1<L_f<N`$ for both constant and fluctuating critical gradient $`z_c`$. Here we consider the dynamics of the more realistic case with random fluctuations in $`z_c`$; the system is robust in that once fluctuations are introduced in the critical gradient, the behavior is essentially insensitive to both their level and spectral properties, see also Ref.. The limit $`L_f=1`$ is the fixed point corresponding to the centrally fueled algorithm of Ref. in one dimension. In the limit $`L_f=N`$ (where $`N`$ is the number of cells in the sandpile) the sandpile is flattened everywhere behind an unstable cell as in Refs.. A real space renormalization group approach shows that the robust scale free dynamics for the limiting case $`L_f=N`$ corresponds to a nontrivial (repulsive) fixed point (see e.g. Ref.). The essential result of Ref. is that different regimes of avalanche statistics emerge, resembling a transition from regular to intermittent dynamics reminiscent of deterministic chaos. The control parameter is the normalized redistribution scalelength $`L_f/N`$ which specifies whether the system is close to the nontrivial $`L_f=N`$ fixed point.

Height profiles for the sandpile with 512 cells, time averaged over many thousands of avalanches, are shown in Fig.1 for three different values of the fluidization length $`L_f`$ in the range $`50<L_f<250`$. The sandpile profile shape, stored gravitational potential energy, and edge structure (smooth decline or pedestal) correlate with each other and with $`L_f`$. As $`L_f`$ is reduced, the edge pedestal steepens and the time averaged stored energy rises; multiple “barriers” (regions of steep gradient) are visible in trace (a) and to some extent trace (b) of Fig.1. Time evolution of the sandpile for $`L_f=50,150`$, and 250 respectively is quantified in Figs.2-4. The top traces show total stored energy; the middle traces show the position of the edge of the sandpile (the last occupied cell); and the bottom traces show the magnitude and occurrence times of mass loss events (hereafter MLEs) in which sand is lost from the system by being transferred beyond the 512th cell. Time is normalized to the mean inter-avalanche time $`\mathrm{\Delta }\tau `$ (proportional to the fueling rate). The sandpile is fueled only at the first cell, so that the great majority of avalanches terminate before reaching the 512th cell (these are classified as internal). While internal avalanches result in energy dissipation (recorded in the upper traces of Figs.2-4), and may alter the position of the edge of the sandpile, they do not result in an MLE; there are corresponding periods of quiescence in the middle and lower traces of Figs.2-4. Conversely the MLEs are associated with sudden inward movement of the sandpile edge, and in this important sense appear to be edge localised. However, MLEs and the associated inward edge movement are in fact the result of systemwide avalanches triggered at the sandpile center (cell $`n=1`$). The character of the MLEs changes with $`L_f`$. In Fig.2, where the mean and peak stored energy are greatest, the MLEs are similar to each other and occur with some regularity. The regularity of MLE occurrence in Fig.3 is less marked, the magnitude of the largest MLEs is greater than in Fig.2, and there is greater spread in MLE size. This trend continues in Fig.4, which also has the lowest stored energy. These effects correlate with the underlying dynamics of the sandpile. Figure 5 plots the relation between average stored energy and $`L_f`$ for the $`N=512`$ system and much larger $`N=4096`$ and $`8192`$ systems (normalized to the system size $`N`$). The curves coincide, demonstrating invariance with respect to system size, with an inverse power law with slope close to $`2`$ for $`L_f/N<1/4`$, and a break at $`L_f/N1/4`$. These two regimes yield the quasi-regular and quasi-intermittent dynamics in Figs.2-4 (see also the plot of avalanche length distribution against $`L_f`$ in Fig.8 of Ref.). The parameter $`L_f/N`$ is a measure of proximity of this high dimensional system to the $`L_f=N`$ nontrivial fixed point. This determines both the apparent complexity of the timeseries in Figs.2-4 and the underlying statistical simplicity described below, which is also invariant with respect to system size.

There is systematic correlation between time averaged stored energy $`<E>`$ and MLE frequency $`f_{MLE}`$, as shown in Fig.6. To obtain these curves, which are again normalized to system size, we have derived MLE frequencies using a standard algorithm previously used to assign frequencies to ELMs observed in tokamak plasmas in the Joint European Torus (JET). Since the sandpile often generates bursts of mass loss with structure down to the smallest timescales, which might not be resolvable under experimental conditions, we have followed Ref. in applying a (relatively narrow) measurement window of width $`450\mathrm{\Delta }\tau `$ to obtain $`f_{MLE}`$. The correlation between $`<E>`$ and $`f_{MLE}`$ is a noteworthy emergent property, furthermore Fig.6’s characteristic curve is very similar to that of Fig.6 of Ref., which relates measured energy confinement to ELM frequency in JET. Energy confinement time $`\tau _c`$ can be defined for the sandpile by dividing the time averaged stored energy $`<E>`$ by the time averaged energy dissipation rate $`<\mathrm{\Delta }E>`$ (where $`\mathrm{\Delta }E`$ is the energy dissipated in a single avalanche). The embedded plot of Fig. 6 shows $`\tau _c`$ against MLE frequency $`f_{MLE}`$.

Finally, we explore the situation where there is a secular change in the redistribution algorithm: in Fig.7, $`L_f`$ decreases slowly, continuously, and linearly with time, from one constant value to another over a period encompassing many tens of thousands of avalanches. There is a corresponding time evolution of the energy confinement properties of the sandpile and of the character of the MLEs. Figure 7(top) shows total stored energy as a function of time as $`L_f`$ changes from $`250`$ at $`t=4\times 10^4`$ to $`50`$ at $`t=1.15\times 10^5`$, while $`10^5`$ avalanches occur: over a period of time corresponding to a few tens of MLEs, the system smoothly evolves from low to high confinement. This is accompanied by a gradual change in character of the time variation in the sandpile edge (position of last occupied cell, Fig.7(middle)) and of the MLEs (Fig.7(lower)), from large amplitude to small and from irregular to regular. Figure 7 can perhaps be regarded as this sandpile’s analogue of, for example, Fig.2 of Ref. or Fig.2 of. The essential point here is that the sandpile apparently freely explores phase space with changing control parameter $`L_f/N`$. Characteristic properties of the dynamics (whether quasi-regular or quasi-intermittent) and correspondingly, confinement properties (such as stored energy and MLE characteristics) smoothly follow changes in this parameter rather than exhibiting a sudden phase transition or catastrophe.

By varying a single control parameter in the sandpile algorithm, we have shown correlations between: stored energy, confinement times, sandpile profile, sandpile edge structure, and the amplitude, frequency, and dynamical character of mass loss events. We have also seen how slow secular change in the control parameter produces a smooth evolution in confinement properties. If a single control parameter analogous to $`L_f/N`$ exists for tokamaks, it can in principle be found from experimental data by examining scaling with respect to system size as above.

The existence of such extensive tokamak-like phenomenology, emergent from a very simple system, is a novel discovery. Insofar as the phenomenological resemblance is close, there is more to be learnt. A minimalist interpretation starts from the premise that the sandpile algorithm of Ref.  provides a simple one-parameter model for studying generic nonlocal transport, conditioned by a critical gradient, in a macroscopic confinement system. Changing the value of the single control parameter $`L_f`$ then corresponds to altering the range in configuration space over which the transport process operates. It then follows from the results in the present paper that this may be the minimum requirement to generate those aspects of tokamak-like confinement phenomenology described. This is a significant conclusion, but one can consider a more far-reaching one. A possible maximalist interpretation attaches greater weight to recent observations of avalanching transport in tokamaks and in largescale numerical simulations thereof, and therefore regards the avalanching transport that is built into sandpile algorithms as an additional point of contact with magnetically confined plasmas. One would then infer from the present results that tokamak observations of avalanching transport are deeply linked to the existence of enhanced confinement and ELMs.

###### Acknowledgements.

We are grateful to Jack Connor, George Rowlands, David Ward and Nick Watkins for comments and suggestions. SCC was supported by a PPARC lecturer fellowship, ROD by Euratom and the UK DTI, and BH by HEFCE.

Captions

FIG.1. Time averaged height profiles of the 512 cell sandpile for $`L_f=`$ (a)50, (b)150 and (c)250. Inset: edge structure.

FIG.2. Time evolution of the 512 cell sandpile with $`L_f=50`$: (top) stored energy, (middle) position of last occupied cell, (lower) magnitude and occurence of mass loss events.

FIG.3. As Fig.2, for $`L_f=150`$.

FIG.4. As Fig.2, for $`L_f=250`$.

FIG.5. Average stored energy versus $`L_f/N`$ for sandpiles of $`N=512,4096,8192`$. Energy is normalized to the $`L_f=1`$ case (effectively to $`N^2`$).

FIG.6. Average stored energy versus MLE frequency, and (inset) $`\tau _c`$ versus MLE frequency for sandpiles of $`N=512,4096,8192`$. Energy and MLE frequency are normalized as in Fig 6.

FIG.7. Time evolution of (top) stored energy, (middle) sandpile edge position and (lower) MLEs, as $`L_f`$ changes slowly and linearly from 250 to 50.