text:
# Ferromagnetic-spin glass transition in four-dimensional random-bond Ising model

## Introduction

The study of how order arises in nature is in the center of interest in several areas of physics, especially in thermodynamics and statistical physics. A special focus is on the question, how order can emerge, even in or due to the presence of quenched disorder. One type of system, which has attracted a lot of attention during the last decades, are spin glasses.

In this work, the transition from one ordered phase, the ferromagnetic phase, to another ordered phase, the spin-glass phase, is studied for a four-dimensional random-bond Ising model. The model is studied by means of ground-state calculations, using the genetic cluster-exact approximation (CEA) method. This approach has the advantage that one does not encounter ergodicity problems or critical slowing down like when using algorithms which are based on Monte-Carlo methods. To the author’s knowledge, there are no numerical studies of the four-dimensional random-bond Ising model (except at $`p=0.5`$). However, for two and three dimensions, the random-bond Ising model has been studied numerically at $`T=0`$ already, resulting in critical values $`p_c^{\mathrm{sq}}=0.10`$ for square lattices respectively $`p_c^{cub}=0.22`$ for cubic lattices.

The model treated here consists of $`N=L^4`$ Ising spins $`\sigma _i=\pm 1`$ on a hypercubic lattice. The Hamiltonian of this system is given by

$$H\underset{i,j}{}J_{ij}\sigma _i\sigma _j$$
(1)
where $`\mathrm{}`$ denotes a sum over pair of nearest neighbors. Systems with quenched disorder of the bonds $`J_{ij}=\pm 1`$ are studied. The ferromagnetic-spin glass transition is driven by increasing the fraction $`p`$ of antiferromagnetic ($`J_{ij}=1`$) (AF) bonds. For low concentration $`p`$, the system is ferromagnetically ordered, while for intermediate values of $`p`$, spin-glass order arises. For very high concentration, the system shows AF ordering.

In previous work, the model has been studied using high-temperature series expansions. In Fig. 1 the schematic phase diagram is presented. The line connecting the multicritical point $`(p^{},T^{})`$ with $`(p=p_c,T=0)`$ starts vertically at the multicritical point. Therefore, $`p_cp^{}`$ can be expected. The values previously obtained for the concentration at the multicritical point are $`p^{}0.3`$ (Ref. ), $`p^{}=0.28(1)`$ (Ref. ) and $`p^{}=0.290(1)`$ (Ref. ).

The rest of the paper is organized as follows: Next, the algorithm to generate spin-glass ground states in four dimensions is briefly explained. Then, the results describing the ferromagnetic-spin glass transition are presented and finally a summary is given.

## algorithm

From the computational point of view, the calculation of spin-glass ground states is very demanding, because it belongs to the class of the NP-hard problems . This means that only algorithms are available, for which the running time on a computer increases exponentially with the system size. Only for the special case of a planar system without magnetic field, e.g. a square lattice with periodic boundary conditions in at most one direction, exist efficient polynomial-time “matching” algorithms.

For higher dimensions, one must rely on exact methods like branch-and-bound  or branch-and-cut, which are able to treat only small systems. For that reason, recently some heuristics have been introduced. By using a hierarchical approach one can calculate true ground states in four dimensions up to size $`L=5`$, while for $`L=6`$ it is not clear whether true ground states were found.

The method applied here, is able to calculate ground-states in four dimensions up to size $`L=7`$. The technique is based on a special genetic algorithm  and on cluster-exact approximation  which is an optimization method designed especially for spin glasses. Now a brief description of the techniques are given.

Genetic algorithms are biologically motivated. An optimal solution is found by treating many instances of the problem in parallel, keeping only better instances and replacing bad ones by new ones (survival of the fittest). The genetic algorithm starts with an initial population of $`M_i`$ randomly initialized spin configurations (= individuals), which are linearly arranged in a ring. Then $`\nu \times M_i`$ times two neighbors from the population are taken (called parents) and two offspring are created using the so called triadic crossover . Then a mutation with a rate of $`p_m`$ is applied to each offspring, i.e. a fraction $`p_m`$ of the spins is reversed. Next, for both offspring the energy is reduced by applying the CEA. This algorithm constructs iteratively and randomly a non-frustrated cluster of spins, whereas spins with many unsatisfied bonds are more likely to be added to the cluster. The non-cluster spins act like local magnetic fields on the cluster spins. For the spins of the cluster, an energetic minimum state can be calculated in polynomial time by using graph-theoretical methods : an equivalent network is constructed , the maximum flow is calculated  and the spins of the cluster are set to the orientations leading to a minimum in energy. This minimization step is performed $`n_{\mathrm{min}}`$ times for each offspring.

Afterwards each offspring is compared with one of its parents. The pairs are chosen in the way that the sum of the phenotypic differences between them is minimal. The phenotypic difference is defined here as the number of spins where the two configurations differ. Each parent is replaced if its energy is not lower (i.e. better) than the corresponding offspring.

After this creation of offspring has been performed $`\nu \times M_i`$ times, the population is halved: From each pair of neighbors the configuration which has the higher energy is eliminated. If not more than 4 individuals remain, the process is stopped and the best individual is taken as result of the calculation.

The whole algorithm is performed $`n_R`$ times and all final configurations which exhibit the lowest energy are stored, resulting in $`n_g`$ statistical independent ground-state configurations. By comparison with an exact method, it was shown that genetic CEA, with an appropriate choice of the parameters ($`M_i`$, $`\nu `$, $`n_{\mathrm{min}}`$, $`p_m`$), indeed calculates true ground states.

The probability that a certain ground-state configuration is found by this method is not equal for all ground states , i.e. the algorithm imposes a bias. In this work, the magnetization is the main quantity of interest. To test, how large the influence of this bias is, for $`L=6`$, $`p=0.27`$, where fluctuations of the magnetization are the largest, 64 realizations were considered. For each realization, (biased) ground states were generated using 200 independent runs. Next, for each realization a (thermodynamically correct) set of ground states was generated in such a way that each configuration contributes with the same probability. For both sets of ground states the magnetization was evaluated. For the biased set an average value of $`m=0.501(11)`$ was obtained, while the correct thermodynamic result is $`m=0.504(12)`$. This shows that the influence of the bias on the result of the magnetization is small. Please note that that is in contrast to other quantities like the distribution $`P(q)`$ of overlaps: The overlap $`q^{\alpha \beta }`$ between two independent ground states $`\{\sigma _i^\alpha \}`$, $`\{\sigma _i^\beta \}`$ is given by $`q=\frac{1}{N}_i\sigma _i^\alpha \sigma _i^\beta `$. When measuring the fraction $`x(0.5)=_{0.5}^{0.5}P(q)𝑑q`$ of small overlaps, one obtains an average value $`x(0.5)=0.051(7)`$ for the biased set of states, while the from the correct thermodynamic treatment one obtains quite a different result $`0.014(8)`$.

Hence, when restricting the measurement to quantities depending solely on the magnetization, it is sufficient to use the biased data. This allows much shorter running times, because only few ground states per realization are needed.

## Results

With the genetic cluster-exact approximation, it it possible to obtain ground-states for four-dimensional systems up to $`L=7`$. Since the study of the transition involves the calculation of a Binder cumulant, it is necessary to average over many samples of the disorder. For that reason, only sizes $`L6`$ were considered in this work. A mutation rate of $`p_m=0.05`$ was used. The other simulation parameters are shown in Tab. I.

Systems with concentrations of the antiferromagnetic bonds in the range $`p[0,0.32]`$ were treated. For each realization $`n_R=5`$ independent runs were performed. The resulting ground-state energy as a function of the parameter $`p`$ is shown in Fig. 2 for different system sizes. For small concentrations $`p`$, the ground state is mainly ferromagnetic. It follows that all the AF bonds are not satisfied, so the ground state energy increases linearly with $`p`$ like $`e(p)4+8p`$. For larger concentrations, the ground-state energy approaches the $`p=0.5`$ limit. With increasing $`L`$ the ground state energy decreases, because the periodic boundary conditions impose less constraints on the system. For $`L\mathrm{}`$ and $`p=0.5`$ a ground state energy of $`e_{\mathrm{}}(0.5)=2.095(1)`$ was found.

From Fig. 2 is clear that the energy as function of concentration is not a suitable quantity for determining the critical concentration $`p_c`$, where the ferromagnetic behavior disappears. For this purpose the Binder cumulant

$$g(p,L)\frac{1}{2}\left(3\frac{[M^4]_J}{[M^2]_J^2}\right)$$
(2)
for the magnetization $`M\frac{1}{N}_i\sigma _i`$ is evaluated. The average $`\mathrm{}`$ denotes the average over different ground states of a realization, while $`[\mathrm{}]_J`$ is the average over the disorder. In Fig. 3 the Binder cumulant is shown for $`L=3,4,5,6`$. To keep the figure clear, only the largest error bar is shown. All curves intersect near $`p_c=0.28\pm 0.01`$. Because of the huge numerical effort, the simulations are restricted to small system sizes, which prevents a more accurate result.

Only few studies concerning the four-dimensional random-bond model have been performed before. Using high-temperature series expansions for the location of the multicritical point $`(p^{},T^{})`$, values of $`p^{}0.3`$ (Ref. ), $`p^{}=0.28(1)`$ (Ref. ) and $`p^{}=0.290(1)`$ (Ref. ) were determined. These results compare well with the value of $`p_c`$ obtained here, indicating that the the line connecting $`(p_c,0)`$ with $`(p^{},T^{})`$ is vertical, or nearly so, as expected.

For the Binder cumulant the following finite-size scaling relation is assumed

$$g(p,L)=\stackrel{~}{g}(L^{1/\nu }(pp_c)).$$
(3)
By plotting $`g(p,L)`$ against $`L^{1/\nu }(pp_c)`$ with correct parameters $`p_c,\nu `$ the datapoints for different system sizes should collapse onto a single curve near $`(pp_c)=0`$. The best results were obtained for $`p_c=0.28`$ and $`1/\nu =1.0`$. In Fig. 4 the resulting scaling plot is shown. It is possible to change the value of $`\nu `$ in a wide range without large effects on the scaling plot. Thus, here a value of $`\nu =1.0(2)`$ is estimated. This is compatible with the value $`\nu =0.8(1)`$ found by a high-temperature series expansion at the multicritical point. Although there is no direct correspondence between the $`T=0`$ ferromagnetic-spin glass transition and the finite temperature ($`p=0.5`$) paramagnetic-spin glass transition at $`T_c=2.03(3)`$, it is remarkable that by Monte-Carlo simulations of four-dimensional spin glasses near $`T_c`$ a similar value of $`\nu =1.0(1)`$ was found. Finally, it should be pointed out, that the result given above is not compatible with the corresponding exponent for four-dimensional bond percolation of $`\nu =0.678(50)`$.

In Fig. 5 the behavior of the average magnetization $`m[M]_J`$ is shown as a function of $`p`$ for different system sizes.

This quantity has the standard finite-size scaling form

$$m(p,L)=L^{\beta /\nu }\stackrel{~}{m}(L^{1/\nu }(pp_c)).$$
(4)
By plotting $`L^{\beta /\nu }m(p,L)`$ against $`L^{1/\nu }(pp_c)`$ with correct parameters $`p_c,\beta ,\nu `$ the datapoints for different system sizes should collapse onto a single curve near $`(pp_c)=0`$. The best result was obtained using $`p_c=0.28`$, $`1/\nu =1.0`$ and $`\beta /\nu =0.4`$. It is shown in Fig. 6 for $`L=3,4,5,6`$. From variations of the value $`\beta /\nu `$, the value of the exponent $`\beta =0.4(1)`$ is estimated. Again, this values significantly differs from the exponent $`\beta =0.693(20)`$ found for the order parameter of the four-dimensional bond-percolation.

