text:
# Chaotic front dynamics in semiconductor superlattices

## I Introduction

Electronic transport in semiconductor superlattices (SL) is known to show strongly nonlinear spatio-temporal dynamics. Either self–sustained current oscillations KAS95 ; HOF96 ; KAS97 ; PAT98 ; STE00a  or a sawtooth–like current–voltage characteristic with many branches associated with static field domains ESA74 ; KAW86 ; GRA91 ; SCH00  have been found. For a recent review see WAC01 . Under time-dependent external voltage conditions, superlattices exhibit a rich menagerie of complex behavior including ac driven chaos and switching scenarios between multistable states. This was studied recently experimentally ZHA96 ; KAS96 ; LUO98 ; SHI97a ; ROG01  and theoretically BUL95 ; AMA01 .

In this paper we present simulations of dynamic scenarios for superlattices under fixed time-independent external voltage in the regime where self-sustained dipole waves are spontaneously generated at the emitter. The dipole waves are associated with traveling field domains, and consist of electron accumulation and depletion fronts which in general travel at different velocities and may merge and annihilate. We find that depending on the applied voltage and the contact conductivity, this gives rise to various oscillation modes and self-synchronization effects as well as different routes to chaotic behavior.

A similar scenario of merging fronts was recently found in the context of a spatially continuous model describing bulk impurity impact ionization breakdown CAN01 . It is also reminiscent of patterns of temperature pulses in globally coupled heterogeneous catalytic systems, e.g. GRA93a .

## II The Model

Weakly coupled superlattices are successfully described by a one dimensional sequential tunneling model for electrons PRE94 ; BON94 ; KAS97 . In the framework of this model electrons are assumed to be localized at one particular well and only weakly coupled to the neighboring wells. The tunneling rate to the next well is lower than the typical relaxation rate between the different energy levels within one well. The electrons within one well are then in quasi–equilibrium and transport through the barrier is incoherent. The resulting tunneling current density $`J_{mm+1}(F_m,n_m,n_{m+1})`$ from well $`m`$ to well $`m+1`$ depends only on the electric field $`F_m`$ between both wells and the electron densities $`n_m`$ and $`n_{m+1}`$ in the wells (in units of $`cm^2`$). For details of the microscopic calculation of $`J_{mm+1}`$ we refer to the literature AMA01 ; WAC01 . A typical result for the current density vs electric field characteristic is depicted in Fig. 1 in the spatially homogeneous case, i.e. $`n_m=n_{m+1}=N_D`$, with donor density $`N_D`$.

In the following we will adopt the total number of electrons in one well as the dynamic variables of the system. The dynamic equations are then given by the continuity equation

$$e\frac{\text{d}n_m}{\text{d}t}=J_{m1m}J_{mm+1}\text{for }m=1,\mathrm{}N$$
(1)
where $`N`$ is the number of wells in the superlattice.

The electron densities and the electric fields are coupled by the following discrete version of Gauss’s law

$$ϵ_rϵ_0(F_mF_{m1})=e(n_mN_D)\text{for }m=1,\mathrm{}N,$$
(2)
where $`ϵ_r`$ and $`ϵ_0`$ are the relative and absolute permittivities, $`e<0`$ is the electron charge, and $`F_0`$ and $`F_N`$ are the fields at the emitter and collector barrier, respectively.

The applied voltage between emitter and collector gives rise to a global constraint

$$U=\underset{m=0}{\overset{N}{}}F_md,$$
(3)
where $`d`$ is the superlattice period.

The current densities at the contacts are chosen such that dipole waves are generated at the emitter. For this purpose it is sufficient to choose Ohmic boundary conditions:

$`J_{01}`$ $`=`$ $`\sigma F_0`$ (4)
$`J_{NN+1}`$ $`=`$ $`\sigma F_N{\displaystyle \frac{n_N}{N_D}}`$ (5)
where $`\sigma `$ is the Ohmic contact conductivity, and the factor $`n_N/N_D`$ is introduced in order to avoid negative electron densities at the collector. Here we make the physical assumption that the current from the last well to the collector is proportional to the electron density in the last well. It is in principle possible to calculate the boundary conditions using microscopic considerations AGU97 ; BON00 , but the qualitative behavior is not changed.

## III Numerical results

In our computer simulations we use an $`N=100`$ superlattice with Al0.3Ga0.7As barriers of width $`b=5\text{nm}`$ and GaAs quantum wells of width $`w=8\text{nm}`$, doping density $`N_D=1.0\times 10^{11}\text{cm}^2`$ and scattering induced broadening $`\mathrm{\Gamma }=8\text{meV}`$ at $`T=20\text{K}`$. The contact conductivity $`\sigma `$ is choosen such that the intersection point with the homogeneous current density vs field characteristic in Fig. 1 is at a current value at which no stationary field domain boundaries exists. By this configuration, accumulation and depletion fronts are generated at the emitter. For large values of $`\sigma `$, e.g. $`\sigma =1.3\mathrm{\Omega }^1\text{m}^1`$, we find that those fronts form a dipole, i.e. a traveling field domain, with leading electron depletion front and trailing accumulation front (cf Fig. 2 (a)). The dipole traverses the sample at almost constant velocity and constant current which shows up as a plateau in the time trace of the current. The current $`J_D`$ and the front velocity in this dipole propagation mode are given by the intersection point of the velocities of accumulation and depletion front (cf. inset of Fig. 1). As the leading depletion front reaches the collector, there is no speed constraint on the remaining accumulation front. As it accelerates, the current rises, and a new dipole is generated at the emitter. Now for a short time there are two accumulation and one depletion fronts in the sample. In order to fulfill eq.(3) for fixed $`U`$ the depletion front has to assume twice the velocity of the accumulation fronts. This constraint fixes the current to $`J_T`$ (cf. inset of Fig. 1) during this tripole regime. Note that in the current time trace the fast small-amplitude oscillations (due to well-to-well hopping of depletion and accumulation fronts in our discrete model) in the dipole and tripole regime are not resolved temporally.

At lower contact conductivity $`\sigma =0.6\mathrm{\Omega }^1\text{m}^1`$, we find that instead of a dipole only a depletion front is generated at the emitter, as the old depletion front reaches the collector (cf. Fig. 2 (b)). Now for a short time a dipole with leading accumulation front exists. The current is fixed again by the constraint of equal velocities of accumulation and depletion front, as for the dipole with leading depletion front. At the time the old accumulation front reaches the collector a new accumulation front is generated at the emitter. This process is accompanied by a dip in the current time trace.

For even lower $`\sigma =0.55\mathrm{\Omega }^1\text{m}^1`$ the fronts at the emitter are generated as dipoles with leading accumulation and trailing depletion fronts (cf. Fig. 2 (c)). The velocity of the fronts is again determined by the current and the number of fronts in the sample. Since the two types of fronts in general move at different velocities, merging and annihilation of accumulation and depletion fronts may occur, which may lead to complicated behaviour including chaos as shown in Fig. 2 (c).

In order to study the bifurcation scenario leading to chaos we now fix the boundary conductivity to $`\sigma =0.5\mathrm{\Omega }^1\text{m}^1`$. Experimentally, low contact conductivities may be realized by changing the doping and barrier structure of the emitter. The temperature dependence of the emitter current may also be exploited. In Fig. 3 (a) a density plot of the positions (well numbers) at which two fronts annihilate is shown as a function of the voltage. We see that for low voltage the annihilation takes place at one definite position in the superlattice with a variation of only a few wells. This distribution broadens for increasing voltage in characteristic bifurcation scenarios reminiscent of period doubling, leading to chaotic regimes. We note that in the chaotic region periodic windows exist. Since we are dealing with a discrete system, the position of the merging of two fronts is also discrete. For more accurate analyses it is convenient to use a continuous variable such as the time difference between two maxima in the electron density in a specific well. The corresponding bifurcation diagram for well no. 20 is shown in Fig. 3 (b). It exhibits a complex structure.

The transition from periodic to chaotic oscillations is enlightened by considering the space-time plots for the evolution of the electron densities at different voltages $`U`$ (Fig. 4). While for $`U=0.5\text{V}`$ the oscillations are regular, we see that at $`U=0.7\text{V}`$ long and short front patterns alternate, which is characteristic for a period doubling bifurcation. Every second pair of fronts travels farther into the superlattice before they annihilate. At $`U=0.79\text{V}`$ a further period doubling occurs, as can be most clearly seen at the alternating length of the longer front patterns. This yields a period-four cycle. Further increase of the voltage finally leads to well developed chaos, with front patterns of different lengths. For $`U1.8\text{V}`$ we find that fronts may even traverse the sample, while the dynamics remains chaotic.

## IV Summary

We have shown that complex chaotic spatio-temporal scenarios may arise in weakly coupled superlattices under time-independent external voltage conditions, whereas previous studies have focussed on ac-driven chaos.

This work was partially supported by DFG in the framework of Sfb 555. We greatfully acknowledge discussion with L. L. Bonilla, M. Rogozia and E. Schomburg.


format: text

text:
# Universal Pulse Shape Scaling Function and Exponents: A Critical test for Avalanche Models applied to Barkhausen Noise

## I  Introduction

