format: text

text:
# Full Aging in Spin Glasses

## Abstract

The discovery of dynamic memory effects in the magnetization decays of spin glasses in 1983 marked a turning point in the study of the highly disordered spin glass state. Detailed studies of the memory effects have led to much progress in understanding the qualitative features of the phase space. Even so, the exact nature of the magnetization decay functions have remained elusive, causing confusion. In this letter, we report strong evidence that the Thermoremanent Magnetization (TRM) decays scale with the waiting time, $`t_w`$. By employing a series of cooling protocols, we demonstrate that the rate at which the sample is cooled to the measuring temperature plays a major role in the determination of scaling. As the effective cooling time, $`t_c^{eff}`$, decreases, $`\frac{t}{t_w}`$ scaling improves and for $`t_c^{eff}<20s`$ we find almost perfect $`\frac{t}{t_w}`$ scaling, i.e full aging.

preprint: APS/123-QED

Since the discovery of aging effects in spin glasses approximately twenty years agoCham83 Lund83 , much effort has gone into determining the exact time dependence of the memory decay functions. In particular, memory effects show up in the Thermoremanent Magnetization (TRM) (or complementary Zero-Field Cooled (ZFC) magnetization) where the sample is cooled through its spin glass transition temperature in a small magnetic field (zero field) and held in that particular field and temperature configuration for a waiting time, $`t_w`$. At time $`t_w`$, a change in the magnetic field produces a very long time decay in the magnetization. The decay is dependent on the waiting time. Hence, the system has a memory of the time it spent in the magnetic field. A rather persuasive argumentBou92  suggests that for systems with infinite equilibration times, the decays must scale with the only relevant time scale in the experiment, $`t_w`$. This would imply that plotting the magnetization on a t/$`t_w`$ axis would collapse the different waiting time curves onto each other. This effect has not been observed.

What has been observedAlba86  is that the experimentally determined magnetization decays will scale with a modified waiting time, $`(t_w)^\mu `$. Where $`\mu `$ is a fitting parameter. For $`\mu <1`$ the system is said to have subaged. A $`\mu >1`$ is called superaging and $`\mu =1`$ corresponds to full aging. For TRM experiments a $`\mu `$ of approximately 0.9 is found for different types of spin glassesAlba86 ; Ocio85 , over a wide range of reduced temperatures indicating subaging. At very low temperatures and temperatures approaching the transition temperature, $`\mu `$ is observed to decrease from the usual 0.9 value. Superaging has been observed in Monte Carlo simulations of spin glassesSibani . This has led to confusion as to the exact nature of scaling.

Zotev et al.Zov02  have suggested that the departures from full $`\frac{t}{t_w}`$ scaling, observed in aging experiments, are mainly due to cooling effects. In a real experimental environment, the situation is complicated by the time it takes for the sample to cool to its measuring temperature. An effect due to the cooling rate at which the sample temperature approaches the measuring temperature has been knownNord87 ; Nord00 . This effect is not trivial, it does not contribute a constant time to $`t_w`$.