## Summary

The ferromagnetic-spin glass transition of the four-dimensional random-bond Ising model was studied at $`T=0`$ using ground-state calculations. The genetic cluster-exact approximation method was applied. Because of the high computational effort, only small systems of size $`N6^4`$ could be studied. Since the ground-state problem is NP-hard, it is very unlikely that significantly larger sizes can be studied in the near future.

By evaluating the Binder cumulant, a critical concentration of $`p_c=0.28(1)`$ was found, which is comparable to values which were obtained by high-temperature series expansions.

Using a finite-size scaling analysis, the critical exponents for the divergence of the correlation length $`\nu =1.0(1)`$ and for the magnetization $`\beta =0.4(1)`$ were determined. Hence, the transition is clearly in a different universality class than four-dimensional bond percolation.

## I Acknowledgements

The author thanks A.P. Young for critically reading the manuscript and various other support and R. Fisch for helpful communications. The work was supported by the Graduiertenkolleg “Modellierung und Wissenschaftliches Rechnen in Mathematik und Naturwissenschaften” at the Interdisziplinäres Zentrum für Wissenschaftliches Rechnen in Heidelberg and the Paderborn Center for Parallel Computing by the allocation of computer time. The author acknowledges financial support from the DFG (Deutsche Forschungsgemeinschaft) under grant Ha 3169/1-1.


format: text

text:
# Neutrino cross-section measurement with neutrinos from muon decay

## 1 The neutrino cross-section measurements

The motivations to measure the neutrino cross-sections at low energies (below 5 GeV) are very strong. For the current generation of large detectors with high accumulated statistics, the systematic errors  start to be a problem, dominating over the statistical errors. At the present time, relatively few measurements of $`\nu _\mu `$ cross-sections at low energies exist (see Figure 1.left) showing, in general poor precisions ($`\mathrm{\Delta }(\sigma _\nu /E_\nu )>`$ 20%) . No measurements at all on $`\nu _e`$ cross-sections have been published up to now.

Several physics analysis could profit from a better knowledge of $`\sigma _\nu `$

* Atmospheric neutrino analysis. The present atmospheric neutrino measurements suffer from a poor knowledge of $`\sigma _\nu `$, which forces the use of double-ratios in order to cancel the systematics (the absolute rates are not reliable). The precision on the measurement of $`\mathrm{\Delta }m_{23}`$ and $`\theta _{23}`$ is limited by systematics on the flux and on the cross-sections. Also, the study whether the atmospheric neutrino anomaly is due to $`\nu _\mu \nu _\tau `$ or $`\nu _\mu \nu _s`$ oscillations is limited by cross-sections uncertainties, and $`\nu _\tau `$ appearance search requires detailed understanding of final state topologies.
* Proton decay searches. A good knowledge of the $`\nu `$ cross-section is required to evaluate in a reliable way the expected backgrounds. An accurate modeling of the background is needed for these channels that require background subtraction.
* Long Baseline neutrinos. One of the main goals of the Long Baseline neutrino programs is the precise measurement of the mixing parameters using appearance experiments, for which it is mandatory a careful understanding of the expected backgrounds.
* Astrophysics. The neutrino-nucleus cross-sections are also required for calculating certain astrophysical processes involving interactions in the outer shells of stars undergoing supernova explosions, which cause further nucleosynthesis to occur, changing the nuclear composition of the star.

In addition, the nuclear targets vastly complicate the description of the neutrino reactions. Effects like the target motion, the binding energy or the final state interactions (specially important in proton decay searches) can depend on the neutrino helicity and on the target type. Even if there is a number of theoretical models trying to account for nuclear effects in neutrino interactions, only a few comparisons with experimental data are available.

For all these reasons, it is clear that a new and precise measurement of the neutrino cross-sections is needed.

## 2 Measurement of $`\sigma _\nu `$ : The idea and the assumptions

The idea proposed here is based on the use of a low energy muon beam as source of neutrinos. This has two clear advantages compared to the traditional pion beams: in one hand, it is possible to measure, at the same time, $`\nu _\mu `$ and $`\nu _e`$ cross-sections and both helicities (selecting the muon sign, $`\mu ^{}\nu _\mu \overline{\nu }_ee^{}`$ or $`\mu ^+\overline{\nu }_\mu \nu _ee^+`$); on the other, the measurement is done in a low contamination environment and one profits from the precise knowledge of the muon flux.

The analysis presented in this paper assumes a low energy ($`E_\mu =`$ 2 GeV) and low intensity $`\mu ^{}`$ beam of $`10^{15}`$ useful222muons decaying on the direction of the detector muons/year. The detector is considered to be a cylinder of 1 m radius with a moderate mass of 100 tons (a narrow and long detector clearly benefits from the beam shape characteristics). Because of the limited geometrical acceptance of the detector, it is clear that placing it as close as possible to the muon storage ring will increase the flux of neutrinos traversing the cylinder: the shorter the baseline, the higher the statistics and the smaller the error. Therefore, the baseline is fixed to 10 meters. The length of the straight section of the muon storage ring is fixed to 30 meters.

As shown in section 3, the above conditions allow a measurement of $`\sigma _\nu /E_\nu `$ below 2 GeV with $``$10% error and a bin size of 100 MeV. Many channels can be studied under this assumptions (and this counts for both, $`\nu _\mu `$ and $`\nu _e`$):

1. Charge Current (CC) interactions ($`\nu NlX`$):

* Quasi-elastic (QE) interaction ($`\nu NlN^{}`$)
    * Single/Multi pion production ($`\nu NlN^{}+n\pi `$, $`n>0`$)
    * Resonant ($`\nu NlN^{}`$ or $`\nu N\nu N^{}`$)
    * Deep-inelastic scattering ($`\nu qlq^{}`$ or $`\nu q\nu q`$)
2. Neutral Current (NC) interactions ($`\nu N\nu N+n\pi ^0`$, $`n0`$)
3. Elastic electron scattering ($`\nu e^{}\nu e^{}`$):

The CC (NC) events are characterized by the presence (absence) of a lepton in the final state. Specially important is the measurement of the QE channels, since the $`\nu _\mu `$ QE interaction is the dominant atmospheric reaction in a water Čherenkov detectors like SuperK ($``$60% of the total CC events). A good particle identification and energy resolutions are desired to separate between the charged ($``$27%) and the neutral ($``$10%) single pion production channels. The study of the Resonant interactions is specially relevant since it can be an important source of background for proton decay events. Finally, if the detector energy threshold for electrons is low enough ($`<`$10 MeV), the measurement of the elastic electron scattering cross-sections can be extremely useful for solar neutrino studies. No particle miss-identification is assumed in this analysis.

## 3  Results

As described in section 2, we assume a flux of intensity $`10^{15}`$ useful muons/year of 2 GeV energy. The detector is considered to be a cylinder of 1 m radius with a total mass of 100 tons. The baseline is 10 meters. Under these conditions, the quoted numbers of expected $`\nu _\mu `$ and $`\overline{\nu }_e`$ events/year from charge-current, quasi-elastic and resonant events can be found in Table 1.

The proposed set-up allows a clean measurement of the different neutrino cross-sections with a small error and a high number of data points, as shown in Figure 1. The four plots on the left part of the figure give the relative errors on $`\sigma _\nu /E_\nu `$ as a function of the neutrino energy ($`\overline{\nu }_e`$ and $`\nu _\mu `$) for charge-current and quasi-elastic interactions. The bin size is 100 MeV. Above 0.5 GeV, the largest part of the points exhibit an error smaller than 10%. Figure 1 (right) shows the final obtained $`\sigma _\nu /E_\nu `$ spectrum for the three $`\nu _\mu `$ interactions. The dotted lines are the theoretical predictions from the NUX model.

One attractive feature of the idea launched on this paper is that with moderate means we can do a relevant measurement. The requirements in terms of intensities and energies are far below the neutrino factory needs, for instance (a factor $`10^5`$ a 10 less, respectively). Also the required detector sizes (a factor $`10^2`$ less) and the baseline are much smaller. Specially interesting should be to perform the cross-section measurements using different types of detectors since some nuclear effects (like final state interactions) depend on the nature of the target (water, Iron, Argon, etc.). Moreover, this can allow a test of the full detector (measuring efficiencies, backgrounds, etc.) in real conditions and much before the starting of the neutrino factory program.


text:
# Faster Monte Carlo Simulations at Low Temperatures. The Waiting Time Method

## 1 Introduction

The archetype of optimization approaches inspired by physics is arguably Simulated Annealing (SA) . In SA, putative solutions of the problem at hand are organized in an ‘energy landscape’ similar to the configuration space of a physical system. This landscape is searched for the global optimum using the Metropolis algorithm  while the temperature is slowly decreased during the simulation. The theoretical underpinning of SA is equilibrium statistical physics and the theory of Markov processes: at each temperature $`T`$, the Metropolis algorithm asymptotically samples the equilibrium Boltzmann distribution which, at $`T=0`$, is supported on the set of ground states.

It is well known to practitioners of SA that maintaining thermalization at arbitrarily low temperatures is practically impossible in landscapes with many local minima: The trajectories get trapped in subsets of configuration space (valleys or pockets) surrounding local minima, and most of the attempted moves are rejected. Thus, ergodicity is broken, and the final solution is just a ‘good’ local optimum, not the global optimum.

To mitigate the effects of ergodicity breaking one must make sure that many different valleys are thoroughly explored. Ways to do so is to employ ensemble methods  or use non-monotonic temperature schedules, as done in bouncing , thermal cycling , and tempering . Another possibility is to select candidate moves from a broad distribution , occasionally producing the ‘long jumps’ in configuration space which are needed to escape the neighborhood of a local minimum.

In this paper we discuss and test a ‘thermal’ Monte Carlo scheme, which we have chosen to call the *Waiting Time Method* or WTM. The WTM belongs to a familiy of rejectionless or ‘event driven’ algorithms  which developed from the so-called n-fold method introduced by Bortz et al. . From an *algorithmic* point of view, the WTM is very close to a recent optimization method, Extremal Optimization (EO) , notwithstanding the fact that the latter is inspired by the Bak-Sneppen model of biological evolution , rather than thermal annealing of a physical system.

Rejectionless methods can offer much improved performance, especially in simulations of glassy systems at low temperature. Yet, they do not seem to have gained widespread acceptance . We hope that the present analysis will contribute to improve the situation.

As a prelude, we show the equivalence of the WTM with the Metropolis algorithm with regard to the equilibrium as well as to the dynamical (relaxation) properties. Secondly, we analyze the performance of the WTM and identify the temperature region where it is computationally advantageous to use it. Finally we compare the performance of Extremal Optimization with the WTM run at constant temperature and with the WTM in combination with the Geman-Geman annealing schedule . The numerical analyses all use a 3D Ising spin glass model with short range Gaussian interactions as a test problem.

## 2 The WTM algorithm

The WTM is well suited for problems where $`N`$ variables contribute to the cost function (energy) in an additive fashion. One such case is the short range Ising spin glass model: $`N`$ spins, $`\sigma _i=\pm 1`$ are placed on a cubic lattice, producing a total of $`2^N`$ configurations. Each configuration has $`N`$ neighbors, which correspond to the $`N`$ possible flips of a single spin. The energy of configuration $`\alpha `$ is given by

$$=\frac{1}{2}\underset{i}{}\sigma _i^\alpha \underset{j}{}J_{ij}\sigma _j^\alpha ,$$
(1)
where the couplings $`J_{ij}`$ are independently drawn from a Gaussian distribution of unit variance and zero mean if $`i`$ and $`j`$ point to adjacent lattice sites. Otherwise $`J_{ij}=0`$.