Real materials, have dirt, or disorder, which often leads to slow (glassy) dynamics due to complex free energy landscapes with diverging energy barriers. On long length scales and practical time scales, thermal effects often become unimportant. Indeed, when such a system is driven by an external field it jumps from one local metastable free energy minimum to the next, and the state of the system depends on its history.

In magnetic materials jumps from one local minimum to the next involve a collective process whereby clusters of magnetic domains change the direction of their magnetization, in an avalanche. These avalanches can be triggered by a slowly but continuously increasing homogeneous external magnetic field $`H`$ (taken from $`\mathrm{}`$ to $`+\mathrm{}`$). These avalanches produce so-called Barkhausen noise, which can be observed experimentally as a voltage signal induced in a pickup coil wound around the magnet. Experiments show that these avalanches come in all sizes; their sizes are typically distributed according to a power law over several decades. Other systems also exhibit a broad range of avalanche sizes and durations following power laws: superconducting vortex line avalanches, resistance avalanches in a superconducting film, capillary condensation of helium in Nuclepore, acoustic emissions in athermal martensites, earthquakes, and many others.

We examine two variants of the non-equilibrium zero temperature random field Ising model (RFIM), which both predict power law distributions of avalanche sizes and universal non-equilibrium collective behavior. The RFIM is a model for a conventional magnet in which magnetic domains are modeled by “spins” on a lattice that can only point up or down. The first variant basically exhibits single front (domain wall) propagation dynamics in which spins at the edge of an existing front flip when it is (locally) energetically favorable to do so. Spins that are not adjacent to this front are very unlikely to flip on their own due to the presence of an infinite-range demagnetizing field, which is present in addition to nearest neighbor ferromagnetic interactions. We call this the front propagation model. In the second variant, called the nucleation model, spins flip anywhere in the system when it is (locally) energetically favorable to do so. In this case there are many interacting fronts, unlike in the front propagation model. Only nearest neighbor ferromagnet interactions are included in the nucleation model.

In this paper we perform a detailed quantitative comparison of the universal avalanche pulse shapes and exponents obtained from: front propagation dynamics, domain nucleation dynamics, and experiment. Our analysis constitutes a test of whether the non-equilibrium zero-temperature RFIM (either variant) is in the same universality class as experimental systems exhibiting Barkhausen noise.

## II  The Model

The RFIM consists of a (hypercubic) lattice of $`N`$ spins ($`s_i=\pm 1`$), which may point up ($`s_i=+1`$) or down ($`s_i=1`$). Spins are coupled to nearest neighbors (through a ferromagnetic exchange interaction $`J`$), and to an external field $`H(t)`$ which is increased adiabatically slowly. To model dirt in the material, we assign a random field, $`h_i`$, to each spin, chosen from a distribution $`P(h_i)=exp(h_i^2/2R^2)/\sqrt{2\pi }R`$, where $`R`$, the disorder, determines the width of the Gaussian probability distribution and therefore gives a measure of the amount of quenched disorder for the system.

The Hamiltonian for the system at a time $`t`$ is given by: $`H=\underset{<ij>}{}Js_is_j_i(H(t)+h_iJ_{inf}M)s_i`$, where $`J_{inf}`$ is the strength of the infinite range demagnetizing field ($`J_{inf}=0`$ for the nucleation model), $`M=\frac{1}{N}\underset{i}{}s_i`$ is the magnetization of the system, and $`<ij>`$ stands for nearest neighbor pairs of spins. Initially, $`H(\mathrm{})=\mathrm{}`$ and all the spins are pointing down. Each spin is always aligned with its local effective field $`h_i^{eff}=J_{<ij>}s_j+H(t)+h_iJ_{inf}M`$. The external field $`H(t)`$ is adiabatically slowly increased from $`\mathrm{}`$ until the local field, $`h_i^{eff}`$, of any spin $`s_i`$ changes sign, causing the spin to flip. It takes some microscopic time $`\mathrm{\Delta }t`$ for a spin to flip ($`\mathrm{\Delta }t1`$ for our simulation). The spin flip changes the local field of the nearest neighbors and may cause them to flip as well, etc. This avalanche process continues until no more spin flips are triggered. Each step of the avalanche, that is, each $`\mathrm{\Delta }t`$ in which a set of spins simultaneously flip, is called a shell. The number of spins that flip in a shell is directly proportional to the voltage $`V(t)`$ during the interval $`\mathrm{\Delta }t`$ that an experimentalist would measure in a pick-up coil wound around the sample. In our simulations we therefore denote the number of spins flipped in a shell at a time $`t`$ by $`V(t)`$. The first shell of an avalanche (one spin flip) is triggered by the external field $`H(t)`$, while each subsequent shell within the avalanche is triggered only by the previous shell, since $`H(t)`$ is kept constant while the avalanche is propagating. $`H(t)`$ is only increased when the current avalanche has stopped, and is increased only until the next avalanche is triggered (i.e. $`\frac{dH}{dt}0`$). The number of shells in an avalanche times $`\mathrm{\Delta }t`$ defines the pulse duration, $`T`$, or the time it took for the entire avalanche to flip. In this paper we will be interested in looking at $`V(t,T)`$ for $`0<tT`$, that is, the voltage as a function of time for an avalanche of a given duration $`T`$.

The front propagation model exhibits self-organized criticality (SOC). This means that as $`H`$ is increased the model always operates at the critical depinning point, and no parameters need to be tuned to exhibit critical scaling behavior (except $`\frac{dH}{dt}0`$). The nucleation model, on the other hand, is a plain old critical system with a continuous second order phase transition with disorder as the tuning parameter. The continuous non-equilibrium phase transition can be understood as follows. For zero disorder, the random fields of each spin will be the same, so that when one spin is flipped the entire lattice of spins will flip. This results in a rectangular hysteresis curve with a macroscopic jump in the magnetization when the external field overcomes the interaction with the neighbors. On the other hand, if the disorder is infinite, each spin will have a very different random field, so as the external field is raised each spin will essentially flip independently, triggering no other spins to flip. This will result in a smooth “hysteresis” curve with an approximately constant slope ($`M(H)H`$ over a wide range of $`H`$). In between these two phases of behavior there is a continuous phase transition at some critical value $`R=R_c`$. At the transition each branch of the hysteresis loop has a single point with infinite slope ($`\frac{dM}{dH}_{\pm H_c}\mathrm{}`$) at a critical field $`H=\pm H_c`$. $`R_c`$ and $`H_c`$ are nonuniversal. In three dimensions they are: $`R_c=2.16`$ and $`H_c=1.43`$ (in units of $`J`$). The transition is characterized by a number of universal critical exponents , and scaling functions. In this paper we specifically focus on only two exponents: $`1/\sigma \nu z`$ and $`\tau `$, and one scaling function. The significance of these exponents is given later. The scaling function we examine can be obtained through experiment as well as simulation. Details of the simulation algorithm are given elsewhere.

## III  The Experiments

In addition to examining results obtained from simulation we study results from three different experiments. We performed one of these experiments, and the results for the other two experiments were obtained from already published results.

Our own experiment was performed on an (unstressed) amorphous alloy, $`Fe_{21}Co_{64}B_{15}`$. The data we present from Durin et al is from an experiment on an amorphous alloy with a different composition, $`Fe_{64}Co_{21}B_{15}`$, under a tensile stress. Studies of the effect of tensile stress on samples of this type indicate that for low tensile stress, the domain structure is a complicated pattern of maze domains, dominated by quenched-in stresses. On the other hand, when such materials undergo tensile stress, the uniaxial anisotropy gives way to a simpler domain structure with a few parallel domains in the direction of the stress. Related to this change in domain structure is a change in the dominant interaction in the material. In amorphous alloys under stress, surface tension effects are thought to be more important than dipolar interactions, while dipolar interactions dominate for polycrystals and materials with small grains .

Durin et al’s sample was under stress so as to enhance stress-induced anisotropy so much that the long range dipolar interactions can be neglected, placing their experiment into the universality class of the front propagation model. Spasojević et al’s sample was a quasi-two-dimensional metal glass, more precisely a commercial VITROVAC 6025 X. Based on the scaling exponents obtained from Spasojević et al’s experiment (given later), their experiment does not seem to fall into any universality class discussed in this paper. Our experiment seems to be in a crossover regime between two universality classes; details are given later in this paper. Further details about Durin et al’s, and Spasojević et al’s experiment can be obtained elsewhere .

## IV  Extracting Pulse Shapes

Nucleation Model: We simulate four realizations of a $`1200^3`$ system near $`R=R_c`$ ($`R=2.2`$) and record the time series $`V(t,T)`$ of avalanches from an $`H`$ window near $`H_c`$ ($`1.42<H<1.43`$). We average avalanches of a fixed pulse duration $`T`$ (within the interval $`[T,1.05T]`$), for various values of $`T`$. In each case we average over $`1000`$ to $`2000`$ avalanches to ensure strong fluctuations have been averaged out. We check finite size effects by performing simulations at $`800^3`$ and $`1000^3`$ and verifying that identical avalanche shapes (within small fluctuations) are obtained for all system sizes.

