text:
# Evidence for the double degeneracy of the ground-state in the 3D ±𝐽 spin glass

## I Introduction

Although a quarter century has passed since the pioneering work on spin glasses by Edwards and Anderson EA75 , many fundamental problems remain unsolved even for the simplest models of these materials Binder-review . In the present paper, we report the results of Monte Carlo simulations of the three-dimensional $`\pm J`$ Ising spin glass, focusing on the nature of the low-temperature phase Hatano-SG1 ; Gubernatis00 . We believe we have found important evidence of a doubly-degenerate ground-state. With the same density of states used to compute this evidence, we found this character changing as the temperature is increased towards the spin glass transition temperature. Just below this temperature, the behavior is consistent with the standard mean-field picture that has an infinitely degenerate ground state.

The spin-glass state is characterized by randomly quenched exchange interactions with frustration. Often systems exhibiting this state are modeled by the $`\pm J`$ Ising model. In three dimensions this model is defined by the Hamiltonian EA75

$$\underset{i,j}{}J_{ij}\sigma _i\sigma _j,$$
(1)
where the Ising spins $`\sigma _i`$ fluctuate thermodynamically, while each exchange interaction $`J_{ij}`$ is quenched to $`\pm J`$ randomly. The summation runs over the nearest-neighbor pair of sites on a cubic lattice.

The frustration of local spin configurations by the exchange interactions generates many local minima in the free energy landscape, each minimum representing a seemingly random configuration of the spins. At low temperatures, the spins, however, may freeze into some configuration, and this freezing is the essence of the spin-glass “order.” Because of the nature of this order, equilibration of spin glasses in simulations and experiments is often very hard.

### I.1 Nature of the spin-glass phase

It is becoming increasingly accepted that the $`\pm J`$ model has a finite-temperature spin-glass phase transition in three dimensions, particularly after several recent Monte Carlo studies Kawashima ; Hukushima ; pal\_and\_car ; Ballesteros . The nature of the low-temperature phase, however, is still controversial. Historically, the controversy has been mainly between advocates of the mean-field picture and the droplet picture. The mean-field advocates maintain the existence of an infinite number of global minima of the free energy in the low-temperature phase. The condition is called replica symmetry breaking RSB . The breaking of this symmetry is rigorously true for the mean-field (Sherrington-Kirkpatrick) model SK75 ; Parisi79  of spin glasses. The question is whether it is also found in finite-dimensional non-mean-field models. The droplet advocates bray ; droplet-letter ; droplet  on the other hand, assert that the free-energy landscape has only two global minima that are connected through spin inversion symmetry. In the droplet picture, the nature of the three-dimensional ground state is seemingly less exotic than in the mean-field picture.

The difference in the ground state degeneracy between the two pictures becomes quantitative when we use the overlap order parameter of the spin-glass phase EA75 . To define this order parameter, we replicate the random exchange interactions in (1) and change the Hamiltonian to

$$\underset{\alpha =1,2}{}\underset{i,j}{}J_{ij}\sigma _i^{(\alpha )}\sigma _j^{(\alpha )},$$
(2)
where the superscript on the spin variables labels a replica. The overlap order parameter is then defined as the Hamming distance between spin configurations of the two replicas EA75 :

$$q=\frac{1}{L^d}\underset{i=1}{\overset{L^d}{}}\sigma _i^{(1)}\sigma _i^{(2)}.$$
(3)
Here $`L`$ is the linear size of the system and $`d`$ is the dimensionality, which is three throughout the paper.

The spins of one replica and those of the other are thermodynamically independent. Still they can be correlated because of the common random exchange interactions. In the high-temperature limit, the spin configurations are uncorrelated, and hence the overlap order parameter tends to zero as $`L^{d/2}`$. In the low-temperature limit, on the other hand, the spin configurations are frozen into energy minima, and hence the overlap order parameter can have a finite value.

In the droplet picture, the overlap order parameter at low temperatures and in the thermodynamic limit can take only two values equal and opposite to each other, as the frozen spin configuration of one replica is either macroscopically identical to the configuration of the other replica or to its inverted configuration. (We note, however, that the identity needs only to be macroscopic; the configurations can differ locally.) Hence the overlap order parameter $`q`$ takes one of the two values with an equal probability. In the mean-field picture, on the other hand, the overlap order parameter can take various values. Because of the many free-energy minima, the frozen spin configurations of the two replicas can be macroscopically different.

More explicitly, we define the order-parameter distribution as

$$P(q)\left[\frac{1}{𝒵}D(E,q)e^{\beta E}𝑑E\right]_{\mathrm{av}.},$$
(4)
where the partition function is given by

$$𝒵\underset{\{\sigma \}}{}e^{\beta E(\{\sigma \})}=𝑑q𝑑ED(E,q)e^{\beta E}.$$
(5)
In these two equations $`D(E,q)`$ is the normalized density of states for the energy $`E`$ and the overlap order parameter $`q`$. The energy in Eqs. (4) and (5) and other energies hereafter (including those in the Monte Carlo simulations) are the energies of the replica Hamiltonians (2). The square brackets in Eq. (4) denote the random average over various samples $`\{J_{ij}\}`$. The distribution $`P(q)`$ is normalized so that

$$P(q)𝑑q=\frac{2}{L3}\underset{q=1}{\overset{q=1}{}}P(q)=1.$$
(6)
Physical quantities that are functions of $`q`$ only are calculated as

$$\left[f(q)\right]_{\mathrm{av}.}=f(q)P(q)𝑑q,$$
(7)
where the angular brackets denote the thermodynamic average; e.g., the spin-glass susceptibility

$$\chi _{\mathrm{sg}}=L^d\beta 2\left[q2q2\right]_{\mathrm{av}.}$$
(8)
and the Binder parameter

$$g_{\mathrm{sg}}(T,L)\frac{3}{2}\left(1\frac{\left[q4\right]_{\mathrm{av}.}}{3\left[q2\right]_{\mathrm{av}.}^2}\right).$$
(9)
The droplet picture and the mean-field pictures have the functional forms of $`P(q)`$ as shown in Fig. 1.

In the droplet picture, the order-parameter distribution in the thermodynamic limit has two delta-function peaks, indicating two possible thermodynamic states (Fig. 1a). In the mean-field picture the order-parameter distribution is continuous between these peaks (Fig. 1b).

In practice, we can only simulate finite systems for which the order-parameter distribution looks like Fig. 1c. We need to see whether the distribution in Fig. 1c converges to Fig. 1a or Fig. 1b as $`L\mathrm{}`$. If $`lim_L\mathrm{}P(0)0`$, then significant doubt is cast upon the mean-field picture.

While some simulations appear to support infinite degeneracy Marinari96 ; Inigues97 ; Marinari99 , others do not Hartmann00 ; Palassini99 ; Palassini00 ; Palassini01 ; komori , and a recent analytic treatment argues against it Newman96 ; Newman98 . We note that the droplet picture is a theory concerning the zero-temperature fixed point, while most studies suggesting the validity of the mean-field picture are based on numerical estimates of $`P(0)`$ at $`T`$ only as low as $`0.7T_c`$, where the glass transition temperature $`T_c`$ is approximately 1.0. If the doubly degenerate ground state predicted by droplet picture is correct, it should at least be seen at temperatures lower than 0.7$`T_c`$, the closer to $`T=0`$ the better. Unfortunately, the slow dynamics of spin glasses has hindered numerical simulations from exploring the vicinity of the zero-temperature fixed point. We comment that a recent numerical study Moore98  demonstrated that results of the Migdal-Kadanoff approximation appear to support the mean-field picture near and below the glass transition temperature, but eventually support the droplet picture as $`T0`$. Uncertainty, however, remains because the Migdal-Kadanoff approximation favors the droplet picture, even for the mean-field model comment1 ; comment2 .

During the course of our work and just after, the possibility that the behavior of the spin glass does not unequivocally fit into one of these two pictures over the entire temperature range began to be appreciated. For example, Krzakala and Martin Krzakala01  proposed the novel perspective that the entropy fluctuations lead to a trivial $`P(0)`$ at zero temperatures, even if there are zero-energy large-scale excitations (complex energy landscape). They further proposed that such a situation should arise in the the three-dimensional $`\pm J`$ Ising spin glass, and argued that if the energy landscape is complex with a finite number of ground state families, then replica symmetry breaking reappears at finite temperatures. This perspective contested the standard picture. Palassini and Young Palassini99 ; Palassini00 ; Palassini01  studied this scenario and concluded that the size dependence of $`P(q)`$ around $`q=0`$ is trivial and does not support the ultrametric picture.

Katzgraber, Palassini, and Young Katzgraber01a ; Katzgraber01b  addressed the issue of Gaussian versus bimodal distribution by simulating at finite temperatures the Gaussian exchange model and concluded that as in the bimodal exchange model $`P(0)`$ is trivial in the thermodynamic limit and suggested the existence of low finite energy excitations that cost finite energy and whose surface has fractal dimension less than the spatial dimension of the system. However, the sizes they studied were quite small and so they concluded there might be a crossover at larger sizes to a different behavior, such as a droplet or replica symmetry breaking picture.

Two related papers by Hed et al. domany1 ; domany2  suggest that the spin-glass phase possesses some characteristics of the mean-field description as a non-trivial $`P(0)`$ and a hierarchical (but not ultrametric) structure of the pure states; nevertheless, they also claim that this phase is consistent with the Fisher-Huse scenario of the droplet picture. Correlated spin domains serve as the cores of zero enegy excitations.

The recent work of Marinari et alMarinari01  adopted the cluster analysis of Hed et alHed01  and found strong continuity among physical features for $`T>0`$ and those found at $`T=0`$, leading to a scenario with emerging mean-field like characteristics that are enhanced in the large volume limit for $`T>0`$. These mean-field like features arise with entropic fluctuations. More recently Lamarcq et alLamarcq01  have studied the fractal dimension of the clusters that are the low-lying excitations of the model.

There are still other papers, for example, Houdayer1 ; Houdayer2 ; stein ; nishimori , all illustrating a shift from simply the droplet versus the mean-field picture to something more subtle, with a consensus still very much evolving. One focus is on the nature and geometrical structure of ground-state excitations.

### I.2 Present Study

In the present study, we significantly reduced the difficulty of the slow dynamics in Monte Carlo simulations by using a bivariate version of Berg and Neuhaus’s multicanonical Monte Carlo method. Multicanonical simulations are performed independent of temperature or of a range of nearby temperatures and estimate the density of states. From it we can in principle calculate expectation values at any temperatures. Our resulting estimates of $`P(0)`$ at low temperatures and as a function of lattice size suggest that at very low temperatures thermodynamic behavior could be consistent with a doubly degenerate ground state, while at higher temperatures, it is consistent with an infinite degenerate ground state. Hence, we found that an intrinsic temperature-independent quantity seemingly exhibited different looking equilibrium behavior at different temperatures.

We present our Monte Carlo results in Sect. II. The order-parameter distribution function $`P(q)`$ at the low temperature of $`T=0.3`$ exhibits features indicative of a doubly degenerate ground state: $`P(0)`$ decreases as the system size is increased. The low-temperature behavior of the Binder parameter also suggests double degeneracy. In Sect. III we discuss further implications of our results for the ground-state degeneracy. In Appendix, we describe our simulation method, namely a bivariate multicanonical Monte Carlo method. Monovariate multicanonical methods Berg91 ; Berg92  have been applied to spin glasses before Berg94 ; Berg96 ; Berg98 ; we found, however, that by using a bivariate version we could reduce the correlation time of the simulation significantly. We show that the autocorrelation time is approximately proportional to the system size.

## II Numerical Results

We have carried out a bivariate multicanonical Monte Carlo simulation of cubic systems with edges $`L=4`$ (1904 samples), $`L=6`$ (2843 samples), $`L=8`$ (1015 samples), and $`L=10`$ (1111 samples). The simulation method directly returns the density of states $`D(E,q)`$, a temperature-independent quantity. We describe the details of the method and its use in Appendix; however, because the method is subtle, relatively new, and quite different from standard methods, we now summarize some of what is discussed in the Appendix.