For any configuration111We will often leave the dependence on $`\alpha `$ understood. $`\alpha `$, the WTM associates to each spin a stochastic variable $`t_i`$ which is the time at which this spin must flip, given that its local field $`_jJ_{ij}\sigma _j`$ stays put. A ‘global time’ value, which is initially set to zero, is stored in the variable $`t_{\mathrm{global}}`$.

Let $`\mathrm{\Delta }_i`$ be the energy change associated with the flip of $`\sigma _i`$. Flipping times $`t_i`$ are initially assigned as

$$t_i=\tau _i\mathrm{log}X_i,$$
(2)
where the $`X_i`$ are independent stochastic variables with a uniform distribution in the unit interval. Thus, each $`t_i`$ has an exponential probability distribution with average waiting time given by

$$\tau _i=\mathrm{max}(1,\mathrm{exp}(\mathrm{\Delta }_i/T)).$$
(3)
After the initialization has been completed, the WTM iterates the following three steps:

1. Flip the spin $`\sigma _m`$ having the lowest flipping time.
2. Update the global time: $`t_{\mathrm{global}}=t_\mathrm{m}`$.
3. Generate from an exponential distribution with average $`\tau _i`$ fresh waiting times $`\delta _i`$ for the spin $`\sigma _m`$ and its lattice neighbors (e.g. six other spins in a cubic lattice). The waiting times are added to the global time in order to obtain the updated flipping times: $`t_i=t_{\mathrm{global}}+\delta _i`$.

In order to identify the next move, one needs to search (and update) a database of size $`N`$ storing the current flipping time of each degree of freedom as well as pointers to the $`z`$ degrees of freedom interacting with it. Ordering in a balanced tree (a heap), with the shortest flipping time placed at the root node, requires a computational effort for the search and update procedure of order $`\mathrm{log}N`$. Since the number of updates per spin flip is equal to $`z+1`$, the computational overhead per spin flip is $`(z+1)\mathrm{log}N`$, independent of the temperature $`T`$.

## 3 Theoretical considerations

For completeness, we discuss below a few basic statistical properties of the WTM. The main results are derived analytically and then checked numerically for illustration purposes. We also briefly describe Extremal Optimization, which is used in the performance comparison of the next section.

### Equilibrium distribution for the WTM

It is simple to establish that the WTM fulfills *detailed balance* . Consider two configurations, $`\alpha `$ and $`\beta _k`$, the latter obtained from the former by flipping spin $`k`$. Given that spin $`k`$ is the next spin to flip, the probability that the system remains in state $`\alpha `$ for at least a time $`t`$ is (by construction)

$$q_k(t)=\mathrm{exp}(t/\tau _k).$$
(4)
The rate $`R_{\beta _k,\alpha }`$ of probability flow from $`\alpha `$ to $`\beta _k`$ is then

$$R_{\beta _k,\alpha }=\frac{1}{\tau _k}=\mathrm{min}(1,\mathrm{exp}(\mathrm{\Delta }_k/T)).$$
(5)
Detailed balance follows since $`R_{\beta _k,\alpha }/R_{\alpha ,\beta _k}=\mathrm{exp}(\mathrm{\Delta }_k/T)=P_{\mathrm{eq}}(\beta _k)/P_{\mathrm{eq}}(\alpha )`$, where $`P_{\mathrm{eq}}(\alpha )`$ is the probability of visiting configuration $`\alpha `$ in equilibrium. As a consequence hereof , the stationary distribution of the WTM is the Boltzmann distribution over the states of the system, as in the case of the Metropolis algorithm.

### Metropolis and WTM dynamics

Since detailed balance is fulfilled, the WTM has the same equilibrium (Boltzmann) distribution as the Metropolis algorithm. To an excellent approximation, there is also a *dynamical* equivalence, in the sense that the transition probabilities from one state to another are very nearly the same222If a number of waiting times $`\tau _i`$, which all equal one, happens to enter the heap in close succession, the WTM has a bias towards flipping the corresponding spins in the same order as the $`\tau _i`$ are submitted. This feature hardly has any practical consequences and is neglected in the formulae. in the two schemes. This equivalence simplifies a performance comparison considerably, since one only needs to compare the computer times spent on average to achieve a flip.

The WTM remains in a state $`\alpha `$ for a time at least equal to $`t`$ with a ‘survival’ probability given by

$$Q_\alpha ^\mathrm{W}(t)=\underset{k=1}{\overset{N}{}}\mathrm{exp}(t/\tau _k)=\mathrm{exp}(t\underset{k=1}{\overset{N}{}}\tau _k^1).$$
(6)
The probability that spin $`k`$ be the first to flip is

$$p_k^\mathrm{W}=_0^{\mathrm{}}\tau _k^1Q_\alpha ^\mathrm{W}(t)𝑑t=\frac{\tau _k^1}{_j\tau _j^1}.$$
(7)
In the Metropolis algorithm, the transition probability $`P(\beta _k,\alpha )`$ from a configuration $`\alpha `$ to one of its $`N`$ neighbors $`\beta _k`$ is

$$P(\beta _k,\alpha )=\frac{1}{N}\mathrm{min}(1,\mathrm{exp}(\mathrm{\Delta }_k/T)).$$
(8)
Hence, the probability that spin $`k`$ be the first to flip is

$$p_k^\mathrm{M}=\frac{P(\beta _k,\alpha )}{_jP(\beta _j,\alpha )}.$$
(9)
Since Eq. (7) is identical to Eq. (9), the WTM and the Metropolis method generate the same Markov chain of spin flips.

We also note that the intrinsic (or ‘physical’) time of the WTM approximately corresponds to the number of MC steps in the Metropolis algorithm. The latter leaves state $`\alpha `$ unchanged for at least $`t`$ MC steps with the probability

$$Q_\alpha ^\mathrm{M}(t)=\left(1\underset{k=1}{\overset{N}{}}P(\beta _k,\alpha )\right)^{tN}.$$
(10)
If $`\alpha `$ is a local minimum all the $`\mathrm{\Delta }_k`$ are positive and the sum in Eq. (10) is close to zero for small $`T`$. A similar conclusion applies to configurations close to a local minimum, which is where the algorithm will spend most of the computer time in a multi-minima landscape at low $`T`$. In this situation one easily obtains

$$Q_\alpha ^\mathrm{M}(t)\mathrm{exp}[t\underset{k=1}{\overset{N}{}}\mathrm{min}(1,\mathrm{exp}(\mathrm{\Delta }_k/T))].$$
(11)
Hence, using the usual Monte Carlo step as time unit makes Eq. (6) and (11) identical. Note however that Eq. (11) is an approximation while Eq. (6) is exact.

Fig. 1 shows, for different temperatures, the density of states of the spin glass problem sampled with the WTM (dashed lines), the Metropolis algorithm (dotted lines) and finally the WTM with a different choice of rates, (henceforth ”choice b”) given in Eq. 12 and discussed below. In all cases we used $`10^6`$ MC steps: At $`T=1.2`$ and $`T=1.0`$ this suffices to equilibrate the system (the discrepancies barely visible at high energies are more pronounced for smaller systems (not shown) and are likely to stem from a finite size effect). At $`T=0.8`$ the WTM and Metropolis are still indistinguishable. The incipient lack of equilibration transpires from the fact that the data sampled by the WTM with choice b of rates now visibly differs from the other two sets at high energies. This difference is enhanced at $`T=0.6`$, while at $`T=0.4`$ the trajectories become trapped in different local minima and the three methods strongly differ due to insufficient statistics.

### Other choices of the average waiting time

Within all possible choices of rates fulfilling detailed balance, Eq. (5) stands out as the one which, besides producing the desired equilibrium distribution, also guarantees dynamical equivalence with the Metropolis algorithm. In systems like spin glasses, where the relaxation is dominated by partial equilibrations in configuration space traps, one nevertheless expects the dynamics to be insensitive to the precise choice of rates. To check this, we considered the choice

$$\tau _i=\mathrm{exp}(\mathrm{\Delta }_i/2T),$$
(12)
whose appeal mainly lies in the fact that the same expression for $`\tau _i`$ applies regardless of the sign of $`\mathrm{\Delta }_i`$.

In Fig. 2 we compare the densities of states sampled at $`T=0.8`$ by the WTM with rates chosen according to Eq. 12 and Eq. 5. The latter choice is, we recall, equivalent to the Metropolis scheme. The data shown are averages over $`10000`$ very *short* runs. All these runs started from the same initial state, with energy per spin $`\epsilon =1.6200`$, and lasted $`100`$ MC steps, which is not even remotely sufficient for equilibration purposes. As anticipated, the densities of states sampled by the two methods are quite similar, even at very short time times.

### Extremal Optimization

The WTM is in many respects *algorithmically* similar to Extremal Optimization . We have, unsuccessfully, tried to map the two methods onto each other, and we also failed to produce an analytical description of the stationary distribution belonging to EO. In short, with optimal choice of parameters for optimization, $`P_{\mathrm{EO}}(\epsilon )`$ is significantly broader than $`P_{\mathrm{WTM}}(\epsilon )`$ .

The quality of the contribution to the overall cost from a degree of freedom is considered in EO as a *fitness* measure. The actual way in which the fitness is assigned is not rigidly specified, but reflects one’s intuition and knowledge of the problem at hand. In our case, we use the local energy contribution $`\sigma _i_jJ_{ij}\sigma _j`$ as the fitness $`f_i`$ of $`\sigma _i`$, in accordance with . The following steps are then iterated:

1. Rank the spins according to their fitness, and pick a spin with rank $`n`$ with probability $`pn^\gamma `$. The exponent $`\gamma `$ is a free parameter which can be tuned for a given class of problems but which is not changed dynamically during the runs.
2. Flip the chosen spin and calculate the ensuing change in the fitness of its neighbors.

If $`\gamma `$ is very large, the least fit spin is chosen every time, and EO resembles the Bak-Sneppen model of evolution .

The major difference between EO and WTM is that only the rank (relative magnitude) of the energy changes is important in the former method, while the absolute magnitudes are important in the latter. Secondly, EO chooses its next move probabilistically from a deterministic ranking of the $`f_i`$. By way of contrast, the WTM chooses the next move deterministically from a list of waiting times produced by a stochastic algorithm. Thirdly, the ranking procedure in EO utilizes the current energy contribution, while the WTM looks ahead and generates the waiting times according to the projected energy change.

## 4 Tests and results

This section presents the salient features of the WTM algorithm, found by simulations of the 3D Gaussian spin glass model with $`N`$ spins. We do not attempt a complete characterization, which would be tedious for those not interested in the particular model used. Likewise, the physical aspects of low temperature relaxation in spin glasses will be discussed elsewhere .

### Acceptance rates and intrinsic WTM time

The ‘intrinsic’ or physical time elapsed in a WTM simulation determines the amount of configuration space explored and how close the system is to thermal equilibrium. Its linear relation to the number of Monte Carlo steps is displayed in the left panel of Fig. 3 for three different temperatures and a number of systems sizes. There is no discernible $`N`$ dependence, except at the lowest temperature, where the *smallest* systems, $`N=5^3`$ and $`7^3`$, are seen to follow slightly steeper lines (longer times) .

The right panel of Fig. 3 shows, for three different system sizes, the empirical temperature dependence of the acceptance rate $`a(T)`$ of the Metropolis algorithm, i.e the ratio of the number of flips to the number of attempted flips. Again, the size dependence is negligible.

The slopes in the left panel are, for each temperature, nearly equal to the reciprocal acceptance rate $`1/a(T)`$ obtained from the right panel. This confirms the correspondence between the WTM time $`t`$ and the number of MC steps discussed in Section 3.