Front Propagation Model: We perform $`100`$ realizations of a $`400^3`$ system and record avalanches from an $`H`$ window within the slated part of the hysteresis loop ($`1.25<H<1.88`$). Even though the front propagation model exhibits SOC due to the infinite range demagnetizing interaction, in order to avoid effects due to initial nucleation of the front (beginning of the hysteresis loop), or when the front encounters the boundaries of the simulation (the end of the hysteresis loop) we must choose avalanches near the middle of the hysteresis loop. We obtain avalanche shapes $`V(t,T)`$ in a manner identical to how they were obtained for the nucleation model. We check finite size effects by performing simulations of $`200^3`$ and $`300^3`$ size systems, and find a consistent avalanche shape for all three system sizes.

Experiment: Measurements were performed on a 21 cm x 1 cm x 30 $`\mu `$m ribbon of $`Fe_{21}Co_{64}B_{15}`$ alloy, a soft amorphous ferromagnet obtained from Gianfranco Durin. The domain walls run parallel to the long axis of the material, with about 50 domains across the width. A solenoid, driven with a triangle wave, applies a magnetic field along the long axis of the sample. Since domain wall motion dominates over other means of magnetization in the linear region of the loop, data were collected in only a selected range of applied fields near the center of the loop. The Barkhausen noise was measured by a small pick-up coil wound around the center of the sample. This voltage signal was amplified, anti-alias filtered and digitized, with care taken to avoid pick-up from ambient fields. Barkhausen noise was collected for both increasing and decreasing fields for $`80`$ cycles of the applied field through a saturation hysteresis loop. The driving frequency was $`0.01`$ Hz; this corresponds to $`c=0.09`$, where c is a dimensionless parameter proportional to the applied field rate and is defined in the Alessandro Beatrice Bertotti Montorsi model (ABBM model) for the Barkhausen effect. In this way, our measurements should be well inside the $`c<1`$ regime identified in the ABBM model, in which we can expect to find separable avalanches rather than continuous domain wall motion.

Due to background noise in real experimental data there is no definitive way to determine when an avalanche begins and when one ends; we set a sensitivity threshold which is high enough to cut out background noise and low enough to capture enough of the avalanche so as not to affect the shape. We check the validity of our threshold by perturbing the threshold by a small amount and noticing that there is no change in the avalanche shape.

## V  Critical Exponents and Data Collapses

In the front propagation model, as well as in the nucleation model near ($`R_c,H_c`$), the voltage $`V(T,t)`$ scales as:

$$V(T,t)=T^{1/\sigma \nu z1}f_{shape}(t/T)$$
(1)
By collapsing average avalanche shapes of various durations $`T`$ we determine the universal scaling function, $`f_{shape}(t/T)`$, and the critical exponent $`1/\sigma \nu z`$. The exponent $`1/\sigma \nu z`$ relates the avalanche size, $`S`$, to the avalanche pulse duration, $`T`$, at criticality by $`ST^{1/\sigma \nu z}`$. We find that the critical exponent $`1/\sigma \nu z`$ for the front propagation and nucleation model, obtained from simulation, is in close agreement with previous theoretical predictions and with experimental values. The collapses are shown in Figs. 1 and 2.

We also determine the avalanche size distribution that scales as $`D(S)S^\tau `$ at criticality. The avalanche size distributions for our simulation and our experiment are given in Fig. 3; the values obtained for $`\tau `$, for the front propagation model and nucleation model, are in close agreement with previously quoted values. Experimentally, the scaling exponent $`\tau `$ weakly depends on $`c`$ in some materials . From this dependence we find that there may be a difference of $`0.02`$ between the value of $`\tau `$ we find in our experiment and the value of $`\tau `$ at $`c=0`$ (zero frequency). However, this difference is within the error bars we give for $`\tau `$. Table I summarizes results for $`\tau `$ and $`1/\sigma \nu z`$ for experiment, simulation, and mean field theory. We see that the critical exponents for the front propagation model agree (within error bars) with exponents from Durin et al’s experiment, as expected.

While Durin et al’s experiment is believed to fall into the front propagation universality class, as discussed above, our experiment is neither in the front propagation nor in the mean field universality class. A priori we would assume that our experiment would be in the mean field universality class since our sample was unstressed, and previous experiments have indicated that unstressed samples will exhibit mean field behavior. The critical exponents found from our experiment indicate that it maybe in a crossover regime between the mean field and the front propagation universality classes. The exponent $`\tau `$ has a value of $`1.46\pm 0.05`$ (see also Table I), which is between the value of $`1.28`$ for the front propagation model (sample with stress) and the mean field value of $`1.5`$ (sample without stress). In addition, the exponent $`\alpha =1.74\pm 0.06`$, determined from the avalanche pulse duration distribution given by $`D(T)T^\alpha `$ at criticality, (see Inset of Fig. 2), is between the value of $`1.5`$ for the front propagation model (with stress), and the mean field value of $`2`$ (without stress). Residual stress on our sample may have resulted in these anomalous exponents.

Although our experiment may be in a crossover regime, we were able to obtain a good collapse of the avalanche pulse shapes (see Fig. 2). In order to reaffirm the validity of the exponent obtained from the collapse, we independently checked the value of $`1/\sigma \nu z`$ from the power spectra ($`P(w)w^{1/\sigma \nu z}`$ at criticality), and found $`1/\sigma \nu z=1.73\pm 0.08`$, which is consistent (within errors bars) with $`1/\sigma \nu z=1.70\pm 0.05`$ obtained from the experimental avalanche pulse collapse.

## VI  Fitting to Orthonormal Polynomials

The mean field shape of avalanches in our model $`V^{MFT}(t,T)`$ is an inverted parabola . In order to study corrections to the mean field shape, in simulation and experiment, we derive a set of orthonormal polynomials $`f_i(t)`$ with $`f_i(L)=f_i(L)=0`$, where $`L=T/2`$ is half the duration of the avalanche. The negative of the first polynomial, $`f_0(t)`$, is proportional to the mean field result. The first four polynomials of the set are:

$`f_0(t)=\sqrt{{\displaystyle \frac{15}{16L^5}}}(t^2L^2)`$ (2)
$`f_1(t)=\sqrt{{\displaystyle \frac{105}{16L^7}}}(t^3L^2t)`$ (3)
$`f_2(t)=\sqrt{{\displaystyle \frac{45}{64L^9}}}(7t^48L^2t^2+L^4)`$ (4)
$`f_3(t)=\sqrt{{\displaystyle \frac{1155}{64L^{11}}}}(3t^54L^2t^3+L^4t)`$ (5)
$`f_4(t)=\sqrt{{\displaystyle \frac{1365}{2048L^{13}}}}(33t^651L^2t^4+19L^4t^2L^6)`$ (6)
We fit a linear combination of the above polynomials to the average avalanche shape obtained from simulation and experiment. $`L`$ is also left as a free parameter in the fits to the average shapes. The pulse duration is then precisely defined as $`T=2L`$. The results of the fits are shown in Fig. 1 and Fig. 2. In Fig. 4 we give the coefficients for the fits, found in simulation and experiment. We also include the coefficients for the fit to an average avalanche shape determined by Durin et al’s experiment , whose avalanche shape is given in Fig. 5. The coefficients are determined from the total fit function:

$$F(t)=a_0f_0(t)+a_1f_1(t)+a_2f_2(t)+a_3f_3(t)+a_4f_4(t)$$
(7)
$`a_0`$, $`a_2`$, and $`a_4`$ are the symmetric  coefficients of the fit, while $`a_1`$, $`a_3`$ are the antisymmetric  coefficients of the fit (i.e. multiplying polynomials that are not symmetric under time reversal $`tt`$). We are particularly interested in the asymmetry of the avalanche shapes.

From inspection of the avalanche shapes we see that the experimental avalanche shapes are strongly asymmetric under time reversal, while the avalanche shapes determined from the simulation of the two models are both very close to symmetric. Quantitatively, the coefficients for the avalanche shapes of the nucleation model are very similar to those of the front propagation model and both are different from the experimental avalanche shapes’ coefficients. While avalanche shapes in both models are slightly asymmetric to the left (i.e. $`V(t,T)`$ increases more slowly than it decreases), the experimental avalanche shapes are strongly asymmetric to the right direction (i.e. $`V(t,T)`$ increases quickly and decreases slowly). The origin of this difference between theory and experiment is not yet understood. It comes as a surprise, since the critical exponents obtained in Durin et al’s experiment are in close agreement with exponents obtained from the front propagation model (see Table I) and the avalanche shapes are expected to be just as universal as the critical exponents. On the other hand, we see in Fig. 5 that different experiments give different pulse shape scaling functions. This difference may be a result of the fact that the experiments are not in the same universality class. Nevertheless, all the experimental pulse shapes do have the same sign of asymmetry.

We did one check on whether the pulse shape function in our model is indeed universal, an assumption often taken for granted: We changed the lattice of our simulation from a simple cubic to a BCC (body centered cubic) lattice and found an identical pulse shape scaling function (within small statistical error), as expected when universality holds.

## VII  discussion and concluding remarks

In performing these analyses we have stumbled upon an interesting observation: after appropriate rescaling of the y-axis, the universal pulse shape function for the front propagation and nucleation model appear to look the same (see Inset for Fig. 1) even though we know  that the models do not belong to the same universality class. However, upon more precise quantitative analysis we see that they are in fact different, their $`a_1`$ coefficients are different by several $`\sigma `$. Furthermore, we tried to collapse two pulse shape functions with $`T52`$ from the two different models and found that they could not be collapsed precisely despite, naturally, scaling by different critical exponents appropriate for the two models. The front propagation pulse scaling function is more asymmetric, as supported by the somewhat larger $`a_1`$ fit coefficient compared to the nucleation model result.

