text:
# Quantum feedback control of a solid-state two-level system

## Abstract

We have studied theoretically the basic operation of a quantum feedback loop designed to maintain the desired phase of quantum coherent oscillations in a two-level system. Such feedback can suppress the dephasing of oscillations due to interaction with environment. Prospective experiments can be realized using metallic single-electron devices or GaAs technology.

The principle of feedback control is used in a wide variety of physical and engineering problems. In particular, it can be applied in a straightforward way to tune the oscillation phase of a harmonic oscillator in order to achieve a desired synchronization with some reference oscillator. An intriguing and fundamental question is whether continuous feedback can be used to control quantum systems; for instance, if it is possible or not to tune the phase of quantum coherent (Rabi) oscillations in a two-level system (TLS).

At first sight the quantum feedback seems to be impossible because according to the “orthodox” collapse postulate  the quantum state is abruptly destroyed by the act of measurement. However, as was shown two decades ago, in particular by Leggett,  in a typical solid-state setup the collapse of a TLS state should be considered as a continuous process rather than as instantaneous event. The reason is typically weak coupling between the quantum system and the detector and also the finite noise of the detector, so that it takes some time until acceptable signal-to-noise ratio is reached and the measurement can be regarded as completed.

While the Leggett’s theory as well as the majority of similar approaches can describe only ensembles of quantum systems, the theory describing the gradual collapse of a single solid-state TLS was developed only recently. (A similar problem in optics was solved much earlier – see, e.g. Refs.  and references in .) Basically, the theory says that the evolution of a single quantum system due to continuous measurement is governed by the information continuously acquired from the detector. Similarly to classical probability, the Bayes formula  which naturally takes into account incomplete information from the detector, can still be applied to the density matrix of the measured quantum system; thus the formalism is called Bayesian.

In case of a poor detector the extra noise acting on the input disturbs the measured system stronger than the limit determined by the uncertainty principle; this leads to gradual decoherence of the measured system. In contrast, when measured with a good (quantum-limited) detector, the quantum system does not loose the coherence (even though the quantum state evolves randomly); moreover, its density matrix can be gradually purified that basically means acquiring as much information about the system as permitted by quantum mechanics.

Since the Bayesian formalism allows us to monitor the continuous evolution of a quantum system in a process of measurement, this naturally gives rise to a possibility of continuous feedback control of a quantum system. In this paper we will study the operation of a feedback loop proposed in Ref.  and designed to maintain a desired phase of quantum coherent oscillations in a solid-state TLS. (Quantum feedback in optics has been proposed and studied earlier – see, e.g., Refs. .) In particular, we will study the dependence of the loop operation on the feedback factor $`F`$, the available bandwidth $`\tau _a^1`$, and the dephasing rate $`\gamma _e`$ due to environment.

As an example of the measurement setup (Fig. 1) we consider a TLS represented by a single electron in a double quantum dot (DQD), the location of which is measured by a quantum point contact (QPC) nearby in a way used in Ref. . If the electron is in the dot 2 (state $`|2`$) which is closer to QPC than dot 1, then the QPC tunnel barrier is higher and so the average current $`I_2`$ through QPC is smaller than the average current $`I_1`$ corresponding to electron in the dot 1 (state $`|1`$). Consequently, from the QPC current one gets information about the electron location. We consider a realistic case of weak response, $`\mathrm{\Delta }II_1I_2I_0(I_1+I_2)/2`$. In this case the measurement time $`S_I/2(\mathrm{\Delta }I)^2`$, which is necessary to achieve signal-to-noise ratio equal to 1 (here $`S_I`$ is the shot noise of the QPC current), is much larger than $`e/I_0`$, so the QPC current $`I(t)`$ is continuous on the measurement timescale and we do not need to consider individual tunneling events in QPC.

The evolution of the TLS density matrix $`\rho `$ during the measurement process is described within the Bayesian formalism by equations

$`\dot{\rho }_{11}=`$ $`\dot{\rho }_{22}=2{\displaystyle \frac{H}{\mathrm{}}}\text{Im}\rho _{12}+\rho _{11}\rho _{22}{\displaystyle \frac{2\mathrm{\Delta }I}{S_I}}[I(t)I_0],`$ (1)
$`\dot{\rho }_{12}=`$ $`\mathrm{i}{\displaystyle \frac{\epsilon }{\mathrm{}}}\rho _{12}+\mathrm{i}{\displaystyle \frac{\mathrm{H}}{\mathrm{}}}(\rho _{11}\rho _{22})`$ (3)
$`(\rho _{11}\rho _{22}){\displaystyle \frac{\mathrm{\Delta }I}{S_I}}[I(t)I_0]\rho _{12}\gamma \rho _{12},`$
where $`\epsilon `$ and $`H`$ are, respectively, the energy asymmetry and tunneling strength of the TLS \[the TLS Hamiltonian is $`_{TLS}=(\epsilon /2)(c_1^{}c_1c_2^{}c_2)+H(c_1^{}c_2+c_2^{}c_1)`$\], and $`\gamma =\gamma _d+\gamma _e`$ is the dephasing rate due to detector nonideality and coupling with environment. Theoretically, $`\gamma _d=0`$ when TLS is measured by a QPC; however, if instead of QPC we use a single-electron transistor, then dephasing $`\gamma _d`$ is always significant.

Notice that the ensemble dephasing rate $`\mathrm{\Gamma }=\gamma +(\mathrm{\Delta }I)^2/4S_I`$ is larger than $`\gamma `$ because of different evolution of the ensemble members due to random $`I(t)`$. Individual realizations can be simulated using the formula

$$I(t)I_0=(\rho _{11}\rho _{22})\mathrm{\Delta }I/2+\xi (t),$$
(4)
where $`\xi (t)`$ is the pure white noise with spectral density $`S_\xi =S_I`$. If Eqs. (1)–(3) are averaged over $`\xi (t)`$ (we use Stratonovich definition for stochastic differential equations), then we get usual ensemble-averaged equations for TLS evolution (terms proportional to $`\mathrm{\Delta }I`$ will disappear and $`\gamma `$ will be replaced by $`\mathrm{\Gamma }`$).

It is natural to characterize the effect of extra dephasing $`\gamma _d`$ by the detector ideality (efficiency) $`\eta 1/[1+\gamma _d4S_I/(\mathrm{\Delta }I)^2]`$. One can show that $`\eta =(\mathrm{}/2ϵ_d)^2`$ where $`ϵ_d`$ is the energy sensitivity of the detector including the effect of backaction but neglecting the correlation between output and backaction noises. So, an ideal case $`\eta =1`$ corresponds to a detector with quantum-limited sensitivity.

To realize a feedback loop (Fig. 1), we can monitor the TLS evolution using the detector current $`I(t)`$ plugged into Eqs. (1)–(3). Then the TLS state is compared with the desired state, and the difference signal is used to control the TLS parameters $`H`$ and/or $`\epsilon `$. In the example studied in this paper the feedback loop is designed to stabilize the quantum oscillations of the state of a symmetric TLS ($`\epsilon =0`$), so the desired evolution is $`\rho _{11}(t)=1\rho _{22}(t)=[1+\mathrm{cos}(\mathrm{\Omega }t)]/2`$, $`\rho _{12}(t)=\rho _{21}^{}(t)=\mathrm{i}\mathrm{sin}(\mathrm{\Omega }\mathrm{t})/2`$, where the frequency is $`\mathrm{\Omega }=(4H^2+\epsilon ^2)^{1/2}/\mathrm{}=2H/\mathrm{}`$. As a difference (“error”) signal we use the phase difference $`\mathrm{\Delta }\varphi `$ between the desired value $`\varphi _0(t)=\mathrm{\Omega }t`$ and the monitored value $`\varphi (t)\mathrm{arctan}\{2\text{Im}\rho _{12}(t)/[\rho _{11}(t)\rho _{22}(t)]\}`$. This difference is used to control the TLS parameter $`H`$ (controlling the barrier height of DQD); we assume a linear dependence: $`H_{fb}=H[1F\times \mathrm{\Delta }\varphi ]`$, where $`F`$ is the dimensionless feedback factor.

In this paper we neglect additional time delay in the feedback network, however, we take into account the finite bandwidth of a line carrying detector current (that is the critical parameter for a possible experiment). More specifically, we average the current $`I(t)`$ with a rectangular window of duration $`\tau _a`$, $`I_a(t)\tau _a^1_{t\tau _a}^tI(t^{})𝑑t^{}`$, before plugging it into Eqs. (1)–(3), so that the “available” density matrix $`\rho _a(t)`$ differs from the “true” density matrix $`\rho (t)`$. Also, to compensate for the corresponding implicit time delay, we use $`\mathrm{\Delta }\varphi =\varphi _a\mathrm{\Omega }(t\kappa \tau _a)`$ with $`\kappa =1/2`$ (we tried various $`\kappa `$ and found that $`\kappa =1/2`$ provides the best operation of the feedback loop).

Let us start with the case of ideal detector, $`\eta =1`$, absence of extra environment, $`\gamma _e=0`$, and infinite bandwidth, $`\tau _a=0`$. Figure 2 shows numerically calculated correlation function $`K_z(\tau )z(t+\tau )z(t)`$ where $`z\rho _{11}\rho _{22}`$, for several feedback factors: $`F=0,`$ 0.05, and 0.5. The curves are obtained using Monte Carlo simulation of the measurement process for moderately weak coupling between TLS and the detector: $`𝒞\mathrm{}(\mathrm{\Delta }I)^2/S_IH=1`$ (notice that the $`Q`$-factor of oscillations is equal to $`8/𝒞`$, so $`𝒞=1`$ is still a weak coupling). In absence of feedback ($`F=0`$) the correlation function decays to zero (Fig. 2) while for finite feedback factor the correlations remain for indefinitely long time (of course, assuming perfect reference oscillator which determines desired evolution). The nondecaying correlations show that the quantum feedback loop really provides the synchronization of quantum oscillations.

The degree of synchronization depends on the feedback factor $`F`$. One can see that for a moderate value of $`F=0.5`$ the synchronization is already very good \[the ideal case would be $`K_z(\tau )=\mathrm{cos}(\mathrm{\Omega }\tau )/2`$\]. For the case of moderate or good synchronization ($`𝒞/16F1`$, $`𝒞/81`$) we have derived the following analytical expression:

$$K_z(\tau )=\frac{\mathrm{cos}\mathrm{\Omega }\tau }{2}\mathrm{exp}\left[\frac{𝒞}{16F}\left(e^{2FH\tau /\mathrm{}}1\right)\right].$$
(5)
The correlation function $`K_I(\tau )I(t+\tau )I(t)`$ of the detector current $`I(t)`$ have somewhat similar dependence, however, it also has the decaying contribution due to correlation $`K_{z\xi }`$ and $`\delta `$-function contribution due to detector noise. The analytical result,

$`K_I(\tau )=`$ $`{\displaystyle \frac{S_I}{2}}\delta (\tau )+{\displaystyle \frac{(\mathrm{\Delta }I)^2}{4}}{\displaystyle \frac{\mathrm{cos}(\mathrm{\Omega }\tau )}{2}}\left(1+e^{2FH\tau /\mathrm{}}\right)`$ (7)
$`\times \mathrm{exp}\left[(𝒞/16F)\left(e^{2FH\tau /\mathrm{}}1\right)\right],`$
agrees well with the Monte Carlo results.