### Tuning of parameters

In ergodic systems, and in particular in finite systems, the best-so-far energy per spin $`\epsilon _{\mathrm{bsf}}(t)`$ seen in a simulation of length $`t`$ must approach the ground state energy as $`t\mathrm{}`$, regardless of the temperature. Nonetheless, for realistic values of $`t`$ and all but the smallest systems, the ground state energy is never reached and one can meaningfully study the parameter dependence of a suitably averaged the best-so-far energy $`\epsilon _{\mathrm{bsf}}(t)`$.

In applying the WTM to the spin glass problem, we averaged over many runs with the same length and different $`J_{ij}`$. We observed a clear temperature dependence of the resulting $`\epsilon _{\mathrm{bsf}}`$ and defined an ‘optimal’ parameter $`T_{\mathrm{opt}}`$ as the temperature minimizing $`\epsilon _{\mathrm{bsf}}`$ after $`10^3N`$ spin flips. Empirically, this $`T_{\mathrm{opt}}`$ is *independent* of the runtime through several orders of magnitude . It is however dependent on the system size in the following way: For $`N=5^3,6^3,8^3,10^3`$ and $`12^3`$, we find $`T_{\mathrm{opt}}=1.0,0.9,0.7,0.65`$ and $`0.6`$, respectively.

While we cannot make any theoretical statements regarding EO, the behavior with respect to $`\epsilon _{\mathrm{bsf}}`$ is very much the same, with $`\gamma `$ in lieu of $`T`$. However, running EO on the same systems as above showed a much weaker $`N`$ dependence: though $`d\gamma _{\mathrm{opt}}/dN<0`$ is observed, $`\gamma =1.1`$ produces good results for $`N12^3`$. Compared to the temperature in the WTM, $`\gamma _{\mathrm{opt}}`$ is insensitive to changes in the system size.

### Runtime considerations

A ball park figure for the temperature $`T_b`$ where the WTM is faster than Metropolis can be obtained by noticing that the generation of waiting times, see Eq. (2), is the most time consuming part of the algorithm. For each performed spin flip the WTM generates $`z+1=7`$ new waiting times, while Metropolis requires $`1/a(T)`$ attempted spin flip operations on average. Equating the two expressions yields $`a(T_b)1/70.15`$ and $`T_b1.5`$. The actually measured value of $`T_b`$ is indicated by the lower corner of the box in Fig. 3 and equals $`1.3`$. This is somewhat lower than the estimate, likely because we neglected the computational cost of order $`(z+1)\mathrm{log}N`$, confirmed in Fig. 4, which is needed in the search and update procedure. From Fig. 3 we see that at $`T=0.5`$ the WTM is approximately three times as fast as standard Metropolis. A closer examination of $`a(T)`$ reveals that at $`T0.25`$ and $`T0.035`$ the WTM is 10 and 100 times faster, respectively, if both methods are to produce the same number of spin flips.

### Performance comparisons

We have compared the WTM with an annealing version and with EO. The result is shown in Fig. 5. In the annealing case we have used the Geman-Geman logarithmic schedule $`T(t)=T_0/(1+\mathrm{log}(1+t))`$. We observed the reasonable fact that $`\epsilon _{\mathrm{bsf}}`$ is minimized when the initial temperature $`T_0`$ is chosen so that the final temperature is slightly below $`T_{\mathrm{opt}}`$ in the WTM.

The WTM is seen to be inferior to EO in terms of spin flips needed to reach the same value of $`\epsilon _{\mathrm{bsf}}`$. However, when implementing EO ad literam333which is usually not done , since an exact ordering is costly and an approximate ordering is assumed to work just as well., at equal runtime EO and the WTM produce similar results for $`N=10^3`$, i.e. they reach approximately the same $`\epsilon _{\mathrm{bsf}}`$. For larger systems EO is slower than the WTM .

When applied to an annealing schedule like Geman-Geman, the WTM is seen to perform better than EO, since varying the temperature in the WTM does not change the runtime notably. We note in passing that while the Geman-Geman schedule is usually dismissed as unrealistically slow, it performs quite satisfactorily in connection with WTM.

## 5 Summary and conclusions

The WTM belongs to a familiy of rejectionless algorithms related to the *n-fold way* of Ref. . These algorithms asymptotically sample the Boltzmann distribution and are mathematically equivalent to the ubiquitous Metropolis algorithm. Compared to Metropolis, they also require a computational and implementation overhead. However they are considerably faster at low temperatures, especially in systems with short range interactions. The n-fold way, in particular, is efficient if $`n`$ is small, i.e. if the possible values of the $`\mathrm{\Delta }_k`$ belong to a small, discrete set. The WTM method does not suffer from the same limitation, and is therefore especially well suited for low $`T`$ simulations of the dynamics of disordered and glassy systems.

Our tests on a spin glass system show that Simulated Annealing with Geman-Geman type of schedule is slightly faster Extremal Optimization, a recent addition to the arsenal of optimization tools. This indicates that ‘thermal’ schemes remain a viable approach to global optimization.

## 6 Acknowledgement

We would like to thank Stefan Boettcher for kindly sharing his work and ideas with us. P.S. thanks the participants of the 1999 Telluride Workshop on Energy Landscapes for valuable input and Statens Naturvidenskabelige Forskningsråd for financial support.


text:
# The Superconducting Gap in MgB2: Electronic Raman Scattering Measurements of Single Crystals

## Abstract

Polarisation-resolved Raman scattering measurements were performed on MgB2 single crystals to determine the magnitude, symmetry and temperature dependence of the superconducting gap. A single sharp peak due to Cooper pair breaking appears in the electronic continuum below $`T_c`$, reaching an average maximum Raman shift of $`105\pm 1`$ cm-1 ($`2\mathrm{\Delta }(0)/k_BT_c=3.96\pm 0.09`$) and showing up to 5 cm-1 anisotropy between polarised and depolarised spectra. The temperature dependence of $`2\mathrm{\Delta }`$ follows that predicted from BCS theory, while the anisotropy decreases with decreasing temperature. It is concluded that the Raman results are consistent with a slightly anisotropic s-wave gap in a conventional BCS superconductor.

The recently discovered superconductor MgB2  shows an unusually high $`T_c`$ for a simple binary compound, immediately raising questions regarding the nature of superconductivity in this material . Experiments on polycrystalline samples generally support a gap of s-like symmetry, but measurements of the gap magnitude vary widely. Some groups report a single superconducting gap, with gap-to-$`T_c`$ ratios ranging from 0.6 to 5 , while other groups find multiple gap-like features with gap-to-$`T_c`$ ratios around 1.1–1.8 and 3.5–4.8 . Due to the difficulties and uncertainties associated with polycrystal measurements, these results strongly motivate further study of the gap feature(s) in single crystals.

Electronic Raman scattering is a useful experimental technique for studying superconductors because it provides a direct probe of the superconducting gap via the breaking of Cooper pairs by the incident light . The frequency, polarisation dependence and temperature dependence of the pair breaking peak reflect the magnitude, symmetry and temperature dependence of the gap . Knowledge of these fundamental properties is essential for the underlying mechanism of superconductivity to be deduced. For polycrystalline samples, the MgB2 Raman spectrum presents a broad phonon feature superimposed on a relatively intense electronic continuum , the latter undergoing a superconductivity-induced renormalization which has been interpreted in terms of two gaps of magnitudes ($`2\mathrm{\Delta }`$) 44 cm-1 and 100 cm-1 . In this letter we report the first Raman measurements of the fundamental gap properties in single crystals of MgB2.

Single crystals of MgB2 were synthesized from a precursor mixture of 99.9% pure magnesium powder and 97% pure amorphous boron. The crystals were grown in a BN container under 4–6 GPa pressure at 1400–1700C for 5 to 60 minutes. Magnetization and resistivity measurements give $`T_c=38.2`$ K with a transition width of 0.3 K . Four circle XRD measurements indicate these are high quality single crystals, with no evidence of bulk impurities. Slab-like crystals with relatively clean and flat in-plane surfaces were selected for measurement. These crystals exhibited broken hexagonal shape of typical dimension 200–300 $`\mu `$m square by 30–40 $`\mu `$m thick. The crystals were mounted for in-plane measurement and oriented with the a-axis horizontal, within $`\pm 5^{}`$, confirmed by x-ray diffraction measurements. Further confirmation was provided by the results of superconducting-state Raman measurements. The incident light was polarised either vertically (V) or horizontally (H), that is, either perpendicular or parallel to the a-axis. Only vertically polarised scattered light was collected. In the $`D_{6h}`$ point group, VV polarisation selects excitations of $`A_{1g}`$ and $`E_{2g}`$ symmetry, while HV polarisation selects only the $`E_{2g}`$ component. Henceforth the VV and HV configurations will be referred to as *polarised* and *depolarised*.

Raman spectra were measured with a Jobin-Yvon T64000 triple monochromator, charge coupled device, and an excitation wavelength of 514.5 nm provided by an Ar-Kr laser. The typical incident power was 0.8 mW, point-focussed to a spot of 100 $`\mu `$m diameter. To determine the degree of laser overheating, spectra were also measured at laser powers higher than typical. The only noticeable effect of higher incident powers was spot heating, and where appropriate subsequently quoted temperatures have been corrected for this effect. Temperature control was provided by a closed-cycle helium refrigerator and at typical laser powers sample overheating was estimated to be less than 2 K. All spectra have been corrected for the Bose thermal contribution.

Polarised and depolarised Raman spectra of a single crystal of MgB2 (crystal A) from 15 to 1300 cm-1 are shown in Fig. 1. The spectral intensity has been normalised at high wavenumbers and the polarised spectra

offset by 0.5 units along the y-axis. Dominating the spectra is a very broad, asymmetric feature located at around 620 cm-1, which is commonly attributed to the Raman-active E2g symmetry phonon in MgB2 . This feature narrows substantially with decreasing temperature, consistent with suggestions that the E2g phonon is very anharmonic , but curiously shows no indication of a large superconductivity-induced hardening below $`T_c`$ predicted by band theoretical calculations . Cursory modelling of this feature with a Fano equation, an unsatisfactory description of the observed spectral lineshape, revealed 7–10 cm-1 hardening between 40 and 21 K. Further experiments must be performed to clarify the character of this mode.

In contrast, a clear superconductivity-induced renormalization of the electronic continuum is visible in the 21 K spectra, where a sharp pair breaking peak forms at around 100 cm-1. The magnitude of the renormalization is greater in depolarised spectra. Accompanying the peak is a decrease of scattering intensity at lower Raman shifts, although some scattering intensity remains directly below the peak. In the case of an isotropic s-wave gap the theoretically calculated electronic Raman spectrum is characterised by a pair breaking peak located at an energy of $`2\mathrm{\Delta }`$ in all scattering channels, below which there is zero scattering intensity . For anisotropic s-wave the scattering threshold remains, located at the minimum gap energy, while pair breaking peaks appear at different Raman shifts in different polarisations, reflecting the gap anisotropy . Since the polarisation dependence of the peak seen in MgB2 is not strongly anisotropic, we directly associate the pair breaking peak energy with $`2\mathrm{\Delta }`$.

Figure 2 shows Raman spectra from crystal A in the vicinity of the pair breaking peak. The spectra have been normalized to unity around 350 cm-1 and successive spectra have been progressively offset by one unit along