By examining the pulse shape scaling functions as a sharper test than merely critical exponents for the universality class of the non-equilibrium zero temperature RFIM, we raised many questions. What accounts for the difference between theory and experiment, and between different experiments? Is the theory incomplete or inaccurate at this level of description? Experimentally, we do not yet know what material features are required to produce universality of the pulse shape function. Differences between experimental results make further experimental tests desirable (see Fig. 5).

Although the three experiments we examined do not fall into the same universality class and do not have very similar pulse shapes, they share one universal aspect not shared by the results from simulation: all the experimental averaged pulse shapes are asymmetric in the leftward direction. This suggests that there may be a phenomena that exists in real materials that has not been accounted for in theory.

ACKNOWLEDGMENTS

We would like to thank M. B. Weissman, G. Durin, A. Travesset, and R. A. White for very useful discussions. We also thank M. Kuntz and J. Carpenter for providing the front propagation model simulation code. K.D. and A.P.M. acknowledge support from NSF via Grant Nos. DMR 99-76550, the Materials Computation Center, through NSF Grant No. 99-72783, and IBM which provided the computers that made the simulation work possible. J.P.S. acknowledges support from NSF via Grant No. DMR 98-73214. A.C.M. acknowledges support from NSF via Grant No. DMR 99-81869. A.P.M. would also like to acknowledge the support provided by UIUC through a University Fellowship, and K.D. gratefully acknowledges support through an A.P. Sloan fellowship.


text:
# Order Parameter to Characterize Valence-Bond-Solid States in Quantum Spin Chains

## Abstract

We propose an order parameter to characterize valence-bond-solid (VBS) states in quantum spin chains, given by the ground-state expectation value of a unitary operator appearing in the Lieb-Schultz-Mattis argument. We show that the order parameter changes the sign according to the number of valence bonds (broken valence bonds) at the boundary for periodic (open) systems. This allows us to determine the phase transition point in between different VBS states. We demonstrate this theory in the successive dimerization transitions of the bond-alternating Heisenberg chains, using the quantum Monte Carlo method.

About two decades ago, Haldane made a striking conjecture that the integer-spin $`S`$ antiferromagnetic Heisenberg chain has a unique disordered ground state with an energy gap, while for a half-odd-integer spin there is no energy gap and the system belongs to the same universality class as the $`S=1/2`$ case Haldane . For $`S=1`$, this conjecture was confirmed by many numerical Nightingale-B  and experimental studies Buyers . Affleck, Kennedy, Lieb, and Tasaki studied an $`S=1`$ isotropic spin chain with special biquadratic interactions which has the exact ground state with a finite gap Affleck-K-L-T . They also proposed the valence-bond-solid (VBS) state for the Haldane gap systems, and concluded that the ground state of the Heisenberg chain is described approximately by the VBS state. Since the spin configurations of the VBS state show the hidden antiferromagnetic order, den Nijs and Rommelse proposed the string order parameter to characterize the $`S=1`$ Haldane phase Nijs-R ,

$$𝒪_{\mathrm{string}}^\alpha =\underset{|kl|\mathrm{}}{lim}\mathrm{\Psi }_0|S_k^\alpha \mathrm{exp}\left[\mathrm{i}\pi \underset{j=k+1}{\overset{l1}{}}S_j^\alpha \right]S_l^\alpha |\mathrm{\Psi }_0,$$
(1)
where $`\alpha =x,y,z`$ and $`|\mathrm{\Psi }_0`$ means the ground state. Thus this order parameter enables us to detect the VBS state indirectly. The string order parameter was generalized to $`S>1`$ cases by Oshikawa Oshikawa .

On the other hand, Haldane’s prediction was also discussed by the Lieb-Schultz-Mattis (LSM) type argument Lieb-S-M . Using this technique, Affleck and Lieb examined Haldane’s conjecture Affleck-L . Later, Oshikawa, Yamanaka, and Affleck generalized this argument, and obtained a necessary condition for a gapped state Oshikawa-Y-A . However, the relation between the LSM argument and the VBS picture including the string order parameter has not been fully understood. In this Letter, we discuss this relation, and show that an overlap integral appearing in the LSM argument \[Eq. (3) below\] plays a role of an order parameter which detects VBS ground states directly. We also show that it can also be applied to determine phase boundaries in between different VBS states. We demonstrate this idea by the quantum Monte Carlo (QMC) simulation for successive dimerization transitions in the bond-alternating Heisenberg spin chains.

First, let us review the LSM argument briefly based on Ref. Oshikawa-Y-A  and introduce the order parameter. We consider a periodic spin chain of length $`L`$ with short range interactions. We assume the translational ($`𝒯:S_j^\alpha S_{j+1}^\alpha `$) and the parity ($`𝒫:S_j^\alpha S_{L+1j}^\alpha `$) invariance. Now we introduce the following “twist operator”

$$U\mathrm{exp}\left[\mathrm{i}\frac{2\pi }{L}\underset{j=1}{\overset{L}{}}jS_j^z\right].$$
(2)
Since this operator rotates all the spins about the $`z`$-axis with relative rotation angle between the neighboring spins $`2\pi /L`$, it creates spin-wave-like excitations. Applying the twist operator $`q`$ times to the unique normalized ground state $`|\mathrm{\Psi }_0`$ generates a set of low-lying excited states $`\{|\mathrm{\Psi }_kU^k|\mathrm{\Psi }_0\}`$ ($`k=1,\mathrm{},q`$). The excitation energy of the state $`|\mathrm{\Psi }_1`$ is evaluated as $`\mathrm{\Delta }E𝒪(L^1)`$. Although $`|\mathrm{\Psi }_1`$ is not necessarily an eigenstate of the Hamiltonian, there exists at least one eigenstate with energy of $`𝒪(L^1)`$, if the state $`|\mathrm{\Psi }_1`$ is orthogonal to the ground state $`|\mathrm{\Psi }_0`$.

In order to consider the orthogonality of these states, we introduce the following overlap integral which plays a central role in this Letter:

$$z_L^{(q)}\mathrm{\Psi }_0|\mathrm{\Psi }_q=\mathrm{\Psi }_0|U^q|\mathrm{\Psi }_0.$$
(3)
The invariance under transformations $`𝒯`$ and $`𝒫`$ yields

$`z_L^{(q)}`$ $`=`$ $`\mathrm{\Psi }_0|𝒯U^q𝒯^1|\mathrm{\Psi }_0=\mathrm{e}^{\mathrm{i2}q(S_1^zm)\pi }z_L^{(q)},`$ (4)
$`=`$ $`\mathrm{\Psi }_0|𝒫U^q𝒫|\mathrm{\Psi }_0=\mathrm{e}^{\mathrm{i2}qm\pi }[z_L^{(q)}]^{},`$ (5)
where $`m`$ is the magnetization per site. Eq. (4) shows that $`z_{\mathrm{}}^{(1)}0`$ is possible only when $`Sm=`$ integer. In this case, the system has a gap without breaking the translational symmetry. On the other hand, for a rational value $`Sm=p/q`$ with $`p`$ being an integer, $`|\mathrm{\Psi }_0`$, $`|\mathrm{\Psi }_1`$,$`\mathrm{}`$, $`|\mathrm{\Psi }_{q1}`$ are mutually orthogonal. This means that the system in the $`L\mathrm{}`$ limit is gapless ($`z_{\mathrm{}}^{(q)}=0`$), otherwise gapped with $`q`$-fold degenerate ground state ($`z_{\mathrm{}}^{(q)}0`$) due to the spontaneous breaking of the translational symmetry. Thus $`q(Sm)=`$ integer is a necessary condition for a gapped state Oshikawa-Y-A . In this case, $`z_L^{(q)}`$ is real or pure imaginary due to Eq. (5).

In the previous works based on the LSM argument Lieb-S-M ; Affleck-L ; Oshikawa-Y-A , however, explicit values of $`z_L`$ have not been calculated. Now let us evaluate $`z_L`$ for various ground states. In the Néel state with two-fold degeneracy, one can obtain $`z_L^{(2)}=(1)^{2S}`$ immediately. In order to calculate $`z_L`$ in VBS states, we introduce the Schwinger boson representation for the spin operators Arovas-A-H :

$$\begin{array}{cc}S_j^+=a_j^{}b_j,& S_j^z=\frac{1}{2}(a_j^{}a_jb_j^{}b_j),\\ S_j^{}=b_j^{}a_j,& S_j=\frac{1}{2}(a_j^{}a_j+b_j^{}b_j),\end{array}$$
(6)
where these bosons satisfy the commutation relation $`[a_i,a_j^{}]=[b_i,b_j^{}]=\delta _{ij}`$ with all the other commutations vanishing. The operator $`a_j^{}`$ ($`b_j^{}`$) increases the number of up (down) $`S=1/2`$ variables under symmetrization. Then, a generalized VBS state discussed in Ref. Guo-K-M  in a periodic system (see Fig. 1(a)) is written as