Another possible explanation for the deviation from full aging comes from the widely held belief that the magnetization decay is an additive combination of a stationary term ($`M_{Stat}=A(\tau _0/t)^\alpha )`$ and an aging term ($`M=f(\frac{t}{t_w})`$ )Cou95 ; Bou95 ; Vin96 . Subtraction of a stationary term, where $`\tau _0`$ is a microscopic spin flipping time, A is a dimensionless constant and $`\alpha `$ is a parameter determined from $`\chi ^{\prime \prime }`$ measurements, was shown to increase $`\mu `$ from 0.9 to 0.97Vin96 .

In this letter we analyze effects of the cooling time through a series of different cooling protocols and we present the first clear and unambiguous experimental evidence that the TRM decays scale as $`\frac{t}{t_w}`$(i.e. full aging).

Three different methods have been regularly employed to understand the scaling of the TRM decays. The first and simplest is to scale the time axis of the magnetization decay with the time the sample has spent in a magnetic field (i.e.$`\frac{t}{t_w}`$)Bou92 . If the decays scale as a function of waiting time it would be expected that the decay curves would overlap. This has not yet been observed.

A second more sophisticated method was initially developed by StruikStu79  for scaling the dynamic mechanical response in glassy polymers and first applied to spin glasses by Ocio et al.Ocio85 . This method plots the log of the reduced magnetization $`M/M_{fc}`$ ($`M_{fc}`$ is the field cooled magnetization), against an effective waiting time $`\xi =\frac{\lambda }{t_w^\mu }`$ where

$`\lambda ={\displaystyle \frac{t_w}{1\mu }}[(1+{\displaystyle \frac{t}{t_w}})^{1\mu }1];\mu <1`$ (1)
or

$`\lambda =t_w\mathrm{log}[1+{\displaystyle \frac{t}{t_w}}];\mu =1`$ (2)
A value of $`\mu `$ = 1 would correspond to perfect $`t/t_w`$ scaling. Previous values of $`\mu `$ obtained on the decays using this method have varied from .7 to .94Alba86  for a temperature range $`.2<T_r<.95`$. A value of $`\mu <1`$ is called subaging.

Finally a peak in the function S(t)

$`S(t)={\displaystyle \frac{1}{H}}{\displaystyle \frac{dM(t)}{d[Log_{10}(t)]}}`$ (3)
as a function of time has been shown to be an approximately linear function of the waiting timeLund83 . This peak occurs at a time slightly larger then the waiting time, again suggesting possible subaging.

In this study we use all three of the above scaling procedures to analyze the data we have produced with different cooling protocols. All measurements in this letter were performed on our homebuilt DC SQUID magnetometer with a $`Cu_{.94}Mn_{.06}`$ sample. The sample is well documentedKen91  and has been used in many other studies. The measurements described in this letter were performed in the following manner: The sample was cooled, in a magnetic field of 20 G, from 35 K through its transition temperature of 31.5 K to a measuring temperature of 26 K. This corresponds to reduced temperature of .83 $`T_g`$. The sample was held at this temperature for a waiting time $`t_w`$, after which time the magnetic field was rapidly decreased to 0G. The resulting magnetization decay is measured 1s after field cutoff to a time greater than or equal to 5$`t_w`$. The only parameters we have varied in this study are $`t_w`$ and the rate and profile at which we cool the sample through the transition temperature to the measuring temperature. The sample is located on the end of the sapphire rod and sits in the upper coil of a second order gradiometer configuration. The temperature measuring thermometer is located 12.5 cm above the sample. Heat is applied to the sample through a heater coil located on the same sapphire rod 17 cm above the sample. Sample cooling occurs by heat transfer with the He bath via a constant amount of He exchange gas, which was previously introduced into each chamber of the double vacuum jacket. We have measured the decay time of our field coil and find that we can quench the field in less then 0.1ms. We have also determined that without a sample, our system has a small reproducible exponential decay, that decays to a constant value less than the system noise within 400 seconds. In order to accurately describe our data we subtract this system decay from all of the data.

In this paper we present TRM data for eight waiting times ($`t_w`$= 50s, 100s, 300s, 630s, 1000s, 3600s, 6310s, and 10000s). The same TRM experiments were performed for six different cooling protocols. In this paper we use four of the cooling protocols. Figure 1 (top row) is a plot of temperature vs. time for four of the cooling protocols. These different cooling protocols, were achieved, by varying applications of heat and by varying the amount of exchange gas in the vacuum jackets. A more detailed description of the cooling protocols will be given in a followup publication. In Figure 1 (bottom row) we plot S(t) (Eq $`\text{3}.)`$ of the ZTRM protocol (i.e. zero waiting time TRM) in order to characterize a time associated with the cooling protocol, $`t_c^{eff}`$. As observed in Figure 1 we have achieved effective cooling times ranging from 406s down to 19 seconds. These times can be compared with commercial magnetometers which have cooling times in the range of 100-400s.

In Figure 2, we plot the data for the TRM decays (first column) with the four cooling protocols. It should be noted that the magnetization (y-axis) is scaled by the field cooled magnetization.

The second column is the same data, as column one, with the time axis(x-axis) normalized by $`t_w`$. It can be observed for $`\frac{t}{t_w}`$ scaling (column 2) that as the effective cooling time decreases the spread in the decays decreases giving almost perfect $`\frac{t}{t_w}`$ scaling for the 19 second cooling protocol.

The last column in Figure 2, is the data scaled,using $`\mu `$ scaling which has previously been described. It has long been known that the rate of cooling affected $`\mu `$ scaling and that $`\mu `$ scaling is only valid in the limit $`t_w>>t_c^{eff}`$. We find this to be true and that the limit is much more rigorous than previously believed. To determine $`\mu `$ scaling, we focused on applying this scaling to the longest waiting time data (i.e., $`t_w`$ = 3600s, 6310s and 10,000s). For the largest effective cooling time data, $`t_c^{eff}=406s`$, we find that we can fit the longest waiting time data with a $`\mu `$ value of .88. This is consistent with previously reported values of $`\mu `$Alba86 . We do find, however, that TRM data with waiting times less then 3600s do not fit on the scaling curve. We find that scaling of the three longest waiting time decays produces $`\mu `$ values which increase as $`t_c^{eff}`$ decreases. We also find that as $`t_c^{eff}`$ decreases, the data with shorter $`t_w`$ begins to fit to the scaling better. It can be observed that at $`t_c^{eff}`$= 19s we obtain almost perfect scaling for all of the data with a value of $`\mu =.999`$. However we find we can reasonably fit the data to a range of $`\mu `$ between .989-1.001. The fitting for the large $`t_w`$ decay curves, $`t_w`$ = 3600s, 6310s and 10,000s, is very very good. Small systematic deviations, as a function of $`t_w`$, occur for $`t_w<3600s`$ with the largest deviations for $`t_w`$= 50s. Even with an effective cooling time two orders of magnitude less then the waiting time, one sees deviations from perfect scaling. We have also scaled the data using Eq. 2. We find no noticeable difference between the quality of this fit and the quality of the $`\mu =.999`$, for $`t_c^{eff}=19s`$, shown in Figure 2. Data with longer cooling times cannot be fit with Eq. 2. We therefore conclude that full aging is observed for the long $`t_w`$ data using the $`t_c^{eff}=19s`$ protocol.

It is clear from Figure 1 (bottom row), that the effect of the cooling time has implications for the decay all the way up to the longest time measured, 10,000 seconds. The form of the S(t) of the ZTRM is very broad. The S(t) function is often thought of as corresponding to a distribution of time scales (or barrier heights) within a system which has infinite equilibration times (barriers). The peak in S(t) is generally associated with the time scale (barrier height) probed in time $`t_w`$. In Figure 1, we observe that for the larger effective cooling times, the waiting times correspond to points on or near the peak of S(t) for the ZTRM. We therefore believe that for the larger effective cooling times there is significant contamination from the cooling protocol over the entire region of $`t_w`$s used in this paper. Only for $`t_c^{eff}=19s`$ cooling protocol do we find that the majority of $`t_ws`$ occur far away from the peak in the S(t).

All the data in figure 3 used the cooling protocol with $`t_c^{eff}=19s`$. In Figure 3a we plot the S(t)(ZTRM) for $`t_c^{eff}=19s`$ with arrows to indicate the waiting times for the TRM measurements. It can be observed that after approximately 1000 seconds the slope of the S(t) function decreases, possibly approaching a horizontal curve, which would correspond to a pure logarithmic decay in M(t). If, on the other hand, the slope is continuously changing this part of the decay may be described by a weak power law. Either way, this region would correspond to aging within a pure non-equilibrated state. We believe that the long waiting time data occurs outside the time regime that has been corrupted by the cooling time and that this is the reason that we have, for the first time, observed full aging.

It has been suggested that subtraction of a stationary component of the magnetization decay will improve scalingVin96 . The very long time magnetization decay is believed to consist of a stationary term that is thought to decay as a power law. We fit a power law, $`M(t)=A(\tau _o/t)^\alpha `$ to the long time decay (1000s-5000s), of the ZTRM for $`t_w`$=19s. Using $`\tau _o=10^{12}s`$, we find $`\alpha =.07`$ and $`A=.27`$. Subtracting this power law form from the magnetization decay destroys scaling. We find that the subtraction of a much smaller power law term with A=.06 and $`\alpha `$=.02 slightly improves scaling at both short and long times. While the $`\alpha `$ values for the two different power law terms we have fit to are quite different, both values fall within the range determined from the decay of $`\chi `$Vin96 . In Figure 3(c-d) and 3(e-f), we plot the two different types of scaling we have performed, with and without the subtraction of the weaker power law term.

We find that even for $`t_c^{eff}=19s`$ the peak in S(t) for $`t_w>1000s`$ occur at a time larger then $`t_w`$(fig 3b). We find that we can fit the effective time associated with the peak in S(t) to $`t_w^{eff}=t_w^{1.1}`$.

In summary, we have performed TRM decays over a wide range of waiting times (50s - 10,000s) for six different cooling protocols. We find that as the time associated with the cooling time decreases, scaling of the TRM curves improves in the $`\frac{t}{t_w}`$ scaling regime and in the $`\mu `$ scaling regime. In $`\mu `$ scaling we find that as the effective cooling time decreases $`\mu `$ increases approaching a value of .999 for $`t_c^{eff}=19s`$. For the $`t_c^{eff}=19s`$ TRM decays, we find that subtraction of a small power law term (A=.06, $`\alpha `$=.02) slightly improves the scaling. It is however likely that the small systematic deviations of the $`t_c^{eff}=19s`$ data as a function of $`t_w`$ are associated with the small but finite cooling rate.

The authors would like to thank V. S. Zotev, E. Vincent and J. M. Hammann for very helpful discussions.


text:
# Energy exponents and corrections to scaling in Ising spin glasses

## I Introduction

Exponents and corrections to scaling play a central role in disordered systems. Consider for instance the dependence of thermodynamic quantities on the size $`L`$ of a system when $`L\mathrm{}`$. In paramagnetic systems, correlations are short range and thus bulk properties converge rapidly to their thermodynamic limits. On the contrary, in systems such as spin glasses below $`T_c`$, correlations are long range and so finite size corrections are “large”. In practice, that means that even if one takes periodic boundary conditions, the disorder average of an intensive quantity converges to its thermodynamic limit slowly: finite size corrections go to zero as an inverse power of $`L`$. This slow convergence is to be contrasted with paramagnetic systems where finite size corrections in disorder averages are exponentially small. One can also consider the sample-to-sample fluctuations of thermodynamic quantities. These fluctuations scale in the thermodynamic limit; an interesting open question is their limiting distribution, if any.

In this work we consider these issues in the context of Ising spin glasses MezardParisi87b ; Young98  at zero temperature and focus on the ground state energy. Given a distribution of disorder realisations (each being represented by $`J`$), the extensive ground state energy $`E_J`$ is a random variable. Denote by $`N`$ the number of spins in the system; $`N=L^d`$ for a $`d`$-dimensional hyper-cubic lattice. We are interested in the probability distribution $`P_N(E_J)`$ and in how the connected moments (cumulants) of this distribution depend on $`N`$ (or $`L`$). Denoting disorder averages by an overline, we have

$$\overline{E_J(L)}=e_0L^d+e_1L^{\mathrm{\Theta }_s}+\mathrm{}$$
(1)
In this expression, $`\mathrm{\Theta }_s`$ is what we call the shift exponent. To leading order, the energy scales with the volume, while $`\mathrm{\Theta }_s`$ gives the leading correction to scaling. The justification for this notation will become clear later; in the mean time, we should compare to the usual notation involving the correction to scaling exponent $`\omega `$; indeed, in the spin glass phase, finite size corrections of intensive quantities decay as a power of $`L`$:

$$\overline{E_J(L)/L^d}=e_0+e_1L^\omega +\mathrm{}$$
(2)
and so $`\omega =d\mathrm{\Theta }_s`$.

We will also look at the width of $`P_N`$ which gives the scaling of the instance to instance fluctuations:

$$\left[\overline{E_J^2(L)}\overline{E_J(L)}^2\right]^{1/2}=\sigma _0L^{\mathrm{\Theta }_f}+\mathrm{}$$
(3)
$`\mathrm{\Theta }_f`$ being the fluctuation exponent. Naturally, one may extend these types of expansions to any cumulant of $`P_N`$, a central question being whether $`P_N`$ has a limiting shape as $`N\mathrm{}`$. One motivation for this is the fact that the ground state energy is an extreme statistic: there are $`2^N`$ configurations of spins ($`S_i=\pm 1`$) and one is interested in the one of minimum energy. (Note that these $`2^N`$ random energies are correlated.) Little is known about such statistics except in a few solvable cases Galambos87 ; CarpentierLedoussal01 . The purpose of this work is to determine $`P_N(E_J)`$ numerically for a variety of spin glass models. From our measurements we extract $`\mathrm{\Theta }_s`$, $`\mathrm{\Theta }_f`$, etc… and compare these results to theory.

The outline of this paper is as follows. We first consider expectations arising from analogies with other systems as well as some known results. In section III we describe the different types of spin glasses used in this study. Then we give in section IV our estimates of the fluctuation exponent $`\mathrm{\Theta }_f`$. This is followed by a study of the probability distribution of the ground state energy in section V; for most models we find that it becomes Gaussian in the large system size limit. Then in section VI we present our results for the shift exponent $`\mathrm{\Theta }_s`$. For all of these studies, we compare the $`d`$-dimensional case to the theoretical predictions. In section VII we show that these exponents sometimes depend on the distribution of the disorder variables (the spin-spin couplings). Finally in section VIII we discuss and conclude this work. Some details of the analytic computations are given in the appendix.

## II Clues from theory

What values should be expected for the two exponents $`\mathrm{\Theta }_s`$ and $`\mathrm{\Theta }_f`$? Suppose we start with the random energy model (REM) Derrida81  as a guide. Its ground state energy $`E_J`$ has a Gumbel distribution in the large $`N`$ limit BouchaudMezard97  with a variance of $`O(1)`$. Also, the disorder averaged ground state energy grows as $`e_0N+O(\mathrm{ln}N)`$ at large $`N`$. These properties show that $`\mathrm{\Theta }_s=\mathrm{\Theta }_f=0`$ and thus in the REM finite size corrections are “small”, though clearly much larger than in a paramagnetic system.

Another model that can guide us is the directed polymer in a random medium (DPRM) Halpin-HealyZhang86 . The problem on a tree is very similar to the REM, although subtle differences appear CarpentierLedoussal01 . The two dimensional case (with one space and one “time” dimension) can also be solved in much detail. In these cases, perhaps surprisingly, a single exponent $`\theta `$ describes the scaling of three a priori unrelated quantities: (1) the corrections to scaling of the disorder averaged ground state energy $`E_J`$; (2) the size of sample-to-sample fluctuations of $`E_J`$; (3) the typical excitation energy of the lowest “system-size” excitation, that is an excitation that is macroscopically different from the ground state. These properties lead to the remarkable relation $`\mathrm{\Theta }_s=\mathrm{\Theta }_f=\theta `$. The DPRM is thus described by a “one-parameter” scaling theory.

It is plausible that an extension of the DPRM scaling theory may apply to Ising spin glasses. At the heart of such a theory, initiated by McMillan McMillan84 , developped by Bray and Moore BrayMoore86  and extended by Fisher and Huse FisherHuse86 , is the exponent $`\theta _{DW}`$. This exponent is analogous to the $`\theta `$ of the DPRM, and gives the scaling of domain wall energies, $`E_{DW}\mathrm{{\rm Y}}L^{\theta _{DW}}`$; numerical estimates give $`\theta _{DW}0.20`$ in $`d=3`$ and $`\theta _{DW}0.28`$ in $`d=2`$. (More generally, a number of other $`\theta `$ exponents have been introduced for spin glasses; they are all associated with excitation energies; this is to be contrasted with our $`\mathrm{\Theta }`$ exponents that are associated with ground state energies.) If we follow the correspondence with the DPRM, we expect that $`\mathrm{\Theta }_s=\theta _{DW}`$, justifying our use of a “theta” notation for $`\mathrm{\Theta }_s`$. Physically, this equality corresponds to the fact that these systems are sensitive to boundary conditions; for some samples, these conditions are such that a “domain wall” must be present in the ground state. However the analogy with the DPRM certainly breaks down for $`\mathrm{\Theta }_f`$: for any short range spin glass in dimension $`d`$, Wehr and Aizenman WehrAizenman90  proved that $`\mathrm{\Theta }_f=d/2`$. This shows that the REM and the DPRM are not good guides for finite dimensional spin glasses.

To have more realistic theoretical predictions, and in particular to preserve $`\mathrm{\Theta }_f=d/2`$, it seems necessary to work with models having a microscopic Hamiltonian defined over configurations of $`N`$ spins. One approach is to use hierarchical (Migdal-Kadanoff) lattices; there, analytical computations as well as powerful numerical methods are possible. We shall also consider mean field spin glasses where spins are coupled amongst one-another at random so there is no geometry to speak of. The $`P_N(E_J)`$ in such models can be referred to as the mean field prediction. Strangely enough, little is know about these systems so we will have to determine their behavior numerically. In some cases these models lead to surprizes as we shall soon see.

## III Models and methods

We focus on three families of Ising spin glass models so that the effects of geometry and dimension are apparent. The Hamiltonian for these models is

$$H=\underset{<ij>}{}J_{ij}S_iS_j$$
(4)
where $`S_i=\pm 1`$ and the $`J_{ij}`$ are quenched Gaussian random variables of zero mean and variance $`1`$ (except in the case of the Sherrington-Kirkpatrick model which we discuss later). The sum $`_{<ij>}`$ is over nearest neighbor spins on a given graph having $`N`$ nodes; the different models we consider vary simply by the nature of that graph.

Our first family of models are of the Edwards-Anderson (EA) type EdwardsAnderson75 : the graphs are square or cubic lattices of linear size $`L`$ (thus $`N=L^d`$ when the dimension is $`d`$), and the edges connect nearest neighbors only; we take periodic boundary conditions in all directions.

Our second family of models come from the Migdal-Kadanoff (MK) approach SouthernYoung77  where one performs a bond-moving real-space renormalization group. This procedure effectively amounts to computing quantities on hierarchical (MK) lattices defined by an iteration process (see Fig. 1). The iteration takes one bond (that is an edge of the current graph) into $`b`$ paths, each made of $`2`$ segments (that is edges); if $`r`$ is the iteration number (beginning with $`r=0`$), the “linear” lattice size $`L`$ grows as $`2^r`$ and the volume (actually the number of edges and thus the number of terms contributing to the energy) grows as $`(2b)^r`$. When using $`\mathrm{}`$ segments instead of $`2`$ in each path, we have

$$L=\mathrm{}^r\mathrm{and}N(\mathrm{}b)^r$$
(5)
so that the dimension is $`d=\mathrm{ln}N/\mathrm{ln}L=1+\mathrm{ln}b/\mathrm{ln}\mathrm{}`$. The usual choice to obtain $`d=3`$ is $`\mathrm{}=2`$ and $`b=4`$, while $`\mathrm{}=b=2`$ corresponds to $`d=2`$ as in Fig. 1.

Our third family of models come from mean field, and here we have considered two types of graphs. First, we use complete graphs where all vertices are connected, corresponding to the Sherrington-Kirkpatrick (SK) model SherringtonKirkpatrick75 . To have an extensive energy, one takes the variance of the $`J_{ij}`$ to be $`1/N`$. Second, we also use diluted models for which the connectivity is fixed and identical for all the vertices of the graph DominicisGoldschmidt89 . The disorder ensemble then consists of the uniform distribution over all graphs satisfying that constraint in addition to the disorder ensemble in the $`J_{ij}`$. Such an ensemble can be used to “model” the Euclidean case by setting the coordination to that of the lattice of interest. Thus to model the $`d`$-dimensional EA model on the hypercubic lattice, we set the coordination to $`2d`$.

For these three families of spin glass models, we will determine the distribution of the ground state energies. The hierarchical lattices allow one to write a recursion for $`P_N(E_J)`$; because of that, there is no need to sample the disorder variables, we can perform the disorder average exactly. On the contrary, for the EA and the mean field cases, we must compute the actual ground state energies for a large number of disorder samples. For that we rely on a previously tested HoudayerMartin01  algorithmic procedure where given enough computational ressources, the ground state is obtained with a very high probability for both Euclidean lattices and for random graphs as long as $`N`$ is not too large. For our runs, we used several months of CPU time on Pentium III personal computers running at 666 MHz. With this amount of CPU, we obtained high statistics for lattices of sizes up to $`L=10`$ (in both $`d=2`$ and $`3`$) and for mean field graphs with $`N`$ up to $`300`$.

## IV Sample to sample fluctuations and the exponent $`\mathrm{\Theta }_f`$

### IV.1 Migdal-Kadanoff lattices

Let us begin with the Migdal-Kadanoff approach for which the important prediction $`\mathrm{\Theta }_f=d/2`$ can be derived. To understand how this relation comes about on the hierarchical lattices, we construct these by “aggregation”, i.e., by recursively (and hierarchically) joining sub-lattices together. (This procedure is to be distinguished from the top-down iteration used in Fig. 1.) We work with the distribution of ground state energies $`E_p`$ (respectively $`E_a`$) subject to fixed boundary conditions, the spins on the ends of the lattice being forced to be parallel (respectively antiparallel.) Let $`\sigma ^2`$ be the variance of the ground state energy at some level of the hierarchical construction. To go to the next level for which $`L`$ will be $`\mathrm{}`$ times larger, first find the ground state energy in one of the $`b`$ paths. Clearly, if we take the ground state configuration in each of the segments of that path, we will have built the ground state for unconstrained end spins. If the result does not give the imposed values for the end spins (parallel or antiparallel), one must add a “correction” term equal to the smallest domain wall energy of the $`\mathrm{}`$ segments of that path. Thus the energy of one path is the sum of the ground state energies of the $`\mathrm{}`$ segments, plus one domain wall energy half of the time. Second, we add up the independent contributions from the $`b`$ different paths, leading to $`E_p`$ and $`E_a`$ at the new level. The ground state energy is then simply $`min(E_p,E_a)`$. If we neglect the “correction” term, then $`\sigma ^2`$ at this new level is just $`\mathrm{}b`$ times larger than at the previous level. At large $`b`$, the correction is in fact small and so it can be neglected. From one level to the next, the volume grows by a factor $`\mathrm{}b`$, just as $`\sigma ^2`$ does, so the variance is linear in the lattice volume and thus $`\mathrm{\Theta }_f=d/2`$. From our numerical study of these hierarchical lattices, we find that this relation holds also for small $`b`$ and for both $`\mathrm{}=2`$ and $`\mathrm{}=3`$. (We did not test for larger $`\mathrm{}`$.)

As mentioned in the introduction, the Wehr-Aizenman theorem WehrAizenman90  shows that $`\mathrm{\Theta }_f=d/2`$ in finite dimensional spin glasses. It is rather conforting that the MK approach also leads to this result, sustaining the belief that it is a useful guide for real (finite dimensional) spin glasses.

### IV.2 Mean field models

Next, consider the mean field prediction for $`\mathrm{\Theta }_f`$. Since there is no geometry in our mean field family of models, we identify $`L`$ with $`N^{1/d}`$ and thus

$$\left[\overline{E_J^2(N)}\overline{E_J(N)}^2\right]^{1/2}=\sigma _0N^{\mathrm{\Theta }_f/d}+\mathrm{}$$
(6)
The issue here is whether mean field also predicts $`\mathrm{\Theta }_f/d=1/2`$ as it should if one believes that this approach correctly describes the large dimensional limit of real spin glasses. To find out, we have performed ground state computations on thousands of samples of the SK model and of fixed connectivity spin glasses, and have extracted for each ensemble its associated $`\mathrm{\Theta }_f/d`$.

First, consider the fixed connectivity spin glasses with the connectivities, $`z=3,4,6`$ and $`10`$. To estimate $`\mathrm{\Theta }_f/d`$, we determine when the root mean square deviation (RMS) of the ground state energy divided by $`N^{\mathrm{\Theta }_f/d}`$ becomes flat as a function of $`N`$. (Our runs were performed for $`50N300`$.) This leads to $`\mathrm{\Theta }_f/d0.5`$ for $`z=3`$, $`\mathrm{\Theta }_f/d0.49`$ for $`z=4`$, $`\mathrm{\Theta }_f/d0.48`$ for $`z=6`$, and $`\mathrm{\Theta }_f/d0.44`$ for $`z=10`$. The drift we observe in $`\mathrm{\Theta }_f/d`$ is most probably an artefact of our procedure, and simply corresponds to the fact that the corrections to the scaling in Eq. 6 are important in our data, especially at large $`z`$. To get better estimates of $`\mathrm{\Theta }_f/d`$, we would need to control these corrections to scaling but our data are not sufficiently precise for that. Nevertheless, it seems very plausible that $`\mathrm{\Theta }_f/d=0.5`$ in all the fixed connectivity models. In direct analogy to what was stated previously, we can also say that this result sustains the belief that mean field models provide a useful guide to real (finite dimensional) spin glasses.

Let us now continue and consider the limit of infinite connectivity, i.e., the SK model. Since there are $`O(N^2)`$ terms contributing to the Hamiltonian, a very simple minded guess would give $`\mathrm{\Theta }_f/d=1`$ and thus larger fluctuations than in the other models. However the opposite happens, revealing that the scaling in the SK is quite subtle. One can get some clues from Kondor’s analytic study Kondor83  performed just below the critical temperature $`T_c`$. In particular, his results have been interpreted by Crisanti et al. CrisantiPaladin92  who argued that the free energy fluctuations should scale as $`N^{1/6}`$, which is very small compared to the naïve estimate. A different estimate was recently proposed by Aspelmeier et al. AspelmeierMoore03b  who argued that free energy fluctuations should scale as $`N^{1/4}`$. Our concern here is the ground state energy; although there are no analytical calculations, it is plausible that the exponent for energy fluctuations at $`T=0`$ is the same as at that for free energy fluctuations at $`0<T<T_c`$. (Such an extrapolation to $`T=0`$ is known to apply to the DPRM where the exponent of the free energy fluctuations at $`T>0`$ is equal to the one of the ground state energy fluctuations.) In that case, we would have $`\mathrm{\Theta }_f/d=1/6`$ according to some authors and $`\mathrm{\Theta }_f/d=1/4`$ according to others.

What do the numerical estimates tell us about this question? We are aware of a study by Cabasino et al. CabasinoMarinari88  who showed beyond any reasonable doubt that $`\mathrm{\Theta }_f/d<0.5`$; in fact their best fit gives $`\mathrm{\Theta }_f/d0.28`$. When we consider our data, motivated by Eq. 6, we find that the ratio of the left and right hand sides is compatible with a constant when we use $`\mathrm{\Theta }_f0.25`$ and $`40N150`$, but the two points at $`N=200`$ and $`300`$ are then below the others as show in Fig. 2. It is difficult to extract an error bar on the value of this exponent, and most likely the terms dropped in Eq. 6 are important just as was the case in the fixed connectivity models. In fact, here the difficulty is even more acute as the standard deviation increases only very slowly with $`N`$. Nevertheless, let us fit the standard deviation to a pure power; then using a $`\chi ^2`$ analysis we obtain $`\mathrm{\Theta }_f/d=0.24\pm 0.005`$ which is very close to the conjectured $`1/4`$ value, a result further supported by other numerical work PalassiniPrivate . Note however that throughout all this paper, the error bars given are statistical only…In our fit, there are $`7`$ degrees of freedom and the resulting $`\chi ^29.5`$ is not bad; nevertheless, a critical examination of the figure leads one to conclude that the actual uncertainty on $`\mathrm{\Theta }_f/d`$ is certainly much greater than $`0.005`$. Note that this estimate is a bit lower that that of CabasinoMarinari88  and is many standard deviations away from the conjecture $`\mathrm{\Theta }_f/d=1/6`$. However, if we impose $`\mathrm{\Theta }_f/d=1/6`$, the data is reasonably straight as a function of $`N^{2/3}`$ as shown in the inset of Fig. 2. Because of this “good” behavior, we cannot rule out $`\mathrm{\Theta }_f/d=1/6`$.

Given that in the dilute spin glasses we expect $`\mathrm{\Theta }_f/d=0.5`$, how does one recover the SK case as $`z\mathrm{}`$? Begin by recalling that, in the SK model, the variance of the $`J_{ij}`$ is taken to be $`1/N`$ to ensure that the ground state energy is extensive. To maintain this property in the diluted spin glasses, we must divide the $`J_{ij}`$ by $`\sqrt{z}`$. The fluctuations for large $`N`$ and $`z`$ in these modified spin glasses scale as

$$\left[\overline{E_J^2(N)}\overline{E_J(N)}^2\right]^{1/2}=\sigma _0(z)N^{\mathrm{\Theta }_f/d}/\sqrt{z}+\mathrm{}$$
(7)
with $`\mathrm{\Theta }_f/d=1/2`$ from what we saw previously. The SK limit corresponds to setting $`z=N`$, so we see that smooth large $`N`$ and $`z`$ limits require $`\sigma _0(z)z^\mu `$ where $`\mu `$ is equal to the SK value of $`\mathrm{\Theta }_f/d`$. We saw that it was difficult to obtain $`\mathrm{\Theta }_f/d`$ in the fixed connectivity models, but obtaining the prefactor $`\sigma _0(z)`$ is even more difficult. Nevertheless, we have extrapolated our data for the root mean square deviation divided by $`N^{0.5}`$ in the different finite connectivity models. This leads to the estimates $`\sigma _0(z=3)0.67`$, $`\sigma _0(z=4)0.70`$, $`\sigma _0(z=6)0.73`$, and $`\sigma _0(z=10)0.76`$. This growth is very slow; it is compatible with the value $`\mu =1/6`$ but much less with the value $`\mu =0.25`$. Given the large uncertainties in our extrapolations however, this has to be considered as only very indirect evidence in favor of $`\mathrm{\Theta }_f/d=1/6`$.

We conclude that understanding the non-trivial size dependence of the fluctuations of the ground state energy in the SK model remains quite a challenge. Let us propose here a simple argument suggesting that $`\mathrm{\Theta }_f/d=1/4`$. Imagine that one changes slowly all the couplings $`J_{ij}=J_{ij}^0/\sqrt{N}`$: $`J_{ij}^0J_{ij}^0+\delta J_{ij}^0`$, where the order of magnitude of $`\delta J_{ij}^0`$ is denoted by $`\epsilon `$. If $`\epsilon `$ is infinitesimal, the ground state remains the same. Up to what value of $`\epsilon `$ will this be true? The change of local field induced by the change of couplings for a fixed configuration $`\{S_i^{}\}`$ of the spins is of order $`_j\delta J_{ij}^0S_j^{}/\sqrt{N}\epsilon `$. Using the fact that the local field distribution vanishes linearly for small fields in the SK model, it is easy to show that the smallest local field is of $`O(1/\sqrt{N})`$. Therefore, the first value of $`\epsilon `$ that will trigger a change of the ground state is $`\epsilon ^{}1/\sqrt{N}`$. Flipping the spin with this low local field should lead to a cascade of flips that lowers the energy and thus to a new ground state that is “substantially” different from the starting one. Furthermore, the variation in the ground state energy when going from $`\epsilon =0`$ to $`\epsilon =\epsilon ^{}`$ is probably $`O(1)`$; indeed, there are excited states of energy $`O(1)`$ above the ground state, and one therefore expects level crossings to occur when $`\epsilon 1/\sqrt{N}`$. Finally, in order to scan the whole distribution of $`J_{ij}^0`$’s, one needs $`\epsilon 1`$. Using this range, the ground state will change $`1/\epsilon ^{}N^{0.5}`$ times. Since between each level crossing, the ground state energy randomly changes by an amount $`O(1)`$, the total expected fluctuation of the ground state energy will be $`O(N^{1/4})`$, so that $`\mathrm{\Theta }_f/d=1/4`$. More analytical work is obviously needed to confirm this speculative result, but note that at the heart of our argument lies the fact that the ground state of the SK model is particularly fragile: a relative change of order $`1/\sqrt{N}`$ of the disorder is enough to substantially change the ground state.

### IV.3 Edwards-Anderson models

In the case of the finite dimensional lattices, we know that $`\mathrm{\Theta }_f=d/2`$ holds exactly because of the Wehr-Aizenman theorem WehrAizenman90 , and recently Aspelmeier and Moore AspelmeierMoore03a  found that this relation holds within replica theory. In spite of these theoretical results, it is instructive to see how this equality transpires numerically. We thus follow the procedure used in the fixed connectivity models where we tested for when the rescaled RMS became size independent. In $`d=3`$, the rescaled data show no obvious trend when $`1.49\mathrm{\Theta }_f1.60`$ while in $`d=2`$ the corresponding range is $`1.00\mathrm{\Theta }_f1.02`$. In particular, in Fig. 3 we show these ratios when $`\mathrm{\Theta }_f`$ is set to $`d/2`$ (the data displayed are from the models with Gaussian $`J_{ij}`$s.)

In principle, it would be interesting to find the corrections to this scaling law. In a renormalization group picture, corrections go as inverse powers of $`L`$. Furthermore, if one believes that the exponent $`\mathrm{\Theta }_s`$ gives the leading corrections to scaling for all extensive quantities, then those for the variance of $`E_J`$ should go as $`L^{\mathrm{\Theta }_s}`$. To consider this possibility, we set $`\mathrm{\Theta }_f=d/2`$ and then ask when our data for the rescaled RMS follow approximately a straight line when plotted as a function of $`L^\omega `$ with the expectation that $`\omega =d\mathrm{\Theta }_s`$. In $`d=3`$, we find $`\omega 1.5`$, to be compared with the “theoretical” value $`\omega 2.8`$ (cf. section VI). Since $`\mathrm{\Theta }_s`$ is small and may in fact be zero, we have performed the analysis for $`\omega =3`$; using the intercept of the curve with the $`y`$ axis, we find that the RMS of $`E_J`$ grows as $`0.765L^{3/2}`$. In $`d=2`$, the rescaled data is too flat and so in practice we cannot give any sensible estimate of $`\omega `$. But we can follow the $`d=3`$ procedure, setting this time $`\omega =2`$ ($`\mathrm{\Theta }_s`$ is also relatively close to zero in $`d=2`$); the corresponding fits give that that the RMS of $`E_J`$ grows as $`0.725L`$. The moral of this story is that corrections to scaling are in general very difficult to determine, even if the leading scaling law is known exactly. Nevertheless, it seems that the same $`\omega `$ may very well describe the dominant corrections to scaling of many observables, as expected from the renormalization group picture.

To put these last numbers in perspective, consider the Mattis model where the couplings are gauge transformed from a ferromagnet having $`J_{ij}`$ chosen randomly on the positive side of a Gaussian. The total ground state energy of such a system is the sum of all these couplings and thus has the expectation value $`dN\sqrt{2/\pi }`$ on the hypercubic lattice of dimension $`d`$. On the other hand its variance is $`dN[12/\pi ]`$. Thus in this Mattis model, the RMS of the ground state energy grows as $`1.04L^{3/2}`$ in $`d=3`$, and as $`0.852L`$ in $`d=2`$. As expected, the Mattis model has larger absolute fluctuations than the EA model. It is also appropriate to compare the relative fluctuations $`\sigma _r`$, that is the RMS of the ground state energy divided by its mean. For the Mattis model, we find $`\sigma _r=0.434/L^{3/2}`$ in $`d=3`$ and $`\sigma _r=0.534/L`$ in $`d=2`$. These should be compared to the values we find in the EA model: $`\sigma _r=0.450/L^{3/2}`$ in $`d=3`$ and $`\sigma _r=0.551/L`$ in $`d=2`$. This shows that the relative fluctuations are slightly smaller in the Mattis model than in the EA model. Although this is in line with what frustration should do, note that the size of the effect is about $`3\%`$ which is very very small.

To summarize our study of $`\mathrm{\Theta }_f`$ for the different spin glass models, we have found that all the models considered seem to satisfy $`\mathrm{\Theta }_f=d/2`$. The notable exception is the SK model for which we made the case that $`\mathrm{\Theta }_f=d/4`$ was likely but $`\mathrm{\Theta }_f=d/6`$ was also possible; no matter what, $`\mathrm{\Theta }_f=d/2`$ is excluded, showing that the $`N\mathrm{}`$ and $`d\mathrm{}`$ limits do not commute. In all other cases, the variance of the ground state energy grows linearly with the system’s volume; this is the scaling expected from a central limit theorem behavior when the different terms contributing to the ground state energy are independent. Our results thus tell us that these terms are only weakly correlated. This feature is completely missed by both the REM and the DPRM; although a one-parameter scaling picture applies to those two models, it cannot apply to spin glasses. Note that this central limit behavior suggests that $`P_N(E_J)`$ tends to a Gaussian; we now turn to see whether this is the case.

## V Probability distribution of ground state energies

If the central limit theorem (CTL) were applicable, not only would we have $`\mathrm{\Theta }_f=d/2`$ for the scaling of the fluctuations, but also the shape of $`P_N(E_J)`$ would become Gaussian at large $`L`$. This behavior indeed arises for the MK lattices, both analytically at large $`b`$ and numerically for all $`b`$. (In our numerical study, we find that the skewness and kurtosis $`P_N(E_J)`$ decrease fast towards zero as $`L`$ grows.) Obviously, the terms contributing to the ground state energy are not independent but their correlations are not strong enough to prevent a CLT large $`L`$ scaling.

The question we address here is whether this simple behavior also holds in the other models. Let us begin with the mean field case. For the fixed connectivity mean field graphs, our data for the skewness and kurtosis of $`P_N(E_J)`$ decrease in magnitude as $`N`$ increases; this decrease is compatible with an extrapolation to zero as $`N\mathrm{}`$ as illustrated in Fig. 4 (These quantities are difficult to measure to high precision, so this should be considered as only suggestive of a Gaussian limit for $`P_N(E_J)`$.) The SK model however is clearly in a different class: not only does it have $`\mathrm{\Theta }_fd/2`$ but also its $`P_N(E_J)`$ is not Gaussian. In particular, its skewness at large $`N`$ stabilizes around $`0.43\pm 0.02`$, while its kurtosis stabilizes around $`0.40\pm 0.03`$. It is instructive to compare this to the values predicted by the REM model where $`P(X_J)`$ is a Gumbel distribution: there, the skewness is $`1.139`$ while the kurtosis is exactly $`2.4`$. Our estimates do not agree either with the values from the Fisher or Weibull universality classes: we thus conclude that the SK ground state energy distribution does not belong to one of the known universality classes of extreme statistics BouchaudMezard97 .

Finally we come to the EA lattices. Given that the MK and diluted mean-field graphs lead to the same conclusion, it will come as no surprise that our data for the EA lattices are also compatible with a Gaussian limiting shape for $`P_N(E_J)`$. Note that this is expected though not proven from the work of Wehr and Aizenman WehrAizenman90 , while it does follow from replica theory calculations AspelmeierMoore03a . An examination of the skewness and kurtosis of $`P_N(E_J)`$ shows that they decay with system size, both in $`d=2`$ and $`d=3`$. Although our measurements lack precision when the number of spins is large, the extrapolations suggest zero limiting values as $`L\mathrm{}`$ as one can see in Fig. 4. Perhaps it is also worth noting that our data are quantitatively similar in the different models; for instance, for $`L=8`$ and $`d=3`$, we find the skewness to be $`0.18`$ in EA and $`0.21`$ in MK.

The overall situation indicates that $`P_N(E_J)`$ is Gaussian. Furthermore we checked whether the convergence to this Gaussian follows the central limit theorem law. Indeed, that law predicts for instance that the skewness scales as $`N^{1/2}`$; thus we have plotted in Fig. 4 the skewness for the different models as a function of that scale. The data is completely compatible with a linear convergence to zero, confirming the CLT scaling.

Only in the case of the SK model are the variables contributing to the ground state energy sufficiently correlated to prevent a central limit theorem behavior. A physical interpretation of this is that as soon as there is a kind of locality which allows one to decompose the sample into quasi-independent subsystems, the central limit theorem behavior will appear, leading to $`\mathrm{\Theta }_f=d/2`$ and a Gaussian $`P_N(E_J)`$. Of course, it is not clear why this should apply to the diluted mean field graphs.

## VI The shift exponent $`\mathrm{\Theta }_s`$

### VI.1 Migdal-Kadanoff lattices

We now move on and study the finite size corrections to the mean energy density. Following Eq. 1, the mean excess of the (extensive) ground state energy is expected to scale as $`L^{\mathrm{\Theta }_s}`$. To have an idea of what this exponent should be, it is again most useful to begin with the hierarchical lattices. The important prediction of that approach is that $`\mathrm{\Theta }_s=\theta _{DW}`$ where $`\theta _{DW}`$ is the domain wall exponent. To see why this is so, reconsider the evolution equation for the energies $`E_p`$ and $`E_a`$ as one applies the recursion. First, along a given path, the energy is the sum of the ground state energies of each of its $`\mathrm{}`$ segments, the sum being sometimes corrected by the energy of the domain wall of the weakest segment in order to satisfy the boundary conditions. This correction shifts the path’s energy by $`O(L^{\theta _{DW}})`$. Second, adding the energies of the different paths does not change the scale of the shifts, so necessarily $`\mathrm{\Theta }_s=\theta _{DW}`$. Naturally, we have confirmed this relation numerically for different values of $`\mathrm{}`$ and $`b`$ by performing fits to moments of $`P(E_p,E_a)`$ (determined with no statistical error for these hierarchical lattices) from which we extracted estimates for $`\mathrm{\Theta }_s`$ and $`\theta _{DW}`$.

When $`b`$ is large, we can derive the analytical expression for $`\theta _{DW}`$ on these MK lattices. Indeed, in this limit, $`P(E_p,E_a)`$ becomes Gaussian and so one can perform a cumulant expansion about this Gaussian, leading to a $`1/b`$ series. This scheme extends the work of Southern and Young SouthernYoung77  who assumed that $`P(E_p,E_a)`$ was Gaussian even for finite $`b`$. For general $`\mathrm{}`$ and $`b`$ we obtain

$$\theta _{DW}(\mathrm{},b)=\frac{d}{2}+a_0(\mathrm{})+\frac{a_1(\mathrm{})}{b}+O(\frac{1}{b^2}).$$
(8)
When $`\mathrm{}=2`$, we find $`a_0=1.2302`$ and $`a_1=0.04573`$ (see the appendix for a derivation). In Fig. 5 we plot the difference between the numerically obtained $`\theta _{DW}(\mathrm{}=2,b)`$ and the terms of the expansion given in Eq. 8. This allows us to determine numerically the next term of the expansion, and we find $`0.045\pm 0.001/b^2`$. This value could be obtained analytically, but we have not pushed the analytical calculation to that order.

Note that the $`1/b`$ expansion corresponds to an expansion in $`e^d`$ where $`d`$ is the dimension of space. This justifies the fact that the $`1/b`$ expansion is quite accurate all the way down to $`d=3`$ (which corresponds to $`b=4`$). Finally, when $`\mathrm{}\mathrm{}`$, we obtain $`a_0=3/2`$. These results show that $`\mathrm{\Theta }_s=\theta _{DW}<\mathrm{\Theta }_f`$, justifying the neglect of the “correction” terms in section IV.1 from which we concluded that $`\mathrm{\Theta }_f=d/2`$.

Before going on to the mean field case, let us remark that the MK value for $`\theta _{DW}`$ is quite close the actual value in the EA model. If we use the standard choice for $`d=3`$, $`\mathrm{}=2`$ with $`b=4`$, the MK prediction is $`\theta _{DW}0.255`$. One can also use the choice $`\mathrm{}=3`$ with $`b=9`$ for which $`\theta _{DW}0.242`$. These values are to be compared to current estimates for $`\theta _{DW}`$ in the EA model, $`\theta _{DW}=0.21\pm 0.02`$ PalassiniYoung99b  and $`\theta _{DW}=0.19\pm 0.02`$ Hartmann99 . A similarly good comparison occurs when $`d=2`$.

### VI.2 Mean-field models

We have seen that the prediction of the MK approach is $`\mathrm{\Theta }_s=\theta _{DW}`$; in $`d=3`$, this gives either $`\mathrm{\Theta }_s=0.25`$ (the MK value) or $`\mathrm{\Theta }_s0.20`$ (if one uses the $`d=3`$ EA model values for $`\theta _{DW}`$). In the case of the mean field models, there is no way to introduce domain walls, and so we simply focus on their prediction for $`\mathrm{\Theta }_s`$.

Consider first the finite size effects in the SK model. Parisi et al. ParisiRitort93b  have computed analytically how various quantities converge to their large $`N`$ limit. An $`N^{2/3}`$ convergence is the general rule, though for the energy density they were able to compute the finite size correction only at the critical temperature $`T_c`$ and on the de Almeida-Thouless line. Nevertheless, the natural extrapolation is that this law should apply to all $`TT_c`$, leading to the prediction $`\mathrm{\Theta }_s/d=1/3`$. (Note that this prediction is very different from that of CrisantiPaladin92  for which $`\mathrm{\Theta }_s/d=\mathrm{\Theta }_f/d=1/6`$.) To our knowledge, the possibility that $`\mathrm{\Theta }_s/d=1/3`$ for the SK ground state energy was first brought up by Palassini PalassiniThesis  in his numerical studies. We can extend his analysis with our data; identifying as before $`L`$ with $`N^{1/d}`$, we perform fits of the mean ground state energy to Eq. 1. We then find $`\mathrm{\Theta }_s/d=0.34\pm 0.02`$, in complete compatibility with $`1/3`$; note that when plotted as a function of $`N^{2/3}`$, the data is very linear starting from $`N=50`$ (see the inset of Fig. 6). In addition, if we perform the fit while forcing $`\mathrm{\Theta }_s/d=1/3`$, we find that the ground state energy density at $`N=\mathrm{}`$ is $`e_0=0.7637\pm 0.0002`$, in very good agreement with the exact MezardParisi87b  result $`e_0=0.7633\mathrm{}`$

In contrast to what happens for $`\mathrm{\Theta }_f/d`$, $`\mathrm{\Theta }_s/d`$ is the same in the SK model and in the fixed connectivity models we have considered; this is illustrated in Fig. 6. For instance, for connectivity $`6`$, a power law fit gives $`\mathrm{\Theta }_s/d=0.35\pm 0.03`$ with $`\chi ^2=9.7`$ for $`7`$ degrees of freedom. In fact, the value $`\mathrm{\Theta }_s/d=1/3`$ works very well for all four connectivities we studied, and we are tempted to consider that this value is the exact exponent. The same conclusion was reached by Boettcher Boettcher03a .

### VI.3 Edwards-Anderson models

Now it is time to compare the MK and mean field “predictions” to our measurements of $`\mathrm{\Theta }_s`$ in the EA models. Let us begin with the case $`d=2`$. We fit our $`4L10`$ data to Eq. 1 where $`e_0`$, $`e_1`$ and $`\mathrm{\Theta }_s`$ are ajustable parameters; that gives us the estimate $`\mathrm{\Theta }_s=0.35\pm 0.01`$. (The associated $`\chi ^2`$ is $`1.9`$ for $`5`$ degrees of freedom.) The resulting fit is displayed along with the data in the inset of Fig. 7. Given our statistical error and perhaps more importantly the systematic effects associated with $`L`$ not being very large, this value is reasonably compatible with the MK prediction $`\mathrm{\Theta }_s=\theta _{DW}`$ since $`\theta _{DW}=0.28`$ RiegerSanten96 ; HartmannYoung01 ; CarterBray02 .

A further argument in favor of $`\mathrm{\Theta }_s=\theta _{DW}`$ is as follows. It is known that whether one modifies the boundary conditions PalassiniYoung99a  or increases the size of the system Middleton99 , the fractal dimension of the surface of the clusters of spins that change is approximately the same. It is thus likely that the same type of excitations are involved in determining $`\mathrm{\Theta }_s`$ and $`\theta _{DW}`$. As a consequence, we believe that the $`d=2`$ EA model is described by a scaling theory with $`\theta _s=\theta _{DW}`$, but also with $`\mathrm{\Theta }_f\theta _{DW}`$ of course. Note that the mean field prediction ($`\mathrm{\Theta }_s=2/3`$) is clearly off; however, one cannot appeal to mean field when $`d=2`$ because one is below the lower critical dimension.

Let us now move on to the $`d=3`$ EA model which is more challenging and has a spin glass transition at $`T_c>0`$. We use the same fitting function (Eq. 1) as before; the best fit then gives a good $`\chi ^2`$ and a mean ground state energy growing as $`1.700L^3+1.9L^{\mathrm{\Theta }_s}`$ with $`\mathrm{\Theta }_s=0.10\pm 0.03`$. This fit is displayed in the main part of Fig. 7. However this value of $`\mathrm{\Theta }_s`$ easily changes by $`0.1`$ when removing some of the data points, and in fact the fit sometimes even leads to negative values for $`\mathrm{\Theta }_s`$. Thus at best we can say that $`\mathrm{\Theta }_s`$ is small, somewhere between $`0.0`$ and $`0.2`$. One can compare this result to “theory”. The mean field value ($`\mathrm{\Theta }_s=1.0`$) is completely ruled out, whereas the MK prediction $`\mathrm{\Theta }_s=\theta _{DW}`$ is quite acceptable since $`\theta _{DW}0.2`$ PalassiniYoung99b ; Hartmann99 . However, another possibility is that the discrepancy we find has a physical origin and that in fact $`\mathrm{\Theta }_s\theta _{DW}`$. Since this issue is important, we push the analysis a bit further as follows. Given a putative value for $`\mathrm{\Theta }_s`$, we adjust $`e_0`$ so that the plot of $`\overline{E_J(L)}e_0L^3`$ versus $`L^{\mathrm{\Theta }_s}`$ passes through the origin. For $`\mathrm{\Theta }_s`$ outside the range $`[0.05,0.2]`$, the data has visible curvature. In the inset of Fig. 8 we show $`\overline{E_J(L)}e_0L^3`$ as a function of $`L`$; we see that there is no clear trend, so even $`\mathrm{\Theta }_s=0`$ seems possible. Such a value could be interpreted from the (mean-field-like) behavior of system-size excitations found in this model KrzakalaMartin00 ; PalassiniYoung00a . The analogous analysis in the $`d=2`$ case is shown in the main part of Fig. 8 when we use the value $`\mathrm{\Theta }_s=0.28`$; there the expected corrections to scaling work quite well.

Of course it would be useful to have a significantly smaller error bar for $`\mathrm{\Theta }_s`$ but we cannot go much beyond what we have done here: the statistical error on $`\overline{E_J}`$ grows as $`L^{d/2}`$, and the amount of computation time grows still faster, so we cannot obtain useful information at large $`L`$.

We have no data for $`d=4`$, but let us briefly consider the published work by Hartmann Hartmann99d  where he used the $`J_{ij}=\pm 1`$ EA model. (One expects its $`\mathrm{\Theta }_s`$ to be the same as in the Gaussian case.) Analyzing his values for the mean ground state energy for $`2L7`$, we find $`\mathrm{\Theta }_s=0.2\pm 0.1`$; this does not compare well with his estimate $`\theta _{DW}=0.65\pm 0.04`$. But, if we remove the $`L=2`$ point from the fit, we find $`\mathrm{\Theta }_s=0.7\pm 0.2`$ which is in good agreement with $`\mathrm{\Theta }_s=\theta _{DW}`$.

We have also analyzed the data of Boettcher and Percus BoettcherPercus01 , and this leads to the same conclusion. In summary, we cannot exclude that $`\mathrm{\Theta }_s<\theta _{DW}`$, but the MK prediction $`\mathrm{\Theta }_s=\theta _{DW}`$ works surprizingly well in the finite dimensional EA models. On the contrary, the mean field prediction is definitely off, and that of the Mattis spin glass is completely wrong since it gives $`\mathrm{\Theta }_s=\mathrm{}`$.

## VII Case of $`+/J`$ couplings

It is widely believed that exponents are universal, i.e., independent of the detailed microscopic nature of the disorder. (Note however that there are longstanding claims of universality violations in spin glasses MariCampbell99 .) In particular, for the $`d=3`$ EA, numerical computations of $`\theta _{DW}`$ confirm this to a large extent: one has $`\theta _{DW}=0.21\pm 0.02`$ for the Gaussian PalassiniYoung99b  and $`\theta _{DW}=0.19\pm 0.02`$ for the $`J_{ij}=\pm 1`$ Hartmann99  models. However, if $`T_c=0`$ as arises in $`d=2`$, one may expect several universality classes and thus some influence of the microscopic properties (i.e., the distribution of the $`J_{ij}`$) on the macroscopic properties (e.g., exponents). Note first that $`\theta _{DW}`$ is known to be different for Gaussian and binary ($`J_{ij}=\pm 1`$) couplings HartmannYoung01  in the $`d=2`$ EA model. A posteriori, that is not so surprising since the associated quantities (for instance domain wall energies) go to zero rather than to infinity; they can thus easily be affected by microscopic details. This issue can be investigated within the framework of MK lattices. We find identical values of $`\theta _{DW}`$ in the Gaussian and the $`J=\pm 1`$ cases whenever $`\theta _{DW}>0`$. However in $`d=2`$ and using $`b=\mathrm{}=2`$, the Gaussian model gives $`\theta _{DW}=0.22`$; on the contrary, the $`J_{ij}=\pm 1`$ case leads to $`\theta _{DW}=\mathrm{}`$, meaning that the domain-wall energies decrease exponentially with the size of the system rather than as a power law. Thus either there is violation of universality (which seems unlikely to us) or there are several universality classes when $`\theta _{DW}0`$; Amoruso et al. AmorusoMarinari03  have given evidence in favor of the latter possibility since the first version of this paper was posted.

Given these amendments to the scope of universality, the microscopic details should indeed be irrelevant for quantities associated with diverging energy scales. The surprising claim we bring forward here is that this expectation is still too strong: a counter example is provided by $`\mathrm{\Theta }_f`$ in the mean field fixed connectivity graphs. Indeed, we saw in the case of Gaussian couplings that $`\mathrm{\Theta }_f/d=1/2`$. Now in Fig. 9 we show the rescaled standard deviation of $`E_J`$ for the case where $`J_{ij}=\pm 1`$ on fixed connectivity graphs; the standard deviations have been divided by $`N^{1/4}`$. If $`\mathrm{\Theta }_f/d=1/2`$, we should see a rapid divergence of the plotted values with increasing $`N`$, but instead the curves are relatively flat and decreasing!

In fact, when performing fits, we find that $`\mathrm{\Theta }_f/d`$ is between $`1/4`$ and $`1/5`$.

To give further evidence that the Gaussian and $`J_{ij}=\pm 1`$ cases scale differently, we plot in the inset of this figure the values of the skewness of the distribution of $`E_J`$ using the same models and values of $`N`$ as in the main figure. Although our data is noisy, the skewness shows no sign of going to zero when $`N\mathrm{}`$. Taken at face value, this means that the distribution of $`E_J`$ is not Gaussian in the large $`N`$ limit when $`J_{ij}=\pm 1`$, in sharp contrast to what happens in the case of Gaussian couplings. To drive this point home further, we find that as $`z`$ increases, the skewness and kurtosis grow in magnitude and seem to approach the values we find in the SK. (But as before, our values become imprecise at large $`N`$.) We thus conjecture that each $`z`$ is associated with a different universality class and that, as $`z\mathrm{}`$, one converges to the class to which the SK model belongs. One way to justify this is to consider that in these models $`z`$ is related to dimension rather than to lattice connectivity; the universality class will then change with $`z`$. To a large extent, all different trends give further credence to the claim that $`\mathrm{\Theta }_f/d=1/4`$ in the SK model if that is the correct value for the fixed connectivity models with binary couplings.

How can one understand this “break-down” of universality when the underlying energies diverge? Recall that for the $`J_{ij}=\pm 1`$ spin glass model on a fixed connectivity graph, the local environment of a spin has no disorder out to finite distances: any sign of the $`J_{ij}`$ can be gauged away so that all the sample-dependent fluctuations arise “at infinity”. On the contrary, in the Gaussian case the $`J_{ij}`$ fluctuations are local, leading to $`O(N^{1/2})`$ fluctuations in the total energy. Taken at face value, this argument also applies to coordination fluctuations; if this is true, the $`\pm 1`$ Viana-Bray model VianaBray85  will have $`\mathrm{\Theta }_f=d/2`$.

To summarize, the exponent $`\mathrm{\Theta }_f`$ depends on the details of the underlying $`J_{ij}`$ distribution even though the energy scale (of fluctuations) diverges. Not surprisingly, we also find that $`P_N(E_J)`$ for that system does not become Gaussian. Finally, in spite of this major change of behavior when going from Gaussian $`J_{ij}`$ to binary values, we find that $`\mathrm{\Theta }_s/d=1/3`$ very precisely in both cases.

## VIII Discussion and conclusion

Let us go over the main results of the present work on spin glasses. First, the sample to sample fluctuations of the ground state energy grow as the square-root of the volume in almost all models so the fluctuation exponent satisfies $`\mathrm{\Theta }_f=d/2`$. Furthermore, the distribution over disorder of the ground state energy probably tends towards a Gaussian in the large volume limit as suggested by Wehr and Aizenman WehrAizenman90 . There are two notable exceptions to this picture: the SK model and the fixed connectivity mean field models having $`J_{ij}=\pm 1`$. In those two cases, the sample to sample fluctuations are much smaller, and we find $`\mathrm{\Theta }_fd/4`$, even though a still smaller value cannot be excluded. (Note that fluctuations which are smaller than square root of the volume $`N`$ are also a characteristic of the directed polymer in a random medium.) On the other hand, for finite dimensional lattices, one expects the finite density of unfrustrated regions to contribute $`O(N^{1/2})`$ fluctuations to the ground state energy. A trivial example where this is the case is the (unfrustrated) Mattis model; there the variance of the ground state energy is obviously maximal and equal to $`zN\left[J^2J^2\right]/2`$ for a connectivity $`z`$.

We also studied the exponent $`\mathrm{\Theta }_s`$ giving the corrections to the scaling of the average (extensive) ground state energy. For the $`d=2`$ and $`3`$ EA models, we find that the equality $`\mathrm{\Theta }_s=\theta _{DW}`$ holds within our limited precision. This means that corrections to scaling are associated with domain walls. However, we were not able to rule out $`\mathrm{\Theta }_s=0`$ in $`d=3`$, leaving the door open to other interpretations. In the context of the MK lattices, the equality $`\mathrm{\Theta }_s=\theta _{DW}`$ holds exactly. Physically, this equality corresponds to the fact that in a finite fraction of the samples, the boundary conditions force a domain wall “defect” to be present in the ground state. (Note that the finite size corrections are indeed always positive for fixed and periodic boundary conditions). Furthermore, we were able to compute $`\theta _{DW}`$ analytically to order $`1/b`$ when the dimension becomes large. In the case of the mean field models (SK or diluted graphs), the situation is quite different: the corrections to scaling grow as $`N^{1/3}`$. Since one expects $`\theta _{DW}d/2`$ in large dimensions, the result $`\mathrm{\Theta }_s/d=1/3`$ for these mean field models differs from the large dimension limit of $`\theta _{DW}/d`$.

Finally, we have exhibited examples where the exponent $`\mathrm{\Theta }_f`$ depends on the distribution of the $`J_{ij}`$, even though energy fluctuations diverge with size, i.e., $`\mathrm{\Theta }_f>0`$.

Our work suggests several paths for further studies. (a) Can one establish the value of $`\mathrm{\Theta }_f/d`$ for the SK model? (b) Since large fluctuations of order $`N^{1/2}`$ are detrimental in a numerical determination of the average ground state energy, is there a way to substract a (computable) contribution of satisfied bonds so as to reduce the variance? (c) Is $`\mathrm{\Theta }_s=\theta _{DW}`$ in finite dimensional spin glasses or is $`\mathrm{\Theta }_s`$ smaller? This second possibility could follow from other types of excitations whose exponents are smaller than $`\theta _{DW}`$.

#### Acknowledgments —

We thank A. Billoire, K. Binder, D. S. Fisher, M. Mézard, M. Palassini, G. Parisi, R. da Silveira and O. White for stimulating discussions. This work was supported in part by the European Community’s Human Potential Programme under the contracts HPRN-CT-2002-00307 and HPRN-CT-2002-00319. F.K. acknowledges financial support from the MENRT. The LPTMS is an Unité de Recherche de l’Université Paris XI associée au CNRS.

## Appendix: Exponents in the Migdal-Kadanoff lattices and the large dimensional limit

In this appendix, we focus on the simplest case $`\mathrm{}=2`$. To write a recursion relation, one needs to keep track of two energies, $`E_p^{(r)}`$ and $`E_a^{(r)}`$ that give the ground state energy of the MK lattice at the $`r`$-th application of the recursion when the two exterior spins are respectively parallel and antiparallel. For a $`b`$-branch lattice with $`\mathrm{}=2`$, the ground state energies at the $`(r+1)`$-th application of the recursion read:

$`E_p^{(r+1)}={\displaystyle \underset{\alpha =1}{\overset{b}{}}}\mathrm{min}`$
$`(E_p^{(r)}(1,\alpha )+E_p^{(r)}(2,\alpha ),E_a^{(r)}(1,\alpha )+E_a^{(r)}(2,\alpha ))`$ (9)
$`E_a^{(r+1)}={\displaystyle \underset{\alpha =1}{\overset{b}{}}}\mathrm{min}`$
$`(E_p^{(r)}(1,\alpha )+E_a^{(r)}(2,\alpha ),E_a^{(r)}(1,\alpha )+E_p^{(r)}(2,\alpha ))`$ (10)
where the index $`1,2`$ refers to the two bonds in the $`\mathrm{}`$ direction. This equation says that along all the $`b`$ branches, one has to choose the orientation of the middle spin (that is decimated) such that the energy contribution is minimized, given the orientations of the external spins. Note that all the random variables that appear in this equation are independent as soon as they live on different bonds.

If one assumes that the distribution $`P^{(r)}(E_p,E_a)`$ of energies at the $`r`$-th generation can be written for large $`r`$ in the scaling form

$$P^{(r)}(E_p,E_a)=1/\sigma _r^2f(\left[E_p_r\right]/\sigma _r,\left[E_a_r\right]/\sigma _r)$$
(11)
then it is immediate to show (using the independence of the branches) that $`_{r+1}=(2b)_r`$ and $`\sigma _{r+1}=\sqrt{2b}\sigma _r`$, so that $`_r(2b)^r`$ and $`\sigma _r(2b)^{r/2}`$. Since the number of spins is given by $`(2b)^r`$, one immediately finds that the ground state energy is extensive and that the fluctuations are described by $`\mathrm{\Theta }_f/d=1/2`$.

Define $`\mathrm{\Delta }^{(r)}`$ as $`E_p^{(r)}E_a^{(r)}`$. It is easy to show that $`\mathrm{\Delta }^{(r)}`$ obeys an autonomous recursion relation SouthernYoung77 :

$$\mathrm{\Delta }^{(r+1)}=\underset{\alpha =1}{\overset{b}{}}ϵ(\alpha )\mathrm{min}(|\mathrm{\Delta }^{(r)}(1,\alpha )|,|\mathrm{\Delta }^{(r)}(2,\alpha )|),$$
(12)
where $`\alpha `$ labels branches and

$$ϵ(\alpha )=\mathrm{sign}\left[\mathrm{\Delta }(1,\alpha )\right]\mathrm{sign}\left[\mathrm{\Delta }(2,\alpha )\right].$$
(13)
Using the independence of the $`\mathrm{\Delta }`$’s corresponding to different branches, one finds that $`\mathrm{\Sigma }_r`$, the RMS of the distribution of the $`\mathrm{\Delta }`$’s, obeys the following equation:

$$\mathrm{\Sigma }_{r+1}^2=4b\mathrm{\Sigma }_r^2_0^{\mathrm{}}\text{d}xg(x)\left[_0^x\text{d}yy^2g(y)+x^2_x^{\mathrm{}}\text{d}yg(y)\right],$$
(14)
where $`g(.)`$ is the distribution of $`\mathrm{\Delta }/\mathrm{\Sigma }`$ that is independent of $`r`$ for large $`r`$. From this relation, one finds that $`\mathrm{\Sigma }_r\lambda ^r`$ with

$$\lambda ^2=4b_0^{\mathrm{}}\text{d}xg(x)\left[_0^x\text{d}yy^2g(y)+x^2_x^{\mathrm{}}\text{d}yg(y)\right]$$
(15)
where of course $`g(.)`$ depends on $`b`$. The energy scale for flipping the relative sign of the exterior spins is $`\mathrm{\Sigma }`$ and so the exponent $`\theta _{DW}`$ is given by

$$\theta _{DW}=\frac{\mathrm{ln}\lambda }{\mathrm{ln}2}.$$
(16)
In the large dimension limit, for which $`b\mathrm{}`$, it is clear that, using the central limit theorem, the distribution of $`\mathrm{\Delta }`$ is Gaussian. Since $`g(.)`$ is then known, $`\lambda `$ can be computed from Eq. 15, giving

$$\lambda ^2=0.36338b,$$
(17)
a result first obtained by Southern and Young.

When $`b`$ is large but not infinite, the first correction to the Gaussian is of order $`1/b`$ (because the $`\mathrm{\Delta }`$ have a symmetric distribution) and reads:

$$g(x)=\frac{1}{\sqrt{2\pi }}\left[1+\frac{\kappa }{24b}\frac{^4}{x^4}+\mathrm{}\right]e^{x^2/2}$$
(18)
where $`\kappa `$ is the kurtosis of the initial variable, i.e., $`ϵ\mathrm{min}(|\mathrm{\Delta }(1)|,|\mathrm{\Delta }(2)|)`$. To first order in $`1/b`$, this kurtosis can be computed by assuming that the $`\mathrm{\Delta }`$ are Gaussian, and one finds $`\kappa =0.434215`$. Injecting the expression of $`g(x)`$ in Eq. (15) then gives $`\lambda `$ to order $`1/b`$: $`\lambda ^2=0.36338b0.023035`$, and finally the result given in the main text after Eq. (8).

This calculation can be extended to next order: $`\kappa `$ will acquire a $`1/b`$ contribution and there will be corrections to $`g(x)`$ of order $`1/b^2`$ coming from the non zero sixth cumulant of $`ϵ\mathrm{min}(|\mathrm{\Delta }(1)|,|\mathrm{\Delta }(2)|)`$, computed as if the $`\mathrm{\Delta }`$’s were Gaussian. One could also, with more work, compute $`\lambda `$ for $`\mathrm{}2`$. In the limit $`\mathrm{}\mathrm{}`$, the problem becomes soluble again using the theory of extreme value statistics for handling the variable $`\mathrm{min}(|\mathrm{\Delta }(1)|,|\mathrm{\Delta }(2)|,\mathrm{},|\mathrm{\Delta }(\mathrm{})|)`$.


text:
# Vortex liquids and vortex quantum Hall states in trapped rotating Bose gases

## 1 Introduction

There has been enormous interest in the last few years in the properties of trapped Bose gases set under rotation, both theoretically and experimentally (a small selection of recent contributions is ). Of particular interest is the behavior of (effectively) two-dimensional gases in highly correlated fractional quantum Hall states (FQHE) . These states lend themselves for an exploration of physics beyond the mean-field picture, where the role of the strong magnetic field in the FQHE of ultrapure 2D electron gases is taken over by a fast rotation of the gas.

The formation of FQHE states in a rapidly rotating, harmonically trapped dilute Bose gas, such that the Landau level mixing due to interactions can be neglected, was studied in Ref. . It was found that, after the vortex lattice melts, at a critical filling factor of the rotational Landau levels the system of particles enters a highly correlated quantum ground state. The many-body wave function corresponding to this state was numerically shown to have large overlap with a bosonic version of the many-body quantum Hall ground state introduced by Laughlin for electrons .

Instead of a quantum Hall fluid of particles, one can equally well consider the vortices in the fluid to be the elementary objects, forming themselves a highly correlated quantum fluid. Previously, a quantum Hall fluid of vortices has been invoked to explain certain features of Hall conductance and magnetization experiments in high-$`T_c`$ superconductors . In the present paper, we propose to investigate vortex liquid phases in a rotating, spatially inhomogeneous, harmonically trapped dilute Bose gas. It is shown, in particular, that there exists a critical dimensionless 2D interaction strength, $`g=g_c0.6`$, above which several phases can coexist (we set $`\mathrm{}=m=1`$) . In the limit of weak interaction, for $`g<g_c`$, and after the vortex lattice has melted, the particles are confined to the lowest rotational Landau level and, for low enough temperatures, the system makes a transition into the particle Laughlin state. For $`g>g_c`$, on the other hand, a transition into a vortex quantum Hall state can occur. Upon partial lattice melting, the system separates into a low-density phase in which the vortex lattice has melted, and an inner high-density region, in which a vortex lattice resides. At a critical rotation velocity, the lattice melts completely, and the inner region consists of a vortex quantum liquid, while the outer shell enters a highly correlated particle quantum Hall state. Finally, at a second critical velocity, the vortex liquid phase disappears and the particle quantum Hall state fills the whole sample.

Below, we first introduce the basic vortex Hamiltonian of superfluid hydrodynamics in the rotating frame, upon which our discussion of the various vortex states will be built. We then give in section III the criterion for the melting of the vortex lattice into a vortex quantum liquid for a general rotating, uncharged superfluid. In section IV, this criterion is specialized to harmonically trapped Bose gases, where it is demonstrated how the phase separation described above is obtained from the inhomogeneous density profile of the gas cloud. In section V, we describe possible experimental means to verify the existence of the vortex quantum liquid. Finally, the conclusions section contains a discussion of how the zero temperature limitation, used in the following sections, can be lifted at practically all temperatures for which the vortex quantum liquid can be observed.

## 2 Vortex Hamiltonian and Plasma Analogy

Consider a two-dimensional (2D) Bose gas at temperatures well below the Kosterlitz-Thouless transition temperature $`T_c^{\mathrm{KT}}=\pi \rho /2`$, where $`\rho `$ is the 2D particle density. In this superfluid state, phase fluctuations are provided by phonons and bound vortex-antivortex pairs. Vorticity is a good quantum number, and the vortices are “fundamental” objects. At sufficiently low temperatures, we can neglect the dissipation which arises from scattering of thermal excitations on vortices, so that their dynamics is conservative. The situation considered then corresponds to the Magnus force dominated limit of vortex motion in very clean superconductors . The vortex motion is governed by a Hamiltonian corresponding to one of point particles, with a charge equal to the quantum of circulation $`\kappa =2\pi `$, interacting with electromagnetic fields. For a large collection of $`N_v1`$ vortices (cf. Fig. 1), in a frame rotating with angular velocity $`\mathrm{\Omega }`$, the vortex Hamiltonian reads :

$`H_v`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_v}{}}}{\displaystyle \frac{\left(𝑷_i\kappa 𝑨_i\right)^2}{2m_v}}\mathrm{\Omega }(𝑿_i\times 𝑷_i)_z`$ (1)
$`{\displaystyle \frac{\rho \kappa ^2}{2\pi }}{\displaystyle \underset{i<j}{\overset{N_v}{}}}\mathrm{ln}\left|{\displaystyle \frac{𝑿_i𝑿_j}{a_c}}\right|.`$
We emphasize that the hydrodynamic vortex Hamiltonian above does in its validity not depend on the (microscopic) form of the many-body wave function; it describes the system accurately if the matter we deal with is a superfluid. The Hamiltonian is also valid if the density varies like in the (harmonically) trapped case of currently experimentally realized Bose-Einstein condensates, as long as this variation takes place on scales much larger than the intervortex distance. This is the case here, because we will consider below large vortex densities, with intervortex separations of order a few times the coherence length.

The first term in (1) is the kinetic energy of vortices, the second stems from a $`𝛀𝑳`$ coupling, the Coriolis force on vortices, and the last term is the repulsive logarithmic Coulomb interaction of point vortices. In the symmetric gauge, the vector potential $`𝑨_i`$, responsible for the vortex velocity part of the Magnus force on the vortex (the analog of the Lorentz force), takes the form $`A_i^a=\frac{1}{2}\rho ϵ^{ab}X_b^i`$ ($`a,b=x,y)`$ . This vector potential corresponds to motion of a charged particle in a “magnetic field” perpendicular to the plane, of magnitude $`\rho `$. The canonical vortex momentum of the vortex $`i`$ is consequently given by the sum of the kinetic and field momenta,

$$𝑷_i=m_v(\dot{𝑿}_i+𝛀\times 𝑿_i)+\kappa 𝑨_i.$$
(2)
The vortex mass is due to the density disturbance a moving vortex causes in the system. For a single central vortex, it is given by $`m_v=[\rho \kappa ^2/(4\pi c_s^2)]\mathrm{ln}[R/a_c]=[\kappa ^2/4\pi g]\mathrm{ln}[R/a_c]`$ , where the core size $`a_c=O(\xi _0)`$ is of order the bulk coherence length $`\xi _0=1/\sqrt{2g\rho }`$, and $`c_s=\sqrt{g\rho }`$ the speed of sound . The interaction strength parameter $`g`$ is related to the compressibility $`\chi `$ of the fluid by $`\chi =1/g\rho ^2=1/\rho c_s^2`$. We will take below the limit of $`m_v0`$, so that the precise value of $`m_v`$, in particular its modification in the presence of a dense vortex lattice, is not of importance to us, but only that $`m_v`$ has a finite value. The latter fact is important for the quantization of the vortex position and momentum degrees of freedom in a straightforward canonical manner, because for $`m_v=0`$ exactly, vortex phase and configuration space coincide , that is, vortex momenta become (gauge dependent) functions of the vortex coordinates alone.

The ground state of the Hamiltonian (1) can be analyzed by setting all the vortex velocities to zero (in the rotating reference frame) and finding the vortex configuration by minimizing the energy:

$$H_v^{^{}}=\frac{1}{2}\rho \kappa \mathrm{\Omega }\underset{i=1}{\overset{N_v}{}}𝑿_i^2\frac{\rho \kappa ^2}{2\pi }\underset{i<j}{\overset{N_v}{}}\mathrm{ln}\left|\frac{𝑿_i𝑿_j}{a_c}\right|.$$
(3)
We neglected the centrifugal potential $`\frac{1}{2}m_v\mathrm{\Omega }^2_{i=1}^{N_v}𝑿_i^2`$ in the above expression, which is justified by taking the small mass limit $`m_v0`$ in (1). In a plasma analogy, the first term in the Hamiltonian $`H_v^{^{}}`$ provides the neutralizing background to the potential created by the vortices. In the rotating frame the system is identical to a system of charges interacting by 2D Coulomb repulsion in a neutralizing background (the total circulation being zero as seen in the rotating frame), such that the ground state of the system is given by the spatially homogeneous vortex density $`n_v=2\mathrm{\Omega }/\kappa =\mathrm{\Omega }/\pi `$.

## 3 Vortex Lattice Melting

The uniform vortex density solution is only a first approximation for the ground state properties of the rotating superfluid. Since the vortices repel each other, the energy of the liquid can be lowered further by forming highly correlated states. At relatively low vortex density, the vortices form a triangular lattice . The excitations above this vortex lattice ground state in a rotating dilute Bose-Einstein condensate have been the subject of recent experimental and theoretical studies . They represent Tkachenko waves , modified by the interaction of the vortex lattice vibrations with sound .

According to the electrodynamical analogy, the Abrikosov–Tkachenko vortex lattice may equivalently be understood as the (bosonic) Wigner crystal of vortices, analogous to the (fermionic) Wigner crystals formed by electrons . On the other hand, for large vortex densities $`n_v`$, the lattice melts, at zero temperature due to quantum fluctuations, and the vortices form a highly correlated quantum liquid. Indeed, long-range logarithmic interactions, experienced by our (bosonic) vortices, have already been shown by Laughlin in his seminal paper on the (fermionic) FQHE to yield excellent overlap with incompressible quantum Hall states (cf. , Table I; see Ref.  for the present bosonic case).

In the limit of small vortex mass $`m_v0`$, we have large quantum fluctuations of the vortices around their equilibrium lattice sites, as well as confinement of the vortices to the lowest Landau level of the Hamiltonian (1). The dimensionless parameter

$$\nu _v=\frac{\mathrm{\Omega }}{\pi \rho }=\left(\frac{d}{\mathrm{}_v}\right)^2,$$
(4)
where $`d\rho ^{1/2}`$ is the interparticle distance and $`\mathrm{}_v=\sqrt{\pi /\mathrm{\Omega }}`$ the mean distance between vortices, then controls the relative strength of quantum fluctuations and therefore the (zero temperature) melting transition into the quantum Hall state. The above parameter is the vortex filling fraction corresponding to the ratio of the number of vortices, associated with $`n_v=\mathrm{\Omega }/\pi `$, and the number of “flux quanta” associated with the “magnetic” field of strength $`\rho `$, which is just the number of particles itself. The filling factor $`\nu _v`$ of the vortex rotational Landau levels is just the inverse of the filling of particle rotational Landau levels, $`\nu =1/\nu _v`$. The Laughlin wave function for the vortices reads, in the usual complex notation for the vortex coordinates $`Z_i=X_i+iY_i`$,

$$\left(\mathrm{\Psi }_v\right)_{\mathrm{Laughlin}}=\underset{i>j}{}\left(Z_iZ_j\right)^{1/\nu _v}\mathrm{exp}\left[\frac{\pi \rho }{2}\underset{i}{}|Z_i|^2\right],$$
(5)
where $`1/\nu _v`$ is even.

The Hamiltonian (1) has been studied in the context of vortex matter phase transitions in layered high-$`T_c`$ superconductors , where the possibility of a melting transition for the vortex lattice is due to the low Cooper-pair superfluid density $`\rho `$ in these materials . In the simplest approach to melting of matter ordered in crystals, one assumes a Lindemann criterion , which simply states that the departure of the vortices from their lattice positions should approximately equal their separation $`\mathrm{}_v`$. In more refined treatments , one takes the vortex lattice and Laughlin state wave functions for the vortices, and evaluates their corresponding energies, equating them at the transition point. In the limit of small vortex mass, the zero temperature transition point from the vortex lattice to the Laughlin state is, using either method, Lindemann or “exact” evaluation of correlation energies, consistently within a factor of order unity found to be at

$$\nu _v=(\nu _v)_m0.1(\text{melting filling factor}).$$
(6)
At absolute zero, for $`\nu _v<(\nu _v)_m`$, the system is in the vortex lattice phase, and for $`\nu _v>(\nu _v)_m`$, it is a vortex quantum liquid, i.e., a “melted” vortex crystal. For small interactions $`g1`$, the resulting state was shown to high accuracy to be a quantum Hall state of particles . In the next section, we argue that for sufficiently large interactions $`g`$, the resulting quantum Hall liquid is, conversely, one of vortices as the elementary constituents of the liquid.

## 4 Vortex Liquid States in Harmonic Traps

### 4.1 Coexistence of Solid and Liquid Phases

An isotropic, harmonically trapped 2D gas possesses, in the Thomas-Fermi limit, an inhomogeneous density profile of the form

$$\rho (r)=\rho _0\left(1\frac{r^2}{R^2}\right).$$
(7)
Here,

$$\rho _0=\left(\frac{N}{\pi g}\right)^{1/2}\left(\omega _{}^2\mathrm{\Omega }^2\right)^{1/2},R^2=2\sqrt{\frac{gN/\pi }{\omega _{}^2\mathrm{\Omega }^2}}$$
(8)
are the central density and the squared Thomas-Fermi radius of the condensate cloud, respectively. In contrast to an essentially homogeneous system like a 2D superconductor, the melting condition (6) then becomes a local criterion due to its dependence on the local density. While the density is a locally dependent quantity, it is important to recognize for the ensuing analysis that the observed vortex arrays possess a striking regularity . Even though the particle density is inhomogeneous, the vortex density is to a very good approximation the constant $`n_v=\mathrm{\Omega }/\pi `$ .

To reach the vortex quantum Hall state in the melted phase, in preference to a particle quantum Hall state, the following two conditions have to be fulfilled. Firstly, for the lowest particle Landau level description to break down, large particle interactions have to mix the first two rotational particle Landau levels , such that $`\omega _{}+\mathrm{\Omega }2\mathrm{\Omega }<g\rho `$. Second, for melting to occur the rotation rate has to be larger than $`\pi \rho (\nu _v)_m`$. These two conditions give

$$\pi \rho (\nu _v)_m<\mathrm{\Omega }<g\rho /2.$$
(9)
We conclude that, for a finite $`\mathrm{\Omega }`$–region to exist in which the vortex quantum Hall state is preferable over other states, $`g>g_c`$ needs to be fulfilled, where

$$g_c2\pi (\nu _v)_m0.6.$$
(10)
In addition, our small mass assumption can only be justified if the Coriolis force on vortices is (much) smaller than the vortex velocity part of the Magnus force, i.e. when $`2m_v\mathrm{\Omega }\rho \kappa `$. This leads to $`\mathrm{\Omega }g\rho /\mathrm{ln}[R/a_c]`$ which is for large systems more stringent than the Landau level mixing criterion $`\mathrm{\Omega }<g\rho /2`$ by a logarithmic factor $`\mathrm{ln}[R/a_c]/2`$. However, for the purpose of the present discussion, in the mesoscopic trapped gases, the difference between the two criteria can be neglected.

We distinguish the following cases, cf. Fig. 2. At sufficiently low angular velocity, $`\mathrm{\Omega }<\mathrm{\Omega }_{c1}`$, with

$$\mathrm{\Delta }\mathrm{\Omega }_{c1}=\omega _{}\mathrm{\Omega }_{c1}=\frac{g\omega _{}}{2\pi N[(\nu _v)_m]^2},$$
(11)
the system separates into three distinct phases, cf. Fig. 2 (a) (it is understood that all expressions for $`\mathrm{\Delta }\mathrm{\Omega }`$ are valid in the large $`N`$ limit, i.e., for $`\mathrm{\Delta }\mathrm{\Omega }/\omega _{}1`$). The center of the trap is occupied by the vortex lattice (region I), surrounded by two vortex liquid phases, in which the vortex lattice has melted. Phase II represents the vortex liquid, and is defined by Eq. (9). The inner and outer radii of the shell are given by

$`{\displaystyle \frac{x_{\mathrm{min}}^2}{R^2}}`$ $`=`$ $`1{\displaystyle \frac{1}{g_c\sqrt{N/g\pi }}}{\displaystyle \frac{2\mathrm{\Omega }}{\sqrt{\omega _{}^2\mathrm{\Omega }^2}}},`$ (12)
$`{\displaystyle \frac{x_{\mathrm{max}}^2}{R^2}}`$ $`=`$ $`1{\displaystyle \frac{1}{\sqrt{gN/\pi }}}{\displaystyle \frac{2\mathrm{\Omega }}{\sqrt{\omega _{}^2\mathrm{\Omega }^2}}},`$ (13)
and the area of the shell $`\mathrm{\Delta }A`$, relative to $`\pi R^2`$, therefore is

$$\frac{\mathrm{\Delta }A}{\pi R^2}=\frac{2\mathrm{\Omega }}{\sqrt{(gN/\pi )(\omega _{}^2\mathrm{\Omega }^2)}}\left(\frac{g}{g_c}1\right).$$
(14)
The number of vortices in the shell is

$`(N_v)_{\mathrm{shell}}`$ $`=`$ $`{\displaystyle \frac{4\mathrm{\Omega }^2}{\omega _{}^2\mathrm{\Omega }^2}}\left({\displaystyle \frac{g}{g_c}}1\right).`$ (15)
For the width of the vortex liquid shell to be a significant fraction of $`R`$ in the limit $`N\mathrm{}`$, such that $`(N_v)_{\mathrm{shell}}N`$, $`\mathrm{\Delta }\mathrm{\Omega }`$ has to be a significant fraction of $`\mathrm{\Delta }\mathrm{\Omega }_{c1}`$. Furthermore, the density must be approximately homogeneous for the “magnetic field” $`\rho `$ to be approximately constant and one well-defined vortex Laughlin state (5) to exist. The two conditions that the region II phases in Fig. 2 be sufficiently thin for a constant density approximation to apply, and $`(N_v)_{\mathrm{shell}}N`$, therefore require that the coupling $`g`$ is sufficiently close to its critical value $`g_c`$. The outermost shell IIIa shown in Fig. 2 (a) is filled by particles confined to their lowest rotational Landau level.

The maximum number of vortices in the vortex liquid shell is reached at $`\mathrm{\Omega }=\mathrm{\Omega }_{c1}`$. For given $`N`$ and $`g`$, this maximum number has the value

$`\mathrm{max}[(N_v)_{\mathrm{shell}}]={\displaystyle \frac{Ng_c}{\pi }}\left(1{\displaystyle \frac{g_c}{g}}\right).`$ (16)
If the rotation velocity grows and reaches $`\mathrm{\Omega }_{c1}`$, the vortex lattice completely melts, see Fig. 2 (b). The gas cloud in the region II remains superfluid. Further increase of $`\mathrm{\Omega }`$ leads to the contraction of region II and eventual disappearance of the vortex liquid phase at $`\mathrm{\Omega }=\mathrm{\Omega }_{c2}`$, where

$$\mathrm{\Delta }\mathrm{\Omega }_{c2}=\omega _{}\mathrm{\Omega }_{c2}=\frac{2\pi \omega _{}}{gN}.$$
(17)
For $`\mathrm{\Omega }_{c1}\mathrm{\Omega }\mathrm{\Omega }_{c2}`$, the number of vortices in the inner melted area II is given by the expression

$`(N_v)_{\mathrm{shell}}`$ $`=`$ $`2\mathrm{\Omega }\left({\displaystyle \frac{gN/\pi }{\omega _{}^2\mathrm{\Omega }^2}}\right)^{1/2}\left(1{\displaystyle \frac{2\sqrt{\pi }\mathrm{\Omega }}{\sqrt{gN(\omega _{}^2\mathrm{\Omega }^2)}}}\right),`$
which vanishes at $`\mathrm{\Omega }=\mathrm{\Omega }_{c2}`$ and coincides with (16) at $`\mathrm{\Omega }=\mathrm{\Omega }_{c1}`$.

The strongly correlated fractional quantum Hall states of particles in a harmonic trap require that the detuning $`\mathrm{\Delta }\mathrm{\Omega }_{\mathrm{particle}}/\omega _{}=1\mathrm{\Omega }_{\mathrm{particle}}/\omega _{}\nu ^2g/N`$ . On the other hand, provided the inequality $`gg_c`$ is fulfilled, the vortex quantum Hall states occurs in an outer shell already at rotation frequencies $`\mathrm{\Omega }<\mathrm{\Omega }_{c1}\mathrm{\Omega }_{\mathrm{particle}}`$ smaller than the rotation frequencies needed for the particle quantum Hall states. Thus, even if only an outer shell of the sample is filled by the vortex quantum Hall state, the number of vortices can be still sufficiently large to observe the characteristics of the vortex Laughlin state to be discussed in the next subsection.

### 4.2 Detection

We now come to discuss experimental procedures to verify that the state of the vortex system is either crystallized or a quantum fluid. We assume in the following that the part of the sample under consideration is occupied either by the vortex lattice or the vortex Laughlin phase.

The melting transition can, for example, be detected by letting the gas expand, and then investigating the resulting spatial configuration of vortices (cf. Fig. 3). The primary difference which reveals itself is that vortex liquids present a disordered structure compared to the regular crystal. The information on the degree of order is encoded in the (assumed isotropic) pair-correlation function $`g(r)`$. This information is obtained from pictures like the ones shown in Fig. 3, by counting the number of vortices inside a shell, of radial size the intervortex spacing, at a distance $`r`$ from a given (central) vortex. For a crystalline ordered structure, there are oscillations in the pair-correlation function around its asymptotic value, whereas these oscillations vanish for the liquid. If the number of vortices is large so that the counting statistics can be made accurate enough to extract the small-distance behavior $`r0`$ of the pair-correlation function, one should find $`g(r)r^{2/\nu _v}`$ for the vortex Laughlin liquid. An alternative method to detect the transition to a melted quantum Hall state was proposed in , by directly verifying a relation connecting condensate fraction and density profile. One may also conceive of using Fourier-space techniques in situ, and analyze Bragg diffraction patterns by scattering light off a collection of vortices. Finally, we mention that a very recent experiment has indeed observed, adding a quartic trapping potential term to achieve $`\mathrm{\Omega }\omega _{}`$, a transition from a vortex lattice to a disordered structure .

## 5 Conclusions

The preceding analysis was carried out at absolute zero. In fact, the limitation to zero temperature can be relaxed. The temperature dependence of the friction force on the vortices can be estimated using the result of ; the 2D friction coefficient $`\eta T^4/c_s^6`$. Magnus force domination corresponds to $`\rho \eta `$, which leads to $`T\rho g^{3/4}`$, provided $`T\rho g`$. On the other hand, since quantum Hall states are incompressible, excitations of the vortex liquid are separated from the ground state, for finite $`m_v`$ by a gap of order of the Landau level spacing $`T_{\mathrm{Landau}}=\kappa \rho /m_v\rho g/\mathrm{ln}[R/a_c]`$ . In the case of large $`g`$, the gap is dominated by the interaction between the vortices, and $`T_{\mathrm{int}}\rho \kappa ^2\mathrm{ln}[R/a_c]`$. Therefore, if $`TT_{\mathrm{Landau}}`$ respectively $`TT_{\mathrm{int}}`$, excitations are frozen out and the state of the system in thermal equilibrium essentially coincides with the Laughlin state. For the interaction strengths under consideration here, the temperatures $`T_{\mathrm{Landau}}`$ and $`T_{\mathrm{int}}`$ are both comparable with the Kosterlitz-Thouless temperature $`T_c^{\mathrm{KT}}=\pi \rho /2`$. The vortex Laughlin state thus forms practically for any temperature well below the transition temperature $`T_c^{\mathrm{KT}}`$ to superfluid behavior.

To realize the vortex quantum Hall liquid, at the low particle densities required, moderately strong particle interactions are necessary, $`g0.6`$, both to accomodate a sufficiently large number of vortices into the rotating sample and to drive the system into the vortex quantum Hall state. The ultimate limit on the maximal value of the particle coupling constant $`g`$ in a rapidly rotating, inhomogeneous 2D low-density Bose gas is not known. However, we may conclude from the study of Bose gases with large scattering lengths in Ref.  that $`gO(1)`$ should be possible before the system enters a solid, uncondensed phase.

Experimental confirmation of the vortex quantum Hall state provides a verification of the Hamiltonian (1), which stems from superfluid hydrodynamics , and the consequent quantization procedure for vortices as massive, “charged” particles. The present investigation should, then, help to shed light on the debate about the physical reality of vortex matter quantum phase transitions, whose occurence has been extensively discussed for the high-$`T_c`$ superconductors. It provides an analogous phenomenon in a pure, uncomplicated system, for which the force fields on the vortices are exactly known.

Finally, the present investigation raises the interesting question if vortices can also exist if only a very small number of particles per vortex, of order $`5\mathrm{}10`$, is present. In an abstract sense, a quantized vortex is a singular hole in the order parameter manifold with a certain residue. The hole effectively acts on the particles, making them spin around the hole in a manner prescribed by the quantization of circulation condition. The question which naturally poses itself, then, is if the vortex vacuum can also exist even in the limit of no particles being present. Formulated differently, one may ask for the lower limit of particle density, if any, for which the entity quantized vortex still remains a well-defined physical object.

Discussions with J. R. Anglin, G. Blatter, J. I. Cirac, F. Dalfovo, E. J. Mueller, and P. Zoller are gratefully acknowledged. P. O. F. has been supported by the Austrian Science Foundation FWF and the Russian Foundation for Basic Research RFRR, U. R. F. by the FWF, and A. R. by the European Union under Contract No. HPRN-CT-2000-00125.

## References