the y-axis. Where appropriate, temperatures have been corrected for sample heating effects. Here, the evolution of the $`2\mathrm{\Delta }`$ peak with temperature is seen and a small anisotropy in peak energy between polarised and depolarised spectra is apparent. In the vicinity of 25 cm-1 the Raman intensity at low temperatures nearly vanishes. To obtain an estimate of the gap energy $`2\mathrm{\Delta }(T)`$, the spectrum in the vicinity of the pair breaking feature was approximated by a Gaussian with a cubic background, and fitted with a least squares algorithm. Subsequently quoted uncertainties in $`2\mathrm{\Delta }(T)`$ are the 95% confidence intervals obtained from the fit. The average anisotropy between the two polarisations thus obtained is $`4\pm 1`$ cm-1, the polarised $`2\mathrm{\Delta }`$ peak always appearing higher in energy than the depolarised peak. Close examination of Fig. 2 suggests that the anisotropy is about 5 cm-1 at temperatures close to $`T_c`$, and falls with decreasing temperature to a little less than 3 cm-1 at the lowest temperatures measured.

A number of single crystals were measured and low-frequency Raman spectra from a second crystal, crystal B, are shown in Fig. 3. In common with crystal A, a pair breaking peak appears in the superconducting state spectrum and $`2\mathrm{\Delta }(T)`$ agrees well between crystals. Anisotropy in $`2\mathrm{\Delta }`$ is also present, amounting to $`2\pm 1`$ cm-1 on average, about half that seen in crystal A. Some additional structure is also apparent. A sharp peak present in most crystal B spectra near 300 cm-1, but not present in crystal A, is attributable to surface contamination. Although the polarised spectra are comparable between crystals, additional scattering intensity is present in the depolarised spectra of crystal B. At low temperatures a

sharp shoulder forms around 30–40 cm-1, which is very similar to that seen in polycrystal results . Residual scattering intensity in the vicinity of 25 cm-1 is non-vanishing in crystal B and the magnitude of the superconductivity-induced renormalization in the depolarised spectra is reduced compared to crystal A. Inspection of crystal B with an optical microscope showed a less regular surface and confirmed the presence of surface impurities. Typically, crystals with flat, clean surfaces showed a sharp $`2\mathrm{\Delta }`$ peak with low residual scattering while crystals with evidence of surface disorder showed higher residual scattering and a shoulder at 30–40 cm-1.