$$|\mathrm{\Psi }_{\mathrm{VBS}}^{(m,n)}\frac{1}{\sqrt{𝒩}}\underset{k=1}{\overset{L/2}{}}(B_{2k1,2k}^{})^m(B_{2k,2k+1}^{})^n|\mathrm{vac},$$
(7)
where $`B_{i,j}^{}a_i^{}b_j^{}b_i^{}a_j^{}`$, $`𝒩`$ is a normalization factor, and $`|\mathrm{vac}`$ is the vacuum with respect to bosons. The integers $`m`$ and $`n`$ satisfy $`m+n=2S`$. Using relations $`Ua_j^{}U^1=a_j^{}\mathrm{e}^{+\mathrm{i}\pi j/L}`$ and $`Ub_j^{}U^1=b_j^{}\mathrm{e}^{\mathrm{i}\pi j/L}`$, a twisted valence bond $`UB_{j,j+1}^{}U^1`$ for $`1jL1`$, and that located at the boundary are calculated as follows,

$`UB_{j,j+1}^{}U^1`$ $`=`$ $`\mathrm{e}^{\mathrm{i}\pi /L}a_j^{}b_{j+1}^{}\mathrm{e}^{\mathrm{i}\pi /L}b_j^{}a_{j+1}^{},`$ (8)
$`UB_{L,1}^{}U^1`$ $`=`$ $`(\mathrm{e}^{\mathrm{i}\pi /L}a_L^{}b_1^{}\mathrm{e}^{\mathrm{i}\pi /L}b_L^{}a_1^{}).`$ (9)
In the latter case, a negative sign appears for each valence bond. This reflects a property of an $`S=1/2`$ spin which needs $`4\pi `$ rotation to return to the original state. Thus the asymptotic form of $`z_L`$ is given by

$$z_L^{(1)}=\mathrm{\Psi }_{\mathrm{VBS}}^{(m,n)}|U|\mathrm{\Psi }_{\mathrm{VBS}}^{(m,n)}=(1)^n[1𝒪(1/L)].$$
(10)
It turns out that $`z_L`$ changes its sign according to the number of valence bonds at the boundary. The value of $`z_L`$ for $`m=n=1`$ was calculated by Totsuka and Suzuki Totsuka-S . A similar relation was also found by Bonesteel who discussed two-dimensional $`S=1/2`$ dimer systems Bonesteel . In Eq. (10), the factor $`(1)^n`$ originates form the relative twist angle $`2\pi `$ at the boundary, so that it does not depend on the detailed structure of the twist operator (2). Therefore, the definition of $`z_L`$ is still meaningful even in cases with broken translational symmetry ($`mn`$) which does not accord with the LSM argument. Especially, in the present definition, $`z_L`$ has a useful symmetry $`z_L(1)^{2S}z_L`$ under interchange of $`m`$ and $`n`$.

In order to calculate $`z_L`$ in more detail, we introduce the matrix product (MP) formalism Fannes-N-W ; Klumper-S-Z . This is useful to calculate the ground-state expectation values and correlation functions. In this formalism, Eq. (7) is expressed as

$$|\mathrm{\Psi }_{\mathrm{VBS}}^{(m,n)}=\frac{1}{\sqrt{𝒩}}\mathrm{Tr}g_1^Ag_2^B\mathrm{}g_{L1}^Ag_L^B.$$
(11)
These matrices for general integer $`S`$ case with $`m=n`$ are given in Ref. Totsuka-S . Based on this result, the matrices for the twisted VBS state $`U^q|\mathrm{\Psi }_{\mathrm{VBS}}^{(m,n)}`$ are given by

$`g_{q,j}^A(s,r)=(1)^{ns+1}\mathrm{e}^{\mathrm{i}(n2s+2)q\pi /L}`$ (12)
$`\times \sqrt{{}_{m}{}^{}\mathrm{C}_{r1}^{}{}_{n}{}^{}\mathrm{C}_{s1}^{}}(a_j^{})^{mr+s}(b_j^{})^{n+rs}|\mathrm{vac}_j,`$
$`g_{q,j}^B(r,s)=(1)^{mr+1}\mathrm{e}^{\mathrm{i}(m2r+2)q\pi /L}`$ (13)
$`\times \sqrt{{}_{m}{}^{}\mathrm{C}_{r1}^{}{}_{n}{}^{}\mathrm{C}_{s1}^{}}(a_j^{})^{mr+s}(b_j^{})^{n+rs}|\mathrm{vac}_j,`$
where $`1rm+1`$, $`1sn+1`$, and $`|\mathrm{vac}_j`$ is the boson vacuum at the $`j`$-th site. By introducing a transfer matrix defined as

$$G_q^A(r_1,r_2;s_1,s_2)g_{0,j}^A(r_1,s_1)g_{q,j}^A(r_2,s_2),$$
(14)
we obtain $`z_L`$ as

$$z_L^{(1)}=(1)^n\frac{\mathrm{Tr}[G_1^AG_1^B]^{L/2}}{\mathrm{Tr}[G_0^AG_0^B]^{L/2}}.$$
(15)
Although the analytic form of $`z_L`$ is complicated in general, one can evaluate numerical values of Eq. (15), which will be shown in Fig. 2 below. In the fully dimerized case $`(m,n)=(2S,0)`$, $`z_L`$ is given by a simple form $`z_L^{(1)}=[_{k=0}^{2S}\mathrm{e}^{\mathrm{i}(2S2k)\pi /L}/(2S+1)]^{L/2}`$.

One can show that $`z_L`$ behaves similarly even in open systems. For even $`L`$ (see Fig. 1(b)), the VBS state of Eq. (7) has $`(n+1)^2`$-fold degeneracy due to the isolated bosons at both ends Oshikawa . Then, $`z_L`$ changes its sign according to the number of isolated bosons at the $`L`$-th site,

$$z_L^{(1)}=(1)^n\mathrm{e}^{\mathrm{i}n\pi /L}[1𝒪(1/L)].$$
(16)
Although $`z_L`$ is complex in this case, the imaginary part vanishes in the large-$`L`$ limit, and the real part of $`z_L`$ behaves in the same way as that of the periodic cases. For open systems with odd $`L`$, the number of isolated bosons at the $`L`$-th site is $`m`$, and the ground state has $`(m+1)(n+1)`$-fold degeneracy. Therefore, $`z_L`$ is simply given by Eq. (16) multiplied by $`(1)^{2S}`$.

Next we consider $`z_L`$ in connection with a low energy effective theory. According to Schulz’s bosonization analysis Schulz , the Lagrangian density of the uniform Heisenberg spin chain is given by the sine-Gordon model,

$$=\frac{1}{2\pi K}\left[\varphi (x,\tau )\right]^2\frac{y_\varphi }{2\pi \alpha ^2}\mathrm{cos}[q\sqrt{2}\varphi (x,\tau )],$$
(17)
where $`\tau `$ is the imaginary time, $`\alpha `$ is a short range cut off, and $`K`$ and $`y_\varphi `$ are the parameters determined phenomenologically. The phase field is related to that used in Ref. Schulz  by $`\varphi =2\sqrt{S}\psi _1`$, and $`q=1`$ ($`q=2`$) for $`S`$ integer ($`S`$ half-odd integer), where $`q`$ corresponds to degeneracy of gapped ground states. In the gapped (gapless) region one has $`y_\varphi (l)\pm \mathrm{}`$ ($`y_\varphi (l)0`$) for $`l\mathrm{}`$ under renormalization $`\alpha \mathrm{e}^l\alpha `$. On the unstable Gaussian fixed line \[$`y_\varphi (0)=0`$ with $`Kq^2<4`$\], a second-order “Gaussian transition” takes place between the two gapped states. In this formalism, and the spin wave excitation created by $`U`$ corresponds to the vertex operator $`\mathrm{exp}(\mathrm{i}\sqrt{2}\varphi )`$, so that $`z_L`$ for $`q(Sm)=`$ integer is related to the ground-state expectation value of the nonlinear term as $`z_L^{(q)}\mathrm{cos}(q\sqrt{2}\varphi )`$ and the three fixed points $`y_\varphi =\pm \mathrm{},0`$ correspond to $`z_{\mathrm{}}=1,0`$, respectively. Thus the Gaussian critical point can be identified by observing $`z_L^{(q)}=0`$ Nakamura-V .

In order to demonstrate the above argument, we consider the bond-alternating Heisenberg chain (BAHC),

$$=J\underset{j=1}{\overset{L}{}}[1\delta (1)^j]𝑺_j𝑺_{j+1}.$$
(18)
For this model, the VBS picture is considered to be realized approximately: The configuration of the valence bonds $`(m,n)`$ changes from $`(0,2S)`$ to $`(2S,0)`$ successively as $`\delta `$ is increased form $`1`$ to $`1`$, meaning the existence of $`2S`$ quantum phase transitions Affleck-H ; Oshikawa . Since the translational symmetry is explicitly broken, the effective model of the BAHC is given by Eq. (17) with $`q=1`$ for any $`S`$ Schulz , and these transitions are described as a Gaussian type.