With the multicanonical method, sometimes called the entropic sampling method, we do not equilibrate the simulation at any value of $`T`$; that is, we never sample the steady-state distribution, generated by a Markov chain, that is supposed to represent the Boltzmann distribution. Instead, we sample from a steady-state distribution, generated by a Markov chain, that is adaptively constructed to be flat on the average. The flatness means we sample all accessible $`(E,q)`$ values with an equal probability. In other words, we sample all thermodynamics states equally. In fact, the sampling emphasizes regions where $`D(E,q)`$ is small and as a consequence are generally difficult to access with many other methods.

Sampling from this flat distribution also allows one to estimate $`D(E,q)`$. We obtain $`D(E,q)`$ for different system sizes. Once we obtain it for a given size, we can in principle obtain the the properties of the system for any temperature; that is, properties of the system at different temperature, no matter how different they may seem, all follow from the same $`D(E,q)`$.

The validity of our low $`T`$ predictions depends on the “flatness” extending to the ground-state energy. As we discuss more fully in the Appendix flatness over the entire range of energy is not essential if not assumed in the evaluation of expectation values. The multicanonical method generally is a good ground-state sampler. This conclusion is supported in part by our to-be-reported rapid convergence of the entropy of the ground state for low $`T`$. As we will show, our estimates of the ground state residual entropy are consistent with the ones which Hartmann obtained by a ground-state counting algorithm Hartmann99 . This agreement is one reason why we are confident our simulations are covering low-energy states properly. At high temperatures, we will also see that the same $`D(E,q)`$ estimates the Binder parameter to within two sigmas with the previous results Kawashima .

### II.1 Order-parameter distribution

From the density of states $`D(E,q)`$, we straightforwardly calculated the order-parameter distribution $`P(q)`$, following Eq. (4). Figure 2 shows the temperature dependence of $`P(q)`$ for $`L=8`$. The function is close to a Gaussian distribution at high temperatures and has a double-peak structure at low temperatures.

The results for $`T=0.3,0.4`$, and $`0.5`$ are shown in Fig. 3 for several values of $`L`$.

We clearly see the decreasing tendency of $`P(0)`$ as $`L\mathrm{}`$. The raw data, however, presumably underestimates the true values of $`P(0)`$: the $`P(q)`$ data for $`L=8`$ noticeably oscillates for $`q0.7`$ with the data point at $`q=0`$ happening to be in a valley of the oscillation. We presume that this oscillation is due to correlations between data at different values of $`q`$, for example, the data point $`P(0)`$ being correlated with $`P(0.01)`$; similar oscillations are seen in the data of others berg94 ; Palassini01 ; Hed01  that was obtained by quite different methods. The present numerical method, as is summarized in Appendix, generates a random walk in the macroscopic phase space. The frequency of access by the random walker may be statistically less in some area of the phase space. The density of states will be underestimated in such area. Thus the correlations between data at different values of $`q`$ can occur.

To obtain proper estimates of $`P(0)`$, we have to smooth out the oscillations. We did this by choosing seven data points at the intervals of $`\mathrm{\Delta }q0.1`$ over the range $`0q<0.7`$ and then least-squares fitting them to the function $`\mathrm{ln}P(q)c_0+c_1q2`$. The fitting parameter $`c_0`$ yields the estimate $`P(0)=\mathrm{exp}(c_0)`$. The smoothed curves and the thus-estimated values of $`P(0)`$ are shown in Fig. 4 for $`T=0.3`$.

We still see the decreasing tendency with increasing $`L`$.

On the other hand, for $`T0.7`$, $`P(0)`$ does not show this tendency, even though its value is calculated with the use of the same $`D(E,q)`$ (Fig. 5).

In fact, $`P(0)`$ appears to converge to a finite value, a behavior which we now argue is spuriously consistent with the low-temperature behavior predicted by the mean-field picture.

Independent of the degeneracy, we know from the scaling ansatz that at the critical point $`T=T_\mathrm{c}1`$, $`P(0)`$ should increase as $`L\mathrm{}`$. On the other hand, at very low temperatures, double degeneracy requires that $`P(0)`$ should decrease and infinite degeneracy requires that it should tend to a constant. What follows from our computed $`D(E,q)`$ for different system sizes is shown in Fig. 6. Near $`T_c`$, $`P(0)`$ does increase with an increasing $`L`$. At lower temperatures, however, it tends to a constant, a behavior supporting infinite degeneracy. At still lower temperatures, it tends to zero, a behavior consistent with double degeneracy.

We note that a crossover scenario from critical, mean-field, to droplet behavior was predicted from the Migdal-Kadanoff approximation Moore98 . In particular we point out the similarity of our Fig. 6 to their Fig. 5. We thus suggest that most previous Monte Carlo studies, claiming to see behavior supporting the mean-field picture (infinite degeneracy), based on results for $`T`$ only as low as $`0.7`$, missed behavior consistent with the droplet picture (double degeneracy) which only appears at much lower temperatures.

Palassini and Young Palassini01  have studied the scaling of the smoothed quantity $`x(1/2)=_{1/2}^{1/2}P(q)𝑑q`$, and found a crossover scaling between $`T=0`$ behavior, where $`x(1/2)`$ becomes trivial for $`L\mathrm{}`$ and finite-temperature behavior, where the non-trivial part of $`P(q)`$ has a much weaker dependence on $`L`$ and is possibly size independent. The crossover is consistent with the qualitative features of our results. In fact, Palassini and Young’s Fig. 4, showing this crossover, is strikingly similar to our Fig. 6. The remark by Palassini and Young that our $`P(0)`$ drops dramatically at low $`T`$ as $`L`$ increase was based on the preliminary analysis of our data where $`P(0)`$ was determined in the absence of the smoothing. Palassini and Young also point to possible different behavior between models with Gaussian distributed exchange interactions, where $`P(0)`$ might be non-trivial as $`L\mathrm{}`$ and bimodal ones, where $`P(0)`$ becomes trivial.

We also note the recent zero temperature work by Hed, Domany, and Hartmann Hed01 . These investigators also found the need to smooth the values of $`P(q)`$, and for better statistics they chose to study $`x^{}=_{0.4}^{0.7}P(q)𝑑q`$. In the ground state they claim that $`x^{}`$ scales to a small non-trivial value. To be more precise, they first separated $`P(q)`$ into a part that comes from the big peaks close to $`q=\pm 1`$ that have an $`L`$-dependent tail at $`q=0`$ which scales to zero and into a part more proper to ground-state excitations whose scaling with $`L`$ is the central issue. For this latter part they claim $`x^{}`$ is non-trivial in the thermodynamic limit.

### II.2 Binder parameter

The Binder parameter for spin glasses $`g_{\mathrm{sg}}`$, defined by Eq. (9), is essentially the kurtosis of the order-parameter distribution $`P(q)`$. Because of the dimensionless combination of the second moment and the fourth moment, this parameter (except for effects due to correction to scaling) is expected to be independent of the system size at fixed points, i.e., $`T=0`$, $`T=T_\mathrm{c}`$, and $`T\mathrm{}`$ . For conventional phase transitions such as ferromagnetic transitions, its temperature dependence for various system sizes has a crossing point at the critical temperature. At the high-temperature $`(T\mathrm{})`$ fixed point, the order-parameter distribution $`P(q)`$ should be Gaussian. Hence we should have

$$q4P(q)𝑑q=3\left(q2P(q)𝑑q\right)2,$$
(10)
or

$$\left[q4\right]_{\mathrm{av}.}=3\left[q2\right]_{\mathrm{av}.}^2.$$
(11)
Thus the high temperature fixed-point value of the Binder parameter is zero, i.e., $`g_{\mathrm{sg}}(\mathrm{},L)=0`$. In the high-temperature phase, the Binder parameter $`g_{\mathrm{sg}}(T,L)`$ is renormalized to be zero from above as $`L\mathrm{}`$. At the low-temperature ($`T=0`$) fixed point, on the other hand, if the order parameter takes only two values $`\pm q_0`$, as in the thermodynamic limit of usual ferromagnets, we have

$$\left[q^4\right]_{\mathrm{av}.}=\left[q^2\right]_{\mathrm{av}.}^2=q_0^4$$
(12)
and hence the fixed-point value of the Binder parameter is unity: $`g_{\mathrm{sg}}(0,L)=1`$. In the low-temperature phase, the Binder parameter $`g_{\mathrm{sg}}(T,L)`$ is renormalized to be unity from below as $`L\mathrm{}`$. At $`T=T_\mathrm{c}`$, the Binder parameter $`g_{\mathrm{sg}}(T_\mathrm{c},L)`$ is expected to have a nontrivial universal value between zero and unity. Thus the crossing point of $`g_{\mathrm{sg}}(T,L)`$ should give the critical temperature $`T_\mathrm{c}`$.

Our Monte Carlo simulation found the crossing point of the Binder parameter as shown in Fig. 7 note2 .

The critical point $`T_c`$ should be in the region $`0.8T_\mathrm{c}1.1`$. For the moment, because of strong corrections to scaling, it is difficult for us to carry out sophisticated scaling analysis and obtain a more accurate estimate of $`T_\mathrm{c}`$. Previous studies Kawashima ; Hukushima ; pal\_and\_car ; Ballesteros  on much larger systems claim $`T_\mathrm{c}1.1`$.

We now offer further evidence for a doubly degenerate ground state, as opposed to an infinitely degenerate one, by reporting our results for the low-temperature behavior of the Binder parameter in Fig. 8.

These results strongly suggest that the Binder parameter tends to unity as $`T0`$ and $`L\mathrm{}`$. This behavior is consistent with Fig. 1a and not with Fig. 1b, as the Binder parameter is less than unity even at $`T=0`$ and $`L\mathrm{}`$. Our Monte Carlo results thus clearly support double degeneracy (see Fig. 9).

### II.3 Residual entropy

We calculated the entropy density from the difference between the energy and the free energy, i.e.

$$s=\frac{\left[E\right]_{\mathrm{av}.}F}{NT},$$
(13)
where $`F`$ is the free energy calculated from the Monte Carlo output $`D(E,q)`$ via

$$F=\frac{1}{\beta }\mathrm{log}Z=\frac{1}{\beta }\mathrm{log}\left(\underset{E,q}{}D(E,q)e^{\beta E}\right).$$
(14)
Figure 10 shows the energy and free-energy densities.

The difference between them is shown in Fig. 11.

It is clear that at $`T=0.1`$ the estimates of the entropy are virtually the residual entropy at $`T=0`$. The residual entropy is plotted in Fig. 12.

The convergence to the thermodynamic limit is quite rapid and the entropy seems to remain finite in the thermodynamic limit. Hartmann Hartmann99  also observed a rapid convergence to a finite entropy from a ground-state search on systems up to $`L=8`$. He obtained a similar estimate of the residual entropy $`s(T=0)=0.051(3)k_\mathrm{B}`$.

Although Hartmann used the existence of the residual entropy as evidence for the mean-field picture, its existence is on the contrary entirely consistent with the droplet picture. The degeneracy of the ground states predicted by these pictures is the degeneracy of thermodynamic (macroscopic) states, while the residual entropy comes from the degeneracy of microscopic states. The distinction is important to note.

Because the energy of a finite-sized $`\pm J`$ model is discrete, there is an inevitable degeneracy of the ground states. The issue is whether the degeneracy arises from these microscopically degenerate states or from many macroscopically different states. To make the distinction clearer, we will now consider a toy model in which we quench the exchange interactions into a periodic configuration with a unit cell of linear size $`l`$. This model has a ground state with a periodic spin configuration. We will assume, however, that every unit cell has one connected cluster of spins such that the spin inversion of the cluster does not change the ground-state energy (See Fig. 13).

We refer to such spin clusters as “zero-energy droplets” bray ; barahona . The number of such droplets in a system of linear size $`L`$ is $`N_{\mathrm{zed}}=(L/l)^d`$, where $`d`$ is the dimensionality. The degeneracy of the ground-state energy is $`2^{N_{\mathrm{zed}}}`$. Therefore, the residual entropy density of the toy model takes a finite value

$$s_{\mathrm{toy}}(T=0)=L^d\times k_\mathrm{B}\mathrm{ln}\left(2^{N_{\mathrm{zed}}}\right)=k_\mathrm{B}l^d\mathrm{ln}2.$$
(15)
This model, on the other hand, produces a $`P(q)`$ consistent with the droplet picture: Consider two replicas of the toy model. Without loss of generality, we can fix the spin configuration of one replica and calculate contributions from different spin configurations of the other replica. The distribution of the order parameter is given by