Within the usual theoretical description of non-resonant electronic Raman scattering, the square of the Raman vertex $`\gamma _𝐤^2`$ weights scattering from different portions of the Fermi surface depending on the incident and scattered light polarisations . Taking the band dispersion given by the nearest neighbour tight binding model  we made a straightforward calculation of $`\gamma _𝐤^2`$ for $`A_{1g}`$ and $`E_{2g}`$ symmetries and found that the former weights excitations from almost the entire hexagonal Brillouin zone, particularly near the vertices, while the latter weights excitations around the edges. Thus, the observation of a smaller gap magnitude in depolarised ($`E_{2g}`$) spectra suggests an anisotropic in-plane gap of the form $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_0[1+ϵ\mathrm{cos}(6\theta )]`$, compatible with s-wave symmetry .

In determining the superconducting gap symmetry, and thus placing constraints on the underlying mechanism of superconductivity, the temperature dependence of the gap provides another convenient test. BCS theory gives an explicit temperature dependence for a s-wave gap with which the experimental data may be compared.

We calculate the reduced gap $`2\mathrm{\Delta }(T)/2\mathrm{\Delta }(0)`$ from our data, extrapolating $`2\mathrm{\Delta }(0)`$ from our lowest temperature measurements. Results of this calculation are shown in Table I.

The magnitude of the superconducting gap energy, averaged between polarisations, is $`2\mathrm{\Delta }(0)=3.96\pm 0.09`$, which is consistent with MgB2 being a moderately strong coupling superconductor. Figure 4 shows $`2\mathrm{\Delta }(T)/2\mathrm{\Delta }(0)`$ for both crystals, which follow the BCS-predicted temperature dependence. Uncertainties in temperature are shown for those data points to which correction was applied. The temperature dependence of the superconducting gap, as directly measured by Raman scattering, indicates that MgB2 is a conventional BCS superconductor.

The presence of Raman scattering intensity below 2$`\mathrm{\Delta }`$ seen in Figs. 2 and 3 initially appears to contravene a simple s-wave interpretation of the data. The theory of Raman scattering in superconductors predicts a zero intensity threshold below the pair breaking peak in an s-wave superconductor  and Raman studies of the Nb3Sn and V3Si conventional superconductors  confirm a sharp decrease in scattering intensity below the 2$`\mathrm{\Delta }`$ peak. In contrast, the high-$`T_c`$ superconductors exhibit broad, anisotropic pair breaking features which are accompanied by significant in-gap scattering intensity in all scattering channels, a consequence of the nodes in the d-wave gap function . Nodes in the gap function of MgB2 might explain the in-gap intensity seen below Tc in crystal A (Fig. 2), however this possibility is discounted by the near-isotropic polarisation dependence of the narrow pair breaking peak. While polycrystal measurements show gap-like features below 100 cm-1 , association of the 40 cm-1 shoulder with a second superconducting gap is problematic since it is not seen consistently in both crystals.

A number of extrinsic factors appear to contribute to the in-gap scattering present in Figs. 2 and 3. Rayleigh scattering, impurity and luminescence components, due to surface contamination, were present to varying degrees during the measurements, and the low-frequency tail of the very broad E2g phonon may also furnish scattering intensity at low frequencies. Significant surface disorder may introduce additional scattering contributions and the similar 40 cm-1 shoulder seen in polycrystal  and crystal B spectra may be related to a c-axis component of the superconducting gap. Such surface effects may also contribute scattering intensity to the below-gap Raman spectrum of crystal A, although to a much lesser degree. Clearly the 40 cm-1 shoulder is not an intrinsic feature of the in-plane Raman spectrum, related to a second pair breaking peak. Rather it appears to be an extrinsic effect related to sample surface quality.

In summary, we have measured the temperature and polarisation dependence of the $`2\mathrm{\Delta }`$ pair breaking peak seen in the in-plane electronic Raman continuum of MgB2 single crystals. This peak shows BCS-like temperature dependence with only slight anisotropy, arguing strongly for conventional s-wave superconductivity. The average gap to $`T_c`$ ratio $`2\mathrm{\Delta }/k_BT_c=3.96\pm 0.09`$ indicates moderately strong electron phonon coupling. We are able to identify only one gap feature in the in-plane Raman spectra, and attribute structure below the pair breaking peak seen in the poorer crystal to surface defects.

This work was supported by the New Energy and Industrial Technology Development Organization (NEDO) as collaborative research and development of fundamental technologies for superconductivity applications.


text:
# Oscillatory decay of a two-component Bose-Einstein condensate

## Abstract

We study the decay of a two-component Bose-Einstein condensate with negative effective interaction energy. With a decreasing atom number due to losses, the atom-atom interaction becomes less important and the system undergoes a transition from a bistable Josephson regime to the monostable Rabi regime, displaying oscillations in phase and number. We study the equations of motion and derive an analytical expression for the oscillation amplitude. A quantum trajectory simulation reveals that the classical description fails for low emission rates, as expected from analytical considerations. Observation of the proposed effect will provide evidence for negative effective interaction.

The experimental realization of Bose-Einstein condensates  offers the possibility of studying quantum matter waves. A paradigmatic effect in such systems is the coherent exchange of particles between weakly connected states, first predicted for superconductors . This so-called Josephson effect is expected to be displayed also by a Bose-Einstein condensate (BEC) in a split trap . It has been argued that an “internal” Josephson effect may also occur in condensates whose atoms can reside in two different hyperfine states coupled by Raman transitions . Depending on the relative strength of the interaction between atoms and the hopping between states, the system may lie in the Rabi regime, where atoms behave as independent, or in the Josephson regime, where interactions are sufficiently strong to create a new collective behavior . The coherent exchange of atoms in the Rabi regime has been realized between two hyperfine states of 87Rb atoms .

It is possible to have two-component condensates in which it is energetically favorable for atoms to be in the same state. This “ferromagnetic” coupling results from an appropriate combination of the three relevant scattering lengths , and does not necessarily require a negative scattering length . On the other hand, atom losses by one- and three-body collisions cause the condensate to decay. In this article, we show that the decay of a two-component condensate with ferromagnetic coupling may display spontaneous oscillations due to a fundamental dependence of the energy on the relative particle number.

We consider Bose condensed atoms with two internal hyperfine levels $`|A`$ and $`|B`$. The atoms interact via $`AA`$, $`BB`$, and $`AB`$ elastic collisions. Additionally, a laser field induces Raman transitions $`|A|B`$, generating an internal Josephson effect . If the orbital wave functions of states $`|A`$ and $`|B`$ may be regarded as rigid, the system can be adequately described by the two-mode Hamiltonian

$$H=\frac{\mathrm{}\omega _R}{2}\left(a^{}b+ab^{}\right)+\frac{E_c}{8}\left(a^{}ab^{}b\right)^2,$$
(1)
where the Rabi frequency $`\omega _R`$ depends on the detailed laser tuning, and where the “charging” energy $`E_c=(u_{AA}+u_{BB}2u_{AB})/4`$ encompasses the various interactions. The effective interaction constants are given by $`u_{ij}=(4\pi \mathrm{}^2a_{ij}/m)d^3x|\phi _i(x)|^2|\phi _j(x)|^2`$, where the mode functions $`\phi _i(x)`$, $`i=A,B`$, are normalized to unity. We have assumed $`u_{AA}=u_{BB}`$, so that a contribution linear in the atom number difference vanishes. In the present work we focus on the case where the effective interaction energy for the internal Josephson dynamics is negative ($`E_c<0`$). This is also the case for a double-well condensate with a negative scattering length.

First we analyze the structure of the classical phase space. For that in the Heisenberg equations of motion, we replace the operators $`a`$,$`b`$ by $`\sqrt{n_A}e^{i\varphi _A}`$, $`\sqrt{n_B}e^{i\varphi _B}`$. For the variables $`z(n_An_B)/N2n/N`$ and $`\varphi (\varphi _A\varphi _B)`$ this yields

$`\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{z}{\sqrt{1z^2}}}\mathrm{cos}\varphi \lambda z,`$ (2)
$`\dot{z}`$ $`=`$ $`\sqrt{1z^2}\mathrm{sin}\varphi .`$ (3)
Hereafter, $`\lambda NE_c/2\mathrm{}\omega _R>0`$ gives the effective interaction strength, and the time is measured in units of $`1/\omega _R`$. Equations (2),(3) can also be obtained from the classical nonrigid pendulum Hamiltonian

$$H(z,\varphi ,\lambda )=\sqrt{1z^2}\mathrm{cos}\varphi \frac{\lambda }{2}z^2,$$
(4)
where $`(z,\varphi )`$ are canonically conjugate coordinates.

The pendulum Hamiltonian (4) obeys the symmetry $`H(z,\varphi ,\lambda )=H(z,\varphi +\pi ,\lambda )`$ and, thus, a sign change of the effective interaction $`\lambda `$ corresponds to a phase shift $`\varphi \varphi +\pi `$ plus time reversal. This means that the $`\pi `$-states which are found for positive $`E_c`$  at maxima of the Hamiltonian, correspond to minima at $`\varphi =0`$ for negative $`E_c`$.

The classical fixed points of the pendulum Hamiltonian can easily be found by linearizing around its extrema. For $`0\lambda <1`$ one finds the stable solution

$$z_0=0,\varphi _0=0,\mathrm{\Omega }_0=\sqrt{1\lambda },$$
(5)
where $`\mathrm{\Omega }_0`$ is the frequency of small oscillations around the potential minimum. For $`\lambda >1`$ its potential curvature becomes negative. Thus at $`\lambda =1`$ there is a bifurcation whereby $`z_0`$ becomes unstable and splits into two symmetry-related fixed points

$$z_\pm =\pm \sqrt{11/\lambda ^2},\varphi _\pm =0,\mathrm{\Omega }_\pm =\sqrt{\lambda ^21}.$$
(6)
As $`\lambda `$ is proportional to the atom number, a decaying condensate will necessarily cross the bifurcation point $`\lambda =1`$, a process which we study below .

The bifurcation process can be appreciated in the inset of Fig. 1, where $`H(z,0,\lambda )`$ is represented for several values of $`\lambda `$. In the conventional case of $`E_c>0`$, it has been noted  that the point where $`\mathrm{\Lambda }E_cN/2\mathrm{}\omega _R=1`$ marks the limit between the Rabi ($`\mathrm{\Lambda }<1`$) and the Josephson ($`\mathrm{\Lambda }>1`$) regimes. Thus we may view the transition from (6) to (5) as the negative $`E_c`$ version of the crossover from an interaction dominated to a hopping dominated regime. Once $`\lambda <1`$, the dynamics is controlled by hopping and becomes insensitive to the sign of the interaction.

The master equation that results from the inclusion of random losses is such that the classical equations of motion (2),(3) remain formally the same, but the effective interaction constant acquires a time dependence,

$$\lambda (t)=\lambda _0e^{\gamma t},$$
(7)
due to the exponential decrease of the total atom number. We integrate the classical equations of motion (2),(3), and (7) with the initial condition $`z(0)=z_+(\lambda _0)`$, $`\varphi (0)=0`$, and $`\lambda _0=2`$, i.e., we start from a minimum of the pendulum Hamiltonian (4) in the bistable regime, where most of the atoms are in the state $`|A`$. Figure 1 shows the corresponding time evolution of the number imbalance. As the total number decreases due to losses, the two minima of the effective potential merge towards a single minimum at $`z=0`$ (see inset). The system evolves adiabatically until $`\lambda `$ approaches and crosses the critical value of 1. We will see that a departure from adiabaticity necessarily occurs when $`\lambda `$ gets sufficiently close to 1. After crossing the bistability threshold the system finds itself distant from the new equilibrium point and begins to oscillate spontaneously around it. These oscillations initiate in the crossover region between the Josephson and the Rabi regimes.

The oscillation amplitude can be estimated analytically if one assumes (i) that outside the transition region $`\lambda 1`$ the effective potential changes sufficiently slowly and (ii) that the crossover region is traversed fast enough for the process to be described as a sudden change in the potential.

For $`\lambda >1`$ the condition for adiabaticity reads

$$1>\left|\frac{d}{dt}\frac{1}{\mathrm{\Omega }_\pm }\right|=\frac{\gamma \lambda ^2}{\mathrm{\Omega }_\pm ^3},$$
(8)
while later, for $`\lambda <1`$, the condition is

$$1>\left|\frac{d}{dt}\frac{1}{\mathrm{\Omega }_0}\right|=\frac{\gamma \lambda }{2\mathrm{\Omega }_0^3}.$$
(9)
For $`\lambda 1`$ the r.h.s. of the inequalities (8) and (9) become $`\gamma /[2(\lambda 1)]^{3/2}`$ and $`\gamma /2(1\lambda )^{3/2}`$, respectively. This indicates that, sufficiently close to the bifurcation point $`\lambda =1`$, adiabaticity breaks down.

Initially, the system follows the potential minimum $`z_+(\lambda (t))`$ adiabatically until the inequality (8) is violated for the first time. This occurs for

$$z=z_+(\lambda )|_{\mathrm{\Omega }_\pm ^3(\lambda )=\gamma \lambda ^2}\gamma ^{1/3},$$
(10)
where $`\lambda 1`$ has been assumed. As $`\lambda `$ crosses the transition region, we may assume that there is a sudden change of the effective potential in the sense that $`z`$ and $`\varphi `$ are effectively constant. Then the adiabatic condition, now reading as in (9), is recovered for $`\lambda `$ sufficiently below 1.

The resulting dynamics can be conveniently studied by noting that, for a harmonic oscillator, the action $`J=2\pi E/\mathrm{\Omega }_0`$ is an adiabatic invariant . The energy $`E=\frac{1}{2}(\varphi ^2+\mathrm{\Omega }_0^2z^2)`$ is estimated from Eq. (10) leading to $`J=2\pi \mathrm{\Omega }_0z_+^2=\pi \gamma `$. Thus the instantaneous amplitude reads $`[\gamma /\mathrm{\Omega }_0(\tau )]^{1/2}`$ and the imbalance oscillates like

$$z(\tau )\sqrt{\gamma /\mathrm{\Omega }_0(\tau )}\mathrm{cos}[\mathrm{\Omega }_0(\tau )\tau ],$$
(11)
where $`\tau `$ is the time elapsed after the oscillation sets in, and $`\mathrm{\Omega }_0(\tau )=(1e^{\gamma \tau })^{1/2}`$ is the instantaneous oscillation frequency. For sufficiently small $`\gamma `$, $`\mathrm{\Omega }_0(\tau )=(\gamma \tau )^{1/2}`$ which yields, for the time $`\tau _k`$ at which the $`k`$-th maximum of the oscillation happens, the condition $`2\pi k=\mathrm{\Omega }_0(\tau _k)\tau _k\gamma ^{1/2}\tau _k^{3/2}`$. The corresponding amplitudes therefore read

$$z_k(\gamma ^2/2\pi k)^{1/6}.$$
(12)
We have numerically checked the values of $`z_k`$ for $`1k5`$ in the range $`10^5<\gamma <10^1`$, and have found good agreement with this analytical prediction.

The classical treatment above requires a smooth potential, such that its curvature is basically constant on the scale of the width $`\mathrm{\Delta }z`$ of the wave packet. While approaching the transition point $`\lambda =1`$, the potential becomes flatter and, thus, the stationary wave packet spreads. At the same time the distance $`|z_+z_{}|`$ between the two potential minima decreases. As long as

$$\mathrm{\Delta }z<|z_+z_{}|,$$
(13)
the classical description should be adequate. When inequality (13) is violated, a quantum description is necessary. For that we take the variables $`z`$ and $`\varphi `$ to obey the commutation relation $`[n,e^{i\varphi }]=ie^{i\varphi }`$, which for small phase uncertainty amounts to $`[z,\varphi ]=2i/N`$. For the Hamiltonian (4) in the harmonic approximation, this yields the number r.m.s. $`(\mathrm{\Delta }z)^2=1/N\mathrm{\Omega }_+`$.

When $`\lambda 1^+`$, the adiabaticity condition (8) reads $`[2(\lambda 1)]^{3/2}>\gamma `$, while (13) amounts to $`[2(\lambda 1)]^{3/2}>1/4N`$. We expect the classical description of the passage through the bifurcation to be correct if, as $`\lambda `$ decreases in the bistable region, adiabaticity breaks down before classicality does  . This is the case if $`N\gamma >1`$, i.e., if losses are sufficiently intense that more than one atom is lost in a typical Rabi cycle. On the contrary, we expect the classical predictions to fail for $`N\gamma 1`$. A similar conclusion is reached for the regime $`\lambda <1`$ if we replace (13) by the condition that the width $`(N\mathrm{\Omega }_0)^{1/2}`$ of the ground state wave function is smaller than the scale $`\sqrt{1\lambda }`$ below which the effective potential behaves as harmonic. Below we see that these expectations are confirmed by the numerical simulation.

To analyze the quantum dynamics of a single experiment, we employ the quantum jump method  where the time evolution over a sufficiently small time interval $`\mathrm{\Delta }t`$ is simulated as follows. If an atom is lost from condensate $`A`$ ($`B`$), the wave function changes as $`|\psi (t+\mathrm{\Delta }t)Q|\psi (t)`$, $`Q=a,b`$. These losses occur with the probabilities $`p_Q(t,\mathrm{\Delta }t)=\gamma \mathrm{\Delta }t\psi (t)|Q^{}Q|\psi (t)`$. If no loss takes place, the system propagates form $`t`$ to $`t+\mathrm{\Delta }t`$ with the effective Hamiltonian $`H_{\mathrm{eff}}=Hi\mathrm{}\gamma (a^{}a+b^{}b)/2`$. Since the propagation with $`H_{\mathrm{eff}}`$ as well as the atom losses are non-unitary, the wave function has to be normalized after each time step. Note that in the case of two identical loss rates the total atom number commutes also with the effective Hamiltonian and terms in $`H_{\mathrm{eff}}`$ which depend only on the total atom number can still be neglected.

We start our simulation in the bistable region with $`\lambda _0=2`$. Thus the initial state is the “local ground state” at the potential minimum $`z_+`$ and satisfies condition (13). Here we make the very reasonable assumption that the parity symmetry is effectively broken due to the large value of $`N_0`$, i.e., that any Schrödinger cat state has evolved into a mixture due to tiny perturbations from the environment.

Figure 2(a) depicts the quantum mechanical time evolution for the same parameters as used in Fig. 1. Both plots agree quite well, which indicates that the system behaves classically. However, for a smaller value of the decay rate \[Fig. 2(b)\] we find substantial differences: (i) the oscillation starts well before the transition to the monostable regime and (ii) its amplitude is larger than in the corresponding classical case (see Fig. 3). This is so because the finite width of the wave packet allows to cross the point $`z=0`$ when there is still a finite barrier.

A comparison of the classical and the quantum mechanical oscillation amplitude (Fig. 3) shows that the quantum dynamics agrees to the classical dynamics for sufficiently large values of the decay rate $`\gamma `$. We find that the quantum mechanical oscillation amplitude is larger than its classical counterpart for $`\gamma 3\times 10^3`$, i.e., roughly at the value for which the classical description was predicted to fail. Then the amplitude is comparable to the phase uncertainty and therefore difficult to observe. A decay rate $`N_0\gamma 10`$ seems optimal for the experimental observation, since then there are still sufficiently many oscillations before the condensate has decayed and the amplitude is still larger than the uncertainty.

Let us assume that the rigid wave functions $`\phi _{A,B}(x)`$ are given by the Thomas-Fermi (TF) approximation for 1000 87Rb atoms in an isotropic trap with frequency $`100\mathrm{Hz}`$. This results in a TF radius of $`4\mu \mathrm{m}`$ and a harmonic oscillator length of $`1\mu \mathrm{m}`$. An effective negative scattering length $`a_{\mathrm{eff}}a_{AA}+a_{BB}2a_{AB}a_{AA}/1000.5\mathrm{nm}`$ yields $`E_c/\mathrm{}0.3\mathrm{s}^1`$. Then, for $`N_0=10^3`$, to have $`\lambda _0=2`$ requires a Rabi frequency $`\omega _R2\pi 10\mathrm{Hz}`$, the “optimum value” for the decay rate corresponds to a condensate life time of $`(\gamma \omega _R)^11.5\mathrm{s}`$.

In summary, we have shown that a two-component condensate with negative interaction energy ($`E_c<0`$) will exhibit spontaneous oscillations in phase and number during its decay due to atom losses. These oscillations have their root in the fundamental dependence of the Josephson coupling energy on the relative boson number. As the atom number decreases, interaction between atoms becomes less important than hopping between states. Thus the system evolves from a Josephson regime, which is bistable due to the negative sign of $`E_c`$, to a monostable Rabi regime dominated by hopping. A classical description of the oscillatory decay dynamics gives a good quantitative account of the decay process for high enough loss rates. The observation of this effect will provide unambiguous evidence of the Josephson effect with negative $`E_c`$, by identifying the existence of a bifurcation that is manifested through the spontaneous appearance of Josephson-Rabi oscillations. By accurately monitoring the decay process, it will be possible to quantitatively check our current understanding of the internal Josephson dynamics of two-component condensates.

We would like to thank J.I. Cirac and P. Zoller for valuable discussions. This work has been supported by the EU TMR Programme under Contract No. MRX-CT96-0042, by the Dirección General de Investigación Científica y Técnica under Grant No. PB96-0080-C02, and by the Acción Integrada No. HU1998-0003.


text:
# Infinite series of magnetization plateaus in the frustrated dimer-plaquette chain

\[

## Abstract

The dimer-plaquette chain undergoes a first order transition driven by frustration to an orthogonal dimer ground state. We analyze the magnetization curve in this orthogonal dimer regime. Besides of the recently found magnetization plateaus at $`m=1/4`$ and the $`m=1/2`$ we find an infinite sequence of plateaus between $`m=1/4`$ and $`m=1/2`$. These additional plateaus belong to magnetizations $`m(n)=n/(2n+2)`$, where $`n`$ is a positive integer number.

\]

The frustrated dimer-plaquette chain (or orthogonal dimer chain) introduced in  has attracted attention because of its first order quantum phase transition to a product dimer ground state and its close relation to the two-dimensional spin model for SrCu2(BO3)2 . Motivated by the recent discovery of magnetization plateaus in SrCu2(BO3)2 Koga et. al.  investigated the magnetization curve for the one-dimensional model. They found magnetization plateaus at $`m=1/4`$ and $`m=1/2`$ of the full moment near the phase transition point using exact diagonalization and DMRG methods. In order to reproduce their results we studied the magnetization curve using large-size exact diagonalization (up to N=40 sites) and analytical considerations based on the product nature of a certain class of eigenstates of the orthogonal dimer chain . We found clear evidence of further magnetization plateaus between $`m=1/4`$ and $`m=1/2`$. In what follows we will explain our results for parameters $`J=1`$ and $`j=J^{}/J=0.7`$ (for notations see Fig. 1 and Refs. ), i.e., in the zero-field regime the ground state is a dimer product state, but the argumentation is more general and can be applied in principle to the (b)-phase of the phase diagram presented in Fig. 9 of .

To show that we have not only a single jump between the $`m=1/4`$ and $`m=1/2`$ plateaus we need two facts.

First we have to know the lowest energy of eigenstates with $`m=1/4`$ and $`m=1/2`$. The energy for $`m=1/4`$ is exactly known as $`4E/N=(E_13/4)/2=1.25595`$ (see table II), because the corresponding eigenstate is a simple product state. For $`m=1/2`$ the ground-state energy can be estimated very accurately to $`4E/N=0.68059(1)`$ by finite-size extrapolation. The table I shows exact diagonalization (ED) results for $`E(m=1/2)`$ up to $`N=40`$ sites. Obviously, the energy per site increases with increasing chain length for an even number of plaquettes, whereas the energy per site decreases with increasing chain length for an odd number of plaquettes. Both energies meet each other at $`N\mathrm{}`$. We note, that for the considered maximal odd and even chain lengths of $`N=36`$ and $`N=40`$ already the first five digits of $`E(m=1/2)`$ coincide.

Second we have to find at least one state with magnetization $`m^{}`$ between $`m=1/4`$ and $`m=1/2`$ having an energy below the straight line connecting the $`m=1/4`$ and $`m=1/2`$ points in the $`E(m)`$ diagram. Then with increasing field the system would occupy next after the state with $`m=1/4`$ the state with $`m^{}`$ and not that one with $`m=1/2>m^{}`$. According to  low-lying states of certain $`m`$ belong to the class of product eigenstates, for which in every $`(k+1)`$-th plaquette the spins $`\stackrel{}{S}_1`$ and $`\stackrel{}{S}_3`$ form a vertical dimer singlet $`(\stackrel{}{S}_1+\stackrel{}{S}_3)^2=0`$ (i.e., the quantum number $`S_{13}`$ in Koga’s notation is zero). The finite strips between two vertical dimer singlets contain $`k`$ vertical dimer triplets ($`S_{13}=1`$). Since these strips are decoupled from each other we call them fragments of length $`k`$ here where such a fragment consist of $`N_k=4k+2`$ spins. The energies $`E`$ of those eigenstates can be expressed by the energies of the finite fragments $`E_k`$ which can be calculated exactly up to a maximum fragment length of $`k=9`$ (see table II). For instance, the energy $`E`$ of states consisting of $`N/(4k+4)`$ fragments of identical length $`k`$ separated by one vertical dimer singlet each is $`4E/N=(E_k3/4)/(k+1)`$.

From the exact numerical data of table II it can be found that all of these states have energies below the straight line connecting the $`m=1/4`$ and $`m=1/2`$ points in the $`E(m)`$ diagram. That clearly shows, that we have not only a single jump between the $`m=1/4`$ and the $`m=1/2`$ plateau.

Let us now present a more detailed analysis of the magnetization process between $`m=1/4`$ and $`m=1/2`$. We start from the $`m=1/4`$ eigenstate, which is a product state of identical fragments of length $`k=1`$ each having total spin $`S_{k=1}=1`$ leading to a total spin of the chain $`S_{chain}=S_kN/(4k+4)=N/8`$. To increase the magnetization we have to insert one fragment of length $`k=2`$ with total spin $`S_{k=2}=2`$. The next eigenstates with higher magnetizations can be generated by replacing further fragments of length $`k=1`$ by fragments of length $`k=2`$. Due to the product nature of these eigenstates the energy necessary to increase the magnetization is constant and defines the critical field for leaving the $`m=1/4`$ plateau. After all fragments of length $`k=1`$ are replaced by fragments of length $`k=2`$ the system has reached a state with $`S_{chain}=S_kN/(4k+4)=N/6`$, i.e., a $`m=1/3`$ plateau. The energy of this $`m=1/3`$ eigenstate is $`4E/N=(E_23/4)/3=1.076056`$ (see table II). This fragmented structure of the relevant eigenstates is also responsible for the finite-size effects, since for finite chains the splitting of the chain in fragments must be compatible with the total length of the chain. To get higher magnetizations all fragments of length $`k=2`$ and $`S_{k=2}=2`$ must now be replaced successively by fragments of length $`k=3`$ and $`S_{k=3}=3`$ and so on. This scenario generates a series of plateaus at $`m=k/(2k+2)`$ where the width of the plateaus decreases rapidly, because the energy difference between fragments of length $`k`$ and $`k+1`$ is decreasing with increasing $`k`$. The critical field between two plateaus can be calculated from the energies of the fragments $`E_k`$

$$h_c(k)=kE_k(k+1)E_{k1}+3/4.$$
(1)
The energies of the fragments $`E_k`$ can be found in table II. Figure 2 gives an overview about the resulting magnetization curve. In contrast to the corresponding Fig .4 in  we get a staircase like infinite series of plateaus. We emphasize, according to the linear programming theorem (see,  and references therein), fragmented states with nonuniform fragment lengths have higher energies than corresponding states built by fragments of identical length. Furthermore our numerical data clearly show that the energy of the above considered fragmented states is significantly below the energies of competing low-lying non-fragmented states (i.e. states where all quantum numbers $`S_{13}=1`$).

We notice, that our result is in accordance with the general rule of Oshikawa et al.  that $`n(s\overline{m})`$ has to be an integer. In our case the period of the ground state $`n`$ is $`4(k+1)`$, the spin $`s`$ is one half and the magnetization per spin $`\overline{m}`$ is $`k/(4k+4)`$ (the magnetization per site $`\overline{m}`$ corresponds to on half of the quantity $`m`$ used in this comment).

Finally we mention, that the magnetic phase diagram on the (j-h) plane presented in Fig. 9 in  has to be extended taking into account the additional plateaus at $`m=k/(2k+2)`$.

In conclusion, we have found an infinite series of magnetization plateaus in a frustrated one-dimensional quantum spin system which is closely related to the possibility of fragmentation of the considered chain. Though the possibility of fragmentation has been observed also for frustrated two-leg spin ladder , in the spin ladder system only very simple ground states seem to be relevant for its magnetization curve. Hence the sequence of plateaus with rational $`m(k)=k/(2k+2)`$ found in this paper seems to be a novel property of quantum spin systems, not found so far.


text:
# Hopping conductivity in the quantum Hall effect – revival of universal scaling

## Abstract

We have measured the temperature dependence of the conductivity $`\sigma _{xx}`$ of a two-dimensional electron system deep into the localized regime of the quantum Hall plateau transition. Using variable-range hopping theory we are able to extract directly the localization length $`\xi `$ from this experiment. We use our results to study the scaling behavior of $`\xi `$ as a function of the filling factor distance $`|\delta \nu |`$ to the critical point of the transition. We find for all samples a power-law behavior $`\xi |\delta \nu |^\gamma `$ with a universal scaling exponent $`\gamma =2.3`$ as proposed theoretically.

The prominent feature of the quantum Hall effect (QHE) is the emergence of quantized plateaus in the Hall resistance of a two-dimensional electron system (2DES) reflecting the localization of states at the Fermi edge. A deeper understanding of this fascinating phenomenon is gained from regarding the transition region between adjacent QHE plateaus Huckestein (1995); Sondhi et al. (1997). In this regime the dissipative transport is governed by the delocalized states in the vicinity of the Landau level center. The degree of localization is expressed by the localization length $`\xi `$ denoting the typical extension of the electron wave function. For any sample with a finite size $`L`$ theory predicts the conductivity tensor in the plateau transition to follow a scaling function $`f\left(L/\xi \right)`$ with a diverging localization length $`\xi =|\delta \nu |^\gamma `$ Pruisken (1988) where $`\delta \nu =\nu \nu _c`$ denotes the filling factor distance to the critical point $`\nu _c`$ of the transition. The critical exponent $`\gamma =2.3`$ was predicted to be universal, its value determined numerically Huckestein and Kramer (1990) and validated by variation of the sample size Koch et al. (1991a). Experimentally the scaling of the conductivity near the critical point was verified with great success for different samples by temperature, current, and frequency dependent measurements of the transition width Wei et al. (1988); Chow et al. (1996); Engel et al. (1993). These new parameters introduce effective lengths $`L_TT^{p/2}`$, $`L_II^{p/(2+p)}`$, and $`L_ff^{1/z}`$ and thereby add additional exponents $`z`$ and $`p`$. These effective lengths replace the physical sample size $`L`$ in the scaling functions $`f(L/\xi )`$. Recent experiments extended the focus to the transition from the quantum Hall state to the Hall insulator Shahar et al. (1995, 1996); Hilke et al. (1998). In these investigations striking similarities to the scaling behavior in the transition between different QH states were observed Shahar et al. (1997) hinting to the same universality class for both types of transitions Pan et al. (1997); van Schaijk et al. (2000).

In spite of the great success of scaling theory in the QHE there are still some experiments not fitting into the picture of universality. Non-universal exponents were observed in the dependence of the transition width on temperature Koch et al. (1991b), current Chow and Wei (1995) and frequency Hohls et al. (2001a). Other experiments seem to contradict scaling theory at all, both for the Hall plateau-insulator transition Shahar et al. (1998) and the transition between QHE plateaus Coleridge (1999); Balaban et al. (1998).

However, before making conclusions on a general failure of scaling theory it has to be considered that nearly all experiments do not measure the localization length $`\xi `$ directly. Therefore an assumption about the functional form and exponents $`z`$ and $`p`$ of the effective length $`L_{\mathrm{eff}}(T,I,f)`$ has to be made. The non-universal scaling exponents deduced from investigating the QHE transition width as a function of $`T`$, $`I`$ and $`f`$ then only reflect non-universal exponents $`p`$ and $`z`$ in $`L_{\mathrm{eff}}`$. Additionally, the measurements mostly focus on the region close to critical points, where $`\xi `$ becomes larger than $`L_{\mathrm{eff}}`$. In this regime, however, the mechanism of electronic transport at nonzero temperature or frequency is not thoroughly understood. Addressing this problem Shimshoni found in a recent theoretical work for the Hall plateau-insulator transition that for $`T>0`$ quantum transport occurs only at some distance to the critical point, namely in the regime of hopping conductivity Shimshoni (1999). Therefore, any lack of universal width scaling with $`T`$ or $`f`$ does not necessarily allow to draw conclusions on the behavior of $`\xi `$.

In order to avoid such complications, we follow a different approach to scaling. We directly evaluate the localization length in the well understood regime of variable-range hopping (VRH) conductivity Polyakov and Shklovskii (1993a). We have recently shown this method to be reliable for frequency dependent measurements Hohls et al. (2001b). VRH dominates the conductivity at low temperatures, when the localization length becomes much smaller than the effective temperature length $`L_T`$. In the QHE regime the VRH conductivity is given as Ebert et al. (1983); Briggs et al. (1983); Polyakov and Shklovskii (1993a)

$$\sigma _{xx}(T)=\sigma _0\mathrm{exp}\left(\sqrt{T_0/T}\right),k_BT_0=C\frac{e^2}{4\pi ϵϵ_0\xi },$$
(1)
with a temperature dependent prefactor $`\sigma _01/T`$. The characteristic temperature $`T_0`$ is determined by the Coulomb energy at a length scale given by the localization length $`\xi `$, the dimensionless constant $`C`$ being in the order of unity. With this direct access to $`\xi 1/T_0`$ it is possible to test for scaling behavior of $`\xi `$ at the edges of the plateau as long as we are careful enough to stay within the localized regime.

Earlier experiments confirmed the expected temperature dependence of $`\sigma _{xx}(T)`$ and extracted $`\xi `$ in the QHE-regime Koch et al. (1995); Furlan (1998), but either did not analyze the scaling behavior or were restricted to a small filling factor range close to the quasi-metallic regime. In a recent experiment on the QH plateau-insulator transition VRH conductivity following Eq. 1 was also established to that regime and was used to determine $`\xi `$, finding a rough agreement with the prediction of universal scaling Murphy et al. (2000).

In this paper we show that when analyzing the conductivity deep into the localized regime of the quantum Hall plateau transition we find a clear universal scaling behavior of the localization length $`\xi |\delta \nu |^\gamma `$ as function of the filling factor distance $`\delta \nu =\nu \nu _c`$ to the critical point at $`\nu _c`$ with a universal $`\gamma =2.3`$. Such a universality is even observed in samples where non-universal exponents extracted from temperature dependent peak-width scaling are found. Going even further, we show that the conductivities $`\sigma _{xx}(T,\delta \nu )`$ in the VRH regime for $`\delta \nu <0.3`$ can be scaled to a single parameter function of $`|\delta \nu |^\gamma /T`$ as predicted by scaling theory for the critical regime.

The samples used in this work are based on modulation doped GaAs/AlGaAs heterostructures with additional scatterers in the active region of the 2DES Ploog (1987). The scatterers are provided by doping the GaAs close to the heterojunction with Si or Be, either as a $`\delta `$-layer or a weak homogeneous background. This results in relatively low mobilities of a few m2/Vs (see table 1). In order to allow a highly sensitive two-point measurement of very low conductivities the samples were patterned into Corbino geometry using contacts fabricated by standard Ni/Au/Ge alloy annealing. Here the conductivity is given by $`\sigma _{xx}=\left(I/2\pi V\right)\mathrm{ln}\left(r_2/r_1\right)`$ where $`r_1=500\mu `$m and $`r_2=550\mu `$m are the inner and the outer radius of the Corbino disk. The samples were mounted onto the cold finger of a dilution refrigerator with base temperature below 20 mK and positioned into the center of a superconducting solenoid.

The sample conductivity $`\sigma _{xx}(B,T)`$ as a function of magnetic field $`B`$ and temperature $`T`$ was extracted from individually measured $`I`$-$`V`$ characteristics by numerical differentiation in the linear regime close to zero bias. For each curve we identified the linear regime which ranged from a few $`\mu `$V around the QHE plateau transition at low temperatures up to a few mV in the plateau center. This method allowed us to measure the conductivity accurately in a huge range between $`10^{13}\mathrm{\Omega }^1`$ in the plateau center up to $`10^5\mathrm{\Omega }^1`$ in the maximum of the plateau transition.

In Fig. 1a $`\sigma _{xx}(B,T)`$ of sample S1 is shown for the transition between the QHE plateaus at $`\nu =1`$ and $`\nu =2`$. We concentrate our studies on this transition where the spin gap is sufficiently large to exclude activated transport. The critical point $`\nu _c`$ of this transition, listed in table 1, is given by the position of the maximum in $`\sigma _{xx}(\nu )`$.

Before analyzing the VRH transport we tested the conventional scaling behavior of the plateau transition width. Data for the full-width at half maximum $`\mathrm{\Delta }\nu `$ of the peak in $`\sigma _{xx}`$ as a function of filling factor $`\nu =hn/eB`$ are shown in Fig. 1b for sample S1. For all investigated samples $`\mathrm{\Delta }\nu `$ follows a power-law behavior $`\mathrm{\Delta }\nu T^\kappa `$ down to temperatures below 60 mK. However, as already reported in earlier work on similar samples Koch et al. (1991b) the critical exponents $`\kappa `$ as presented in table 1 considerably deviate from the proposed universal value $`\kappa =0.43`$ Wei et al. (1988).

Let us now turn to a closer analysis of the VRH regime. As shown in Fig. 1c for sample S1 the data fit well the predicted temperature behavior of Eq. 1 with a prefactor $`\sigma _01/T`$. We also tested for activated behavior, Coulomb gap with constant $`\sigma _0`$, and Mott hopping $`\sigma _{xx}T^m\mathrm{exp}\left((T_0/T)^{1/3}\right)`$ for various $`m`$, all matching our data worse. From these fits to Eq. 1 we are able to extract the characteristic temperature $`T_0`$ in the VRH regime. To stay well inside this localized regime we only take into account filling factors where the low temperature conductivity is at least two orders of magnitude below the critical conductivity $`\sigma _c`$.

In Fig. 2 the characteristic temperature $`T_0`$ is plotted against the distance $`|\delta \nu |=|\nu \nu _c|`$ to the critical point $`\nu _c`$. The data follows a power-law behavior $`T_0|\delta \nu |^\gamma `$ up to distances as large as $`|\delta \nu |0.3`$. This demonstrates that the localization length $`\xi 1/T_0`$ follows a scaling behavior deep into the localized regime.

The most interesting result is that we find an accurate universal scaling of the localization length $`\xi `$ with temperature. Our extracted values for $`\gamma `$ agree within experimental accuracy well with the numerical results $`\gamma =2.3`$ (Huckestein (1995) and references therein) as well as with the exponents determined from size dependent scaling experiments Koch et al. (1991a). This underlines the fact that our direct scaling analysis in the VRH regime suits much better as an access to the localization length and its scaling behavior than temperature dependent peak-width scaling where no universality is observed. In fact, for such a conventional analysis we are still lacking a sufficient knowledge of the temperature dependent transport mechanisms in the metallic regime of the QHE plateau transition.

Until now we used the term scaling in a rather reduced sense as synonym for a power-law behavior of the localization length $`\xi `$. Of course scaling includes much more, namely the existence of a single parameter scaling function $`\sigma _{xx}=f(x)`$ with a parameter $`x(T,\delta \nu )`$. Using our above results in the VRH regime an appropriate definition is $`x=|\delta \nu |^\gamma /TT_0/T`$ with the above deduced exponents $`\gamma `$. The postulation of single parameter scaling together with the finding of $`\sigma _01/T`$ then fixes the prefactor in Eq. 1 to $`\sigma _0(\delta \nu ,T)=\sigma ^{}|\delta \nu |^\gamma /T`$ with a constant $`\sigma ^{}`$. This yields a scaling of the conductivity in the form

$$\sigma _{xx}\left(\frac{|\delta \nu |^\gamma }{T}\right)=\sigma ^{}\frac{|\delta \nu |^\gamma }{T}\mathrm{exp}\left(\sqrt{T^{}\frac{|\delta \nu |^\gamma }{T}}\right)$$
(2)
where $`T^{}`$ is constant.

In Fig. 3 we have plotted all the conductivity data $`\sigma _{xx}(\nu ,T)`$ for $`\nu >\nu _c`$ as function of a single parameter $`|\delta \nu |^\gamma /T`$. The values for $`\gamma `$ are taken from the power-law fits of $`T_0`$ as shown in Fig. 2. Rescaling the axes in a proper way indeed shows that all experimental points fall onto straight lines as represented by the scaling function in Eq. 2. Therefore, the conductivity in the VRH regime for $`\delta \nu <0.3`$ and $`\nu >\nu _c`$ follows the stringent postulation of single parameter scaling and thus demonstrates the large validity range of the universal scaling phenomenon in the transition between QHE plateaus.

For $`\nu <\nu _c`$, i.e. approaching the spin gap induced plateau at $`\nu =1`$, the single parameter rescaling does not look as good as for approaching the Landau gap at $`\nu =2`$, although Fig. 2 showed the same quality of $`T_0`$-scaling for both sides of the transition. This discrepancy probably mirrors a more complicated behavior of the prefactor $`\sigma _0`$ in Eq. 1. Obviously it can no more be written as a simple function $`|\delta \nu |^\gamma /TT_0/T`$. Such a more complicated behavior in this regime is no surprise taking into account the redistribution of electrons between spin states and the emerging of spin textures as a function of $`\nu `$ at this filling range Girvin and MacDonald (1997).

In addition to the temperature dependence of $`\sigma _{xx}`$ in the VRH regime we have also investigated its dependence on the electric field strength $`F`$. In the nonlinear regime at high voltages the effect of $`F`$ can be interpreted as an effective temperature $`T_F=eF\xi /2k_B`$ Polyakov and Shklovskii (1993a). Using Eq. 1 the conductivity is then rewritten to $`\sigma _{xx}=\sigma _0^F(F)\mathrm{exp}\left(\sqrt{F_0/F}\right)`$ with $`F_0T_0/\xi 1/\xi ^2`$. Analogous to the single parameter scaling with temperature in the VRH regime the natural reduced parameter for scaling with electric field is defined by $`y(F,\nu )=|\delta \nu |^{2\gamma }/F`$. Applying such an analysis to our data without any additional fit parameter (i.e. using the exponents $`\gamma `$ as determined from $`T_0(\nu )`$, see Fig. 2) we find for $`\nu >\nu _c`$ for all data in the range $`\sigma _{xx}10^4e^2/h`$ a single parameter scaling; an example for sample S2 is shown in the inset of figure 3.

In conclusion we have investigated the temperature and electric field dependence of the conductivity in the QHE plateau transition for samples where no universal behavior is found in the conventional temperature dependent peak-width scaling experiments. We find a revival of universality in the VRH regime where the localization length $`\xi `$ scales as $`\xi |\delta \nu |^\gamma `$ with an experimentally deduced scaling exponent close to the theoretically expected value $`\gamma =2.3`$. Even further, we have shown that all the data on the Landau gap side of the transition can be rescaled on a single parameter function over more than five orders of magnitude in the conductivity.

We thank F. Evers, B. Huckestein, B. Kramer and D. G. Polyakov for stimulating discussions. The samples were grown by K. Ploog at the Max Planck-Institut für Festkörperforschung.