In Fig. 2, we show $`z_L`$ (with $`q=1`$) of the $`L=64`$ periodic BAHC with $`S=\frac{1}{2}`$, 1, $`\frac{3}{2}`$, and 2 as a function of $`\delta `$. For the calculation of $`z_L`$, we employ the QMC method with the continuous-time loop algorithm Todo-K . We used the multi-cluster variant of the loop algorithm. The QMC steps are $`10^3`$ for thermalization, and $`5\times 10^5`$$`10^6`$ for the measurement, and the inverse temperature $`\beta `$ is taken large enough so that the value of $`z_L`$ can be identified as that of the ground state. The largest inverse temperature used in this calculation is $`\beta J=128`$.

As seen clearly in Fig. 2, the $`\delta `$ dependence of $`z_L`$ agrees qualitatively with the present interpretations based on the VBS picture, which predicts $`z_{\mathrm{}}=(1)^n`$ with a symmetry $`z_L(1)^{2S}z_L`$ for $`\delta \delta `$ ($`mn`$). We have also calculated $`z_L`$ for open chains with even and odd $`L`$’s, and confirmed that the results agree with our predictions including Eq. (16). In Fig. 2, we also present $`z_L`$ for the VBS states calculated by the MP method \[Eq. (15)\] as horizontal lines. The difference between $`z_L`$ for the BAHC and that for the VBS states becomes larger as $`S`$ increases. This means that the VBS picture becomes poor for the BAHC with large $`S`$.

Next, we determine the critical point $`\delta _\mathrm{c}`$ by observing $`z_L=0`$ with $`L`$ up to $`320`$ (Fig. 3). In this calculation, the inverse temperature is taken as $`\beta J=L/2S`$, being assumed the Lorentz invariance. Extrapolation to the $`L\mathrm{}`$ limit has been done by the least-squares fitting by assuming the asymptotic form as $`\delta _\mathrm{c}(L)=\delta _\mathrm{c}(\mathrm{})+A/L^2+B/L^4+C/L^6`$. This polynomial with even powers of $`1/L`$ is justified by the parity symmetry which ensures that $`z_L`$ is real \[Eq. (5)\]. For $`S=1`$, we have obtained the critical value for the transition between the $`(1,1)`$ and $`(2,0)`$ phases as $`\delta _\mathrm{c}=0.25997(3)`$, where $`()`$ denotes $`2\sigma `$. This result is consistent with the previous estimates: $`\delta _\mathrm{c}=0.2595(5)`$ by the QMC calculation for the susceptibility Kohno-T-H  and $`0.2598`$ by the level-crossing method Kitazawa-N . Similarly, we identify the critical point of the $`S=3/2`$ case \[$`(2,1)`$-$`(3,0)`$\] and those of the $`S=2`$ case \[$`(2,2)`$-$`(3,1)`$, $`(3,1)`$-$`(4,0)`$\]. We obtain $`\delta _\mathrm{c}=0.43131(7)`$, 0.1866(7), and 0.5500(1), respectively (see the inset of Fig. 3). They are also consistent with $`\delta _\mathrm{c}=0.4315`$, 0.1830, and 0.5505 obtained by the level-crossing method Kitazawa-N2 , but much more accurate.

Here, we comment on the method used in Refs. Kitazawa-N  and Kitazawa-N2  proposed by Kitazawa who pointed out that the Gaussian transition with $`q=1`$ can be identified by a level crossing of excitation spectra under twisted boundary conditions Kitazawa . Since this method is also explained by the sine-Gordon theory Kitazawa  and the VBS picture Kitazawa-N ; Kitazawa-N2 , the results are considered to be equivalent to ours. Although finite-size corrections in the level-crossing point in VBS states tend to be smaller than those of the $`z_L=0`$ point, the application of the level-crossing method to larger systems is difficult, because it relies on the exact diagonalization to obtain excitation spectra which needs the whole Hilbert space. On the other hand, our approach is based only on the ground state quantity, $`z_L`$, so that various numerical methods such as the present QMC method can be employed. This makes it possible to deal with systems with enormous Hilbert space, such as large-$`S`$ and ladder systems. Note that the density matrix renormalization group method is also suitable for our approach, since $`z_L`$ remains as a meaningful order parameter even in open systems. Furthermore, our method also has an advantage to comprehend the physical picture as already discussed.

Finally, we discuss the relation between $`z_L`$ and the string order parameter. In Ref. Oshikawa , Eq. (1) for $`|\mathrm{\Psi }_{\mathrm{VBS}}^{(m,n)}`$ was calculated as $`𝒪_{\mathrm{string}}^z=[\frac{m+n+2}{2(m+2)}]^2\delta _{n,\mathrm{odd}}`$, where finite and vanishing $`𝒪_{\mathrm{string}}^z`$ correspond to negative and positive $`z_L`$, respectively. However, this calculation is limited only for $`S`$ integer, and the twist angle $`\pi `$ in Eq. (1) is needed to be generalized for other cases Oshikawa . In addition, it is difficult to determine accurate phase transition points by $`𝒪_{\mathrm{string}}^z`$, because it does not change the sign. Thus $`z_L`$ turns out to be more rational order parameter to describe VBS ground states in a unified way.

In summary, we have introduced $`z_L`$ given as the ground-state expectation value of the twist operator (3) as an order parameter to characterize various ground states in quantum spin chains. Especially, $`z_L`$ changes the sign according to the configuration of valence bonds. This property enables us to determine the critical point between different VBS states by observing $`z_L=0`$. We have demonstrated this theory by using the QMC simulation for the successive dimerization transitions of the BAHC, and determined the phase boundaries with quite high accuracy.

M.N. thanks M. Oshikawa and J. Voit for useful discussions, and Y. Oi Nakamura for encouragement. S.T. acknowledges support of the Swiss National Science Foundation. The simulations were performed on SGI 2800 at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.


text:
# 1 Introduction

## 1 Introduction

It is generally agreed that the tremendous densities reached in the centers of neutron stars provide a high pressure environment in which numerous particles processes are likely to compete with each other. These processes range from the generation of hyperons to quark deconfinement to the formation of kaon condensates and H-matter . Another striking possibility concerns the formation of absolutely stable strange quark matter, a configuration of matter even more stable than the most stable atomic nucleus, iron. In the latter event all neutron stars would in fact be strange (quark matter) stars , objects largely composed of pure strange quark matter, eventually enveloped in a thin nuclear crust made up of ordinary, hadronic matter.

There has been much recent progress in our understanding of quark matter, culminating in the discovery that if quark matter exists it will be in a color superconducting state . The phase diagram of such matter appears to be very complex . At asymptotic densities the ground state of QCD with a vanishing strange quark mass is the color-flavor locked (CFL) phase. This phase is electrically neutral in bulk for a significant range of chemical potentials and strange quark masses . If the strange quark mass is heavy enough to be ignored, then up and down quarks may pair in the two-flavor superconducting (2SC) phase. Other possible condensation patters are the recently discovered CFL–$`K^0`$ phase  and the color-spin locked (2SC+s) phase . The magnitude of the gap energy lies between $`50`$ and $`100`$ MeV. Color superconductivity thus modifies the equation of state (eos) at the order $`(\mathrm{\Delta }/\mu )^2`$ level, which is only a few percent. Such small effects can be safely neglected in present determinations of models for the eos of neutron stars and strange quark matter stars. There has been much recent work on how color superconductivity in neutron stars could affect their properties . These studies revealed that possible signatures include the cooling by neutrino emission, the pattern of the arrival times of supernova neutrinos, the evolution of neutron star magnetic fields, rotational (r-mode) instabilities, and glitches in rotation frequencies. In this review I shall complement this list by reviewing several, most intruigung astrophysical implications connected with the possible absolute stability of strange quark matter. (Surface properties of strange matter are discussed in Usov’s paper elsewhere in this volume.) This is followed by a discussion of two astrophysical signals that may point at the existence of quark matter in both isolated neutron stars as well as in neutron stars in low-mass x-ray binaries (LMXBs). We recall that a convincing discovery of quark matter in neutron stars would demonstrate that strange quark matter is not absolutely stable, ruling out the absolute stability of strange quark matter and the existence of strange quark stars, for it is not possible for neutron stars to contain quark matter cores and strange matte quark stars to both be stable .

## 2 Nuclear crusts on strange matter stars

Since stars in their lowest energy state are electrically charge neutral to very high precision, any net positive quark charge must be balanced by leptons. As a general feature, there is only very little need for leptons, since charge neutrality can be achieved essentially among the quarks themselves. (This is specifically the case for superconducting CFL quark matter in the asymptotic limit.) If electrons form a component of absolutely stable strange quark matter, their presence is crucial for the possible existence of a nuclear curst on such matter . The reason being that the electrons, which are bound to strange matter by the Coulomb force rather than the strong force, extend several hundred fermi beyond the surface of strange matter. Associated with this electron displacement is a very strong electric dipole layer which can support, out of contact with the surface of the strange matter, a crust of nuclear material, which it polarizes. The maximal possible density at the base of the crust (inner crust density) is determined by neutron drip, $`ϵ_{\mathrm{drip}}=4.3\times 10^{11}\mathrm{g}/\mathrm{cm}^3`$, at which neutrons begin to drip out of the nuclei and form a free neutron gas. Being electrically charge neutral, the neutrons do not feel the repulsive Coulomb force and hence would gravitate toward the quark matter core, where they become converted into strange matter. Neutron drip thus sets a strict upper limit on the crust’s maximal inner density. The actual value may be slightly smaller though . The somewhat complicated situation of the structure of a strange matter star with crust can be represented by a proper choice of eos , which consists of two parts. At densities below neutron drip it is represented by the low-density eos of charge-neutral nuclear matter, for which we use the Baym-Pethick-Sutherland eos. The star’s strange matter core is described by the bag model.