$$\begin{array}{ccc}|q|=1\hfill & \text{with the probability}\hfill & 2^{N_{\mathrm{zed}}},\hfill \\ |q|=12v/L^d\hfill & \text{with the probability}\hfill & 2^{N_{\mathrm{zed}}}\times N_{\mathrm{zed}},\hfill \\ |q|=14v/L^d\hfill & \text{with the probability}\hfill & 2^{N_{\mathrm{zed}}}\times \frac{1}{2}N_{\mathrm{zed}}(N_{\mathrm{zed}}1),\hfill \\ \mathrm{}\hfill & & \\ |q|=12nv/L^d\hfill & \text{with the probability}\hfill & 2^{N_{\mathrm{zed}}}\times \left(\begin{array}{c}N_{\mathrm{zed}}\\ n\end{array}\right)\hfill \\ \mathrm{},\hfill & & \end{array}$$
(16)
where $`v`$ is the volume of each zero-energy droplet. In the thermodynamic limit $`L\mathrm{}`$ or equivalently $`N_{\mathrm{zed}}\mathrm{}`$, the probability for $`n=N_{\mathrm{zed}}/2`$ overwhelms and hence $`P(q)`$ converges to a delta-function peak at

$$q=1\frac{N_{\mathrm{zed}}}{2}2vL^d=1\frac{v}{l^d}.$$
(17)
and to a second delta-function peak at $`q`$.

## III Ground-state degeneracy of the $`\pm J`$ model

The toy model just introduced is easily converted to a more realistic one for the ground-state degeneracy of the $`\pm J`$ model by allowing unequally-sized droplet volumes $`v_i`$ ($`i=1,2,3,\mathrm{},N_{\mathrm{zed}}(L)`$). In defining the volume of these droplets, we always take the minimal volume; for example, if a zero-energy droplet contains another zero-energy droplet, we decompose them into two droplets as shown in Fig. 14.

In addition, the possible maximum volume of the zero-energy droplet is $`L^d/2`$: if there is a zero-energy droplet larger than the half of the system, we re-define the original ground-state spin configuration by flipping the largest zero-energy droplet.

The ground-state degeneracy is $`2^{N_{\mathrm{zed}}(L)}`$, and hence the residual entropy density is given by

$$s(T=0;L)=L^dk_\mathrm{B}\mathrm{ln}\left(2^{N_{\mathrm{zed}}(L)}\right)=L^dN_{\mathrm{zed}}(L)k_\mathrm{B}\mathrm{ln}2.$$
(18)
We now define the distribution of the zero-energy droplets as $`n_{\mathrm{zed}}(v)`$ so that we have

$$N_{\mathrm{zed}}(L)=L^d_1^{L^d/2}n_{\mathrm{zed}}(v)𝑑v.$$
(19)
and rewrite Eq. (18) as

$$s(T=0;L)=k_\mathrm{B}\mathrm{ln}2\times _1^{L^d/2}n_{\mathrm{zed}}(v)𝑑v.$$
(20)
Next, we consider the ground-state value of $`|q|`$, and again we fix the spin configuration of one replica to the original ground-state spin configuration. Various contributions come from the spin configurations of the other replicas with some of the zero-energy droplets flipped. We extend Eq. (16) and obtain

$$\begin{array}{cc}|q|=1\hfill & \text{with the probability }2^{N_{\mathrm{zed}}}\text{,}\hfill \\ |q|=12v_i/L^d\hfill & \text{with the probability }2^{N_{\mathrm{zed}}}\text{ for }i=1,2,\mathrm{},N_{\mathrm{zed}}\text{,}\hfill \\ |q|=12(v_i+v_j)/L^d\hfill & \text{with the probability }2^{N_{\mathrm{zed}}}\text{ for }i,j=1,2,\mathrm{},N_{\mathrm{zed}}\text{ with }i<j\text{,}\hfill \\ \mathrm{}\hfill & \\ |q|=1(2/L^d)_{m=1}^nv_{i_m}\hfill & \text{with the probability }2^{N_{\mathrm{zed}}}\hfill \\ & \text{for }i_1,i_2,\mathrm{},i_n=1,2,\mathrm{},N_{\mathrm{zed}}\text{ with }i_1<i_2<\mathrm{}<i_n\text{,}\hfill \\ \mathrm{}\hfill & \end{array}$$
(21)
Summing all contributions, we find

$$|q|=1\frac{2}{2^{N_{\mathrm{zed}}}L^d}\left[\underset{i}{}v_i+\underset{i<j}{}(v_i+v_j)+\mathrm{}\right]=1\alpha _{\mathrm{zed}},$$
(22)
where

$$\alpha _{\mathrm{zed}}(L)\frac{1}{L^d}\underset{i=1}{\overset{N_{\mathrm{zed}}}{}}v_i=_1^{L^d/2}vn_{\mathrm{zed}}(v)𝑑v.$$
(23)
We can now compute $`q2`$ in the same way and obtain

$$q2=\frac{1}{2^{N_{\mathrm{zed}}}}\left[\underset{i}{}\left\{1\frac{2}{L^d}v_i\right\}^2+\underset{i<j}{}\left\{1\frac{2}{L^d}(v_i+v_j)\right\}^2+\mathrm{}\right].$$
(24)
which reduces to

$$q2=(1\alpha _{\mathrm{zed}})2+\beta _{\mathrm{zed}},$$
(25)
or

$$q2|q|2=\beta _{\mathrm{zed}},$$
(26)
where

$$\beta _{\mathrm{zed}}(L)\frac{1}{L^{2d}}\underset{i=1}{\overset{N_{\mathrm{zed}}}{}}v_i2=\frac{1}{L^d}_1^{L^d/2}v2n_{\mathrm{zed}}(v)𝑑v.$$
(27)
Similarly we have

$$q43q22+2|q|4=2\gamma _{\mathrm{zed}},$$
(28)
where

$$\gamma _{\mathrm{zed}}(L)\frac{1}{L^{4d}}\underset{i=1}{\overset{N_{\mathrm{zed}}}{}}v_i4=\frac{1}{L^{3d}}_1^{L^d/2}v4n_{\mathrm{zed}}(v)𝑑v.$$
(29)
If the order-parameter distribution $`P(q)`$ has only two delta-function peaks at the $`L\mathrm{}`$ limit, then the left-hand sides of Eq. (26) and Eq. (28) vanish. This condition is satisfied if the density distribution of the zero-energy droplets $`n_{\mathrm{zed}}(v)`$ decays fast enough as $`v\mathrm{}`$. If $`n_{\mathrm{zed}}(v)`$ decays exponentially, the integrals in Eqs. (27) and (29) give finite values; hence $`\beta _{\mathrm{zed}}(L)=O(L^d)`$ and $`\gamma _{\mathrm{zed}}(L)=O(L^{3d})`$. If $`n_{\mathrm{zed}}(v)`$ decays as $`v^x`$ with $`x>2`$, we have $`\beta _{\mathrm{zed}}(L)=O(L^{d\mathrm{min}(1,x2)})`$ and $`\gamma _{\mathrm{zed}}(L)=O(L^{d\mathrm{min}(3,x2)})`$.

On the other hand, if the order-parameter distribution $`P(q)`$ has a non-trivial part as $`L\mathrm{}`$, $`\alpha _{\mathrm{zed}}(L)`$, $`\beta _{\mathrm{zed}}(L)`$, and $`\gamma _{\mathrm{zed}}(L)`$ all give finite values as $`L\mathrm{}`$. These conditions are satisfied if the distribution function decays as $`n_{\mathrm{zed}}(v)v^2`$.

In our numerical results, both $`\beta _{\mathrm{zed}}(L)`$ and $`\gamma _{\mathrm{zed}}(L)`$ appear to be decreasing as $`L`$ is increased (Fig. 15).

The decrease is roughly $`O(L^1)`$, although the statistical errors are too large to make a more definite statement. This behavior is additional evidence for a doubly degenerate ground state. From our numerical results we suggest that the size distribution of the zero-energy droplets decays as $`n_{\mathrm{zed}}(v)v^x`$ with $`x>2`$.

## IV Summary

The results of our multicanonical Monte Carlo calculation suggest a doubly degenerate ground state. Our main findings supporting double degeneracy are: (i) $`P(q)`$ near $`q0`$ decreasing at low temperatures as the system size is increased; (ii) the Binder parameter approaching unity at low temperatures; (iii) the effect of the ground-state degeneracy on moments of the overlap order parameter. In different temperature ranges, however, we further demonstrated, within finite size limitations, that the same density of states gave different scaling behaviors as a function of the system size. Only at very low temperatures is the scaling consistent with a doubly degenerate ground state.

We believe our numerical results are consistent with recent analytic and numerical work and spin-glass themes pointing to a decoupling of the degeneracy of the ground state from the standard droplet and mean-field pictures and placing emphasis on the nature of the low-energy excitations and entropic fluctuations with them. Our numerical work provides information on these excitations only through the density of states. The $`P(0)`$ predicted from this density of states is trivial at low temperatures. For finite systems, our density of states predicts crossover thermodynamic behavior as the temperature is lowered.

###### Acknowledgements.

We thank E. Domany for a helpful conversation. We gratefully acknowledge computational support from the Applied Mathematics program of the Department of Energy for the use of the massively parallel computers at NERSC, from the Center for Nonlinear Studies (CNLS) at Los Alamos National Laboratory for the use of its parallel computer Avalon, and from the Los Alamos Directed Research program at Los Alamos for computer time on the ASCI Blue Mountain computer. The use of the random-number generator library SPRNG of National Center for Supercomputing Applications (NCSA) at University of Illinois is also gratefully acknowledged. A part of the computation was also carried out on a PC-cluster parallel computer ARK of Aoyama Gakuin University Nakata00 .

## Appendix A Multicanonical Monte Carlo Method

We will describe the multicanonical Monte Carlo method from a viewpoint Gubernatis00  slightly different from other authors Berg91 ; Berg92 ; Berg94 ; Berg96 ; Berg98 . We first compare it with random-sampling and canonical methods.

Suppose that we perform a Monte Carlo simulation aimed at sampling from the weight $`W(E)`$. After the simulation has equilibrated, that is, after the Markov chain reaches a steady-state, each microscopic state (spin configuration) is generated at a rate proportional to $`W(E)`$. Also suppose that in “equilibrium” we construct a histogram $`h(E)`$ of the energy $`E`$ values associated with each configuration (state) generated. The rate at which a state with the energy $`E`$ appears is proportional to $`D(E)W(E)`$, where $`D(E)`$ is the density of states.

In a random-sampling simulation, one generates all microscopic states at the same rate. In other words, $`W(E)=\mathrm{const}.`$, and hence the resulting histogram is proportional to the density of states: $`h(E)D(E)`$ (see Fig. 16).

In many cases, the density of states becomes very small in low-energy regions. Hence simulations based on the random sampling generally have difficulty investigating low-temperature properties because the low-energy states, which dominate the thermodynamic average at low temperatures, are generated only infrequently, if at all.

Importance sampling algorithms, like the Metropolis algorithm for canonical ensemble, were developed to overcome this difficulty. Here states are generated with the importance weight $`W(E)e^{\beta E}`$, and hence the energy histogram is $`h(E)e^{\beta E}D(E)`$. Thus, in ideal cases, low-energy states appear at a greater frequency at low temperatures. In spin-glass simulations, however, the simulation tends to get stuck in local minima of the free-energy landscape.

In the multicanonical simulations, the importance weight is set to $`W(E)1/D(E)`$, and hence the energy histogram should be flat. Of course, one cannot set $`W(E)1/D(E)`$ because the density of states is a priori unknown. However, the multicanonical method is a procedure that makes the importance weight converge to the reciprocal of the density of states Smith95 . The aim is to ensure better statistics of low-energy states than in the random-sampling method and, at the same time, generate more high-energy states than in the canonical method so that the system can escape local free-energy minima.

Berg and others have invented Berg91 ; Berg92  and applied Berg94 ; Berg96 ; Berg98  this method to spin glasses, taking the weight of each spin configuration $`\{\sigma \}`$ either as

$$W(\{\sigma \})1/D\left(E(\{\sigma \})\right)$$
(30)
or

$$W(\{\sigma \})1/D\left(q(\{\sigma \})\right),$$
(31)
where $`D(E)`$ is the density of states with respect to the energy and $`D(q)`$ is the one with respect to the overlap order parameter.