The spectral density $`S_I(\omega )`$ of the detector current can be obtained as a Fourier transform of $`K_I(\tau )`$. While in absence of feedback the quantum oscillations in TLS can provide only a moderate peak of $`S_I(\omega )`$ around frequency $`\mathrm{\Omega }`$ (the peak height cannot be larger than 4 times the noise pedestal), the feedback synchronization leads to the appearance of a $`\delta `$-function at the frequency of desired oscillations. (In principle the desired frequency can differ a little from $`\mathrm{\Omega }`$; however, in this case the performance of the feedback loop worsens.)

Besides the correlation function and spectral density, we have studied one more characteristic, $`D`$, of the synchronization degree. We define $`D`$ as the average scalar product of the unity-length vector on the Bloch sphere corresponding to the desired state and the vector corresponding to the actual state of TLS. The equivalent definition is $`D2\text{Tr}\rho \rho _d1`$, where $`\rho _d`$ is the density matrix of the desired pure state. Perfect synchronization corresponds to $`D=1`$. It is simple to show that in the case of weak coupling and unshifted desired frequency, $`D^2/2`$ is equal to the nondecaying amplitude of $`K_z(\tau )`$ dependence.

Upper solid line in Fig. 3 shows the dependence of $`D`$ on the feedback factor $`F`$ for $`𝒞=1`$ and $`\tau _a=0`$. One can see that $`D`$ is proportional to $`F`$ for small $`F`$ (“soft” onset of synchronization) and $`D`$ is asymptotically approaching 1 at large $`F`$. Using Eq. (5) it is simple to obtain expression $`D=\mathrm{exp}(𝒞/32F)`$ which is very close to numerical results for moderate and good synchronization (dashed line in Fig. 3).

Finite available bandwidth of the detector current $`I(t)`$ (finite averaging time $`\tau _a`$ in our formalism) worsens the performance of the quantum feedback loop. The solid lines in Fig. 3 show the dependence of the synchronization degree $`D(F)`$ for $`\tau _a/T=0`$, 1/3, and 2/3, where $`T=2\pi /\mathrm{\Omega }`$ is the oscillation period. Obviously, a significant information loss occurs when $`\tau _a`$ becomes comparable to $`T`$, leading to a decrease of $`D`$. The curves $`D(F)`$ saturate at large $`F`$ allowing us to introduce the dependence $`D_{max}(\tau )`$. Calculations for parameters of Fig. 3 show pretty good synchronization, $`D_{max}=0.993`$, for $`\tau _a=T/30`$, while $`D_{max}=0.98`$ for $`\tau _a=T/10`$, $`D_{max}=0.92`$ for $`\tau _a=T/3`$, and $`D_{max}=0.57`$ for $`\tau _a=2T/3`$.

The main potential practical importance of the quantum feedback is the ability to suppress the effect of TLS dephasing caused by interaction with environment (see Fig. 1). Solid lines in Fig. 4 show the dependence $`D(F)`$ for several magnitudes of the dephasing due to environment, $`d_e=0`$, 0.1, and 0.5, where $`d_e\gamma _e/[(\mathrm{\Delta }I)^2/4S_I]`$ is the ratio between TLS coupling to the environment and to the detector (we still assume an ideal detector). First of all, we see that feedback still maintains the phase coherence of TLS for infinitely long time. However, for finite $`d_e`$ the degree of synchronization $`D`$ saturates at the level less than unity. The dependence $`D_{max}(d_e)`$ is linear at small $`d_e`$, numerically we obtained $`D_{max}10.5d_e`$ for parameters of Fig. 4.

Notice that the solid lines shown in Figs. 3 and 4 are calculated assuming the TLS feedback control $`H_{fb}=H[1F\times \mathrm{\Delta }\varphi ]`$ even when $`H_{fb}`$ becomes negative (this is also an assumption for analytical results). To eliminate this unphysical assumption we have also performed numerical calculations with restriction $`H_{fb}>0`$ and with restriction $`H_{fb}>H/2`$. This leads to rather minor modifications of the presented curves (dashed and dotted lines in Fig. 4 show the results for $`d_e=0`$ and $`\tau _a=0`$). However, important difference is that $`D(F)`$ goes down at large $`F`$, so the optimum $`D_{max}`$ is achieved at some finite value of $`F`$. More detailed study of this problem will be presented elsewhere.

Besides the discussed feedback based on $`\mathrm{\Delta }\varphi `$ calculation, we have also studied much simpler feedback loop in which $`H_{fb}(t)/H1=F\{2[I_a(t)I_0]/\mathrm{\Delta }I\mathrm{cos}[\mathrm{\Omega }(t\tau _a/2)]\}\mathrm{sin}[\mathrm{\Omega }(t\tau _a/2)]`$ (in an experiment it would eliminate the necessity to solve fast the Bayesian equations). Quite surprisingly, such “naive” feedback can also provide a good phase synchronization of quantum oscillations if $`F/𝒞`$ is close to 1/4 (see dotted line in Fig. 3). However, it requires more careful choice of $`F`$ and $`\tau _a`$ than for the Bayesian feedback, and also suffers more significantly from the restriction on $`H_{fb}`$ variation. The results in more detail will be discussed elsewhere.

Experimentally, besides the realization of quantum feedback control of a DQD continuously measured by a QPC, one can also think about the TLS based on a single-Cooper-pair box measured by a single-electron transistor (see discussion in ). This realization can be preferable because of a rapid progress of metallic single-electronics technology. However, the problems are high output impedance of the single-electron transistor and its significant nonideality as a quantum detector. The third potential realization can be based on SQUIDs. For any realization the major problem is bandwidth: the feedback should be faster than the TLS dephasing due to environment. Because of that, the quantum feedback of a TLS should probably be attempted only after the realization of recently proposed Bell-type two-detector correlation experiment.

The authors thank A. J. Leggett and G. J. Milburn for useful discussions. The work was supported by NSA and ARDA under ARO grant DAAD19-01-1-0491.


format: text

text:
# Natural Quintessence and Large Extra Dimensions

## I Introduction

If current experimental indications are to be believed, the universe is composed of islands of luminous material embedded amongst darker baryonic matter, wrapped in the riddle of nonbaryonic dark matter inside the enigma of dark energy. Best fits to cosmological models give the last three of these in the rough relative proportions of few% : 30% : 70% . It is staggering that essentially nothing is known about either of the two dominant components to the universe’s present energy density.

The very existence of dark energy presents serious puzzles. The simplest account of the cosmic acceleration seems to be a non-zero cosmological constant $`\mathrm{\Lambda }`$. However the energy scale of the cosmological constant is problematic, requiring an energy density $`\rho _\mathrm{\Lambda }v^4`$ with $`v10^3`$ eV. The long-standing cosmological constant problem stems from the fact that vacuum fluctuations contribute to the cosmological constant, and simple estimates calculated by assuming a Planck-mass cutoff give $`\mathrm{\Lambda }`$ a value $`120`$ orders of magnitude greater than is be consistent with observations. This discrepancy is reduced somewhat in models where the natural cutoff is the supersymmetry breaking scale, but even in that case there is a very serious problem.