## 3 Complete sequences of strange matter stars

Since the nuclear crusts surrounding the cores of strange stars are bound by the gravitational force rather than confinement, the mass-radius relationship of strange matter stars with crusts is qualitatively similar to the one of purely gravitationally bound stars, neutron stars and white dwarfs, as illustrated in Fig. 1 .

The strange star sequences are computed for the maximal possible inner crust density, $`ϵ_{\mathrm{crust}}=ϵ_{\mathrm{drip}}`$, as well as for an arbitrarily chosen, smaller value of $`ϵ_{\mathrm{crust}}=10^8\mathrm{g}/\mathrm{cm}^3`$, which may serves to demonstrate the influence of less dense crusts on the mass-radius relationship . From the maximum mass star (dot in the upper left corner), the central density decreases monotonically through the sequence in each case. The stars located along the dashed line represent unstable configurations . The fact that strange stars with crusts tend to possess somewhat smaller radii than neutron stars implies smaller mass shedding (Kepler) periods $`P_\mathrm{K}`$ for the former. This is already indicated by the classical expression $`P_\mathrm{K}=2\pi \sqrt{R^3/M}`$ and carries over to the full general relativistic determination of $`P_\mathrm{K}`$ ,

$$P_\mathrm{K}=2\pi \left(\omega +\frac{\omega ^{}}{2\psi ^{}}+e^{\nu \psi }\sqrt{\frac{\nu ^{}}{\psi ^{}}+\left(\frac{\omega ^{}}{2\psi ^{}}e^{\psi \nu }\right)^2}\right)^1.$$
(1)
This expression is to be computed simultaneously in combination with Einstein’s field equations for a rotating compact body ,

$$R^{\kappa \lambda }\frac{1}{2}g^{\kappa \lambda }R=\mathrm{\hspace{0.33em}8}\pi T^{\kappa \lambda }(ϵ,P(ϵ)).$$
(2)
It is found that, due to the smaller radii of strange stars, the complete sequence of such objects, and not just those close to the mass peak as is the case for neutron stars, can sustain extremely rapid rotation . In particular, model calculations indicate for a strange star with a typical pulsar mass of $`1.45M_{}`$ Kepler periods in the range of $`0.55\mathrm{msec}\stackrel{<}{_{}}\mathrm{P}_\mathrm{K}\stackrel{<}{_{}}0.8\mathrm{msec}`$, depending on the thickness of the nuclear curst and the bag constant . This range is to be compared with $`P_\mathrm{K}1\mathrm{msec}`$ obtained for neutron stars of the same mass.

The minimum-mass configurations of the sample strange star sequences in Fig. 1 have masses of about $`0.017M_{}`$ (about 17 Jupiter masses) and $`10^4M_{}`$, depending on the valuer of $`ϵ_{\mathrm{crust}}`$. For inner crust densities smaller than $`10^8\mathrm{g}/\mathrm{cm}^3`$ one obtains stable strange matter stars that can be by orders of magnitudes lighter than Jupiters. If abundant enough, these light strange stars could be seen by the gravitational microlensing searches . Strange stars located to the right of the minimum mass configuration of each sequence consist of small strange cores, typically smaller than about 3 km, surrounded by nuclear crusts (ordinary white dwarf matter) that are thousands of kilometers thick. Such objects are called strange dwarfs. Their cores have shrunk to zero at the crossed points. What is left are ordinary white dwarfs with central densities equal to the inner crust densities of the former strange dwarfs. A stability analysis of strange stars against radial oscillations  shows that all strange dwarf sequences that terminate at stable ordinary white dwarfs are stable against radial oscillations. Strange stars that are located to the left of the mass peak of ordinary white dwarfs (solid dot in upper right corner), however, are unstable against oscillations and thus cannot exist in nature. So, in sharp contrast to neutron stars and white dwarfs, the branches of strange stars and strange dwarfs are stably connected with each other . Finally we would like to stress that strange dwarfs with $`10^9\mathrm{g}/\mathrm{cm}^3<ϵ_{\mathrm{crust}}<4\times 10^{11}\mathrm{g}/\mathrm{cm}^3`$ form entire new classes of stars that contain nuclear material up to $`4\times 10^4`$ times denser than in ordinary white dwarfs of average mass, $`M0.6M_{}`$ (central density $`10^7\mathrm{g}/\mathrm{cm}^3`$). The entire family of such strange stars owes its stability to the strange core. Without the core they would be placed into the unstable region between ordinary white dwarfs and neutron stars .

Until recently, only rather vague tests of the theoretical mass-radius relation of white dwarfs have been possible. This has changed because of the availability of new data emerging from the Hipparcos project . These data allow the first accurate measurements of white dwarf distances and, as a result, establishing the mass-radius relation of such objects empirically. Figure 2 shows a

comparison of several data from the Hipparcos project with the mass-radius relationships of strange dwarfs (solid lines) and ordinary white dwarfs computed for different compositions.

## 4 Post-glitch behavior of strange stars

Ordinary neutron stars older than a few months have crusts made of a crystal lattice or an ordered inhomogeneous medium reaching from the surface down to regions with a density of $`2\times 10^{14}\mathrm{g}/\mathrm{cm}^3`$. This crust contains only a few percent of the total moment of inertia. Strange stars, in contrast, can only support a crust with a density below neutron drip ($`4.3\times 10^{11}\mathrm{g}/\mathrm{cm}^3`$), for reasons discussed in section 2. Such a strange star crust contains at most $`10^5`$ of the total moment of inertia. This is an upper bound, since the strange star may have no crust at all, depending on its prior evolution. The difference in the moment of inertia stored in the crust of neutron stars and strange stars seems to pose significant difficulties for explaining the glitch phenomenon observed in radio pulsars with models based on strange stars. Glitches are observed as sudden speed-ups in the rotation rate of pulsars. The fractional change in rotation rate $`\mathrm{\Omega }`$ is $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }10^610^9`$, and the corresponding fractional change in the spin-down rate $`\dot{\mathrm{\Omega }}`$ is of order $`\mathrm{\Delta }\dot{\mathrm{\Omega }}/\dot{\mathrm{\Omega }}10^210^3`$. In  it has been demonstrated that strange stars with crusts obey

$$\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}\left(10^5\mathrm{to}10^3\right)f,(0<f<1),$$
(3)
where $`f`$ represents the fraction of the crustal moment of inertia that is altered in the quake, $`|\mathrm{\Delta }I|=fI_{\mathrm{crust}}`$. Since observed glitches have relative frequency changes of $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }10^9`$ to $`10^6`$, a change in the crustal moment of inertia by less than 10% would cause a giant glitch ($`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }10^6`$) even in the least favorable case. Of course there remains the question of whether there can be a sufficient build up of stress and also of the recoupling of crust and core which involves the healing of the pulsar period. This is probably a very complicated process that does not simply involve the recoupling of two homogeneous substances. For $`\mathrm{\Delta }\dot{\mathrm{\Omega }}/\dot{\mathrm{\Omega }}`$ we established that

$$\frac{\mathrm{\Delta }\dot{\mathrm{\Omega }}}{\dot{\mathrm{\Omega }}}>(10^1\mathrm{to}10)f,$$
(4)
yielding a small $`f`$ value as before, namely $`f<10^4\mathrm{to}10^1`$. We have used measured values of the ratio $`(\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega })/(\mathrm{\Delta }\dot{\mathrm{\Omega }}/\dot{\mathrm{\Omega }})10^6\mathrm{to}10^4`$ for the Crab and Vela pulsars respectively. So the observed range of the fractional change in $`\dot{\mathrm{\Omega }}`$ is consistent with the crust having the small moment of inertia calculated and the quake involving only a small fraction, $`f`$, of that, just as in Eq. (3). Nevertheless, without undertaking a study of whether the nuclear solid crust on strange stars could sustain a sufficient buildup of stress before cracking to account for such a sudden change in relative moment of inertia, or whether the healing-time and intervals between glitches can be understood, one cannot say definitely that strange stars with a nuclear solid crust can account for any complete set of glitch observations for a particular pulsar.

## 5 Formation of strange dwarfs