In the ideal situation, the multicanonical simulation effectively generates a random walk in the macroscopic phase space. Because all the thermodynamic states should appear with the same probability, one can expect that the autocorrelation time of a thermodynamic variable ($`E`$ in the above case) is proportional to the volume (the number of the sites). In the simulation of the $`\pm J`$ model, the energy change after each spin flip is of the order of $`J`$. On the other hand, the upper and lower bounds of the energy are of the order of $`NJ`$, where $`N=L^d`$ is the number of the spins. Hence the number of spin flips necessary for the random walk to cover the entire energy phase space is on the average of the order of $`N^2`$. Thus if random walk is in energy space, a one-dimensional space, the measured energy should decorrelate after $`N^2`$ Monte Carlo steps or $`N`$ Monte Carlo sweeps. (A sweep is the attempt to flip each Ising spin once on average.)

This ideal situation was not achieved in the previous applications of the multicanonical simulations to spin glasses Berg94 ; Berg96 ; Berg98 . Berg et al. reported $`\tau N^{2.8}`$ for the importance weight (30Berg94  and $`\tau N^{2.42}`$ for the importance weight (31Berg98 . (Here $`\tau `$ is measured in the unit of Monte Carlo sweep.) This suggests that the slow dynamics in the low-energy region was not removed completely in these mono-variate multicanonical simulations. (Note, however, that Berg et al. did not use the conventional definition of the autocorrelation time in measurement of $`\tau `$. This might be a contribution to the above power law behavior.)

In the present study, we carried out the simulation using the bivariate importance weight:

$$W1/D(E,q).$$
(32)
As shown in Fig. 17 (a), we almost achieved the ideal situation $`\tau N`$.

The random walk in the two-dimensional phase space is illustrated in Fig. 17 (b). The autocorrelation plotted in Fig. 17 (a) was computed as follows: We first measured the autocorrelation as the simple Monte Carlo average without any weighting:

$`C_E(t)`$ $``$ $`{\displaystyle \frac{1}{N_{\mathrm{MCS}}}}{\displaystyle \underset{j}{}}E(j)E(j+t)`$ (33)
$`C_q(t)`$ $``$ $`{\displaystyle \frac{1}{N_{\mathrm{MCS}}}}{\displaystyle \underset{j}{}}q(j)q(j+t),`$ (34)
where $`N_{\mathrm{MCS}}`$ is the number of Monte Carlo steps for measurement. Then the autocorrelation time was computed as the integrated time:

$`\tau _E{\displaystyle \frac{1}{C_E(0)}}{\displaystyle \underset{t=0}{\overset{\stackrel{~}{t}_E}{}}}C_E(t),`$ (35)
$`\tau _q{\displaystyle \frac{1}{C_q(0)}}{\displaystyle \underset{t=0}{\overset{\stackrel{~}{t}_q}{}}}C_q(t),`$ (36)
where $`\stackrel{~}{t}_E`$ ($`\stackrel{~}{t}_q`$) is defined as the time when $`C_E(t)`$ ($`C_q(t)`$) first becomes negative.

There is a drawback with the bivariate simulation: large electronic memory is needed to store the data from which the bivariate histogram is constructed. For the Ising model we used the number of all possible discrete states to equal the number of bins in the energy phase space. This number was of the order of $`L^d`$ and was similarly sized in the $`q`$ phase space. Hence an array of the size $`L^{2d}`$ was necessary to store each density of states $`D(E,q)`$, the histogram $`h(E,q)`$, and related work space. The combined storage requirements amounted to nearly 30MB per sample for $`L=10`$ and 400MB per sample for $`L=16`$. Based on direct comparisons with the original mono-variate multicanonical method, we have, however, concluded that for spin-glass simulations the use of both $`E`$ and $`q`$ is essential for improving the slow dynamics Hatano-SG1 . We comment that bivariate multicanonical Monte methods have been previously used in different contexts, e.g., protein folding simulations bivar-MCMC1 ; bivar-MCMC2 .

The actual algorithm of the multicanonical method has the following structure:

1. Make a guess at the density of states, $`D_0(E,q)`$. Start the following loop with $`i=0`$.
2. The outer loop (multicanonical-iteration loop):

1. Set the importance weight to $`W_i=1/D_i`$. Calculate the energy and the order parameter of the initial spin configuration, $`E_0`$ and $`q_0`$. Start the following loop with $`j=0`$.
    2. The inner loop (spin-update loop):

    1. Choose a spin to be updated. Calculate the energy and the order parameter, $`E_j^{}`$ and $`q_j^{}`$, assuming that the spin is actually flipped.
        2. Calculate the spin-flip probability as

        $$P_{\mathrm{flip}}=\mathrm{min}(1,\frac{W_i(E_j^{},q_j^{})}{W_i(E_j,q_j)})$$
        (37)
        3. Draw a random number $`R[0,1]`$. If $`R<P_{\mathrm{flip}}`$, flip the spin and make the substitution $`E_{j+1}=E_j^{}`$ and $`q_{j+1}=q_j^{}`$. Otherwise, keep the current spin configuration and make the substitution $`E_{j+1}=E_j`$ and $`q_{j+1}=q_j`$.
        4. Increment the histogram bin $`h_i(E_{j+1},q_{j+1})`$ by one. This step should be skipped during the burn-in stage (specifically given below), that is, until the simulation reaches “equilibrium” for the importance weight $`W_i`$.
        5. Increment the number of steps $`j`$ by one and go to the step (i) until $`j`$ exceeds a pre-determined number.
    3. Measure physical quantities $`Q(E,q)`$ via

    $$Q_i=\frac{1}{Z_i}\underset{E,q}{}Q(E,q)e^{\beta E}h_i(E,q)D_i(E,q),$$
    (38)
    where the partition-function estimate is given by

    $$Z_i=\underset{E,q}{}e^{\beta E}h_i(E,q)D_i(E,q).$$
    (39)
    In practice, this step was skipped until the histogram $`h_i`$ became reasonably flat (see below for more specifics).
    4. Update the guess of the density of states as follows: Smith95

    $$D_{i+1}(E,q)=\frac{D_i(E,q)\left(h_i(E,q)+1\right)}{\mathrm{Normalizaion}}.$$
    (40)
    The normalization constant is chosen so that the integration of the density of states becomes unity. (An extra count is added to the histogram in (40) in order to avoid zero division in setting the next importance weight $`W_{i+1}=1/D_{i+1}`$ Smith95 .)
    5. Increment the number of multicanonical iterations $`i`$ by one and go to the step (a), until the histogram meets certain conditions of being flat.

The method is an iterative procedure requiring an initial guess $`D_0(E,q)`$ of $`D(E,q)`$. We found our average results to be quite independent of this guess, but used the following particularly convenient and efficient choice: Our choice was a flat one for $`L=4`$, a case where the convergence was very rapid. The results for this system size were then scaled and used for the starting points for larger systems.

For a given iteration, step (iv) should be skipped until the simulation reaches “equilibrium” for the importance weight $`W_i`$. For $`L=4`$, 6, 8, and 10, we used 10 000, 40 000, 50 000, and 100 000 Monte Carlo sweeps for this burn-in stage. The number of the burn-in sweeps were approximately a hundred times longer than the autocorrelation time given in Fig. 17a. The autocorrelation times were determined before we began our production running. After the burn-in we used 1 000 000, 4 000 000, 5 000 000, and 10 000 000 sweeps for the different systems sizes during each iteration step; that is, the number of sweeps in each iteration was 10 000 times the measured autocorrelation time and one hundred times the burn-in time.

Step (c) was skipped until the histogram $`h_i`$ became reasonably flat. The number of iterations required for this varied as function of $`L`$ and from one random bond configuration to another for a given $`L`$. For $`L=4`$, the range was 2 to 15; for $`L=6`$, 2 to 23; for $`L=8`$, 2 to 276; and for $`L=10`$, 4 to 227. True flatness is never achieved nor is it necessary. What is important is the histogram be reasonably flat, extends to the band edges, and lacks spikes at the band edges. In defining a “reasonable” flatness criteria, which is the criterion for the iteration to have converged, we were motivated by the following observations: While converging, the histogram typically is first concentrated in the mid-part of the density of states and expands outward towards the upper and lower bounds of the energy (band edges) and the order parameter until it becomes flat. At the leading edges, spiky structures sometimes appear as the histogram is crossing over from occupied to unoccupied states. When the edges are reached, the spikes disappear. We defined the histogram to be flat, whenever the total count at the minimum energy and that at $`E=0`$ were within one sigma of the average histogram count. In contrast to the usual Ising model, the band width for a bond configuration is not known a priori, but it is easy to estimate from short runs and then set it slightly larger than expected to handle extremal states. In general, we found the multicanonical method to be an efficient estimator of the lowest energy of a given bond configuration.

After convergence, we repeated the outer loop five times after so that we could detect erroneous estimates of the density of states; that is, using the just determine $`D(E,q)`$ as the starting point, we repeated multicanonical procedure until five consecutive iterations maintained convergence. Testing determined that five was an adequate number to avoid false convergence. We also divided the inner loop into ten blocks, from which we calculated the statistical errors. Because the length of each block is about 1000 times longer than the auto-correlation time, the absence of correlations between block averages is expected and standard estimates for statistical error were used.

Because it is an importance estimator, (38) could in principle be applicable before the histogram becomes flat; that is, achieving flatness is not an essential requirement for the accuracy of the method if importance weighted estimators are used. Even with convergence, we used such estimators for the sake of caution. We only collected data for converged results. With sufficient statistics, the histogram $`h_i`$ would be proportional to $`W_i(E,q)D(E,q)=D(E,q)/D_i(E,q)`$, where $`D(E,q)`$ is the true density of states. Thus the estimator (38) would give

$$Q_i\frac{1}{Z_i}\underset{E,q}{}Q(E,q)e^{\beta E}D(E,q),$$
(41)
with the partition function

$$Z=\underset{E,q}{}e^{\beta E}D(E,q).$$
(42)
We comment again that we could have calculated all physical quantities after the simulations were completed, provided we had stored the density of states for each random bond configuration on disk. This procedure would requires a huge amount of disk space because of the need to simulate a large number of random samples.

We also note that Marinari et al. Marinari\_comment  have recently tried to implement our bivariate version of Berg and Neuhaus’s multicanonical method and reported experiences very different from ours. They reported difficulties in ergodicly sampling the importance function, finding it only remained “flat” only by using an order of magnitude larger number of Monte Carlo sweeps than we used. This would correspond needing burn-in times and block sizes 1000 times our estimated auto-correlation time. More importantly the very small values of $`h_i`$ erratically riddling their histograms would translate into very small values of $`h_i(E,q)D_i(E,q)`$, which should have caused our measured values for a given bond configuration to vary widely. We did not observe this, and our estimates of autocorrelation time were consistent for different system sizes. As mentioned above, most of our data analysis was done “on-the-fly” to avoid storage problems so we cannot re-examine the actual data for our reported results. However, when a referee pointed out to us the Marinari et al. difficulties, we immediately executed a number of simulations for the $`L=8`$ system size causing Marinari et al. problems and looked at the same quantity that they reported in their Fig. 3. We never saw the behavior that they claimed occurs. The origin of the discrepancy between their results and ours is unknown to us.


format: text

text:
# Volume 85, Number 4 PHYSICAL REVIEW LETTERS 24 July 2000

Schulman Replies: Casati et al.  have failed to understand the goals of my article . They and I find ourselves in a situation that is common in discussions of the “arrow of time,” namely no disagreement on technical issues and no agreement on basic assumptions. Instead of criticizing their goals (presumably “to clear out the too-long standing confusion around (ir)reversibility of statistical laws” ), which occupy many workers in this area , I will try to make clear what problem I do address.

About 40 years ago Gold  argued that the thermodynamic arrow of time followed the cosmological one. This is a thesis of beauty and scope, but I found  that for logical consistency one needed to introduce a new element into the argument: the use of boundary conditions at two (usually remote) times. This is because the use of macroscopic initial conditions already fixes the thermodynamic arrow. A further benefit of this conceptual point is that if we are in a big-bang/big-crunch universe and if, relieved of the initial conditions perspective, one would give roughly symmetric boundary conditions, one would then have a framework within which to look for physical effects arising from the future big crunch. See , where these arguments are spelled out in greater detail and the relevant literature, particularly due to J. A. Wheeler, is presented.

A consequence of this framework is that one can imagine subsystems that, because of a particular history of isolation, survive the era of maximal expansion with arrow intact. Such a scenario has been the basis for attacks on Gold’s thesis . It is this isolation and survival phenomenon that is captured by the boundary conditions given in my article . The two times given there are some intermediate times, not particularly near the end points, such that at those times one system is dominated by the conditioning at one end, the other at the other end. If one desired a more complete picture, one could (in principle) arrange boundary conditions near the big bang and big crunch such that relatively isolated chunks of matter would pass the point of maximal expansion with arrow intact. (Obviously this requires greater spatial and dynamical richness than the cat map.) If the time interval under study is $`[t_1,t_2]`$, then one system would have low entropy at $`t_1`$, one at $`t_2`$ (and both times would be significantly distant from the end points). I modeled this by restricting each subsystem to a single coarse grain in its phase space at its respective “initial” time. But as part of a larger dynamical scheme (i.e., as a slice out of big bang to big crunch boundary value problem), each system would carry a great deal more information encoded in its microscopic coordinates. This further information ensures that each system possesses its arrow “before” (by its own clock) its encounter with the other. In terms of the remarks of Casati et al. this microscopic information prevents the bi-directional entropy increase envisioned in their criticism.

The goal of many researchers in this area (and one that Casati et al. inaccurately ascribe to me in ) is the elucidation of macroscopic irreversibility in the face of microscopic (near) reversibility. In fact I take for granted that this high purpose can be achieved. My intent is to see whether in our actual world there can be departures from the uniform direction that all (so far) observed systems appear to possess. This requires a formulation, namely two-time conditioning, that does not eliminate the phenomenon by definition (as would happen with initial conditions). Once one has in this way formulated the problem, one can seek actual physical realization. An important aspect of this is the conclusion  that we could actually see the stuff, contrary to the views of Wiener .

Because a (roughly) time-symmetric universe is the scenario in which opposite-arrow material is most likely to occur (frankly, I can’t think of any other reason for it to exist), the ability to identify and characterize such material and the possibility of its observation, have a direct bearing on significant physical questions.

L. S. Schulman
Physics Department
Clarkson University
Potsdam, New York 13699-5820
Received 3 April 2000
PACS numbers: 05.70.Ln, 05.20. y, 95.35. +d


text:
# 1 Introduction

## 1 Introduction

In order to understand the role of impurities and inhomogeneities in real physical systems, many statistical models with quenched randomness have been proposed. The effect of randomness on continuous phase transitions has thus been of great interest for many years. The main question is to determine whether the randomness leaves the critical properties of the pure system unchanged. One prediction concerning models with random bonds comes from the Harris criterion , which states that bond randomness changes the values of the critical exponents only if the specific heat exponent $`\alpha `$ is positive. The effects of bond randomness have been studied intensively in the $`2D`$ Ising model and more recently in the $`q<4`$ states Potts model. In the former case, the $`log`$ singularity of the specific heat is transformed in a $`log(log)`$ singularity , whereas in the latter case a new critical point with new critical exponents is reached . Furthermore, it has been recently realized that a weak bond randomness can also have strong effects on $`2D`$ spin systems possessing a first order transition. Following earlier work of Imry and Wortis , Hui and Berker  have shown, using a phenomenological renormalization group, that bond randomness (even infinitesimal) induces a second order transition (i.e. vanishing latent heat) in a system that would have undergone a first order one. This result has been proved rigorously by Aizenman and Wehr . At this level, a question naturally arises concerning the universality class of this second order transition. In order to test these predictions, the $`q`$ states Potts model ($`q>4`$ is necessary to get a first order transition) with bond randomness has been recently extensively analyzed numerically , the conclusion being that the transition is of second order type due to disorder and the critical behaviour depends on $`q`$. Analytical results on the $`q`$ states Potts model are, as usual, much more difficult to obtain. However, this is not the case for those systems presenting what Cardy called a “weak fluctuation-driven first order transition” . The typical example Cardy provided is a system of $`N`$ critical Ising models coupled by their energy density. Without disorder, this perturbation drives the system in a strong coupling regime where perturbative renormalization group (RG) analysis fails. Yet, for this model, it can be shown non-perturbatively that this runaway flow is indeed associated to a massive regime (with a finite correlation length depending on the coupling). This is a way to define this “weak fluctuation-driven first order transition” phenomenon. The addition of weak bond randomness, even infinitesimal, makes the system flow to $`N`$ decoupled critical Ising models . This result has been extended to the case of $`N`$ coupled Potts models . The result was a non-Ising like second order transition, the universality class depending on the sign of the coupling between the models. This class of examples is interesting because some infinitesimal bond randomness is able to change the “weak fluctuation-driven first order transition” into a continuous one. Concretely speaking, the RG runaway flow is completely turned over by disorder preventing a strong coupling regime. Nevertheless a lot of questions still remain open. What happens if we take a different coupling than the energy density one? Does infinitesimal disorder will also change the order of the transition?

To give some insights to these questions, we would like to study the following system: two Ising planes coupled by a spin-spin interaction where disorder (dilution or bond-randomness) is added in each plane. When the two Ising planes are critical, the spin-spin perturbation which is much stronger that the energy-energy perturbation drives the system in a strong coupling regime. It has been shown using integrable Toda field theories that a finite correlation length is generated , and the full spectrum of the theory can be computed too. In this sense, the situation is analogous to the model treated by Cardy , except that now the coupling is different and much stronger. Therefore, this model can bring some clues to the aforementioned questions. Moreover, it is worth noticing that $`N`$ Ising planes (with $`N1`$) coupled by the spin-spin interaction is a way to reach a $`3D`$ Ising model. $`N=2`$ may be regarded as a kind of first step in this direction. Therefore, we expect Aizenman-Wehr theorem not to apply here.

The plan of the paper is the following. In section 2, we present analytical field theory arguments showing that a finite amount of disorder is necessary to change the “weakly fluctuation-driven order first order transition” into a continuous one. In section 3, we confirm it by Monte Carlo numerical simulations. In section 4, we will summarize the results and discuss the whole phase diagram of the model.

## 2 Analytical results

### 2.1 Pure case

In this section, we first summarize results available in the pure case. The Hamiltonian then reads

$$H=J\underset{}{}(\sigma _i^1\sigma _j^1+\sigma _i^2\sigma _j^2)J^{}\underset{i}{}\sigma _i^1\sigma _i^2,$$
(1)
where $`\sigma _i^{1,2}=\pm 1`$, the first sum is over the nearest neighbours on a 2D square lattice, the same coupling $`J`$ has been chosen for the two models and a coupling $`J^{}`$ has been considered between the two planes. This model can be described in the continuum limit by the following action

$$𝒜=𝒜_{Is}^1+𝒜_{Is}^2+m\mathrm{d}^2x(\epsilon _1(x)+\epsilon _2(x))+g\mathrm{d}^2x\sigma _1(x)\sigma _2(x),$$
(2)
where $`𝒜_{Is}^i`$ denotes the action of the critical Ising model, $`ϵ_i(x)`$ the Ising energy operator, $`\sigma _i`$ the Ising spin operator, $`m(TT_c)`$ and $`gJ^{}`$.

When $`g=0`$, the action correspond to two decoupled massive Ising models which can be described by free massive Majorana fermions or by free Dirac fermions or by a sine-Gordon model. If the strongly relevant coupling $`g`$ is switched on at the $`m=0`$ critical point, the system is known to be driven into a massive regime whose mass spectrum has been computed exactly . When $`g0`$, the bosonic representation is more suitable (in terms of Majorana fermions representation, this interaction is indeed non-local), and the generalized sine-Gordon reads

$$𝒜=\mathrm{d}^2x\frac{1}{2}(\varphi )^2+m\mathrm{cos}(\beta \varphi )+g\mathrm{cos}(\alpha \varphi \pm \delta ),$$
(3)
with $`\beta =\sqrt{4\pi },\alpha =\sqrt{\pi }`$ and $`\delta =\frac{\pi }{2}`$.

In order to obtain this expression, we write the Hamiltonian in terms of Majorana fermions and then use standard bosonization formula (see for example  for a review). It is worth noticing that the addition of a four spin interaction like $`K\underset{<i,j>}{}\sigma _i^1\sigma _j^1\sigma _i^2\sigma _j^2`$ in the original Hamiltonian (1) (which therefore defines a generalized Ashkin-Teller model) would not change the form of the action (3) except the substitution $`(\sqrt{4\pi },\sqrt{\pi })(\beta ,\frac{1}{2}\beta )`$ with $`\beta =\beta (K)`$.

The action (2) is a priori hard to study since it contains two strongly relevant perturbations and it is unfortunately non integrable. Notice that $`g`$ is much stronger that $`m`$ and scales around the critical point like $`g|m|^{7/4}|TT_c|^{7/4}`$. A perturbed conformal field theory approach is unappropriate and beyond feasibility because the flow is quickly driven in a strong coupling regime.

The action (3) which is one particular case of the double frequency sine-Gordon model has been qualitatively studied in the general case by Delfino and Mussardo . It has been argued by looking at the solitonic structure in both weak and strong coupling regime ($`mg`$ and $`mg`$) that when the ratio $`\frac{\alpha }{\beta }\frac{n^{}}{n}`$ is a rational number and $`\delta =\frac{\pi }{n}`$, a phase transition should occur for $`m>0`$ (i.e. $`T>T_c`$). Concerning our particular case, we can easily convince ourselves about this transition by looking at the evolution of the shape of the potential $`𝒱=m\mathrm{cos}(\beta \varphi )g\mathrm{sin}(\beta \varphi /2)`$ which passes from a periodic two degenerate minima to a periodic one absolute minimum situation when increasing the ratio $`g/m`$ (see Figure 1). From a Ginzburg-Landau point of view, it suggests that the transition is in the universality class of $`\mathrm{\Phi }^4`$, namely is of Ising type.

Let us come back to the original Hamiltonian (1). In the limit $`g\mathrm{}`$, the spins $`\sigma _i^1`$ and $`\sigma _i^2`$ are locked to the same value and the model reduces to only one Ising model with a critical temperature $`2J_c`$. For a generic value of $`g`$ the same scenario should occur for a critical temperature $`J_c<T_c(g)<2J_c`$ since $`g`$ tends to order the system in contrast to temperature. Therefore, there should be a critical line in the $`(T,g)`$ plane interpolating between the fixed points $`(J_c,0)`$ and $`(2J_c,\mathrm{})`$. In the $`gJ`$ limit, we can rewrite the Hamiltonian (1) in the basis $`\sigma _i=\sigma _i^1;\tau _i=\sigma _i^1\sigma _i^2`$ and integrate over the $`\tau `$ variable using a mean field approximation (we refer to  for a more detailed analysis in a different context). It remains one Ising model with approximate effective coupling $`2J(1c\frac{J^2}{g^2})2J_{\mathrm{eff}}`$, with $`c`$ a positive number satisfying a self-consistent equation. The Ising criticality is reached when $`J_{\mathrm{eff}}=J_c`$. This analysis is valid only close to the $`(2J_c,\mathrm{})`$ Ising fixed point, namely in the strong coupling limit.

According to Zamolodchikov c-theorem  (without disorder the theory is still unitary), the flow goes from the $`(J_c,0)`$ fixed point with central charge $`c=1`$ toward the ending Ising fixed point $`(2J_c,\mathrm{})`$ with central charge $`c=\frac{1}{2}`$. From the continuum limit point of view, it is worth noticing that the elementary excitations around this fixed point have nothing to do with the ones around the $`c=1`$ fixed point (the correspondence is indeed non local explaining why a standard perturbation theory around the $`c=1`$ fixed point is hopeless).

### 2.2 Disordered case

Let us now add bond randomness or dilution in each Ising planes. In the two limits $`g=0`$ and $`g=\mathrm{}`$ the problem reduces to the well known disordered Ising model and there are no new fixed points. Can we say something in the general situation? The phase diagram is parametrized by three variables, $`T,g`$ and $`\mathrm{\Delta }`$ the strength of the disorder. In the planes $`\mathrm{\Delta }=0`$ and $`g=0`$, the shape of the transition separating the ferromagnetic phase from the paramagnetic one is known respectively from results concerning the disordered Ising model and from the pure case analysis described above. Therefore, by a continuity arguments, we expect a critical surface joining the two Ising critical lines. Let us focus on the $`T=T_c`$ plane. In the continuum limit, the effect of disorder is to change $`m`$ into $`m(x)`$ in the action (3) with $`\overline{m(x)}=0`$, $`\overline{m(x)m(x^{})}=\mathrm{\Delta }\delta (xx^{})`$ for a Gaussian disorder. This action is very hard to study with usual tools and we can only give qualitative arguments. At fixed value for $`g`$, when $`m(x)>m^{}(g)=g/4`$, it implies a local high temperature region surrounded by low temperature regions. In the limit $`\mathrm{\Delta }m^{}(g)`$ (still at fixed $`g`$), the phenomenology of the critical disordered Ising model is recovered with a mixture of low and high temperature regions. This suggests that there should exist a finite critical value $`\mathrm{\Delta }^{}(g)`$ of the disorder strength able to disorder the system. In fact, the situation resembles the one of a disordered Ising model with $`\overline{m(x)}=m_0>0`$. This is clearly a non-perturbative phenomenon. We conjecture this transition to be of Ising type governed by the Ising fixed point (with $`c=\frac{1}{2}`$) at $`(2J_c,\mathrm{},\mathrm{\Delta }=0)`$ with a $`log(log)`$ behavior for the specific heat due to the disorder. Another non trivial fixed point at $`\mathrm{\Delta }^{}>0`$ cannot a priori be excluded but is very unprobable according to the symmetries of the model under consideration. Note that contrary to the examples treated in  where some infinitesimal disorder changes the weakly driven first order transition into a continuous one, here a finite value of the disorder is needed!

## 3 Numerical analysis

To put on firmer grounds these qualitative arguments, we have performed Monte Carlo simulations at $`T=T_c`$. Let us present our numerical results. We have considered two $`2D`$ Ising models, with bond dilution $`\mathrm{\Delta }`$ in both planes, i.e. $`P(J)=\mathrm{\Delta }\delta (J)+(1\mathrm{\Delta })\delta (J1)`$. The disorder is independent in each plane. We let the planes interact via a spin-spin coupling term, whose strength is $`g`$.

We have performed numerical simulations with 3 different lattice sizes ($`L=20,40,60`$) and for three values of the disorder: $`\mathrm{\Delta }=0.05,0.15,0.25`$. Lattice sizes may seem small with respect to previous studies on the disordered Ising model , however they seem to be in the scaling regime (see below) and so the extrapolation to large sizes is safe, at least up to the precision we are interested in. For any disorder value we have found clear evidences that, increasing the value of the coupling $`g`$, the system undergoes a phase transition from a paramagnetic to a ferromagnetic phase. The critical coupling values can be well fitted by a power law: $`g_c(\mathrm{\Delta })\mathrm{\Delta }^{1.8}`$ (see left inset in Fig. 2, where the errors are smaller than the symbol size). Note that the value of this exponent is very close to the theoretical one obtained by scaling arguments close to the $`c=1`$ fixed point, which is $`7/4`$.

For every $`\mathrm{\Delta }`$, the critical coupling $`g_c`$ has been determined as the crossing point of the Binder cumulants. This crossing point has almost no finite size effects. Two more estimations of $`g_c`$ can be obtained from the $`g`$ values where the specific heat and the susceptibility have a peak. In fact these values converge to $`g_c`$ in the large $`L`$ limit. However they show much larger finite size effects than the Binder cumulants crossing point. Then, in order to have a better estimation of $`g_c`$, we have used the data for every lattice size and we have considered finite size corrections through the formula $`g_c(L)=g_cAL^{1/\nu }`$, with $`\nu =1`$. We have found that our data are compatible with that formula for any disorder value and the $`g_c`$ estimation coincides with the one extracted from the Binder cumulants (see right inset in Fig. 2). The effect of increasing $`\mathrm{\Delta }`$ simply reflects in a larger value for the non-universal coefficient $`A`$.

In order to verify that the transition is governed by the pure $`2D`$ Ising model fixed point, we have analyzed our numerical data using finite size scaling and the critical exponents of the pure $`2D`$ Ising model, i.e. $`\nu =1`$ and $`\gamma =7/4`$. We have found very good data collapse for any disorder value $`\mathrm{\Delta }`$. In particular in Fig. 2 we show the collapse of the rescaled susceptibility for the larger disorder value $`\mathrm{\Delta }=0.25`$. For smaller disorder amplitudes the scaling is even better and the scaling function is more peaked signaling that the critical region is narrower.

## 4 Discussion and conclusions

In this paper we have presented analytical and numerical evidences that in two coupled 2D Ising models the “weak fluctuation-driven first order transition” changes to a continuous one only when a finite amount of disorder is added. Moreover, we have seen that even for finite and large disorder amplitudes, the phase transition is in the $`2D`$ disordered Ising universality class (Ising transition modified by logarithm corrections). The results we get can be summarized in the schematic $`3D`$ phase diagram depicted in Fig. 3. Our analysis strongly suggests that almost all the critical surface shown in Fig. 3 is governed by the $`2D`$ Ising fixed point, except the curve in the $`g=0`$ plane (governed by the $`c=1`$ fixed point) and the one in the $`T=0`$ plane (governed by a percolation fixed point111Note that the $`T=0`$ critical line is discontinuous at $`g=0`$, as can been seen in Fig. 3. In fact for every $`g0`$ the bond percolation threshold is $`1/\sqrt{2}`$, while for $`g=0`$ is $`1/2`$). Because of the known flow lines on this surface (the ones with an arrow in Fig. 3), we believe that any transition across the surface is governed by the $`c=1/2`$ fixed point located in $`(2J_c,\mathrm{})`$ (the rightmost big dot in Fig. 3) because the $`c=1`$ fixed point is strongly repulsive in the $`g`$ direction (even if it is attractive on the $`g=0`$ plane ). For an even clearer numerical check, a numerical estimation of the central charge would be wellcome.

As a conclusion, we would like to mention that several directions remain open. First of all, a theoretical framework to describe quantitatively such non-pertubative transition is needed. Secondly it would be nice extending this analysis for $`N>2`$ coupled Ising models and finally it is worth performing a similar analysis for correlated disorder.

Acknowledgements One of us (PS) would like to acknowledge interesting discussions with Vl. S. Dotsenko, M.-A. Lewis, G. Mussardo and A. A. Nersesyan. This work has been partially supported by the EC TMR Programme Integrability, non–perturbative effects and symmetry in Quantum Field Theories, grant FMRX-CT96-0012.


text:
# 1 Introduction

## 1 Introduction

Quantizations of triangular Lie bialgebras $`𝐋`$ with antisymmetric classical $`r`$-matrices $`r=r_{21}`$ are defined by a twisting element $`=f_{\left(1\right)}f_{\left(2\right)}𝒜𝒜`$ which satisfies the twist equations :

$$\begin{array}{c}\left(\right)_{12}\left(\mathrm{\Delta }\mathrm{id}\right)=\left(\right)_{23}\left(\mathrm{id}\mathrm{\Delta }\right),\hfill \\ \left(ϵ\mathrm{id}\right)=\left(\mathrm{id}ϵ\right)=1.\hfill \end{array}$$
(1)
Explicit form of the twisting element is quite important in applications because it provides explicit expressions for the quantum $``$-matrix $`_{}=_{21}^1`$ and for the twisted coproduct $`\mathrm{\Delta }_{}\left(\right)=\mathrm{\Delta }_{}\left(\right)^1`$.

The first nontrivial explicitly written twisting elements $``$ were given in the papers , ,  and . These twists can be defined on the following  carrier algebra $`L`$:

$$\begin{array}{cccc}[H,E]=E,& & & [H^{},E]=\gamma ^{}E,\\ [H,A]=\alpha A,& & & [H^{},A]=\alpha ^{}A,\\ [H,B]=\beta B,& & & [H^{},B]=\beta ^{}B,\\ [E,A]=[E,B]=0,& & & [A,B]=E\\ \alpha +\beta =1,& & & \alpha ^{}+\beta ^{}=\gamma ^{}.\end{array}$$

Explicit expressions for their twisting elements are

$$\begin{array}{c}\begin{array}{c}\mathrm{\Phi }_{}=\mathrm{e}^{\mathrm{H}\mathrm{H}^{}},\hfill \\ \mathrm{\Phi }_𝒥=\mathrm{e}^{\mathrm{H}\sigma },\hfill \\ \mathrm{\Phi }_𝒥=\mathrm{\Phi }_{}\mathrm{\Phi }_𝒥=\mathrm{e}^{\mathrm{A}\mathrm{B}\mathrm{e}^{\beta \sigma }}\mathrm{e}^{\mathrm{H}\sigma },\hfill \end{array}\begin{array}{c}r_{}=HH^{},\hfill \\ r_𝒥=HE,\hfill \\ r_𝒥=HE+AB,\hfill \end{array}\hfill \\ \sigma =\mathrm{ln}\left(1+E\right).\hfill \end{array}$$
(2)
Here $`r_{},r_𝒥,r_𝒥`$ are the corresponding classical $`r`$-matrices.

Carrier subalgebras $`𝐋`$ can be found in any simple Lie algebra $`g`$ of rank greater than 1.

It was demonstrated in  that these twists can be composed into chains. They are based on the sequences of regular injections constructed for the initial Lie algebra

$$g_pg_{p1}\mathrm{}g_1g_0=g.$$

To form the chain one must choose an initial root $`\lambda _0`$ in the root system $`\mathrm{\Lambda }\left(\mathrm{g}\right)`$, consider the set $`\pi `$ of its constituent roots

$$\begin{array}{c}\pi =\{\lambda ^{},\lambda ^{\prime \prime }\lambda ^{}+\lambda ^{\prime \prime }=\lambda _0;\lambda ^{}+\lambda _0,\lambda ^{\prime \prime }+\lambda _0\mathrm{\Lambda }\left(\mathrm{g}\right)\}\hfill \\ \pi =\pi ^{}\pi ^{\prime \prime };\pi ^{}=\left\{\lambda ^{}\right\},\pi ^{\prime \prime }=\left\{\lambda ^{\prime \prime }\right\}.\hfill \end{array}$$

and the subset $`\mathrm{\Lambda }_{\mathrm{\Lambda }_0}^{}`$ of roots orthogonal to $`\lambda _0`$ (the corresponding subalgebra in $`g`$ will be denoted by $`g_{\lambda _0}^{}`$ ).

It was shown that for the classical Lie algebras $`g`$ one can always find in $`g_{\lambda _0}^{}`$ a subalgebra $`g_1g_{\lambda _0}^{}g_0=g`$ whose generators become primitive after the extended twist $`\mathrm{\Phi }_𝒥`$. Such primitivization of $`g_kg_{k1}`$ (called the matreshka effect ) provides the possibility to compose chains of extended twists of the type $`\mathrm{\Phi }_𝒥`$,

$$\begin{array}{c}_{_{0p}}=_{k=0}^p\mathrm{\Phi }__\mathrm{k}\mathrm{\Phi }_{𝒥_\mathrm{k}},\hfill \\ \mathrm{\Phi }__\mathrm{k}\mathrm{\Phi }_{𝒥_\mathrm{k}}=_{\lambda ^{}\pi _\mathrm{k}^{}}\mathrm{exp}\left\{\mathrm{E}_\lambda ^{}\mathrm{E}_{\lambda _0^\mathrm{k}\lambda ^{}}\mathrm{e}^{\frac{1}{2}\sigma _{\lambda _0^\mathrm{k}}}\right\}\mathrm{exp}\{\mathrm{H}_{\lambda _0^\mathrm{k}}\sigma _{\lambda _0^\mathrm{k}}\}.\hfill \end{array}$$
(3)
Chains of twists quantize a large variety of $`r`$-matrices corresponding to Frobenius subalgebras in simple Lie algebras .

## 2 Construction of a full chain of twists

The main point in the construction of a chain is the invariance of $`g_{k+1}`$ with respect to $`\mathrm{\Phi }_{_\mathrm{k}𝒥_\mathrm{k}}`$. When these subalgebras are proper the canonical chains have only a part of $`B^+\left(g\right)`$ as the twist carrier subalgebra:

$$\begin{array}{ccccccc}\mathrm{}& g_{\lambda _0^k}^{}& g_{\lambda _0^{k1}}^{}& \mathrm{}& g_{\lambda _0^1}^{}& g_{\lambda _0^0}^{}& g\\ & & & & & & \\ \mathrm{}& g_{k+1}& g_k& \mathrm{}& g_2& g_1& g_0\end{array}$$
(4)
We would like to demonstrate that the effect of primitivization is universal and extends to the whole subalgebra $`g_{\lambda _0^k}^{}`$. It was shown in  that the invariance of a subalgebra in $`g_{\lambda _0^k}^{}`$ is only one of the forms of the primitivization . In general this is the existence (in the twisted Hopf algebra $`U_{_k𝒥_k}\left(g_{\lambda _0^k}^{}\right)`$ ) of a primitive subspace $`V_G^{k+1}`$ with the algebraic structure isomorphic to $`g_{\lambda _0^k}^{}`$. On this subspace the subalgebra $`g_{\lambda _0^k}^{}`$ is realized nonlinearly so $`V_G^{k+1}`$ is called deformed carrier space .

In this context the situation with the twists for $`U\left(sl(N)\right)`$ is degenerate: the subalgebra $`\left(sl(N)\right)_{\lambda _0^k}^{}`$ coincides with $`\left(sl(N)\right)_{k+1}`$ , i.e. $`V_G^{k+1}=V_{\left(sl(N)\right)_{\lambda _0^k}^{}}`$ .

In the case of $`U\left(so(M)\right)`$ the situation is different. Let the root system $`\mathrm{\Lambda }\left(\mathrm{s}\mathrm{o}(\mathrm{M})\right)`$ be

$$\{\pm e_i\pm e_ji,j=1,2,\mathrm{}M/2;ij\}$$

for even $`M`$ and

$$\{\pm e_i\pm e_j;\pm e_ki,j,k=1,2,\mathrm{}\left(M1\right)/2;ij\}$$

for odd $`M`$. Take $`e_1+e_2`$ as the initial root. Here the subalgebras $`g_{\lambda _0^{k1}}^{}`$ and $`g_k`$ in (4) are related as follows,

$$g_{\lambda _0^{k1}}^{}=g_kso^{\left(k\right)}(3)=so\left(M4k\right)so^{\left(k\right)}(3).$$

Consider the invariants of the vector fundamental representations of $`g_{k+1}=so(M4\left(k+1\right))`$ acting on $`g_k`$ :

$$\begin{array}{c}I_{2N+1}^a=\frac{1}{2}E_a^2+_{l=3}^N\left(E_{a+l}E_{al}\right),\hfill \\ I_{2N+1}^{ab}=E_aE_b+_{l=3}^N\left(E_{a+l}E_{bl}+E_{al}E_{b+l}\right),\hfill \end{array}$$
(5)
$$\begin{array}{c}I_{2N}^a=_{l=3}^N\left(E_{a+l}E_{al}\right),\hfill \\ I_{2N}^{ab}=_{l=3}^N\left(E_{a+l}E_{bl}+E_{al}E_{b+l}\right),\hfill \end{array}$$
(6)
The $`so^{\left(k\right)}(3)`$ summands are non-trivially deformed by $`\mathrm{\Phi }_{_{\mathrm{k}1}𝒥_{\mathrm{k}1}}`$:

$$\begin{array}{cc}\mathrm{\Delta }_{_{\mathrm{k}1}𝒥_{\mathrm{k}1}}\left(\mathrm{E}_{12}^\mathrm{k}\right)=& E_{12}^k1+1E_{12}^k+\left(1e^{\frac{1}{2}\sigma _{1+2}^{k1}}\right)I_{M4k}^{11}\\ & +I_{M4k}^1\left(e^{\sigma _{1+2}^{k1}}1\right),\\ \mathrm{\Delta }_{_{\mathrm{k}1}𝒥_{\mathrm{k}1}}\left(\mathrm{E}_{21}^\mathrm{k}\right)=& E_{21}^k1+1E_{21}^k+\left(e^{\sigma _{1+2}^{k1}}1\right)I_{M4k}^2e^{\sigma _{1+2}^{k1}}\\ & +\left(1e^{\frac{1}{2}\sigma _{1+2}^{k1}}\right)I_{M4k}^{22}.\end{array}$$

According to the main principle formulated above (despite the deformed co-structure of $`V_{g_{\lambda _0^{k1}}^{}}`$) the primitivization is realized on its isomorphic image $`V_G^{k+1}`$ contained in $`U_{_{k1}𝒥_{k1}}\left(g_{\lambda _0^{k1}}^{}\right)`$. To find this deformed carrier subspace $`V_G^{k+1}`$ it is sufficient to inspect the coproducts of invariants (5) and (6),

$$\begin{array}{c}\mathrm{\Delta }_{_\mathrm{k}𝒥_\mathrm{k}}\left(\mathrm{I}_{\mathrm{M}4\mathrm{k}}^1\right)=\hfill \\ =I_{M4k}^1e^{\sigma _{1+2}^k}+1I_{M4k}^1+I_{M4k}^{11}\left(1e^{\frac{1}{2}\sigma _{1+2}^k}\right),\hfill \end{array}$$

$$\begin{array}{c}\mathrm{\Delta }_{_\mathrm{k}𝒥_\mathrm{k}}\left(\mathrm{I}_{\mathrm{M}4\mathrm{k}}^2\mathrm{e}^{\sigma _{1+2}^\mathrm{k}}\right)=\hfill \\ =I_{M4k}^2e^{\sigma _{1+2}^k}1+e^{\sigma _{1+2}^k}I_{M4k}^2e^{\sigma _{1+2}^k}+I_{M4k}^{22}\left(1e^{\frac{1}{2}\sigma _{1+2}^k}\right).\hfill \end{array}$$

Now one can construct the following nonlinear primitive generators

$$\begin{array}{c}G_{12}^{k+1}=E_{12}^kI_{M4k}^1,\hfill \\ G_{21}^{k+1}=E_{21}^kI_{M4k}^2e^{\sigma _{1+2}},\hfill \\ H_{12}^{k1},\hfill \end{array}\begin{array}{c}\mathrm{\Delta }_{_\mathrm{k}𝒥_\mathrm{k}}\left(\mathrm{G}_{12}^{\mathrm{k}+1}\right)=\mathrm{G}_{12}^{\mathrm{k}+1}1+1\mathrm{G}_{12}^{\mathrm{k}+1},\hfill \\ \mathrm{\Delta }_{_\mathrm{k}𝒥_\mathrm{k}}\left(\mathrm{G}_{21}^{\mathrm{k}+1}\right)=\mathrm{G}_{21}^{\mathrm{k}+1}1+1\mathrm{G}_{21}^{\mathrm{k}+1},\hfill \\ \mathrm{\Delta }_{_\mathrm{k}𝒥_\mathrm{k}}\left(\mathrm{H}_{12}^{\mathrm{k}1}\right)=\mathrm{H}_{12}^{\mathrm{k}1}1+1\mathrm{H}_{12}^{\mathrm{k}1}.\hfill \end{array}$$

The subspace spanned by $`\{H_{12}^k,G_{12}^{k+1},G_{1+2}^{k+1}\}`$ forms the algebra $`so_G^{\left(k+1\right)}(3)`$ $`so^{\left(k+1\right)}\left(3\right)`$ :

$$\begin{array}{c}[H_{12}^k,G_{12}^{k+1}]=G_{12}^{k+1},\hfill \\ [H_{12}^k,G_{21}^{k+1}]=G_{21}^{k+1},\hfill \\ [G_{12}^{k+1},G_{21}^{k+1}]=2H_{12}^k.\hfill \end{array}$$

Therefore we obtain the deformed primitive space

$$V_G^{k+1}\left(g_{\lambda _0^k}^{}\right)=V\left(g_{k+1}\right)V\left(so_G^{\left(k+1\right)}(3)\right),$$

that can be considered as a carrier for the twists (2). The next extended Jordanian twist in the chain (that is defined on $`g_{k+1}`$) does not touch the space $`V\left(so_G^{\left(k+1\right)}(3)\right)`$. Consequently after all the steps of the chain we will still have a primitive subalgebra

$$𝒟=\underset{k=0}{\overset{p}{}}{}_{}{}^{}so_G^{\left(k+1\right)}(3)$$

defined on the sum of deformed spaces $`V\left(so_G^{\left(k+1\right)}(3)\right)`$.

Thus in the twisted Hopf algebra $`U_{_{0p}}\left(so\left(M\right)\right)`$ one can perform further twist deformations with the carrier subalgebra in $`𝒟`$ . The most interesting among them are the Jordanian twists defined by

$$\mathrm{\Phi }_{𝒥_\mathrm{k}}^\mathrm{G}=\mathrm{exp}\left(\mathrm{H}_{12}^\mathrm{k}\sigma __\mathrm{G}^\mathrm{k}\right)\mathrm{with}\sigma __\mathrm{G}^\mathrm{k}\mathrm{ln}\left(1+\mathrm{G}_{12}^{\mathrm{k}+1}\right)$$

This means that in the general expression for the twisting element $`_{_{0p}}`$ one can insert in the appropriate $`k0`$ places the Jordanian twisting factors defined on the deformed carrier spaces, i.e. to perform a substitution

$$\begin{array}{c}\mathrm{\Phi }__\mathrm{k}\mathrm{\Phi }_{𝒥_\mathrm{k}}\mathrm{\Phi }_{𝒥_\mathrm{k}}^\mathrm{G}\mathrm{\Phi }__\mathrm{k}\mathrm{\Phi }_{𝒥_\mathrm{k}}\mathrm{\Phi }_{𝒢_\mathrm{k}}\hfill \\ \mathrm{exp}\left\{I_{M4k}^{12}\left(1e^{\frac{1}{2}\sigma _{1+2}^k}\right)\right\}\mathrm{exp}\{H_{1+2}^k\sigma _{1+2}^k\}\hfill \\ \mathrm{exp}\left(H_{12}^k\sigma __G^k\right)\mathrm{exp}\left\{I_{M4k}^{12}\left(1e^{\frac{1}{2}\sigma _{1+2}^k}\right)\right\}\mathrm{exp}\left(H_{1+2}^k\sigma _{1+2}^k\right)\hfill \end{array}$$

This gives the full chain in the following form

$$\begin{array}{c}_{𝒢_{0p}}=_{k=p}^0\mathrm{\Phi }_{𝒢_\mathrm{k}}=\hfill \\ _{k=p}^0\left(\mathrm{exp}\left(H_{12}^k\sigma __G^k\right)\mathrm{exp}\left\{I_{M4k}^{12}\left(1e^{\frac{1}{2}\sigma _{1+2}^k}\right)\right\}\mathrm{exp}\{H_{1+2}^k\sigma _{1+2}^k\}\right).\hfill \end{array}$$
(7)
Obviously the additional twistings by $`\mathrm{\Phi }_{𝒥_\mathrm{k}}^\mathrm{G}`$ cannot be performed before the deformation of the corresponding spaces $`V_G^{k+1}`$ by the extended Jordanian twists $`\mathrm{\Phi }__\mathrm{k}\mathrm{\Phi }_{𝒥_\mathrm{k}}`$.

## 3 Applications

The previous result means that we have constructed explicit quantizations

$$_{𝒢_{0p}}=\left(_{𝒢_{0p}}\right)_{21}\left(_{𝒢_{0p}}\right)^1$$

of the following set of classical $`r`$-matrices:

$$r_{𝒢_{0p}}=\underset{k=0}{\overset{p}{}}\eta _k\left(H_{1+2}^kE_{1+2}^k+\xi _kH_{12}^kE_{12}^k+I_{M4k}^{12}\right)$$

Here all the parameters are independent.

The dimensions of the nilpotent subalgebras $`N^+\left(so\left(M\right)\right)`$ in the sequence $`g_{\lambda _0^p}^{}g_{\lambda _0^{p1}}^{}\mathrm{}g_{\lambda _0^0}^{}g`$ are subject to the simple relation:

$$dim\left(N^+\left(so\left(M\right)\right)\right)dim\left(N^+\left(so\left(M4\right)\right)\right)=2\left(dimd_{so\left(M4\right)}^v+1\right).$$

Taking this into account we see that the chains (7) are full in the sense that for $`p=p^{\mathrm{max}}=\left[M/4\right]+\left[\left(M+1\right)/4\right]`$ their carrier spaces contain all the generators of $`N^+\left(so\left(M\right)\right)`$. When $`M`$ is even-even or odd the total number of Jordanian twists in a maximal full chain $`_{𝒢_{0p\mathrm{max}}}`$ is equal to the rank of $`so\left(M\right)`$. Thus in the latter case the carrier subalgebra is equal to $`B^+\left(so\left(M\right)\right)`$.

It was demonstrated in  how to construct new Yangians using the explicit form of the twisting element. These new Yangians are defined by the corresponding rational solution of the matrix quantum Yang-Baxter equation (YBE). In particular, for the orthogonal classical Lie algebras $`so\left(M\right)`$ one needs the twisting element $``$ in the defining (vector) representation $`d^v`$ and the auxiliary operators: the flip $`P:vwwv`$ ($`P\mathrm{Mat}\left(M\right)\mathrm{Mat}\left(M\right)`$) and the operator $`K`$, which is obtained from $`P`$ by transposing its first tensor factor. The following expression gives the corresponding deformed rational solution of the YBE:

$$ud^v\left(_{21}^1\right)+P\frac{u}{u1+M/2}d^v\left(_{21}\right)Kd^v\left(^1\right)$$

Here $`u`$ is a spectral parameter. In  such deformed solutions were obtained in the explicit form for the canonical chains $`=_{_{0p}}`$ .

All the calculations can be reproduced for the twisting elements $`=_{𝒢_{0p}}`$ of the full chains. This will lead to a new set of so called deformed Yangians .

This work was partially supported by the Russian Foundation for Basic Research under the grant 00-01-00500 (VDL) and 98-01-00310 (PPK).


text:
# ADONIS observations of the HD 100546 circumstellar dust disk 11footnote 1based on observations collected on the ESO 3.6m telescope at La Silla, proposal 63.H-0239

## 1 Introduction

The discovery of infrared excesses around an important fraction of main-sequence stars (Aumann et al. 1984; Plets & Vynckier 1999) and attributed to the presence of cool dust grains orbiting the star and geometrically arranged in a disk, has triggered a lot of studies because these disks may be linked to planetary formation. For a long time the disk around the star $`\beta `$ Pictoris was the only example to be resolved both in the visible (Smith & Terrile 1984) and in the mid-infrared range (Lagage & Pantin 1994). Recent discoveries of young “debris” circumstellar disks around relatively “old” and isolated (i.e. not associated to any star forming region) stars have shown that this phenomenon extends towards Pre-Main-Sequence stars (cf the photometric survey by Malfait et al. 1998a). The detection and resolution of disk around stars such as HD 141569 (Weinberger et al. 1999; Augereau et al. 1999) or HR 4796A thanks to high-resolution observations in the visible/near-infrared range (Schneider et al. 1999; Augereau et al. 1999) or mid-infrared images (Koerner et al. 1998; Jayawardhana et al. 1998) have shown the possibility to observe and study precursors to main-sequence dust disks. These so-called “baby-$`\beta `$ Pic” dust disks are the denser precursors to main-sequence debris disks. IRAS and ISO/SWS observations have shown that they usually produce a huge infrared excess (typically 250 times the infrared excess produced by the $`\beta `$ Pic disk); some of them, as HD 100546 for instance, showing prominent signatures of crystalline water ice (Waelkens et al. 1996, Malfait et al. 1998b), are particularly interesting targets in the visible/near-infrared because of high particle albedo. Observing this class of disks at various stages of evolution will help in finding a comprehensive scenario for the origin, the evolution, and the lifetime of the Vega phenomenon, i.e. how these disks form, how long they last, and how they disappear.
We report in this paper the first images of the disk around HD 100546, an isolated young main-sequence star. Using a model of scattering by dust grains, we derive the morphology of the disk outwards 10 AU and show evidence for a density maximum around 40 AU. The structure of the disk is compared to the best known example of the Vega Phenomenon, the disk around $`\beta `$ Pic.

## 2 Observations

We observed the HD 100546 star (m$`_\text{V}`$=6.7) using the SHARP II+ camera coupled with the adaptative optics system ADONIS, and mounted on the ESO 3.6m telescope at La Silla, Chile. A pre-focal optics coronograph (Beuzit et al. 1997) was used to reject the direct starlight and increase the integration time in each elementary exposure. The mask size we used has a diameter of 0.74″. J (1.25 $`\mu `$m) and Ks (2.15 $`\mu `$m) band exposures were obtained on nights 23 and 24 of June, 1999. The seeing was quite variable the first night, but had decreased to a value of $`0.7^{\prime \prime }`$ at the time of the J and Ks observations. It was stable at a value around $`0.85^{\prime \prime }`$ during the second night. We spent a total observing time of 1550 s in the Ks band and 2600 s in the J band. The Point Spread Function (PSF) was frequently monitored thanks to interlaced observations of two reference stars, HD 97218 and HD 101713. These reference stars were used later in the reduction process to remove the wings of the stellar PSF and for photometric calibration. Observations of empty fields were also performed to estimate and remove the background flux which can be important in the Ks band.

## 3 Reduction procedure

First, standard reduction techniques were applied to the data : bias subtraction, flat-field correction. For each filter, we obtained a set of HD 100546 observations and corresponding PSFs. In spite of the use of a coronograph mask, the disk surface brightness is still dominated by the stellar emission at any distance from the star. In order to retrieve the disk image, one has to remove numerically the starlight wings. The rough substraction of a scaled PSF to HD 100546 images gives poor results because of slight shifts in position on the array between the reference and the object, uncertainties on the fluxes (given by the literature), and a residual background (ADONIS bench emission, different airmasses). We developed a specific method to get an optimum subtraction. For each couple of HD 100546 image (Obj) and corresponding Psf, we have to find 3 parameters : a shift ($`\delta x`$,$`\delta y`$) between the two images, a scaling factor $`R`$, and a residual background $`Bg`$. These parameters are estimated by minimizing the following error functional :

$$J=|\text{Obj}S(\text{Psf}/R,\delta x,\delta y)Bg|$$
(1)
where the S function stands for image shift. The sum is performed on a set of pixels (typically 1000) located in a region of the images where no disk emission is expected, and from which non accurate ones are excluded (bad pixels, pixels belonging to area contaminated by diffracted light from the coronograph support). The functional minimum is found using a 0-order minimization algorithm called Powell method (Press 1996). The method was first checked on simulated data whose input parameters (shifts, scaling factor) were recovered with an accuracy better than 5 %. In order to evaluate the errors in the subtraction process, test the stability of the PSF, and make sure that the disk obtained is not an artifact created by the process, the same substraction process was applied to the two different PSFs. Figure 1 shows the resulting images in the J and Ks filters of the disk and the corresponding residuals between the two reference stars. The position of the star under the mask, and therefore the disk center is also carefully determined in this procedure thanks to offset observations (out of the coronograph mask) of the reference stars and the object.

## 4 Results

Since the dataset in J band is more accurate (because of a larger integration time, better observing conditions, and a better sensitivity of the whole system) than in Ks band, the geometrical parameters and the disk model are based on the J band image. We checked in a second step that the results are compatible with the Ks band data.

### 4.1 Disk geometrical parameters

Given a disk inclination with respect to the line of sight which is $``$ 45, we expect to observe some deviations from the perfect ellipse case. Indeed, particles relatively large with respect to the wavelength produce non-isotropic scattering. One simple way to reduce this effect is to symmetrize the disk i.e. construct an isotropic-like virtual disk by averaging the original disk image and its symmetrical image with respect to the major axis. Best ellipse fitting on this symmetrized images gives the following results:

* disk position angle on the sky of $`37\pm 5`$ degrees with respect to east direction.
* ellipse major to minor axis ratio of 1.2 implying a disk inclination of $`50\pm 5`$ degrees with respect to the line of sight (or from the edge-on case).
* disk extension : $``$ 2.0 arcsec from the star corresponding to a distance to the star around 200 AU (assuming a distance to the Earth of 103 parsecs (van den Ancker et al. 1998) and a limiting sensitivity of 5 mJy/arcsec2.

### 4.2 Photometry

Using the stars HD 101713 (B9V) and HD 97218 (G8) as references whose V magnitude is known and using appropriate V-J and V-Ks color indices, we derived the following photometry for the disk in J and Ks band: a disk flux in J band of $`0.32\pm 0.06`$ Jy outside the coronograph mask and a maximum surface brightness of $`0.5\pm 0.1`$ Jy/arcsec2. From images of HD 100546 and corresponding reference star images (HD 101713 and HD97218) obtained when shifting the targets outside of the mask, we can deduce a total of scattered flux in the J band around 1 Jy; note that we get also a disk-shaped emission when removing a scaled reference to the HD 100546 image obtained outside the mask, but its signal to noise is very low. This result is compatible with the infrared excess computed considering a Kurucz model of the star (B9V spectral type). Consequently, the disk flux outside the coronograph mask represents roughly 30% of the total disk emission.
In the Ks band, we find an emission of $`0.25\pm 0.1`$ Jy outside the coronograph mask and a maximum surface brightness of $`0.25\pm 0.1`$ Jy/arcsec2

### 4.3 Model and observation fit

In order to derive physical parameters of the disk, we built a numerical model of starlight scattering by an optically thin disk made of dust particles. Assuming that no multiple scattering is occurring in the disk, and since the exact physical properties of the particles are poorly known (composition, shape, size distribution), we chose to use a global, single-parameter, scattering phase function. The most practical one, although not realistic, is the Henyey-Greenstein phase function (Henyey & Greenstein 1941). The disk is assumed to be seen with an inclination derived above, and to have a proper thickness similar to the one of $`\beta `$ Pic parameterized following Artymowicz et al. 1989. The disk midplane density is described by a scattering area which follows a series of broken power laws. We start with a first guess deduced from the radial surface brightness, and we iterate on the parameters until we obtain a satisfactory fit along the four directions defined on figure 2.

### 4.4 Fit results

The best fit model provides the following parameters (one must keep in mind that they are model-dependent) :

* The Henyey-Greenstein phase function parameter found is $`g=0.2`$ corresponding to dominant spherical particles with size around 0.1 micron.
* A radial normal optical thickness fitted by a series of broken power laws whose exponents are : 0.6 in the range AU, $``$0.4 in the range AU, $``$1.35 in the range AU, $``$5 outwards from 100 AU, see figure 3.
* An exponentially decreasing vertical profile, and an opening angle around 0.1 radian.
* A total scattering area of 1029 cm2 and a dust mass of 0.02 M.
* An approximate disk surface density at maximum around $`3\times 10^{14}`$ particles/m2 assuming single-sized particles of 0.1 $`\mu `$m, or an equivalent maximum scattering area of the order of the unity/m2.

## 5 Discussion

The radial surface brightness distribution shows that the disk density peaks at a distance of 40 AU from the star corresponding to the edge of the planetary zone in our solar system. The normal optical thickness found is more than 100 times larger than in the case of the $`\beta `$ Pic disk. Simultaneously, the presence of sub-micronic particles sizes, well below the radiation pressure cut-off in sizes, suggests that huge quantities of very small particles are produced continuously by evaporation or collisions of planetesimals. We are probably witnessing a analogue of the young Kuiper belt with planetesimals undergoing violent agitation. If the dust particles are the result of collisions, a high rate is required to produce the observed disk. This implies that gravitational perturbations occur frequently in this disk, maybe because of giant planets. If the disk is replenished by planetesimal evaporation as suggested by Lecavelier des Etangs (1998) in the case of the $`\beta `$ Pic disk, 109 planetesimals (40 km size) whose orbits were recently perturbed by planet migration are required to maintain the disk. Some of them which orbits are gravitationally perturbed might be the origin of the casual redshifted absorption profiles observed in some of the spectra of HD 100546 (Grady et al. 1997), a phenomenon highly similar to what is observed in the system of $`\beta `$ Pictoris (Lagrange-Henry et al. 1987, Beust et al. 1996).

## 6 Conclusion

Thanks to coronograph adaptive optics observations we could detect the disk surrounding the PMS star HD 100546. This detection of a disk which is believed to represent a class of precursors to $`\beta `$ Pic-like disks adds a new clue in the understanding of the formation and evolution of these debris disks. Other imaging observations of disks at various stage of evolution are necessary to sketch more closely the possible scenario of such a disk creation and evolution. This should be achieved by next generation instruments such as the VLT/NAOS (adaptive optics) and VLT/VISIR (mid-infrared) instruments or the NGST.

###### Acknowledgements.

We would like to express our warm gratitude to the persons who helped in achieving these observations : D. Le Mignant, F. Marchis, O. Marco, V. Meri$`\stackrel{~}{n}`$o, E. Wenderoth, and the whole 3.6m team.