Many physicists are attracted by the idea that the cosmological constant problem could eventually be solved by some symmetry argument that would set $`\mathrm{\Lambda }`$ precisely to zero (although very interesting alternative views do exist. The observed cosmic acceleration creates problems with that line of reasoning, since it appears to indicate that $`\mathrm{\Lambda }`$ actually is non-zero. One way out is to assume that $`\mathrm{\Lambda }`$ really is zero, and that something else is causing the cosmic acceleration. This “something else” is usually called “dark energy” or “quintessence”.

For quintessence, typically one proposes the existence of a scalar field whose mass is small enough for it to be evolving in a cosmologically interesting way at present. So far, essentially all quintessence models present the following challenges or puzzles. P1: Why don’t radiative corrections destabilize the fantastic hierarchy between $`v`$, the scalar masses which are required ($`mv^2/M_p10^{33}`$ eV) and the other known scales of physics? (Only the first quintessence model attempts to address this problem) P2: Why isn’t the very long-range force mediated by the new light scalar detected in precision tests of gravity within the solar system?

We here report on a class of quintessence models which is motivated by recent efforts to understand the hierarchy problem, in which all observed particles except gravitons (and their supersymmetric partners) are trapped on a (3+1)-dimensional surface, or brane, in a higher-dimensional ‘bulk’ space . It turns out that within this framework the above two puzzles turn out to be resolved in a natural way. We believe ours is the only quintessence models to successfully address both these particular challenges*** Ref.  discusses another type of quantum correction which corrects classical rolling of the quintessence field off a local maximum. Those calculations do not apply to the type of potentials considered here..

A great virtue of the class of models which we explore is that they generally arise in the low-energy limit of theories which were originally proposed to solve the gauge hierarchy problem, without any cosmological applications in mind. As we shall see, the generic low-energy limit of these models looks like scalar-tensor gravity , with the four-dimensional scalars of relevance to cosmology coming from some of the extra-dimensional modes of the six-dimensional metric. We show here that although these models were not proposed with cosmology in mind, the couplings they predict for these scalars have following properties:

1. Good News: They can predict a late-time cosmology which provides a quintessence-style description of the dark energy. And they do so with scalars whose masses are extremely small, yet technically natural.
2. Good News: They can evade the constraints on the existence of very long-ranged forces because the couplings of ordinary matter to the the very light scalars are predicted to be field dependent, and so evolve over cosmological timescales. As has been observed elsewhere , once free to evolve, such couplings often like to evolve towards zero, implying acceptably small deviations from General Relativity within the solar system at the present epoch.
3. Good News: The expectation values of the scalar fields whose cosmology we explore determine the value of Newton’s constant, and the success of Big-Bang nucleosynthesis (BBN) implies strong limits on how much this can have changed between then and now. We find the evolution of Newton’s constant changes very little over cosmologically long times, for a reasonably wide range of initial conditions when we enter the BBN epoch: Newton’s constant tends not to change once the motion of the scalar field is either kinetic- or potential-dominated.
4. Bad News: The previous two cosmological successes are not completely generic, in the following sense. First, they depend on dimensionless couplings taking values which are of order $`ϵ1/50`$ or so. (Although couplings of this order of magnitude naturally do arise in the models of interest, we require some modest coincidences in their values to ensure sufficiently small evolution in Newton’s constant.) Second, although many initial conditions do produce an acceptable cosmology, this success is not completely generic. Attractor solutions exist which would cause too much variation in Newton’s constant if they described the cosmology between BBN and today. Furthermore, obtaining the present-day value of Newton’s constant requires adjusting the features of the scalar potential.

We present our results here as an example of a viable cosmology, based on using the brane-world’s extra-dimensional ‘radion’ mode as a quintessence field. Our ability to do so may itself be considered something of a surprise, since extra-dimensional aficionados had earlier examined the radion as a quintessence candidate, and found it wanting. Our conclusion here differs, because of our recognition that the radion couplings tend to be field dependent, in a way we make precise below.

We regard this work as a first exploration of the conditions which a successful cosmology must assume. We believe our results to be sufficiently promising to justify further exploration of this class of models, to see to what extent the unattractive features can be improved upon, and to ascertain what observations can test whether the universe indeed proceeds as we describe here.

Our presentation is organized as follows. In the next section we briefly describe the brane models we use, and their motivations and problems. The low-energy action relevant for cosmology is also presented in this section. Section 3 describes some of the cosmological solutions of the theory, focusing in particular on analytic expressions on which intuition can be based. More complicated, and realistic, solutions are then described numerically. We close in Section 4 with a discussion of our results.

## II The Model

Our point of departure is the brane-world scenario in which all observed particles except the graviton are trapped on a (3+1)-dimensional surface, or brane, in a higher-dimensional ‘bulk’ space. In particular, we imagine the scale of physics on the brane is $`M_b1`$ TeV and the bulk space is approximately six-dimensional, comprising the usual four, plus two which are compact but with radius as large as $`1/r10^3`$ eV . The large ratio $`M_br10^{15}`$ is phenomenologically required to reproduce an acceptably weak gravitational coupling in four dimensions, $`M_p=(8\pi G)^{1/2}M_b^2r10^{18}`$ GeV.

These kinds of models face two main difficulties, one theoretical and one phenomenological. Theoretically, these models trade the large hierarchy $`M_b/M_p=10^{15}`$ for the hierarchy $`M_br=10^{15}`$, but an explanation of this last value ultimately must be provided. Recently progress towards providing this explanation has been made, by proposing energetic reasons why the extra-dimensional radii might want to prefer to take large values . Our quintessence model does not rely on the details of these mechanisms, but our choice for the scalar field couplings is strongly motivated by the mechanism of ref. .

Phenomenologically, there are strong upper limits on how large $`M_br`$ can be. These limits come from supernova  and cosmological observations  which preclude the existence of Kaluza-Klein (KK) partners of the bulk-space graviton and its superpartners. We believe the present-day observational constraints on KK modes can be accommodated if only two dimensions have $`r`$ as large as $`(10^3\text{eV})^1`$, but prevent $`r`$ from being taken any larger. Although the cosmological constraints generally impose stronger upper limits on $`r`$ they depend more on the details of the assumed cosmology, and can be evaded under certain conditions (which we describe in more detail later).

### A The Low-Energy Action

We start by sketching the brane world scenario, and its low energy action. Our discussion closely follows that of ref. .

The couplings of the brane-world model relevant to cosmology are those of the six-dimensional metric, $`𝒢_{MN}`$. (Six dimensions are required in order to allow the extra dimensions to have radii as large as we need.) The leading terms in the derivative expansion for the six-dimensional action have the form $`S=S_B+S_b`$, where the Bulk action describing the metric degrees of freedom is:

$$S_B=\frac{M_b^4}{2}d^4xd^2y\sqrt{𝒢}+\mathrm{},$$
(1)
where $``$ denotes the scalar curvature built from the six-dimensional metric, and we write no cosmological constant for the bulk.This can be ensured by symmetries, such as in a supersymmetric compactification . Like everyone else we assume the vanishing of the effective four-dimensional cosmological constant, and our ideas do not add to the understanding of why this should be so.

The dependence of the brane action on its various matter fields is not crucial for our purposes, but what is important is that it depends only on the metric evaluated at the brane’s position:

$$S_b=d^4x\sqrt{g}_b(g_{\mu \nu },\mathrm{}),$$
(2)
where the ellipses denote all of the other fields on which $`_b`$ could depend.

The tree-level dimensional reduction of this action using the metric

$$𝒢_{MN}=\left(\begin{array}{cc}\widehat{g}_{\mu \nu }(x)& 0\\ 0& \rho ^2(x)h_{mn}(y)\end{array}\right),$$
(3)
with $`\rho =M_br`$, gives the leading contribution to the effective four-dimensional Lagrangian:

$$_{\mathrm{kin}}=\frac{M_b^4r^2}{2}\sqrt{\widehat{g}}\left[R(\widehat{g})2\left(\frac{r}{r}\right)^2+\frac{R(h)}{M_b^2r^2}\right],$$
(4)
where we have adopted the conventional normalization $`M_b^2d^2y\sqrt{h}=1`$. Here $`R(\widehat{g})`$ and $`R(h)`$ denote the curvature scalars computed from the metrics $`\widehat{g}_{\mu \nu }`$ and $`h_{mn}`$, respectively.

Quantum corrections modify this action in several important ways. First, the kinetic terms can become renormalized, introducing the possibility that their coefficients can acquire a dependence on $`r`$. Second, terms not written in eq. (4) can be generated, including in particular a no-derivative potential, $`U(r)`$. At large $`r`$, $`U(r)`$ can typically be expanded in powers of $`1/r`$, with the leading term falling off like $`1/r^p`$, for some $`p0`$. For instance, for toroidal compactifications, the leading contribution can arise due to the Casimir energy, which produces the potential $`U(r)=U_0/r^4+\mathrm{}`$.

We are in this way led to consider a scalar-tensor theory , for the couplings of $`r`$ and the four-dimensional metric, which we write as follows:

$`{\displaystyle \frac{}{\sqrt{\widehat{g}}}}`$ $`=`$ $`{\displaystyle \frac{M_b^4r^2}{2}}\left[A(r)R(\widehat{g})2B(r)\left({\displaystyle \frac{r}{r}}\right)^2\right]`$ (6)
$`U(r){\displaystyle \frac{_m(\widehat{g})}{\sqrt{\widehat{g}}}}.`$
As just described, the tree-level dimensional reduction of the 6d Einstein-Hilbert action predicts $`A=B=1`$ and $`U=0`$, and with these choices eq. (6) describes a Jordan-Brans-Dicke theory with coupling parameter $`\omega =1/2`$ . Notice that the matter Lagrangian, $`_m`$, describes the couplings of all of the known particles, and only depends on $`\widehat{g}_{\mu \nu }`$ but not on $`r`$. (Ordinary matter does not couple directly to $`r`$ in the brane-world picture because it is trapped on the brane, and does not carry stress energy in the extra dimensions.)

Notice that if the scale of $`U(r)`$ is completely set by $`r`$, as is true for the cases discussed in , then the mass which is predicted for this scalar is incredibly small: of order $`mv^2/M_p10^{33}`$ eV. It is the existence of such a small mass which allows $`r`$ to potentially play a cosmologically interesting role right up to the present epoch, and so to be a candidate for a quintessence field.

Since the predictions for $`A,B`$ and $`U`$ depend crucially on the approximation made in their derivation, in general we expect them to be given by more complicated functions. This turns out to be true, in particular, within the hierarchy-generating mechanism described in ref. . Since we use this proposal to motivate our choices for these functions, we briefly pause here to summarize its implications.

The proposal supposes that the spectrum of bulk states includes a six-dimensional scalar, $`\varphi `$, which counts amongst its interactions a renormalizable cubic self-coupling, $`g\varphi ^3`$. Since the coupling $`g`$ of this interaction is dimensionless, it can introduce corrections to the low-energy action which are logarithmic in $`r`$. Since all other couplings have negative mass dimension in 6d, their quantum contributions are suppressed by powers of $`1/r`$, implying the coupling $`g`$ can dominate at large $`r`$.

The bottom line of this scenario is that the dominant corrections to $`A,B`$ and $`U`$ are of the form $`1+c_1\alpha (r)+c_2\alpha ^2(r)+\mathrm{}`$, where the running loop-counting parameter is given by $`\alpha (r)=\alpha _0+\alpha _0^2\mathrm{log}(r/r_0)+\mathrm{}`$, with $`\alpha _0=g_0^2/(4\pi )^3`$. In this scenario the potential $`U(r)`$ is plausibly minimized for $`\alpha (r)1`$, implying a minimum for which there is the desired large hierarchy: $`M_br\mathrm{exp}[+1/\alpha _0]1`$. It is generic in this scenario that the quantum corrections to the other low-energy functions, $`A`$ and$`B`$, are also likely to be large.

For illustrative purposes when exploring cosmology, we choose the functions $`A,B`$ and $`U`$ which are suggested by this picture in the perturbative regime, for which $`A1+a\mathrm{log}(M_br)`$, $`B1+b\mathrm{log}(M_br)`$ and $`U\left(\frac{U_0}{r^4}\right)\left[1+c\mathrm{log}(M_br)\right]`$. Here $`a,b`$ and $`c`$ are small, since they are proportional to the dimensionless coupling constant, $`ϵg^2/(4\pi )^3`$, and higher orders in $`ϵ`$ give higher powers of $`\mathrm{log}(M_br)`$ although these have not been written explicitly. In the spirit of this picture we suppose $`ϵ`$ to be of order $`1/50`$, as required in order to ensure that $`U`$ is minimized when $`M_br10^{15}`$. For cosmological purposes, we shall be interested in the regime for which $`ϵ\mathrm{log}(M_br)<1`$ is negligible, but $`ϵ\mathrm{log}^2(M_br)`$ is not. We take $`U_0`$ to be generic in size: $`U_0O(M_b^4)`$.

### B Naturalness

Since naturalness is one of our motivations, we briefly summarize the arguments presented in ref.  as to why radiative corrections do not destabilize the scalar mass. The argument takes a different form for the contribution due to the integrating out of scales above, or below, $`v1/r10^3`$ eV.

For scales below $`v`$, the effective theory consists of a four-dimensional scalar-tensor theory coupled to itself, and to ordinary matter with couplings, $`\kappa `$, of gravitational strength: $`\kappa 1/M_p`$. Radiative corrections to the scalar mass within this effective theory are then generically acceptably small, being of order $`\delta m^2\kappa ^2\mathrm{\Lambda }^4`$, with ultraviolet cutoff $`\mathrm{\Lambda }v`$. These are therefore the same size as the tree-level mass, as is required if they are to be technically natural.

Radiative corrections to the masses coming from higher scales, between $`v`$ and $`M_b`$, are the dangerous ones. Generically, these are too large and it is the integrating out of these scales which introduces the naturalness problem to most quintessence models. Within the framework of ref. , the contributions from these scales are suppressed by supersymmetry, which is proposed to be broken at the scale $`M_b`$ on the brane on which observable particles live. (This ensures acceptably large splittings between the masses of ordinary particles and their superpartners.)

Supersymmetry suppresses corrections to the radion mass because it relates $`r`$ to the four dimensional graviton, which is massless. Indeed, above the scale $`v`$ the bulk space is really six-dimensional, and the radion is simply a component of the six-dimensional metric. As such, its properties are tied to those of the four-dimensional metric by the supersymmetry of the bulk space. Since the metric and its superpartners only couple to the supersymmetry breaking with gravitational strength, supersymmetry-breaking effects in the bulk sector are suppressed to be of order $`M_b^2/M_p1/rv`$.

Since the scalar potential for the radion cannot arise until supersymmetry breaks, it is naturally of order $`v^4F(\phi )`$ where $`\phi `$ is the dimensionless component of the six-dimensional metric which defines the radion. $`\phi `$’s kinetic term, on the other hand, comes from the 6d Einstein-Hilbert action, and so is allowed by unbroken supersymmetry. At tree level it is of order $`M_p^2(\phi )^2`$. and the coefficient $`M_p^2`$ is much larger than the contributions produced by integrating out physics at scales smaller than $`M_bM_p`$. Again we are led on dimensional grounds to a mass which is of order $`mv^2/M_p`$.

### C Cosmological Equations

The cosmology of the radion is most easily described by passing to the Einstein frame, for which the radion and metric are canonically normalized. For the choices of $`A=1+a\mathrm{log}(M_br)`$ and $`B=1+b\mathrm{log}(M_br)`$ the canonically normalized variables are:

$$\widehat{g}_{\mu \nu }=g_{\mu \nu }/[(M_br)^2A^2],\text{and}\mathrm{log}(M_br)=\widehat{a}\widehat{\chi }(1\widehat{b}\widehat{\chi }),$$
(7)
with $`\widehat{a}=(13a/4)/2`$, $`\widehat{b}=(ab)/16`$ and $`\widehat{\chi }=\chi /M_b`$. This choice of variables is called the Einstein frame, to distinguish it from the Jordan frame, as defined by eq. (6).A note on numerics: In the Jordan frame $`M_pM_b^2r`$, ordinary particle masses are $`m_pM_b`$ and Kaluza-Klein masses are $`m_{\mathrm{kk}}1/r`$, while in the Einstein frame $`M_pM_b`$, $`m_p1/r`$ and $`m_{\mathrm{kk}}1/(M_br^2)`$. Clearly, ratios of masses like $`m_{\mathrm{kk}}/m_p=1/(M_br)`$ are the same in either frame, but $`M_b=1`$ corresponds to using TeV units in the Jordan frame, but Planck units in the Einstein frame.

In the Einstein frame the action, eq. (6), becomes

$$\frac{}{\sqrt{g}}=\frac{M_b^2}{2}R\frac{1}{2}(\chi )^2V(\chi )\frac{_m(g,\chi )}{\sqrt{g}},$$
(8)
where

$`V(\chi )`$ $`=`$ $`{\displaystyle \frac{U[r(\chi )]}{(M_br)^4A^2}}`$ (9)
$`=`$ $`V_0\left[1+(v2a)\widehat{a}\widehat{\chi }+\mathrm{}\right]\mathrm{exp}[\lambda _0\widehat{\chi }(1\widehat{b}\widehat{\chi })],`$ (10)
and $`\lambda _0=8\widehat{a}`$. (Even though $`\widehat{b}=O(ϵ)`$ we do not expand the second term in the exponent in this expression because such an expansion presupposes $`ϵ\widehat{\chi }^21`$, which is not true for our applications.)

For cosmological applications we take $`ds^2=dt^2+a^2d\stackrel{}{x}^2`$ in the Einstein frame, in which case the Lagrangian (8) implies the usual Friedman equations:

$$3H^2M_b^2=\rho ,2\dot{H}M_b^2=(\rho +p),$$
(11)
where both the energy density, $`\rho `$, and pressure, $`p`$, receive contributions from $`\chi `$ and the radiation and matter sectors – $`\rho =\rho _r+\rho _m+\rho _\chi `$ and $`p=p_r+p_m+p_\chi `$, where $`p_r\frac{1}{3}\rho _r`$ and $`p_m0`$.

In the Einstein frame the matter contributions also depend explicitly on $`\chi `$:

$$\rho _m=\frac{\rho _{m0}}{a^3(M_br)},\rho _r=\frac{\rho _{r0}}{a^4},$$
(12)
where the quantities $`\rho _{m0}`$ and $`\rho _{r0}`$, may be related to the density of matter and radiation at the present epoch (which we choose to be $`a=1`$).

Eqs. (11) must be supplemented with the evolution equation for $`\chi `$, which is:

$$\ddot{\chi }+3H\dot{\chi }+V^{}(\chi )+\frac{\eta (\chi )}{M_b}T_{}^{\mu }{}_{\mu }{}^{}=0,$$
(13)
where $`T_{\mu \nu }`$ is the stress energy of the radiation and matter, and $`\eta (\chi )=\widehat{a}(1+a/22\widehat{b}\widehat{\chi })`$. The successful comparison of General Relativity with precision tests of gravity within the solar system require $`\eta ^2(\chi _0)<10^3`$ when evaluated at the present epoch, $`\chi =\chi _0`$ . Notice that since radiation satisfies $`T_{}^{\mu }{}_{\mu }{}^{}=0`$, only matter contributes to this last coupling, which is proportional to $`T_{}^{\mu }{}_{\mu }{}^{}=\rho _m`$.

## III Cosmology

We now show that cosmologically reasonable solutions may be found to eqs. (11) and (13), which satisfy all of the bounds to which the theory must be subject. Since our goal is to demonstrate the existence and illustrate some of the features of these solutions, our analysis is not systematic and should be regarded as only preliminary explorations of models of this type.

### A Constraints

Our main assumption is that the universe enters into an effectively four-dimensional radiation-dominated evolution at some temperature $`T_{}`$ higher than a few MeV, and so before Big Bang Nucleosynthesis (BBN). This assumption is a prerequisite for describing cosmology in terms of four-dimensional fields, and seems likely to be required if the success of standard BBN is not to be forsaken.

There is a limit to how high $`T_{}`$ may be chosen, since the four-dimensional picture requires the extra-dimensional KK modes to be cold. If $`T_{}`$ is too high, then these KK modes can be thermally excited too efficiently through their couplings with ordinary particles. If these modes are too abundant they can overclose the universe or decay too frequently into photons after recombination, and so contribute unacceptably to the diffuse gamma ray background . If the KK modes all decay into observable daughters, then these constraints require $`T_{}`$ not to be larger than the TeV scale, and can preclude having $`1/r10^3`$ eV.

As the authors of refs.  themselves point out, these constraints may be evaded, depending on the cosmology which is assumed at earlier times in the universe’s history. For instance, contributions to the diffuse gamma ray background is suppressed if the KK modes decay quickly enough into unobserved modes to be no longer present after the recombination epoch, when photons decouple from matter. Such fast decays can be arranged, for instance, if the KK modes decay more frequently into particles on other branes. Such decays may, in fact, be likely in cosmological scenarios for which the very early universe is dominated by a gas of branes and antibranes . (Another escape is possible if the extra-dimensional space is hyperbolic instead of toroidal, with its volume much greater than the appropriate power of its radius of curvature: $`Vr_c^2`$ . In this case KK modes can be very heavy, $`m_{\mathrm{kk}}1/r_c`$, even though the hierarchy problem is solved by having $`V`$ large.) Although our present interest is for $`T_{}<1`$ TeV, we have no difficulty imagining some variation on these themes being invoked if scenarios requiring $`T_{}>1`$ TeV were of interest.

Our second major constraint is the requirement that the quintessence field, $`\chi `$, not destroy the successes of BBN. This constraint arises because the success of BBN may be regarded as indicating that the hierarchy $`M_w/M_p(G_N/G_F)^{1/2}`$ is not much different at the epoch of nucleosynthesis than it is today. A sufficient condition for this to be true is to ensure that $`M_br`$ not to have evolved by more than about 10% of its present value between BBN and the present epoch. This constraint is much stronger than what is required by garden-variety quintessence models. It is very strong because it requires that the field $`r`$ cannot have rolled significantly over cosmologically long times. A special property of the cosmological solutions we shall examine is that they often have the property that $`r`$ becomes fixed over cosmologically large timescales. Motivated by the BBN constraints we are particularly interested in those cosmological solutions for which $`r`$ does not strongly vary after BBN. We shall see that such solutions can arise for reasonably generic initial conditions for the functions $`V(\chi )`$ and $`\eta (\chi )`$ we are using.

Next, we must demand agreement with present-day observations, which come in several important types. First, precision measurements of gravitational forces indicate that the coupling $`\eta (\chi )`$ must be $`<0.03`$ in order not to be in conflict with the success of General Relativity. Second, if $`\chi `$ is to provide an explanation of the dark energy, then its energy density must now be beginning to dominate that of matter, with $`w=p_\chi /\rho _\chi <\frac{1}{3}`$ if this energy density is to be currently accelerating the universe.

Finally, we require $`\chi `$ not to ruin the success of scale-invariant fluctuations in describing present-day observations of the cosmic microwave background radiation (CMB). Since a detailed analysis of the CMB in these models is beyond the scope of the present paper, we instead apply a poor-man’s bound by requiring cosmological properties not to be too different from the standard picture near the epochs of radiation/matter crossover and recombination.

### B Approximate Cosmological Solutions

Before laying out a cosmology which satisfies the above constraints, it is worth building intuition by outlining several approximate analytical solutions to the cosmological equations which turn out to control the features which are seen in the numerical results. Since the microphysics motivates using the potential: $`V(\chi )=V_pe^{\lambda \chi }`$, with $`V_p(\chi )`$, $`\lambda (\chi )`$ and the coupling $`\eta (\chi )`$ varying much more slowly as (polynomial) functions of $`\chi `$ than does the exponential $`e^{\lambda \chi }`$, we explore solutions to the equations when $`\lambda ,V_p`$ and $`\eta `$ are all constants. The cosmological evolution of such fields in the limit $`\eta 0`$ has been much studied , and our results go over to these in the appropriate limits.

There are two types of solutions in this limit which are of principal interest. These solutions may be classified according to which the terms dominate in the $`\chi `$ equation, eq. (13).

#### 1 Scaling Solutions

The first class of approximate solutions we consider are the scaling, or tracker, solutions, of which there are two different types. These solutions are characterized by having all energy densities scaling as a power of the scale factor, and either: ($`i`$) $`\ddot{\chi }3H\dot{\chi }V^{}\eta \rho _m`$ or, ($`ii`$) $`\ddot{\chi }3H\dot{\chi }\eta \rho _mV^{}`$. Both types of solutions can arise regardless of whether it is the scalar, radiation or matter which is dominating the evolution of the universe.

Assume, therefore, the scale factor scales with time as $`a/a_0=(t/t_0)^k`$, and so $`H=\dot{a}/a=k/t`$. From the Friedman equation, $`H^2\rho `$, we see that if the dominant energy density scales as $`\rho a^n`$, then we must have $`k=2/n`$.

($`i`$) $`\ddot{\chi }3H\dot{\chi }V^{}\eta \rho _m`$: In this regime we assume

$$Ve^{\lambda \widehat{\chi }}\frac{1}{r^8}\frac{1}{a^s},$$
(14)
and we find from the $`\chi `$ field equation, eq. (13), that $`s=n`$.

If radiation or matter dominate the energy density, then $`n=3`$ or $`n=4`$. In this case we find solutions for which $`K=\frac{1}{2}\dot{\chi }^2V(\chi )1/a^n`$, and so $`ra^{n/8}`$. Eq. (13) also implies that these scaling solutions are trackers in the sense that $`\rho _\chi `$ tracks at precisely the fraction $`n/\lambda ^2`$ of the dominant energy density as the universe expands. Solutions of this type only are possible if $`\lambda ^2>n`$. For $`\lambda ^2<n`$ the scalar energy density tends towards scalar domination of the universal expansion.

If it is $`\chi `$ itself which dominates the energy density, then $`k=2/\lambda ^2`$ and a scaling solution of the type we seek exists so long as $`\lambda ^2<6`$. In this case both $`K`$ and $`V`$ scale as $`1/a^{\lambda ^2}`$, and so $`ra^{\lambda ^2/8}`$.

These are the tracker solutions which are familiar from standard discussions of scalars evolving with exponential potentials . Clearly these tracker solutions imply a relatively rapid variation of $`r`$ as the universe evolves.

($`ii`$) $`\ddot{\chi }3H\dot{\chi }\eta \rho _mV^{}`$: We next consider the case where the $`\eta T_{}^{\mu }{}_{\mu }{}^{}`$ term dominates the scalar potential term in eq. (13). Notice that this can hold even if both $`\rho _r`$ and $`\rho _\chi `$ are much larger than $`\rho _m`$, so the tracker solution we here obtain can be relevant even when nonrelativistic matter makes up only a very small part of the universe’s total energy budget.

Under these assumptions, because the last term of eq. (13) dominates the potential term, we ask $`\rho _m`$ to scale like $`\ddot{\chi }`$ and $`H\dot{\chi }`$:

$$\rho _m\frac{1}{a^3r}\frac{1}{a^s},$$
(15)
which, when inserted into eq. (13) implies $`s=n`$ and so $`K\rho _m\rho 1/a^n`$. As usual, if the universe is radiation or matter dominated then $`n=4`$ or $`n=3`$. (Interestingly, both matter and radiation fall as $`1/a^4`$ within this kind of tracker solution in a radiation-dominated universe.) Alternatively, $`K\rho _m\rho 1/a^{\lambda ^2}`$, if the universe is $`\chi `$-dominated.

Although the scaling of scalar kinetic energy is the same for tracker solutions of type $`(i)`$ and type $`(ii)`$, they differ markedly in the scaling which they imply for $`r`$. In the present solution the scaling of $`r`$ vs $`a`$ follows from the known scaling of $`\rho _m`$. For radiation/matter domination, where $`\rho 1/a^n`$ we find $`ra^3a^n`$ and so $`ra^{n3}`$. Clearly $`ra`$ for radiation domination, but $`r`$ only varies logarithmically during matter domination! During a scalar-dominated phase, on the other hand, we have $`ra^{\lambda ^23}`$. Again, $`r`$ does not vary if $`\lambda ^2=3`$.

What is striking about these solutions where $`r`$ changes so slowly is that this happens in spite of its not sitting at the minimum of the potential. It happens, even though the radion’s kinetic energy dominates its potential energy, because the amount of roll is limited by the Hubble damping generated by the dominant energy component of the universe.

#### 2 Transient Solutions

What is perhaps even more striking than the trackers for which $`r`$ slowly varies is the behavior of $`r`$ during the transient evolution as the tracker – which is an attractor – is approached. Since the approach to the tracker solution may be approached from initial conditions for which $`KV`$ or $`KV`$, we consider each of these in turn.

$`KV`$ (Slow Roll): If $`KV`$ then it is a good approximation to neglect $`\ddot{\chi }`$ in eq. (13) relative to the other terms. If $`\eta \rho _m`$ is also negligible relative to $`V^{}`$, then the evolution is described by $`3H\dot{\chi }+V^{}=0`$, which has as solutions:

$$V(a)=\frac{2V_{\mathrm{exit}}a_{\mathrm{exit}}^p}{a^p+a_{\mathrm{exit}}^p},$$
(16)
where $`p=n`$ for radiation ($`n=4`$) and matter ($`n=3`$) domination, and $`p=\lambda ^2`$ for scalar domination. Along this trajectory the kinetic energy grows as $`K/V=\lambda ^2V/(6\rho )`$.

This solution implies $`VV_{\mathrm{exit}}`$ is a constant (and so also must be $`r`$) until $`aa_{\mathrm{exit}}`$, after which point the solution joins the tracker solution, for which $`\rho _\chi 1/a^p`$ and $`\rho _\chi /\rho =p/\lambda ^2`$.

$`KV`$ (Kinetic Domination): If $`K`$ dominates the scalar evolution, then $`V^{}`$ and $`\eta \rho _m`$ may be dropped relative to $`\ddot{\chi }`$ and $`3H\dot{\chi }`$. In this case the $`\chi `$ equation may again be solved analytically, with solution: $`\chi =C_1C_2\left(a_0/a\right)^{31/p}`$, where $`C_1,C_2`$ are integration constants with $`C_2`$ related to the initial velocity by $`C_2=(31/p)a_0(d\chi /da)_0`$. Here, as usual, $`p=2/n`$ for radiation/matter domination, and $`p=2/\lambda ^2`$ for scalar domination.

This solution implies the kinetic energy, $`K=\frac{1}{2}\dot{\chi }^2`$ varies as:

$$K(a)=K_0\left(\frac{a_0}{a}\right)^6,$$
(17)
regardless of the value of $`p`$.

The evolution of the potential energy is more interesting in these solutions. Consider first the case of radiation/matter domination, for which $`p=2/n`$ so

$$V(a)=V_0\mathrm{exp}\left\{c\left[\left(\frac{a_0}{a}\right)^{3n/2}1\right]\right\},$$
(18)
with $`c=(31/p)\lambda a_0(d\chi /da)_0`$.

As the universe expands according to this solution, $`\chi `$ carries an ever-smaller fraction of the energy density since $`K(a)1/a^6`$ falls faster than does the dominant matter or radiation, $`\rho 1/a^n`$. Eventually $`K`$ becomes comparable in size to $`V(a)`$, which for large $`a`$ asymptotes to a constant value: $`V(a)V_{\mathrm{}}=V_0e^c`$. At this point the approximation of dropping the scalar potential in eq. (13) breaks down, and we match onto the slow-roll solution discussed above. The amount of expansion which occurs before $`V`$ and $`K`$ are comparable is determined by $`a/a_0=(e^cK_0/V_0)^{1/6}`$.

What is remarkable about this solution is that the constant $`c`$ cannot be made arbitrarily large. To see this, notice that $`c`$ can be related to the initial kinetic energy by using $`K=\frac{1}{2}H^2a^2(d\chi /da)^2`$. Using the upper limit $`K3H^2`$, which is dictated by the inequality $`K\rho `$ and the Friedman equation, we see that $`c\sqrt{6}\lambda (31/p)`$. For $`\lambda =4`$ and $`p=1/2`$ (radiation domination) we have $`c4\sqrt{6}9.6`$.

This upper limit on $`c`$ implies $`V`$ cannot fall by more than a factor of $`e^ce^{\lambda \sqrt{6}}`$ during a radiation dominated epoch, despite the fact that $`r`$’s evolution is being dominated by its kinetic energy. Since $`V1/r^8`$ we see that $`r/r_0<e^{\lambda \sqrt{6}/8}=e^{0.3\lambda }`$ varies hardly at all during the entire kinetic roll. Physically the paradoxical upper limit on $`r`$ during a kinetic dominated roll arises due to the Hubble damping which increases with the kinetic energy, if the kinetic energy gets large. This observation has also been made in another context in ref. .

A similar result holds for the kinetic roll in the scalar-dominated phase, although with results which differ strongly depending on the value of $`\lambda `$. If $`\lambda ^2<6`$, the above analysis applies directly, with $`p=2/\lambda ^2`$. As before $`K1/a^6`$ and

$$V(a)=V_0\mathrm{exp}\left\{c\left[\left(\frac{a_0}{a}\right)^{3\lambda ^2/2}1\right]\right\}.$$
(19)
Provided $`\lambda ^2<6`$, $`V(a)`$ again asymptotes to $`V_0e^c`$, with $`c`$ once more subject to a very general upper bound. Once $`K`$ falls to be of order $`V`$, this solution matches onto a slow roll solution, and approaches the tracker solution, which exists for $`\lambda ^2<6`$.

If $`\lambda ^2>6`$, however, the large-$`a`$ limit of $`V(a)`$ is very different. In this case $`V(a)\mathrm{exp}[|c|(a/a_0)^{(\lambda ^26)/2}]`$, implying $`V(a)`$ falls very quickly to zero as the roll proceeds. This solution describes the attractor for the case $`\lambda ^2>6`$, which is a kinetic-dominated roll rather than a tracker solution.

#### 3 Realistic Complications

Although these solutions describe well most features of the numerical evolution we describe below, there are also important deviations. These deviations arise because $`\eta `$ and the quantities $`\lambda `$ and $`V_p`$ (of the potential $`V=V_pe^{\lambda \chi }`$) used in the numerical evolution are all polynomials in $`\chi `$. Motivated by the microphysical discussion given in section II we choose:

$$\eta =\eta _0+\eta _1\widehat{\chi };\lambda (\chi )=\lambda _0+\lambda _1\widehat{\chi };$$
(20)
$$\text{and}V_p=V_0+V_1\widehat{\chi }+\frac{V_2}{2}\widehat{\chi }^2,$$
(21)
with $`\eta _0`$, $`\lambda _0`$ and $`V_0`$ all taken to be $`O(1)`$; $`\eta _1,\lambda _1`$ and $`V_1`$ assumed to be $`O(ϵ)`$ and $`V_2`$ taken $`O(ϵ^2)`$, where $`ϵ1/50`$. (The second-order quantity $`V_2`$ is included here in order to ensure $`V`$ does not change sign for any value of $`\chi `$.)

Because these quantities vary more slowly with $`\chi `$ than does the exponential in $`V`$, their evolution over cosmological timescales is also slow, and the above approximate solutions provide a good description for most of the universe’s history. We have found that their $`\chi `$-dependence can nevertheless play an important role when constructing viable cosmologies, by allowing important deviations from the approximate solutions, including:

* The $`\chi `$-dependence of $`V_p`$ becomes important for those $`\chi `$ satisfying $`V_p^{}/V_p\lambda `$, since the potential $`V`$ can acquire stationary points near such points. These kinds of stationary points can slow or trap the cosmological $`\chi `$ roll, in much the same way as was exploited in ref. .
* The $`\chi `$-dependence of $`\lambda `$ can cause transitions between the domains of attraction of the various attractor solutions. For instance, in a scalar dominated era, the evolution of $`\chi `$ can drive $`\lambda `$ across the dividing point, $`\lambda ^2=6`$, forcing a transition between the scaling and kinetic-dominated attractor solutions.
* The $`\chi `$-dependence of $`\eta `$ is crucial for evading the bounds on long-range forces during the present epoch, because it allows $`\eta `$ to vanish for some choices of $`\chi `$. The bounds can be evaded provided the current value of $`\chi `$ is sufficiently close to such a zero of $`\eta `$. Furthermore, as described in detail in ref. , in many circumstances $`\chi `$’s cosmological evolution is naturally attracted to the zeros of $`\eta `$ as the universe expands.

Because of the natural hierarchy satisfied by the coefficients of these functions, the consequences of $`\chi `$-dependence generically occur for $`\chi O(1/ϵ)`$, corresponding to $`r`$ close to its present value.

### C A Realistic Cosmological Evolution

We next describe a realistic cosmology which exploits some of these properties. We present this model as an existence proof that cosmologically interesting evolution can really be found using the scalar properties which are suggested by the microphysics. We find that the main constraint which governs the construction of such solutions is the requirement that $`r`$ not vary appreciably between BBN and the present epoch.

In order to fix ideas, we imagine starting the universe off at $`T_{}>1MeV`$, perhaps after an epoch of earlier inflation, in a kinetic-dominated roll. We note that such a state is the generic endpoint if we start with $`\lambda 4>\sqrt{6}`$, as is motivated by the tree-level expression for the radion potential, and assume the scalar dominates the energy density of the universe. (As we will discuss briefly later, one might chose to make other assumptions about the initial conditions.)

With this assumption the scalar energy density must eventually fall below the density of radiation, thus initiating the radiation-dominated era. We assume that this radiation-domination begins shortly before BBN, so that nucleosynthesis occurs during radiation domination.

Our second assumption is dictated by the requirement that $`r`$ does not change appreciably between BBN and now. Although the kinetic-dominated roll before radiation domination will automatically ensure that $`VK`$, we ask that $`V`$ be within an order of magnitude of to its present value as radiation domination begins. This ensures that $`r`$ will have an initial value at the onset of radiation domination which is very close to its present-day value. We make this (fairly ugly) assumption to demonstrate the feasibility of models along the lines we are presenting, but hope to be able to eventually identify a more natural frameworks in which such a condition could be ensured by the system’s dynamics rather than as an initial condition.

For $`\lambda 4`$, once the universe becomes radiation dominated the scalar in the domain of validity of the tracker solutions described above. With the given initial conditions the scalar enters radiation domination in a kinetic-dominated roll, and so tries to approach a tracker solution by first damping its kinetic energy in a continued kinetic-dominated roll, and then performing a slow-roll at fixed $`r`$. We have seen that the universe expands by an amount $`(a/a_0)^2K_0/V_0`$ during the transient phase, before the tracker solution is reached, so given the low values of $`V_0`$ assumed at the onset of radiation domination, $`K_0/V_0`$ is very large, and so the transient evolution to the tracker can easily take the entire time from BBN to the present epoch.

We have seen on general grounds that the value of $`r`$ does not change by more than a factor of $`O(1)`$ during this transient phase, including both the transient’s kinetic-dominated and slow-rolling parts. The challenge in making this scenario realistic is to have this transient evolution survive right up to the present day, without first being intercepted by a tracker solution.

We have been able to find initial conditions and plausible couplings for which this is accomplished. One such is illustrated in Figures (1) through (4), which plot the evolution of various quantities from the onset of radiation domination to the present. The evolution described by this simulation uses the values suggested by the lowest-order microscopic action, $`\eta _0=1`$, $`\lambda _0=4`$ and $`V_0=1`$, plus the following correction terms: $`\eta _1=0.015`$, $`\lambda _1=0`$, $`V_1=0.030`$ and $`V_2=0.00046`$.

Figure (1) shows how the energy density in the universe is distributed between radiation, matter and the kinetic and potential energy of the scalar field, as a function of the universal scale factor, $`a`$, in the Einstein frame. Figure (2) plots the same information, but normalized as a fraction of the total universal energy density, and also plotting the scalar equation-of-state parameter, $`w=p_\chi /\rho _\chi `$. Figure (3) gives a plot of $`r`$ vs scale factor, showing that $`r`$ does not vary appreciably between BBN and the present epoch. Finally, Figure (4) shows the evolution of $`\eta `$ against scale factor, showing that $`\eta `$ evolves to sufficiently small values at the present epoch.

As is evident from the figures, this solution satisfies all of the cosmological and present-day bounds we were listed earlier. This was accomplished by choosing appropriately the constants governing the scalar interactions, by ensuring $`V`$ to have a minimum relatively near to a zero of $`\eta `$. This was required because although $`r`$ naturally does not evolve significantly for most of the universe’s history (after the onset of radiation domination), for pure exponential potentials it generically tends to join a tracker solution near the onset of matter domination. Because $`r`$ evolves too far once in this tracker solution, it must be avoided. The solution illustrated in the figure does so by having $`\chi `$ become snagged by the minimum of $`V`$ before entering into the matter-dominated phase.

## IV Discussion

We have shown that viable cosmologies may be built using the radion as the quintessence field. We use radion interactions which were suggested by a microphysical model which was recently proposed to naturally solve the hierarchy problem within a brane-world framework. The radion in this model is naturally light enough to play a role in late-time cosmology, and this small mass is technically natural in that it is not destabilized by radiative corrections. The model also satisfies all current bounds which constrain the existence of long-range scalar-mediated forces.

The success of the radion in filling the quintessence role is somewhat of a surprise, since radion models have long been believed to have phenomenologically unacceptable couplings. The model evades these problems because the couplings are predicted to be weakly field dependent, allowing the scalar couplings to evolve over cosmological timescales. The model evades all present-day bounds on new forces because the relevant couplings evolve towards small values at late times.

Although most of the parameters and initial conditions chosen were natural in size, the model does have two features which we believe need further improvement. First, although all scalar couplings were chosen with natural sizes, the precise values chosen were adjusted to arrange the scalar potential to have a minimum close to a zero of $`\eta `$. We would prefer to find an attractor solution, à la Damour and Nordtvedt , which more naturally draws $`\eta `$ to small values in the present epoch.

A second unsatisfactory feature of the cosmology presented was its reliance on the entering of radiation domination with $`r`$ having close to its present value. Although large values of $`r`$ are naturally generated by kinetic-dominated rolls during scalar dominated epochs, there is no natural reason why they should have precisely the current value. This problem might be resolved by a more developed understanding of radion initial conditions. Perhaps one could argue that there are no a-priori constraints on the initial conditions for $`r`$. In that picture, it may be the case that the largest amount of phase space that matches the observation that Newton’s constant is constant and non-zero today would be one where $`r`$ has the current value throughout the cosmic evolution. (In that reasoning, the constancy of $`G`$ is not a prediction, but an input.)

A natural way to alleviate both of these problems would be to consider more moduli than just the radion, since other moduli would be less constrained by the requirement that Newton’s constant not appreciably change since nucleosynthesis. We believe work along this lines to be worthwhile, given the model’s successful addressing of the naturalness issues of quintessence models, and the promising cosmology to which it leads despite its being devised to solve purely microphysical problems such as the hierarchy problem.

## V Acknowledgments

We would like to acknowledge useful conversations with A. Aguirre, S.Basu and S. Gratton, and the Aspen Center for Physics for providing the pleasant environs in which this work started. CB’s research is partially funded by NSERC (Canada), FCAR (Québec) and the Ambrose Monell Foundation. AA and CS are supported by DOE grant DE-FG03-91ER40674 and U.C. Davis.


text:
# Ground-State Phase Diagram of Quantum Heisenberg Antiferromagnets on the Anisotropic Dimerized Square Lattice

## I Introduction

Recently many low-dimensional antiferromagnets with excitation modes separated from a ground state by a finite energy gap have been synthesized, and effects of impurities or magnetic fields on those materials have been investigated experimentally in relation to the impurity-induced long-range order (LRO) and magnetic-field-induced LRO. For example, there are the $`S=1/2`$ quasi-one-dimensional (Q1D) Heisenberg antiferromagnet (HAF) with bond dimerization, CuGeO3 $`S=1`$ Q1D HAF’s, NENP,  NDMAZ,  NDMAP,  and PbNi2V2O8 and $`S=1`$ Q1D HAF’s with bond alternation, NTEAP  and NTENP.  Those materials have attracted our interest since they reveal various aspects of the quantum phase transition between the quantum-disordered spin-gapped phase and the classical (Néel) long-range-ordered phase.

Intrachain spin interaction with or without bond alternation and interchain interaction are considered to be the most basic ingredients to understand the quantum phase transitions mentioned above. More explicitly, they are expected to be modeled effectively by the quantum HAF on the anisotropic dimerized square lattice, which is described by the Hamiltonian:

$``$ $`=`$ $`{\displaystyle \underset{i,j}{}}𝐒_{2i,j}𝐒_{2i+1,j}+\alpha {\displaystyle \underset{i,j}{}}𝐒_{2i+1,j}𝐒_{2i+2,j}`$ (1)
$`+`$ $`J^{}{\displaystyle \underset{i,j}{}}𝐒_{i,j}𝐒_{i,j+1}.`$ (2)
Here $`𝐒_{i,j}`$ is the quantum spin operator at site $`(i,j)`$ on the square lattice. The first two terms in r.h.s. represent the one-dimensional antiferromagnetic (AF) Heisenberg chains with alternating coupling constants, 1 and $`\alpha `$ $`(0\alpha 1)`$, and the last term does the AF interchain exchange interaction ($`J^{}0`$). We choose the $`x`$-axis as being along the chain direction and the $`y`$-axis as in the perpendicular one. The bond arrangement of this model is shown in Fig. 1.

The ground state of decoupled chains, i.e., $`J^{}=0`$, has been well established. In particular, that of the $`S=1/2`$ chain  is the dimer state with a finite spin gap except for the uniform case ($`\alpha =1`$), which has a critical ground state. On the other hand, for the $`S=1`$ chain there exist two spin-gapped phases: the Haldane phase  at $`\alpha >\alpha _\mathrm{c}`$ and the dimer phase at $`\alpha <\alpha _\mathrm{c}`$. At the critical point $`\alpha =\alpha _\mathrm{c}`$ between these two phases the gap vanishes. The value of $`\alpha _\mathrm{c}`$ has been estimated to be 0.5879(6). In both cases, the critical point is considered to belong to the Gaussian universality class.

For the AF LRO to appear the higher-dimensionality effect, i.e., the interchain interaction $`J^{}`$, is indispensable. In most of the numerical works reported so far, the effect of interchain coupling has been examined in certain approximated or perturbed ways. For example, Sakai and Takahashi  estimated the critical strength, $`J_\mathrm{c}^{}`$, for the uniform case ($`\alpha =1`$) by the exact diagonalization method for the intrachain interactions combined with the mean-field approximation for the interchain interaction, and obtained $`J_\mathrm{c}^{(S=1/2)}=0`$ and $`J_\mathrm{c}^{(S=1)}0.025`$. More recently, Koga and Kawakami  investigated the $`S=1`$ model by the cluster-expansion method, and obtained $`J_\mathrm{c}^{}=0.056(1)`$ for $`\alpha =1`$. However, there have been only a very limited number of numerical works, in which both of the interchain and intrachain interactions are treated on an equal footing.  Such numerical analyses are certainly required, since the mean-field-like approximation is not necessarily appropriate even in the Q1D regime.

In the present paper, we report the results of quantum Monte Carlo (QMC) simulations by using the continuous-imaginary-time loop algorithm  on the $`S=1/2`$ and $`S=1`$ HAF model described by Eq. (1). The present paper is organized as follows. In Sec. II, the method of our numerical analyses is explained. In Secs. III and IV, the ground-state phase diagram parameterized by the strength of the bond alternation, $`\alpha `$, and that of the interchain coupling, $`J^{}`$, is determined precisely for $`S=1/2`$ and $`S=1`$, respectively. Especially, for the $`S=1`$ system with $`\alpha =1`$, we obtain $`J_\mathrm{c}^{}=0.043648(8)`$, which is consistent with $`0.040(5)`$ suggested by the recent QMC work, but is much more accurate. Furthermore, both in the $`S=1/2`$ and $`S=1`$ systems, the quantum phase transitions between the spin-gapped phases and the AF-LRO phase are confirmed to belong to the same universality class with that of the 3D classical Heisenberg model: the exponent of the correlation lengths is $`\nu =0.71(3)`$ for $`S=1/2`$ and $`\nu =0.70(1)`$ for $`S=1`$, which coincides fairly well with that of the latter model, $`\nu =0.7048(30)`$ We also show the results on the correlation length and the gap in the spin-gapped phase. In Sec. V, the topology of the phase diagram is discussed in detail based on the result of the present QMC calculation. We show that all the spin-gapped phases, such as the Haldane and dimer phases, are adiabatically connected with each other in the extended phase parameter space. This is in a sharp constrast to the strict 1D case, in which the spin-gapped phases are classified into different classes in terms of the so-called string-order parameter. It is of interest that the 1D spin-gap phases, which have different hidden symmetry, are connected without encountering any singularity in the 2D phase diagram. The final section is devoted to the concluding remarks.

## II Method

We consider the system descrived by the Hamiltonian (1) with $`S=1/2`$ and $`S=1`$. The real-space size is $`L_x\times L_y`$ and the inverse temperature, i.e., the imaginary-time size, is $`\beta =1/T`$. Periodic boundary conditions are imposed in the $`x`$\- and $`y`$-directions. We use the continuous-imaginary-time loop algorithm with multi-cluster update. Especially, for the $`S=1`$ system we adopt the subspin-representation technique, in which the $`S=1`$ system is represented by an $`S=1/2`$ system with special boundary conditions in the imaginary-time direction. By using these techniques, we can perform the simulatation up to $`L_x\times L_y=336\times 48`$ with $`\beta =100`$ for the $`S=1`$ case without encountering any difficulty.

The imaginary-time dynamical structure factor is defined by

$`S_\mathrm{d}(q_x,q_y,\omega )`$ (3)
$`={\displaystyle \frac{1}{L_xL_y\beta }}{\displaystyle \underset{i,j}{}}{\displaystyle _0^\beta }𝑑t𝑑t^{}\mathrm{e}^{i𝐪(𝐫_i𝐫_j)i\omega (tt^{})}S_i^z(t)S_j^z(t^{}),`$ (4)
where $`𝐪=(q_x,q_y)`$ is the wave-number vector, $`S_i^z(t)`$ is the $`z`$-component of the spin on site $`i`$ at imaginary time $`\tau `$, and $`\mathrm{}`$ is the thermal average. By using $`S_\mathrm{d}(q_x,q_y,\omega )`$, the staggered correlation length along the $`x`$-direction, $`\xi _x`$, is then evaluated by the second-moment method,

$$\xi _x=\frac{L_x}{2\pi }\sqrt{\frac{S_\mathrm{d}(\pi ,\pi ,0)}{S_\mathrm{d}(\pi +2\pi /L_x,\pi ,0)}1}.$$
(5)
The correlation length in the $`y`$-direction, $`\xi _y`$, and that in the imaginary-time ($`\tau `$) direction, $`\xi _\tau `$, which is related to the gap $`\mathrm{\Delta }`$ by $`\mathrm{\Delta }=1/\xi _\tau `$, are calculated similarly. Finally the staggered susceptibility, $`\chi _\mathrm{s}`$, is evaluated by

$`\chi _\mathrm{s}`$ $`=`$ $`S_\mathrm{d}(\pi ,\pi ,0)`$ (6)
$`=`$ $`{\displaystyle \frac{1}{L_xL_y\beta }}{\displaystyle \underset{i,j}{}}{\displaystyle _0^\beta }𝑑t𝑑t^{}e^{i\pi (𝐫_i𝐫_j)}S_i^z(t)S_j^z(t^{}).`$ (7)
All the structure factors are calculated by using the improved estimators. The period of $`10^2`$$`10^3`$ Monte Carlo steps (MCS) is used for thermalization and that of $`10^3`$$`10^5`$ MCS for the evaluation of physical quantities.

Near the critical point $`(\alpha _\mathrm{c},J_\mathrm{c}^{})`$ of the ground-state transition, the correlation lengths diverge as

$`\xi _x,\xi _yt^\nu `$ (8)
$`\xi _\tau t^{z\nu }=t^\nu ,`$ (9)
where $`t`$ is the distance from the critical point and $`\nu `$ is the critical exponent for the correlation length. Here we have put $`z=1`$ assuming the Lorenz invariance. Furthermore, the following finite-size-scaling (FSS) formula holds near $`(\alpha _\mathrm{c},J_\mathrm{c}^{})`$ and $`T=0`$ for systems with the fixed ratio $`L_x:L_y:\beta `$,

$$\xi _x/L_xf(tL_x^{1/\nu },L_x^zT)=f(tL_x^{1/\nu }),$$
(10)
and similar ones for $`\xi _y`$ and $`\xi _\tau `$, and

$$\chi _\mathrm{s}L_x^{\gamma /\nu }g(tL_x^{1/\nu },L_x^zT)=L_x^{\gamma /\nu }g(tL_x^{1/\nu }).$$
(11)
Here $`f`$ and $`g`$ are scaling functions and $`\gamma `$ the exponent for $`\chi _\mathrm{s}(t^\gamma )`$. Note that $`L_xT`$ is put constant in the above equations. We assume a polynomial up to the second order for the scaling functions. By using least-squares fitting, we obtain the critical point $`(\alpha _\mathrm{c},J_\mathrm{c}^{})`$ and the associated critical exponents $`\nu `$ and $`\gamma `$.

In addition, at some points in the spin-gapped phase, we explicitly evaluate the correlation lengths, $`\xi _x`$ and $`\xi _y`$, and the gap, $`\mathrm{\Delta }`$, at $`T=0`$ in the thermodynamic limit $`L_x,L_y\mathrm{}`$. For this purpose we extrapolate the simulated data first to the ground state $`T0`$ and then to the thermodynamic limit $`L_x,L_y\mathrm{}`$.

## III Results for $`𝐒=\mathrm{𝟏}/\mathrm{𝟐}`$

### A Ground-state phase diagram

In Fig. 2 we show the ground-state phase diagram of the $`S=1/2`$ system obtained by the FSS analysis explained in the previous section. As an example of the FSS analysis, we show in Fig. 3 that of $`\xi _\tau `$ against $`\alpha `$ for $`J^{}=1`$ (dotted line in Fig. 2). The aspect ratio of the $`(2+1)`$-dimensional system is taken as $`L_x:L_y:\beta =1:1:1`$. By the least-squares fitting, the exponent $`\nu `$ and the critical coupling $`\alpha _\mathrm{c}`$ are estimated as 0.71(1) and 0.31407(5), respectively. Here, the figure in parentheses denotes the statistical error ($`1\sigma `$) in the last digit. We also perform the same analyses on other lines in the $`\alpha `$-$`J^{}`$ plane, whose results are presented by the solid circles in Fig. 2. For example, on the line $`\alpha =J^{}`$ (dashed line in Fig. 2), we obtain $`\alpha _\mathrm{c}=J_\mathrm{c}^{}=0.52337(3)`$ and $`\nu =0.71(3)`$. In the phase diagram we can see that the ground state of the chain ($`J^{}=0`$) is the dimer state with a spin gap except for $`\alpha =1`$ and that the region of the AF-LRO phase enlarges monotonically as $`J^{}`$ increases.

Our phase diagram is qualitatively the same as that of Katoh and Imada (KI),  but not quantitatively. In particular, we obtain the critical point, $`\alpha _\mathrm{c}=0.31407(5)`$, on the line $`J^{}=1`$ (dotted line in Fig. 2), which is significantly smaller than their estimate $`\alpha _\mathrm{c}=0.398`$. More importantly, in the present simulation, the critical exponent $`\nu `$ on the transition points is evaluated as $`\nu =0.71(1)`$, which is consistent with $`\nu =0.7048(30)`$ for the 3D classical Heisenberg model.  The similar results have been obtained for the 2D $`1/5`$-depleted HAF model. On the other hand, KI concluded $`\nu =1`$. The reason of the these discrepancies might be due to the smallness of the system sizes and the inverse temperature used in the study by KI.

### B Correlation lengths and the gap

We also evaluate explicitly the ground-state correlation lengths and the gap on some points in the dimer phase by using the dynamic structure factors. Unless the points are very close to the critical line, these quantities in each systems with $`L`$ ($`=L_x=L_y`$) saturate to the ground-state values at temperatures we have simulated. For example, the $`T`$-dependences of these quantities are not to be discernible at $`T=0.05`$ and $`0.01`$ for $`\alpha =J^{}=0.4`$ and 0.5, respectively. On the other hand, the $`L`$-dependence still remains in sizes we have calculated. The $`L`$-dependences of the ground-state spatial correlation lengths and the gap are shown in Fig. 4 for $`\alpha =J^{}=0.5`$, which is close to the critical point $`\alpha _\mathrm{c}=J_\mathrm{c}^{}=0.52337(3)`$. Their values in the thermodynamic limit are estimated by fitting $`\xi _k(L)`$ to $`\xi _k(L)=\xi b\mathrm{exp}(cL)`$, where $`k=x`$, $`y`$, or $`\tau `$, $`\xi `$ is the value in the thermodynamic limit, and $`b`$ and $`c`$ are fitting parameters. As a result, we obtain $`\xi _x=3.0089(9)`$, $`\xi _y=2.2097(6)`$ and $`\mathrm{\Delta }=0.32261(4)`$ for $`\alpha =J^{}=0.4`$ and $`\xi _x=11.998(9)`$, $`\xi _y=9.312(10)`$, and $`\mathrm{\Delta }=0.0913(2)`$ for $`\alpha =J^{}=0.5`$. As $`\alpha `$ ($`=J^{}`$) becomes smaller, i.e., the system becomes more distant from the critical point, $`\mathrm{\Delta }`$ becomes larger and $`\xi `$ smaller.

## IV Results for $`𝐒=\mathrm{𝟏}`$

### A Overview on the phase diagram

Before going into the QMC analysis, let us here summarize the ground-state phase diagram of the $`S=1`$ system argued so far, which is shown in Fig. 5. For some points in the phase diagram the ground state is well understood by the previous theoretical and numerical studies: (1) $`(\alpha ,J^{})=(0,0)`$: The system consists of a set of the isolated antiferromagnetically-coupled spin pairs. The ground state is a trivial tensor product of dimer singlets sitting on each bond. (2) $`(\alpha ,J^{})=(1,0)`$: The system consists of isolated $`x`$-parallel Haldane chains. (3) $`(\alpha ,J^{})=(1,1)`$: The system is a uniform and isotropic 2D HAF. There exists an AF LRO in the ground state. (4) $`J^{}=\mathrm{}`$: The system consists of $`y`$-parallel Haldane chains. Note that in this limit the value of $`\alpha `$ becomes irrelevant (see also discussions in Sec. V).

In their analysis by the cluster expansion method, Koga and Kawakami derived three phases in the Q1D region, and they called the regions which includes point (1), (2) and (3) the dimer phase, the Haldane phase, and the AF-LRO phase, respectively. The region near point (4) is another Haldane phase. Therefore we call here the region which includes point (2) the Haldane I (H-I) phase and the one which includes the line (4) the Haldane II (H-II) phase. For the uniform systems with $`\alpha =1`$ the H-I and H-II phases are equivalent when we exchange the roles of the $`x`$-axis and the $`y`$-axis, and of $`J^{}`$ and $`1/J^{}`$.

### B Haldane-AF phase transition in the non-dimerized system

To demonstrate our FSS analyses, let us begin with critical behavior near $`(1,J_\mathrm{c}^{})`$, which separates the H-I and AF phases ($`J_\mathrm{c}^{}0.04`$ due to Ref. ). We sweep $`J^{}`$ near supposed $`J_\mathrm{c}^{}`$ with $`\alpha `$ fixed to unity. In Fig. 6, $`\xi _x/L_x`$, $`\xi _y/L_y`$, and $`\xi _\tau /\beta `$ of the systems with $`L_x=L_y=\beta L`$ ($`=16`$, 24,$`\mathrm{}`$, 64) are plotted. As one sees immediately, the data suffer from quite large corrections to scaling, i.e., the crossing point of the scaled correlation lengths with two different $`L`$’s clearly shifts to larger $`J^{}`$ as the system size increases. We attribute these large corrections to the strong spatial anisotropy in the coupling constants ($`J^{}1`$). Indeed, the value of $`\xi _y`$ at $`J^{}=0.0435`$ is $`9.29(7)`$ (see Fig. 6 (b)) even for the largest system size $`L=64`$, which is quite smaller than those in the other directions ($`\xi _x=60.4(3)`$ and $`\xi _\tau =24.7(2)`$). Among the three correlation lengths, $`\xi _x`$ is the largest: the growth of the correlation is dominated only by the system size in the $`x`$-direction. This indicates that we need larger lattices, especially in the $`x`$-direction, in order to perform a precise FSS analysis.

To simulate larger lattices with minimal costs, we therefore optimize the aspect ratio, $`L_x:L_y:\beta `$, as explained below. Expecting that the scaled correlation lengths, $`\xi _x/L_x`$, $`\xi _y/L_y`$, and $`\xi _\tau /\beta `$, become nearly equal with each other at the critical point, we set the aspect ratio $`L_x:L_y:\beta `$ as $`7:1:2`$ based on the data presented in Fig. 6. With this ratio we simulate systems with $`L_x=168`$, 224, 280, and 336 and perform the FSS analyses. The raw data of $`\xi _x/L_x`$, $`\xi _y/L_y`$, and $`\xi _\tau /\beta `$ with this aspect ratio are shown in Fig. 7. Now $`\xi _y=12.02(5)`$ even for the smallest system size ($`L_x=168`$ and $`L_y=24`$) at $`J^{}=0.0435`$. The ratios $`\xi _x/L_x`$, $`\xi _y/L_y`$ and $`\xi _\tau /\beta `$ in a common range of $`J^{}`$ become nearly equal, and corrections to scaling become much smaller than in the previous ones as we expected. The FSS plot for $`\xi _x`$ is shown in Fig. 8. The resultant $`J_\mathrm{c}^{}`$ and $`\nu `$ are as follows: $`(J_\mathrm{c}^{},\nu )=(0.043648(9),0.69(1))`$ from $`\xi _x`$, $`(0.043649(8),0.71(1))`$ from $`\xi _y`$, and $`(0.043648(7),0.69(1))`$ from $`\xi _\tau `$.

Averaging these three values we conclude with

$`J_\mathrm{c}^{}`$ $`=`$ $`0.043648(8)`$ (12)
and

$`\nu `$ $`=`$ $`0.70(1).`$ (13)
Fixing the value of $`J_\mathrm{c}^{}`$ and $`\nu `$ thus determined, we next perform the FSS analysis on the staggered susceptibility $`\chi _\mathrm{s}`$. The raw data of $`\chi _\mathrm{s}`$ vs. $`J^{}`$ are shown in Fig. 9 and the fitting result is shown in Fig. 10. The latter yields

$$\gamma =1.373(3).$$
(14)
The exponents $`\nu `$ and $`\gamma `$ we have obtained for the $`S=1`$ system again agree with those of the 3D classical Heisenberg model.

The critical point obtained just above is consistent with the previous result by the method involving the mean-field approximation, $`J_\mathrm{c}^{}0.025`$, and also with that by the recent QMC method, $`J_\mathrm{c}^{}=0.040(5)`$. On the other hand, the present result significantly differs from that obtained by the cluster-expansion method, $`J_\mathrm{c}^{}=0.56(1)`$. The present FSS analyses on our extensive QMC results make us possible to obtain the critical value with much higher accuracy than the other methods.

### C Ground-state phase diagram

In a similar way as described in the previous subsection, we obtain other critical points on the $`\alpha `$-$`J^{}`$ phase diagram as shown in Fig. 11. First, the scaled correlation lengths, $`\xi _x/L_x`$, $`\xi _y/L_y`$, and $`\xi _\tau /\beta `$, are calculated up to $`L_x=64`$ with either $`\alpha `$ or $`J^{}`$ fixed. Sweeping $`J^{}`$ or $`\alpha `$ with sufficiently high resolution, we regard a crossing point of these $`\xi `$’s as the critical point. Note that the optimal aspect ratio depends strongly on the value of $`\alpha `$ and $`J^{}`$. However, we adopt $`L_x:L_y:\beta =1:1:1`$ for simplicity. Although the results thus obtained suffer from relatively larger systematic corrections than those presented in the last subsection, the absolute magnitude of the systematic error in the estimates should be still smaller enough than the symbol size in Fig. 11.

For some critical points the FSS analysis as in the previous subsection is also carried out. We obtain the exponents $`\nu `$ and $`\gamma `$ which are consistent with Eqs. (13) and (14), respectively. This supports that the quantum critical phenomena in the $`S=1`$ system also belong to the same universality class as that of the 3D classical Heisenberg model. An exception is the 1D critical point located at $`(\alpha _\mathrm{c},J_\mathrm{c}^{})=(0.5879(6),0),`$ which separates the dimer phase from the Haldane phase. The apparent value of $`\nu `$ starts to deviate from (13) when $`\alpha `$ becomes closer to $`\alpha _\mathrm{c}`$. This is attributed to the crossover to the critical phenomena belonging to the Gaussian universality class. We confirm that the AF-LRO phase exists between the two spin-gap phases at least down to $`J^{}=0.01`$ at $`\alpha =\alpha _\mathrm{c}`$. Although in the present simulation it is quite difficult to prove the existence of the AF-LRO phase at smaller $`J^{}`$, we believe that the the point $`(\alpha _\mathrm{c},0)`$ is tricritical: the 1D critical point is unstable against an infinitesimal interchain coupling and the AF LRO immediately appears as the same as in the $`S=1/2`$ uniform chain.

Interestingly, the H-II and D phases are adiabatically connected with each other. The gapless AF-LRO phase does not touch the line of $`\alpha =0`$ as seen in Fig. 12, where the part of the whole phase diagram (Fig. 11) near $`\alpha =0`$ is magnified. Indeed, on the $`\alpha =0`$ line, which corresponds to the $`S=1`$ two-leg ladder, it is shown that there exists no critical point by the recent QMC study. Thus, the D phase can be identified with the H-II phase, and there are only two distinct spin-gap phases, H-I and H-II in Fig. 11. The closeness of the critical line to the $`\alpha =0`$ line is due to the strong AF fluctuations, which already exist in the two-leg ladder system.

### D Correlation lengths and the gap

We obtain the explicit values of $`\xi _x`$, $`\xi _y`$, and $`\mathrm{\Delta }`$ in the ground state at $`(\alpha ,J^{})=(1,0.04)`$. They are calculated for systems with sizes $`L_x=168`$, 224, 280, and 336 and with the aspect ratio $`L_x:L_y=7:1`$ at $`T`$ regarded as zero temperature. Their $`T`$-dependences are negligible at $`T=0.01`$. We extrapolate the finite-size data to the thermodynamic limit in the same way as explained for the $`S=1/2`$ system. We obtain $`\xi _x=39.2(1)`$, $`\xi _y=5.67(1)`$, and $`\mathrm{\Delta }=0.0632(2)`$. As $`J^{}`$ becomes smaller, $`\xi _x`$ becomes smaller and $`\mathrm{\Delta }`$ larger to reach at $`J^{}=0`$ the single chain values $`\xi _x=6.0153(3)`$ and $`\mathrm{\Delta }=0.41048(6)`$.

## V Discussions

The analyses presented in the preceding section revealed that the ground-state phase diagram of the $`S=1`$ system has a rather complicated topology, i.e., the H-II and D phases are adiabatically connected with each other, though the channel between them is quite narrow (Fig. 11). On the other hand, as for the H-I and D phases, in the 1D system ($`J^{}=0`$) these two spin-gapped phases are distinctively separated by the critical point at $`\alpha _\mathrm{c}=0.5879(6)`$, and they are distinguished by the string-order parameter, which is zero in the former phase and finite in the latter one. The transition can be viewed as a rearrangement of dimer-singlet pattern between the (1,1)- and (2,0)-valence-bond-solid (VBS) states. We emphasize that once $`J^{}`$ is introduced, however, the string-order parameter should vanish even in the H-I phase, being similar to the $`S=1`$ ladder. Still one may consider that the two phases essentially differ with each other since they are separated by the AF-LRO phase. If, however, we introduce the bond alternation also in the $`y`$-direction, the two phases can be connected without passing the gapless state as explained below.

Let us consider the 2D HAF model defined in the extended parameter space:

$``$ $`=`$ $`J_x\left\{{\displaystyle \underset{i,j}{}}𝐒_{2i,j}𝐒_{2i+1,j}+\alpha _x{\displaystyle \underset{i,j}{}}𝐒_{2i+1,j}𝐒_{2i+2,j}\right\}`$ (15)
$`+`$ $`J_y\left\{{\displaystyle \underset{i,j}{}}𝐒_{i,2j}𝐒_{i,2j+1}+\alpha _y{\displaystyle \underset{i,j}{}}𝐒_{i,2j+1}𝐒_{i,2j+2}\right\}.`$ (16)
The original Hamiltonian (1) corresponds to the case with $`J_x=1`$, $`\alpha _x=\alpha `$, $`J_y=J^{}`$, and $`\alpha _y=1`$. The Hamiltonian (15) is invariant under the exchange between $`(J_x,\alpha _x)`$ and $`(J_y,\alpha _y)`$.

To draw the phase diagram in this extended parameter space, it is convenient to introduce a parameter, $`RJ_y/(J_x+J_y)`$. The limits $`J_y0`$ and $`J_y\mathrm{}`$ correspond to $`R=0`$ and 1, respectively. Since $`\alpha _x`$ ($`\alpha _y`$) becomes irrelevant in the limit $`J_x0`$ ($`J_x\mathrm{}`$), the whole phase diagram in the three-dimensional parameter space is shaped as a tetrahedron. In Fig. 13, we present the ground-state phase diagram parametrized by $`R`$, $`\alpha _x(1R)`$, and $`\alpha _yR`$ for $`S=1/2`$ (a) and $`S=1`$ (b). It should be noted that the phase diagram should be invariant under the transformation, $`(R,\alpha _x(1R),\alpha _yR)(1R,\alpha _yR,\alpha _x(1R))`$, reflecting the symmetry in the Hamiltonian explained above. In the phase diagram, the edge AB (CD) corresponds to isolated $`x`$-parallel ($`y`$-parallel) decoupled chains, the edge AD isolated four-spin plaquettes, and the edge AC (BD) the two-leg ladders in $`y`$-direction ($`x`$-direction).

The face ABC (and also CDB) in Fig. 13 corresponds to the original phase diagram shown in Figs. 2 and 11, though the $`y`$-parallel-chain limit $`J^{}\mathrm{}`$ in the original diagram is represented by one vertex C (B) in the new ones. In the extended phase diagram the shaded (unshaded) area represents the AF-LRO (spin-gapped) phase on the ABC- and CDB-faces.

It should be emphasized that on the ACD- and DBA-faces there is no AF-LRO phase, since the system is one-dimensional dimerized two-leg ladder. There exist only the 1D critical points discussed already. Especially, there is no critical point on the edge AD. Therefore, in the $`S=1`$ case, the three spin-gapped phases, H-I, D, and H-II, are connected by the path C$``$A$``$D$``$B. Similarly, in the $`S=1/2`$ case, the two dimer phases, which correspond to the vertex A and D, respectively, are connected directly by the path A$``$D. Thus, in both cases, there are only two phases, namely, the spin-gapped phase and the AF-LRO one.

## VI Concluding Remarks

In this paper, we have investigated the ground-state phase diagram of $`S=1/2`$ and $`S=1`$ HAF on the anisotropic dimerized square lattice by means of the extensive QMC simulation with the continuous-imaginary-time loop algorithm and the FSS analyses. It is confirmed that, for both $`S=`$1/2 and 1, the quantum critical phenomena in the model belong to the same universality class as that of the 3D classical Heisenberg model, except for the 1D critical points, which belong to the same universality class as that of the Gaussian model. We have also demonstrated that the spin-gapped phases of the 1D chain are connected when we introduce the interchain couplings with bond alternation. In the 2D system, only one spin-gapped phase exists in both of the $`S=1/2`$ and $`S=1`$ systems.

The results obtained in the present work are considered to be the proper basis for investigation of peculiar phenomena observed in the Q1D HAF’s mentioned at the beginning of this paper. We have already reported the QMC analysis on the site-dilution-induced AF LRO in these materials based on the Hamiltonian (1).  In order to discuss various experimental results quantitatively, it is certainly necessary to take into account other ingredients than in Eq. (1), such as the next-nearest-neighbor intrachain interaction and the single-ion anisotropy. They are beyond the scope of the present work. The present results, however, demonstrate the role of the higher dimensionality which has been overlooked so far.

Most of the numerical calculations in the present work have been performed on the DEC Alpha, SGI ORIGIN 2000, SGI 2800, and RANDOM at the Materials Design and Characterization Laboratory, Institute for Solid State Physics, University of Tokyo and on the Hitachi SR-2201 at the Supercomputer Center, University of Tokyo. The program used in the present simulation was based on the library ‘Looper version 2’ developed by S.T. and K. Kato and also on the ‘PARAPACK version 2’ by S.T. The present work was supported by the “Research for the Future Program” (JSPS-RFTF97P01103) of the Japan Society for the Promotion of Science. S.T’s work was partly supported by the Swiss National Science Foundation.