At present there is neither a well-studied model for the formation of hypothetical strange dwarfs, nor exists a study that determines their abundance in the universe. One possible scenario would be the formation of strange dwarfs from main sequence progenitors that have been contaminated with strange nuggets over their lifetimes. We recall that the capture of strange matter nuggets by main sequence stars is an inevitable consequence if strange matter were more stable than hadronic matter  because then the Galaxy would be filled with a flux of strange nuggets which would be acquired by every object they come into contact with, i.e. planets, neutron stars, white dwarfs, and main sequence stars. Naturally, due to the large radii of the latter, they arise as ideal large-surface long integration time detectors for the strange matter flux. Nuggets that are accreted onto neutron stars and white dwarfs, however, never reach their centers, where the gravitational potential is largest, because they are stopped in the lattice close to the surface due to the large structural energy density there. This prevents such stars from building up a cores of strange matter. The situation is different for main sequence stars which are diffuse in comparison with neutron stars and white dwarfs. In this case the accreted nuggets may gravitate to the star’s core, accumulate there and form a strange matter core that grows with time until the star’s demise as a main sequence star. Another plausible mechanism has to do with primordial strange matter bodies. Such bodies of masses between $`10^2`$ and $`1M_{}`$ may have been formed in the early universe and survived to the present epoch . Such objects will occasionally be captured by a main sequence star and form a significant core in a single and singular event. The core’s baryon number, however, cannot be significantly larger than $`5\times 10^{31}(M/M_{})^{1.8}`$ where $`M`$ is the star’s mass. Otherwise a main sequence star is not capable of capturing the strange matter core . Finally we mention that in the very early evolution of the universe lumps of hot strange matter will evaporate nucleons which are plausibly gravitationally bound to the lump. The evaporation will continue until the quark matter has cooled sufficiently. Depending on the original baryon number of the quark lump, a strange star or dwarf, both with nuclear crusts will have been formed.

## 6 Quark matter in isolated neutron stars

Whether or not quark deconfinement occurs in neutron stars makes only very little difference to their static properties, such as the range of possible masses and radii, which renders the detection of quark matter in such objects extremely complicated. This turns out to be strikingly different for rotating neutron stars (i.e. pulsars) which develop quark matter cores in the course of spin-down. The reason being that as such stars spin down, because of the emission of magnetic dipole radiation and a wind of electron-positron pairs, they become more and more compressed. For some rotating neutron stars the

mass and initial rotational frequency may be just such that the central density rises from below to above the critical density for dissolution of baryons into their quark constituents. This effects the star’s moment of inertia dramatically , as shown in figure 4. Depending on the ratio at which quark and normal matter change with frequency, the moment of inertia can decrease very anomalously, and could even introduce an era of stellar spin-up lasting for $`10^8`$ years . Since the dipole age of millisecond pulsars is about $`10^9`$ years, we may roughly estimate that about 10% of the $`25`$ solitary millisecond pulsars presently known could be in the quark transition epoch and thus could be signaling the ongoing process of quark deconfinement! Changes in the moment of inertia reflect themselves in the braking index, $`n`$, of a rotating neutron star, as can be seen from $`(I^{}\mathrm{d}I/\mathrm{d}\mathrm{\Omega },I^{\prime \prime }\mathrm{d}^2I/\mathrm{d}\mathrm{\Omega }^2)`$

$$n(\mathrm{\Omega })\frac{\mathrm{\Omega }\ddot{\mathrm{\Omega }}}{\dot{\mathrm{\Omega }}^2}=3\frac{3I^{}\mathrm{\Omega }+I^{\prime \prime }\mathrm{\Omega }^2}{2I+I^{}\mathrm{\Omega }}.$$
(5)
The right-hand-side of this expression reduces to the well-known canonical constant $`n=3`$ if $`I`$ is independent of frequency. Evidently, this is not the case for rapidly rotating neutron stars, and it fails completely for stars that experience pronounced internal changes (phase transitions) which alter the moment of inertia significantly. Figure 4 illustrates this for the moment of inertia of the neutron star of Fig. 4. Because of the changes in $`I`$, caused

by the gradual transition of hadronic matter into quark matter, the braking index deviates dramatically from 3 at the transition frequency, when pure quark matter is generated. Such dramatic anomalies in $`n(\mathrm{\Omega })`$ are not known for conventional neutron stars, because their moments of inertia appear to vary smoothly with $`\mathrm{\Omega }`$ . The future astrophysical observation of such anomalies in the braking behavior of pulsars may thus be interpreted as a signal for quark deconfinement in neutron stars.

## 7 Quark matter in accreting x-ray neutron stars

Accreting x-ray neutron stars provide a very interesting contrast to the spin-down of isolated neutron stars discussed in sect. 6. These x-ray neutron stars are being spun up by the accretion of matter from a lower-mass ($`M\stackrel{<}{_{}}0.4M_{}`$), less-dense companion. If the critical deconfinement density falls within that of canonical pulsars, quark matter will already exist in them but will be “spun out” of x-ray stars as their frequency increases during accretion. This scenario has been modeled in  and will be discusses next. The spin-up torque experienced by a neutron star causes a change in the stars’ angular momentum that is described by the relation

$`{\displaystyle \frac{dJ}{dt}}=\dot{M}\stackrel{~}{l}(r_\mathrm{m})N(r_\mathrm{c}),`$ (6)
where $`\dot{M}`$ denotes the accretion rate and $`\stackrel{~}{l}(r_\mathrm{m})=\sqrt{Mr_\mathrm{m}}`$ is the angular momentum added to the star per unit mass of accreted matter. The quantity $`N(r_\mathrm{c})=\kappa \mu ^2r_\mathrm{c}^3`$ stands for the magnetic plus viscous torque terms, with $`\mu R^3B`$ the star’s magnetic moment. The quantities $`r_\mathrm{m}`$ and $`r_\mathrm{c}`$ denote the radius of the inner edge of the accretion disk and the co-rotating radius, respectively, and are given by $`(\xi 1)`$ $`r_\mathrm{m}=\xi r_\mathrm{A}`$ and $`r_\mathrm{c}=\left(M\mathrm{\Omega }^2\right)^{1/3}`$. Accretion will be inhibited by a centrifugal barrier if the neutron star’s magnetosphere rotates faster than the Kepler frequency at the magnetosphere. Hence $`r_\mathrm{m}<r_\mathrm{c}`$, otherwise accretion onto the star will cease. The Alfén radius $`r_\mathrm{A}`$, where the magnetic energy density equals the total kinetic energy of the accreting matter is defined by $`r_\mathrm{A}=\mu ^{4/7}(2M\dot{M}^2)^{1/7}`$. The rate of change of a star’s angular frequency $`\mathrm{\Omega }=2\pi /\nu `$ $`(=J/I)`$ then follows from Eq. (6) as

$$I(t)\frac{d\nu (t)}{dt}=\frac{\dot{M}\stackrel{~}{l}(t)}{2\pi }\nu (t)\frac{dI(t)}{dt}\kappa \mu (t)^2r_\mathrm{c}(t)^3,$$
(7)
with the explicit time dependences as indicated. Evidently, the second term on the right-hand-side of Eq. (7) depends linearly on $`\mathrm{\Omega }`$ while the third terms grows quadratically with $`\mathrm{\Omega }`$. The temporal change of the moment of inertia of accreting neutron stars which undergo phase transitions is crucial . This also renders the calculation of the moment of inertia, given by

$`I(t)=2\pi {\displaystyle _0^\pi }𝑑\theta {\displaystyle _0^{R(\theta ;t)}}𝑑re^{\lambda +\mu +\nu +\psi }{\displaystyle \frac{ϵ+P}{e^{2\nu 2\psi }(\omega \mathrm{\Omega })^2}}{\displaystyle \frac{\mathrm{\Omega }\omega }{\mathrm{\Omega }}},`$ (8)
very cumbersome, for each quantity on the right-hand-side varies accordingly during stellar spin-up. (We assume rigid body rotation, which renders the star’s frequency $`\mathrm{\Omega }`$ constant throughout the star.) The solution of equation (7) is shown in Fig. 6. The result is most striking. One sees that quark matter remains relatively dormant in the stellar core until the star has been spun up to frequencies at which the central density is about to drop below the threshold density at which quark matter exists. As known from Fig. 4, this manifests itself in a significant increase of the star’s moment of inertia. The angular momentum added to a neutron star during this phase of evolution is therefore consumed by the star’s expansion, inhibiting a further spin-up until the quark matter has been converted into a mixed phase of matter made up of hadrons and quarks. Such accreters, therefore, tend to spend a greater length of time in the critical frequencies than otherwise. There will be an anomalous number of accreters that appear at or near the same frequency, as shown in Fig. 6. This is what was found recently with the Rossi x-ray Timing Explorer (shaded area in Fig. 6). Quark deconfinement constitutes an attractive explanation for this anomaly , though alternative explanations were suggested too .

Acknowledgments: I am grateful to the organizers of Compact Stars in the QCD Phase Diagram, R. Ouyed and F. Sannino.

Discussion

S. Balberg (Hebrew University): The “anomalous” behavior in spin down (up) of an isolated (pair in a LMXB) neutron star attributed to the quark phase is really a feature of the equation of state, not of quark matter in itself. In other words, we need a relatively sharp feature in the equation of state to get a significant change in the moment of inertia, which would, in principle be due to other phenomena as well.

Weber: One needs a pronounced feature in the equation of state in order to get an anomaly in the breaking index (spin-up). Hyperons and meson condensates lead to features in the equation of state that seem to be way to week to cause such a striking anomaly . A softening alone, as caused by hyperons and/or bosons, is not sufficient. What is really required is a softening followed by a gradual stiffening at higher densities which is naturally obtained for hadron-quark matter gradually compressed to higher densities.

H. Heiselberg (NORDITA): How does $`M(R)`$ look for your quark matter star?

Weber: Glendenning and Kettner demonstrated the existence of non-identical neutron star twins, depending on whether or not the low-density sequence of neutron stars terminates at central densities that fall close to the end of the mixed phase of quarks and hadrons .