text:
# Critical behavior of a one-dimensional fixed-energy stochastic sandpile

## I Introduction

Sandpile models are the prime example of self-organized criticality (SOC), or scale-invariance in the apparent absence of tuning parameters . SOC in a slowly driven sandpile has an analog in the absorbing-state phase transition of the corresponding nondriven or fixed-energy sandpile (FES) . While most studies of sandpiles have probed the slow-driving limit (addition and loss of sand grains at an infinitesimal rate), there is great interest in understanding the scaling properties of FES models as well . One may hope that these systems, which are isolated and translation-invariant, will prove easier to characterize than driven sandpiles, which have eluded a rather massive theoretical and computational effort to obtain precise critical exponents or assign universality classes. In this paper we present extensive numerical results for a particularly simple one-dimensional FES. Analogous studies of higher dimensional systems can be found in  (A review of these results and related models can be found in Refs. ).

One of the central features of sandpile models is the presence of a conserved field, which may be thought of as the density of sand grains, or, in our nomenclature, the energy density. This field is conserved locally as well as globally, and couples to the activity density, which is the order parameter. The configuration evolves through a series of “toppling” events, i.e., redistributions of energy from an active site to its neighbors, which may be either deterministic or stochastic. The well known Bak-Tang-Wiesenfeld (BTW) sandpile has a deterministic toppling rule, and is also Abelian (the final absorbing configuration is independent of the order of toppling), which allows many stationary-state properties of the driven sandpile to be found exactly . A less desirable aspect of the deterministic dynamics, on the other hand, is that in the steady-state only a small subset of the possible configurations (determined by the initial state) are visited . This leads to strong nonergodic effects in the FES version of the BTW automaton . Here we study a one-dimensional stochastic FES sandpile to elucidate the consequences of stochastic rules.

In our model, a variant of the Manna sandpile , the pair of particles liberated when a site topples move independently, and with equal probability, to the nearest neighbors on the left or right. We may therefore think of the particles as independent random walkers. (There is no restriction on the number of walkers at a given site.) While the walkers are independent as far as their hopping direction is concerned, their mobility does involve an interaction: isolated particles cannot move, but if two or more particles occupy the same site, the site may topple. Thus for densities $`\zeta =N/L1`$ ($`N`$ is the number of walkers on a ring of $`L`$ sites), absorbing configurations exist, in which all walkers are paralyzed. But since the fraction of absorbing configurations vanishes as $`\zeta 1`$, it is reasonable to expect a phase transition from an absorbing to an active phase at some $`\zeta _c<1`$. Simulations bear this out, and show that there is a continuous transition at $`\zeta _c0.9488`$.

Our goal in the present paper is to determine the critical behavior of the model. To this end, we analyze the static behavior and the approach to the steady state in the language of absorbing-state phase transitions, studying the order parameter, its fluctuations, and temporal correlations. We also investigate the cumulative local activity, whose dynamics obeys a linear diffusion equation in the presence of an effective quenched noise term that reflects the stochastic redistribution rules, the Linear Interface Model (LIM) . Our main finding is that the one-dimensional Manna model defines a universality class different both from that of directed percolation, and of the linear interface model. The existence of long-range correlations that arise from a conservation law is believed to be at the basis of this new universality class .

The paper is structured as follows. In Sec. II we define the model and our simulation algorithm. Numerical results are analyzed, in the contexts of absorbing-state phase transitions and of driven interfaces, in Sec. III. In Sec. IV we summarize and discuss our findings.

## II Model

The model is defined on a one-dimensional lattice with periodic boundaries. The configuration is specified by the energy (or number of walkers) $`z_i=0,1,2,\mathrm{}`$ at each site; sites with $`z_i2`$ are said to be active. A Markovian dynamics is defined by the toppling rate, which is unity for all active sites, and zero for sites with $`z_i<2`$. When a site $`i`$ topples, $`z_iz_i2`$ and the two particles liberated move independently to randomly chosen nearest neighbors $`j`$ and $`j^{}`$ ($`j`$, $`j^{}\{i+1,i1\}`$). (Thus $`j=j^{}`$ with probability $`1/2`$.) In most cases we use random sequential dynamics: the next site to topple is chosen at random from a list of active sites, which must be updated following each toppling event. The time increment associated with each toppling is $`\mathrm{\Delta }t=1/N_A`$, where $`N_A`$ is the number of active sites just prior to the toppling. In this way, $`N_A`$ sites on average topple per unit time, just as they would in a simultaneously updated version of the model, in which all active sites at time $`t_n`$ are toppled simultaneously and $`\mathrm{\Delta }t1`$. We expect the two dynamics to be equivalent insofar as scaling properties are concerned, and note that the latter choice was used in some of the interface representation studies discussed below.

In most of our simulations, the initial condition is generated by distributing $`\zeta L`$ particles randomly among the $`L`$ sites, yielding an initial (product) distribution that is spatially homogeneous and uncorrelated. Once the particles have been placed, the dynamics begins (we verified that allowing some toppling events during the insertion phase has no effect on the stationary properties). The condition of having at least one active site in the initial configuration is trivially satisfied on large lattices, for the $`\zeta `$ values of interest, i.e., close to the critical value. In fact, for large $`L`$, the initial height at a given site is essentially a Poisson random variable, and the probability of having no active sites decreases exponentially with the lattice size. It is worth remarking that while the initial conditions are statistically homogeneous, the energy density is not perfectly smooth. For $`1lL`$, the energy density on a set of $`l^d`$ sites is essentially a Gaussian random variable with mean $`\zeta `$ and variance $`l^d`$. Such initial fluctuations relax on a time scale $`l^{z_D}`$ with $`z_D2`$. That is, the exponent governing diffusive relaxation must be at least as large as for independent random walkers. (Recall that in our model regions of lower density will exhibit less activity, hence a slower redistribution of particles.) The initial value of the critical-site density, $`\rho _c`$ (the density of sites with $`z_i=1`$), moreover, is generally far from its stationary value, complicating relaxation to the steady state. These observations are pertinent to the studies of modified initial conditions reported in Sec. IIId.

If after some time the system falls into a configuration with no active sites, the dynamics is permanently frozen, i.e., the system has reached an absorbing configuration. We shall see that as we vary $`\zeta `$, the model exhibits a phase transition separating an absorbing phase (in which the system always falls into an absorbing configuration), from an active phase possessing sustained activity.

## III Simulation Results

### A Absorbing-state phase transition

The presence of an absorbing-state phase transition is a general feature of FES models. In analogy with standard analyses of such transitions, it is possible to derive a set of Langevin equations describing their critical behavior in high dimensions. Following Ref., we can write the dynamical equations for the density of active sites $`\rho _a`$ and the energy density $`\zeta `$ as

$$_t\rho _a=D^2\rho _a+\mu \rho _ab\rho _a^2+w\zeta \rho _a+\rho _a^{1/2}\eta $$
(1)
$$_t\zeta =\lambda ^2\rho _a,$$
(2)
where $`\mu ,b,w`$ and $`\lambda `$ are coupling constants and $`\eta `$ is a Gaussian noise. Eqs. (1-2) closely resemble the Langevin field theory of directed percolation (DP), apart from an additional coupling between the activity and the energy field, which is conserved by the dynamics. Due to this coupling, the FES critical behavior is not in the DP universality class, but belongs to a new universality class embracing all systems with multiple absorbing states in which the order parameter is coupled to a static conserved field . (By ‘static’ we mean that this field is frozen in regions where the order parameter is null.) The field theory for a reaction-diffusion model in this universality class has been derived exactly; it reproduces in some limits Eqs. (1-2) . It is thus interesting to study the critical behavior in $`d=1`$, to see if it is possible to generalize the known results for fixed-energy sandpiles  to one-dimensional systems.

We simulated the model using system sizes ranging from $`L=100`$ to about $`10^4`$ sites. (Since $`\zeta =N/L`$ with $`N`$ an integer, we are obliged to use different sets of $`L`$ values to study different energy densities $`\zeta `$.) In stationary-state simulations, we collect data over an interval of $`t_m`$ time units, following a relaxation period of $`t_r`$. For small system sizes, $`t_m`$ and $`t_r`$ are of the order of 103, but for our largest systems ($`L10^4`$) we used $`t_r5\times 10^6`$ and $`t_m=2.5\times 10^6`$. We verified that our results show no systematic variation with time for $`t>t_r`$. A run consists of $`N_s`$ independent trials, each with a different initial configuration, with $`N_s`$ ranging from $`2\times 10^5`$ for $`L=100`$, to 500 or 1000 for $`L10^4`$. In practice $`t_m`$ is limited because for $`\zeta \zeta _c`$, the survival probability decays sensibly over this time scale; in some cases only about 25% of the trials survive to time $`t_r+t_m`$.

Fig. 1 shows the overall dependence of the stationary active-site density as a function of $`\zeta `$; the points represent extrapolations of results for $`L=100`$ \- $`5000`$ to the $`L\mathrm{}`$ limit. The data indicate a continuous transition from an absorbing state ($`\rho _a`$ = 0) to an active one at $`\zeta _c`$ in the vicinity of 0.95.

Our first task is to locate the critical value $`\zeta _c`$. We studied the stationary active-site density, $`\rho _a`$, and its second moment $`\rho _a^2`$, anticipating that as in other absorbing-state phase transitions, the active-site density (i.e., the order parameter) will obey finite-size scaling :

$$\overline{\rho _a}(\mathrm{\Delta },L)=L^{\beta /\nu _{}}(L^{1/\nu _{}}\mathrm{\Delta }),$$
(3)
where $`\mathrm{\Delta }\zeta \zeta _c`$, and $``$ is a scaling function with $`(x)x^\beta `$ for large $`x`$, since for $`L\xi \mathrm{\Delta }^\nu _{}`$ we expect $`\overline{\rho _a}\mathrm{\Delta }^\beta `$ (here $`\xi `$ is the correlation length). When $`\mathrm{\Delta }=0`$ we have that $`\overline{\rho _a}(0,L)L^{\beta /\nu _{}}`$. For $`\mathrm{\Delta }>0`$, by contrast, $`\overline{\rho _a}`$ approaches a stationary value, while for $`\mathrm{\Delta }<0`$ it falls off as $`L^d`$. Thus in a double-logarithmic plot of $`\overline{\rho _a}`$ versus $`L`$ (see Fig. 2), supercritical values ($`\mathrm{\Delta }>0`$) are characterized by an upward curvature, while for $`\mathrm{\Delta }<0`$ the graph curves downward. Using this criterion (specifically, zero curvature in the data for $`L1000`$), we find $`\zeta _c=0.94887(7)`$, with the uncertainty reflecting the scatter in our numerical results for the curvature (see Fig. 2, inset). The associated exponent ratio is $`\beta /\nu _{}=0.235(11)`$. A similar analysis of the data for $`\rho _a^2`$ yields $`\zeta _c=0.94883(5)`$ with an exponent of $`2\beta /\nu _{}=0.483(18)`$. We therefore adopt the estimates $`\zeta _c=0.94885(7)`$ and $`\beta /\nu _{}=0.239(11)`$.

We measured the autocorrelation function for the number of active sites $`N_a`$,

$$C(t)=\frac{N_a(t_0+t)N_a(t_0)N_a^2}{N_a^2N_a^2}$$
(4)
in the stationary state. To obtain clean results for $`C(t)`$ we study surviving trials in relatively long runs (from $`t_m=2\times 10^5`$ for $`L=625`$, to $`5\times 10^6`$ for $`L=10^4`$; this obliges us to reduce our sample to 200 surviving trials for $`L2500`$ and 100 surviving trials for $`L5000`$). Results for $`\zeta =0.9488`$ are shown in Fig. 3: C(t) decreases monotonically, but does not follow a simple exponential decay. To study the dependence of the relaxation time on system size, we determine the temporal rescaling factor $`r`$ required to obtain a data collapse between $`C(t;L/2)`$ and $`C(t/r;L)`$. A good collapse is possible (see Fig. 3), but the rescaling factor depends on $`L`$; for $`L=2^n625`$, we use $`t^{}=t/r^n`$ with $`r=2.93`$, 2.91, and 2.89 for $`n=1`$, 2, and 3, respectively. The rescaling factor $`r`$ appears to approach a limiting ($`L\mathrm{}`$) value of 2.80(5), corresponding to a relaxation time that scales as $`\tau L^{\nu _{||}/\nu _{}}`$ with $`\nu _{||}/\nu _{}z=\mathrm{ln}r/\mathrm{ln}2=1.49(3)`$. We also studied the half-life $`\tau _H`$ (the time for the survival probability to decrease by a factor of 1/2 in the stationary state) at $`\zeta _c`$. This relaxation time exhibits a clean power-law dependence on system size, $`\tau _HL^z`$, with $`z=1.42(1)`$.

Next we examine the stationary scaling of the order parameter away from the critical point. We determined the stationary active-site density $`\overline{\rho _a}(\zeta ,L)`$, for $`\zeta `$ in the vicinity of $`\zeta _c=0.94885`$, and for system sizes $`L=100`$ \- 5000. We analyze these data using the finite-size scaling form of Eq. (3), which implies that a plot of $`L^{\beta /\nu _{}}\overline{\rho _a}(\mathrm{\Delta },L)`$ versus $`L^{1/\nu _{}}\mathrm{\Delta }`$ should exhibit a data collapse. We shift each data set (in a log-log plot of $`L^{\beta /\nu _{}}\rho _a`$ versus $`L^{1/\nu _{}}\mathrm{\Delta }`$), vertically by $`(\beta /\nu _{})\mathrm{ln}L`$, using $`\beta /\nu _{}=0.239`$ as found above, and determine the horizontal shift $`S(L)`$ required for data collapse. These are found to follow the relation $`S(L)=\nu _{}^1\mathrm{ln}L`$, with $`\nu _{}^1=0.553(3)`$. That these values yield an excellent data collapse is evident from Fig. 4. The slope of the scaling plot (linear regression using the 25 points with $`\mathrm{ln}(L^{1/\nu _{}}\mathrm{\Delta })>0.5`$), yields $`\beta =0.410(4)`$. This is somewhat smaller than, but consistent with, the estimate $`\beta =0.43(2)`$ obtained by combining $`\nu _{}^1=0.553(3)`$ and $`\beta /\nu _{}=0.239(11)`$. We adopt $`\beta =0.42(2)`$ as our final estimate.

### B Interface Representation

Interface representations are useful for studying both sandpiles  and particle systems such as the contact process . Using such a representation, the activity history of the sandpile can be described in terms of an equation of motion that is an exact analog of a driven interface advancing in a quenched disordered medium. Our one-dimensional case is an example of a larger class of models, for which a mapping can be provided from the sandpile to an interface field . Such a description can be useful, for example, in understanding the role of microscopic rules and in defining novel quantities for study.

The mapping is constructed as follows. Let the interface height variable $`H_i(t)`$ count the number of topplings at site $`i`$ up to time $`t`$. The dynamics of $`H_i(t)`$ follows from the observation that the toppling of a site depends (through its height $`z_i`$) on the number of particles it has received, and on the number of times it has toppled. Each time a neighbor of site $`i`$ topples, site $`i`$ gains one, two, or no particles, with probabilities of 1/2, 1/4 and 1/4, respectively. This allows a ‘projection trick’, in which the dynamics of $`H_i(t)`$ is constructed from particle counting . Topplings at the sites neighboring $`i`$ give rise to an average incoming flux $`\overline{n}_i^{in}`$ at site $`i`$, plus a fluctuating part $`\delta n_i^{in}`$, which records, in an exact manner, the random choices made for each particle. The average flux minus the number of topplings at $`i`$ gives rise to a discrete Laplacian, while the fluctuating part can be mapped to a noise term $`\tau (i,H)`$, which can be evaluated for each toppling. Following these prescriptions, we find a discrete interface equation with quenched noise, which reproduces the sandpile activity :

$$\frac{H_i}{t}=\mathrm{\Theta }[^2H_i+F_i1+\tau (i,H)]\mathrm{\Theta }(f_i)$$
(5)
where $`\mathrm{\Theta }`$ is the step function, so that $`H_i`$ is nondecreasing. Since $`H_i(t)`$ is restricted to integer values the time derivative should be interpreted as the rate of transitions from $`H`$ to $`H+1`$. The ‘force’ $`f_i`$ at site $`i`$ consists of the following contributions. The (discrete) Laplacian term arises, as noted above, from the balance between the average number of particles site $`i`$ gains from, and loses to, its neighbors in toppling events. The ‘columnar’ noise term $`F_i`$ represents the number of particles initiallly at site $`i`$, while $`\tau (i,H)`$ is the fluctuation part of the particle flux associated with the $`H`$-th toppling of site $`i`$. The step-function character of the interface equation is discussed in greater detail in Ref. .

The noise-term $`\tau (i,H)`$ is a quenched variable that arises from a mapping of the ‘annealed’ disorder of the original sandpile. Note that one could also consider the random choices of the particles’ destinations at each toppling as a quenched, random landscape that is chosen in advance for each sandpile configuration. It is then easier to see that this maps into a quenched disorder term for the associated interface equation. This noise term is particle-conserving, and it has at each $`i`$ a random-walk-like correlator $`[\tau (H+\mathrm{\Delta }H,i)\tau (H,i)]^2\mathrm{\Delta }H`$. If these correlations could be neglected, the scaling of the interface would follow that of the linear interface model in a quenched random field, as appears to be the case in two- and three-dimensional systems . One should note, however, that long range-correlations in the noise can change the universality class, as has been shown numerically for another interface model . The possible LIM universality classes are still under debate .

In the interface description the model undergoes a depinning transition at the critical force $`F_c\zeta _c`$ . At this point power-law correlations develop in the history of topplings. The interface behavior, assuming simple scaling, is described by two exponents: the roughness exponent $`\alpha `$, and the early-time exponent $`\beta _W`$. Introducing the width $`W`$ as usual,

$$W^2(t,L)=[H_i(t)\overline{H}(t)]^2,$$
(6)
(here $`\overline{H}(t)`$ is the mean height), these exponents are defined via

$$W^2(t,L)\{\begin{array}{cc}t^{2\beta _W}\hfill & tt_\times \hfill \\ L^{2\alpha }\hfill & tt_\times \hfill \end{array},$$
(7)
where the crossover time $`t_\times L^z`$. Assuming that simple scaling holds and that the correlations in the interface can be described by a single length scale, we have the exponent relation $`\beta _Wz=\alpha `$.

This scaling picture, familiar from the study of surface growth, was recently shown to apply in the case of a simple absorning-state phase transition . For the one-dimensional Manna model the situation is complicated by several factors. First, the noise appearing in the interface description has two components, $`\tau (i,H)`$ and $`F_i`$. The interface behavior will therefore exhibit a crossover from a regime dominated by the initial configuration (reflected in $`F_i`$) to a randomness-dominated regime. This effect also appears in higher dimensions, but in $`d=1`$, due to the meager phase space, relaxation is much slower and transient effects may be much more severe. Note that the $`H`$-independent terms on the r.h.s. of the interface equation can be understood as an initial height profile, which then relaxes to the asymptotic rough state . It is easy to see that by varying such an initial condition different transients can be generated.

A problem with one-dimensional models is that the two-point correlation function of the surface roughness scales with a different exponent, $`\alpha _{loc}`$, which for fundamental reasons is limited to $`\alpha _{loc}1`$. The exponent $`\alpha `$ can attain a larger value, for example $`\alpha =1.25`$ for the one-dimensional linear interface model . These exponents are related via $`\alpha =\alpha _{loc}+\kappa `$, where $`\kappa `$ measures the divergence of the height difference between neighboring sites with $`L`$. In the corresponding 1-d linear interface model, this is termed ‘anomalous scaling’ . This means that as $`t\mathrm{}`$, the typical height difference between neighboring sites increases without limit. Since larger systems have a greater lifetime, this has implications for the roughness as measured by $`W_{sat}^2lim_{tt_c}W^2(t,L)`$, with $`t_c`$ the time at which the absorbing state with no activity is reached, i.e., the lifetime. The saturation width $`W_{sat}`$ scales as $`L^\alpha `$, with $`\alpha `$ related to $`\beta _W`$ and $`z`$ as above. In models exhibiting an absorbing state, such as the contact process or a FES, the width saturates only because all activity eventually ceases; the width in surviving trials does not saturate. This is in marked contrast to the behavior of interface models, in which the width saturates due to the Laplacian term representing surface tension. (Continuum descriptions of absorbing-state phase transitions, and their associated interface representations, likewisw include a Laplacian term, but the noise driving the interface grows without limit, so saturation is not required while there is activity .)

Finally, in absorbing-state models, interface scaling appears to be strongly linked to the approach to the stationary state. In a model with simple scaling (i.e., unique diverging length and time scales, and no conserved quantities), the growth exponent $`\beta _W`$ is related to the critical exponent $`\theta `$ governing the initial decay of activity via $`\beta _W+\theta =1`$ . In the present case relaxation is complicated by effects that may mimic (for a certain time) columnar disorder, and slowly relaxing perturbations (the long-wavelength density fluctuations mentioned in Sec. II).

In simulations we first studied the time-dependent mean interface width $`W^2(t,L)`$ in systems of size $`L=1253`$, 2506, 5012, 10024, and 20048, at the critical point $`\zeta =0.94892`$. The dependence of the saturation width on system system size yields $`\alpha =1.42(1)`$. We may then attempt to collapse the data for $`W^2(t,L)`$ using this exponent, and varying $`z`$ to obtain the best collapse; in this way we find $`z=1.65(2)`$. The resulting scaling plot (Fig. 5), of $`\stackrel{~}{W}^2W^2/L^{2\alpha }`$ versus $`\stackrel{~}{t}t/L^z`$, shows a good collapse, and an apparent power-law growth in the roughness, following an initial transient. From the scaling relation $`\beta _W=\alpha /z`$ we obtain $`\beta _W=0.863(13)`$, while a direct fit to the time-dependent width data yields $`\beta _W=0.87(2)`$.

For comparison, an independent series of studies at $`\zeta =0.9490`$ were performed to determine the lifetime $`t_c`$ and saturation width $`W_c`$ for $`L=400`$ to $`L=6400`$. Power-law fits to these data yield essentially consistent results, i.e., $`\alpha =1.48(2)`$, $`\beta _W=0.86(2)`$ and $`z=1.70(3)`$. Fig. 6 shows a clear power-law dependence of the saturation width on the lifetime in individual runs. This indicates that a normal Family-Vicsek-style scaling plot of $`W(t,L)`$ as in Fig. 5 is problematic, since as noted above the width in any surviving run does not saturate as would be the case in normal interface models, and the saturation width $`Ww_c`$ will therefore depend on how long the particular configuration stays alive.

Similarly to the case of the linear interface model, we find that there is an independent local roughness exponent $`\alpha _{loc}`$ that describes the two-point $`k`$-th order height-height correlation function $`G_k(r)=|H_{i+r}H_i|^kr^{k\alpha _{loc}}`$ for $`r<\xi (t)t^{1/z}`$. We find $`\alpha _{loc}=0.59(3)`$ for $`k=1\mathrm{}8`$. This shows that the interface is simply self-affine, not multifractal. We note, however, that the interface exponents of the Manna sandpile are not those of the one-dimensional LIM. This is most likely due to the fact that, perhaps differently from the two- and higher-dimensional cases, here it is important that the noise term increases in strength with the propagation of the interface, or with sustained activity. This arises from the fact that the fluctuations in the number of grains received at each site, reflected in $`\tau (H,i)`$, follow random-walk characteristics. Thus the anomaly exponent $`\kappa `$ indicates an even stronger dependence of the step height on $`L`$ at saturation than in the LIM. The same is also true if the step height is considered as a function of time for $`t<t_c`$. We find that $`\kappa _{Manna}0.82\kappa _{LIM}+0.5`$. Due to this time- or height-dependence of $`\tau `$ there is a direct coupling of $`t_c`$ and $`w_c`$ with the final average height, $`H_c`$ (see Fig. 7). As discussed above, this result follows, in the interface picture, from the fact that the noise strength increases with time. In the LIM, $`t_c`$ and $`w_c`$ are not correlated on a sample-to-sample basis.

### C Initial Relaxation

A comparison of our results for stationary correlations ($`z1.45`$) with those for the interface ($`z1.7`$) suggests that the relaxation to the steady state and relaxation of fluctuations within the stationary regime show distinct scaling properties. This is supported by our results for the relaxation of the active-site density. At the critical point of a simple absorbing-state model such as the contact process, the activity density $`\rho _a`$ exhibits an initial power-law decay, $`\rho _at^\theta `$, followed by a crossover to the quasi-stationary value $`\rho _aL^{\beta /\nu _{}}`$ . As noted above the growth exponent is related to the activity-decay exponent via $`\theta +\beta _W=1`$ if only one timescale is present . A plot of $`L^{\beta /\nu _{}}\rho _a(t)`$ versus $`t/L^z`$ yields a data collapse to a scaling function that is independent of $`L`$. In the present case (Fig. 8), we see that the collapse is imperfect, and that the form of $`\rho _a(t)`$ changes with $`L`$. Here $`z`$ was chosen so as to optimize the collapse at long times, yielding $`z=1.75(3)`$. For large systems, the active-site density exhibits three distinct regimes before reaching the quasi-stationary state: an initial power-law decay (I), followed by a crossover to a slower power-law regime (II), and finally a rapid decay (III) to the stationary state. For $`L=20048`$, the exponents associated with regimes I and II are 0.163 and 0.144, respectively. While the latter exponent is in reasonable agreement with the scaling relation $`\beta _W+\theta =1`$, it is clear that relaxation to the stationary state is more complicated for the sandpile than for, say, the contact process, which presents a unique power-law regime. A qualitative explanation may be found in the interface representation: the initial dynamics is dominated by relaxation of the initial grain profile $`z_i`$, which in the interface language means that at short times the columnar noise, $`F_i`$, dominates.

Thus another facet of the relaxation process is the approach of the mean height to its global value, $`\zeta `$, at a site with initial height $`z(0)`$. In Fig. 9 we plot the mean height $`z(t)|z(0)`$ in at system of 1400 sites at $`\zeta _c`$, averaged over 2000 trials. The inset shows that the asymptotic approach to $`\zeta `$ is approximately power-law, $`|z(t)|z(0)\zeta |t^\varphi `$, with $`\varphi =0.46`$, 0.45, 0.47, 0.50, and 0.53 for $`z(0)=0`$, 1, 2, 3 and 4, respectively. All of these exponents are close to $`\varphi =1/2`$, the value expected for uncorrelated diffusion.

### D Effects of initial preparation

In light of the complicated pattern of relaxation noted in the preceding subsection, it is well to recall that initial-decay studies of simple models such as the contact process employ a uniform, featureless initial configuration, far from the absorbing state. This is usually arranged by setting the the activity density to unity, an option that is not available in the present case, if we want to fix the particle density to its critical value $`\zeta _c`$. Since the height is a discrete variable, some local variations in density and activity are inevitable. All the results described up to now were obtained using the random deposition (RD) initial preparation described in Sec. II. In this subsection we report on studies using initial configurations in which the $`N=\zeta L`$ particles are distributed so as to reduce density fluctuations. We again employ the sequence $`L=2^n1253`$, $`N=2^n1189`$ ($`n=0,\mathrm{},4`$) corresponding to $`\zeta =0.94892`$.

One way of suppressing density fluctuations is via restricted random deposition (RRD): we divide the system into blocks, each consisting of 25 (or in some cases 26) sites, and depositing at random, within each block, a fixed number (24 or 25) of particles. In this way we generate a height distribution very similar to that of standard RD, but with virtually no density fluctuations on scales $`>50`$ sites. Our conclusions from RRD simulations are as follows. First, the stationary activity density is identical, to within uncertainty, to that found previously. Thus we may continue to use $`\beta /\nu _{}=0.239`$ in the scaling analysis. Figure 10, a scaling plot of the activity density (as in Fig. 8), shows a reasonable data collapse at long times, whereas the initial phase does not collapse. (It does collapse in an unrescaled plot, showing that the early stage is a size-independent, presumably purely local process.) Once again $`z`$ has been chosen to optimize the collapse at long times, which in this case gives $`z=1.47(6)`$. The final approach to the stationary state appears to follow a power law, $`\rho _at^\theta `$, with $`\theta =0.16(1)`$. Analysis of the interface representation yields $`z=1.57(3)`$, $`\alpha =1.13(2)`$ and $`\beta _W=0.77(3)`$, consistent with the scaling relation $`\beta _W=\alpha /z`$. (The growth of the survival time with $`L`$ suggests yet another value, $`z1.3`$.) Thus for RRD, all of the exponents describing transient properties are different from those found in RD studies. The exponent $`z`$ has decreased to a value consistent with stationary relaxation, perhaps due to suppression of slowly relaxing, long-wavelength density variations. The sum $`\theta +\beta _W`$, however, is now somewhat smaller than unity.

A second kind of initial configuration with reduced density fluctuations may be termed maximum activity (MA). In this case all sites are either empty or have height 2 (except for a single site, in case $`N`$ is odd), and the spacing between sites with $`z=2`$ is fixed at 1 or 2 sites, in a regular pattern, so the average density is constant on a scale of 50 or so sites. While the density is again uniform on large scales, the number of active sites is as large as possible (as in usual initial-decay studies), and consequently the height distribution is quite different from that in the RD or RRD configurations. For MA IC’s we find $`\alpha =1.12(4)`$, $`\beta _W=0.78(1)`$, and $`\theta =0.168(6)`$. From the collapse of $`\rho _a(t)`$ and of $`W^2(t)`$, we obtain, respectively, $`z=1.47(4)`$ and $`z=1.49(2)`$, while the survival time data again yield $`z1.3`$. Thus our results for the two kinds of initial conditions, RRD and MA, with suppressed density fluctuations are in good agreement. This suggests that the precise form of the initial height distribution does not affect asymptotic scaling behavior.

A general conclusion from these studies is that suppressing density fluctuations reduces the apparent value of the exponent $`z`$, suggesting that the unusually slow relaxation observed using random deposition initial conditions is due to long-wavelength density variations. The exponents $`\alpha `$ and $`\beta _W`$ decrease as well, indicating that such density variations also provide a source of roughness for the interface dynamics. Further studies are required to determine whether the interface-growth exponents reported here are typical of an entire class of initial conditions (i.e., with reduced density fluctuations), or if their values may depend on other details of the initial preparation. Until such studies are performed, we cannot be sure that the exponents describing roughness and initial relaxation are well defined, describing asymptotic evolution, or merely represent apparent power laws over some limited range of times and system sizes. In none of the cases studied do we observe the simple initial relaxation of $`\rho _a`$ typical of the contact process. This may be due to the impossibility of preparing a completely uniform, high-activity initial state, and/or to the conservation of particles.

### E Ergodicity

In a recent study we found strong nonergodic effects in the BTW sandpile . Given the randomness in the Manna sandpile rules, we do not expect such effects here. To test this hypothesis we plot (see Fig. 11) the active-site density $`\rho _a`$ versus time in four individual trials in a system of 1000 sites near $`\zeta _c`$. The curves for $`\rho _a`$ in the various trials are completely intermixed; their ranges of fluctuation about a common mean are roughly equal. This is in complete contrast to the BTW case, where a similar plot shows distinct mean values, and very different ranges of variation of $`\rho _a`$, in each trial . Thus the Manna sandpile appears to be ergodic, in the sense that a time-average of $`\rho _a`$ in a single realization is independent of the initial condition (after a short transient), and yields the same value as the stationary average of $`\rho _a`$ over many trials. Nonergodicity in the BTW model has been associated with “hidden” conservation laws, i.e., a memory of the initial configuration that persists over the entire trial. The Manna sandpile seems to be free of such conservation laws: while we can not rigorously exclude their presence, their effect must in any event be much weaker than in the BTW case. This is of course consistent with the fact that the BTW dynamics is deterministic, while the redistribution of particles is stochastic in the present case. (As a simple example of this, recall that a nearest-neighbor pair with both sites empty ($`z=0`$) is forbidden in the stationary state of the BTW sandpile, whereas no such restriction exists for stochastic toppling rules.)

## IV Discussion

We studied the scaling behavior of a one-dimensional fixed-energy sandpile with the same local dynamics as the Manna model. The model exhibits a continuous phase transition between an absorbing state and an active one at a critical particle density $`\zeta _c=0.94885(7)`$. In the interface (integrated activity) picture this translates into a depinning transition with a fixed force (fixed-energy model) or a force that is ramped up at an infinitesimal rate (driven case) .

The phase transition in the one-dimensional stochastic sandpile is characterized by the critical exponents $`\beta =0.42(2)`$ and $`\nu _{}=1.81(1)`$, which differ significantly from those associated with directed percolation ($`\beta =0.2765`$, $`\nu _{}=1.0968`$) and linear-interface depinning ($`\beta =0.25(3)`$, $`\nu _{}1.3`$). While absorbing-state phase transitions are expected to fall generically in the directed percolation universality class , it is reasonable to exempt the Manna model from this rule, due to local conservation of particles; this conservation law is expected to alter the universality class. In fact, studies of various models with the same local conservation law as the Manna sandpile, in dimensions $`d>1`$, indicate a new, common universality class for models sharing this feature .

Studying the interface representation of the model, we obtain the roughness exponent $`\alpha =1.48(2)`$ and growth exponent $`\beta _W=0.86(2)`$, which should be compared with $`\alpha =1.33(1)`$, $`\beta _W=0.839(1)`$ for DP and $`\alpha =1.25(1)`$, $`\beta _W=0.88(2)`$ for LIM. Study of the height-height correlation function yields the local roughness exponent $`\alpha _{loc}=0.59(3)`$; the corresponding DP value is $`0.63(3)`$ . Changing the initial condition to suppress long-wavelength density fluctuations, we obtain $`\alpha =1.13(2)`$ and $`\beta _W=0.77(2)`$, which again differ from the values associated with DP and LIM. Comparing these results with the apparent agreement of the Manna exponents with those of the LIM in two dimensions implies that either the rough numerical equivalence is fortuitous, or that the noise $`\tau `$ in the interface equation has a fundamentally different structure depending on the dimension. It is worth remarking that the measured roughness exponent is rather close to $`\alpha =3/2`$, which is the value one expects if only the columnar component of the noise is relevant . The other exponents, however, seem to be far from the columnar disorder universality class (i.e., $`\nu _{}=2`$, $`z=2`$, $`\beta =1`$, $`\beta _w=3/4`$) .

In linear interface models, translational invariance of the noise can be used to derive the scaling relation $`(2\alpha )\nu _{}=1`$ . This relation does not appear to be verified by our numerical results. Notice that in the context of interface depinning $`z`$ can be linked to the other exponents by a scaling relation  which reads $`z=\beta /\nu _{}+\alpha `$. Inserting the values of $`\beta `$, $`\nu _{}`$ and $`\alpha `$ measured in simulations, we obtain $`z=1.7`$, which is consistent with the exponent obtained analyzing the scaling of the width, for RD initial conditions. Our results for initial conditions with reduced fluctuations do not, however, follow this scaling relation. Understanding this in the interface picture by considering the height-dependent quenched noise term is a challenge.

Our results for the dynamic exponent $`z`$ are conflicting. We have studied various definitions of the relaxation or correlation time: that for the decay of the activity autocorrelation in the stationary state (yielding $`z=1.49(3)`$), one associated with the survival time ($`z=1.42(1)`$ for RD initial conditions, $`z1.3`$ for RRD and MA IC’s), and those associated with the surface roughness and the decay of activity ($`z=1.70(5)`$ for RD IC’s, $`z=1.52(5)`$ for initial configurations with suppressed density fluctuations). At a simple absorbing-state phase transition, all of these relaxation times are of course governed by the same exponent. Whether the unusually large $`z`$ associated with time-dependent behavior using random deposition derives, as we have suggested, from slow relaxation of long-wavelength density fluctuations in the initial configuration, is a subject for further investigation. Our studies of modified initial states support this notion, but to reach a definitive conclusion other classes of initial conditions will need to be examined. Indeed, we have noted that, given their apparent dependence on the initial state, it is not clear whether the exponents $`z`$, $`\alpha `$, $`\beta _W`$ and $`\theta `$ are well defined, i.e., whether the asymptotic behavior follows simple power laws.

We note that the present model does not exhibit the strong nonergodic effects observed in the fixed-energy version of the BTW sandpile. The relaxation of the mean height $`z_i(t)`$ from its initial value to the average, $`\zeta `$, follows a power law with an exponent $`1/2`$. We find good evidence for finite-size scaling, in contrast with most driven sandpiles . In summary, we have identified a one-dimensional sandpile model that exhibits an absorbing-state phase transition as the relevant temperature-like parameter (the energy density) is varied. It appears to be the “minimal model” for absorbing-state phase transitions belonging to a new universality class associated with a conserved density. Preliminary studies indicate that the driven version of the model exhibits scale-invariant avalanche statistics . We may therefore hope that analysis of the driven model, and of spreading of activity at the critical point of the fixed-energy system, will permit us to establish detailed connections between scale invariance under driving and the underlying absorbing-state phase transition.

Let us finally stress that while in higher dimensions the linear interface universality class and that of systems with absorbing states in the presence of a conserved static field seem to coincide, the results of this paper show that this equivalence appears to be broken in $`d=1`$. It will be interesting to study other one-dimensional systems with absorbing states and an order parameter coupled to a static conserved field  in order to compare the critical exponents and anomalies with those reported in this paper.

Acknowledgments

We thank R. Pastor-Satorras and F. van Wijland for helpful comments and discussions. RD thanks CNPq for financial support and CAPES and FAPEMG for support of computer facilities. We acknowledge partial support from the European Network contract ERBFMRXCT980183. MAM acknowledges support from the Spanish DGESIC project PB97-0842, and Junta de Andalucía project FQM-165.

FIGURE CAPTIONS

FIG. 1. Stationary active-site density versus energy density $`\zeta `$. The points represent extrapolations of data for $`L=100`$ \- 5000 to the $`L\mathrm{}`$ limit.

FIG. 2. Stationary active-site density versus system size. From bottom to top, $`\zeta `$ = 0.948, 0.94857, 0.94864, 0.94874, 0.9488, 0.94892, 0.949, and 0.95. The inset shows the curvature $`b`$ of the log-log plot as a function of $`\zeta `$, for $`L1000`$. The straight line is a least-squares linear fit.

FIG. 3. Inset: stationary autocorrelation function $`C(t)`$ versus $`t`$ at $`\zeta =0.9488`$, for (left to right) $`L=625`$, 1250, 2500, 5000, and $`10^4`$. In the main graph the data are plotted versus a rescaled time $`t^{}`$ defined in the text.

FIG. 4. Scaling plot of the active-site density versus $`\mathrm{\Delta }\zeta \zeta _c`$. Symbols: $`+`$: $`L=100`$; $``$: $`L=200`$; $`\times `$: $`L=500`$; $``$: $`L=1000`$; $`\mathrm{}`$: $`L=2000`$; $`\mathrm{}`$: $`L=5000`$.

FIG. 5. Scaling plot of $`\stackrel{~}{W}^2W^2/L^{2\alpha }`$ versus $`\stackrel{~}{t}t/L^z`$, for $`\zeta =0.9489`$, using $`\alpha =1.41`$, and $`z=1.645(13)`$. System sizes (top to bottom) $`L=1253`$, 2506, 5012, 10024, and 20048.

FIG. 6. Scaling plot of saturation width versus lifetime in individual trials, $`\zeta =0.9490`$, $`L=400`$-6400.

FIG. 7. Saturation height versus lifetime in individual trials, $`\zeta =0.9490`$, $`L=400`$-6400.

FIG. 8. Scaled active-site density $`\rho ^{}=L^{\beta /\nu _{}}\rho _a`$ versus scaled time $`t^{}=t/L^z`$ for $`\zeta =0.9489`$, system sizes $`L=1253`$,…, 20048 as indicated.

FIG. 9. Mean height $`z(t)|z(0)`$ of sites with initial height $`z(0)`$, in at system of 1400 sites at $`\zeta _c`$, averaged over 2000 trials. From bottom to top: $`z(0)`$ = 0, 1, 2, 3, and 4. The inset is a plot of $`\mathrm{ln}|z(t)|z(0)\zeta |`$ versus $`\mathrm{ln}t`$.

FIG. 10. Scaled active-site density versus scaled time for $`\zeta =0.9489`$ as in Fig. 8, but for RRD initial conditions.

FIG. 11. Active-site density $`\rho _a(t)`$ in four different trials at $`\zeta =0.95`$ in a system of 1000 sites.


format: text

text:
# Cosmological Adiabatic Geometric Phase of a Scalar Field in a Bianchi Spacetime11footnote 1Presented in the IX Regional Conference on Mathematical Physics, held in the Feza Gursey Institute, Istanbul, August 1999.

## I Introduction

In Ref.  a two-component formulation of the Klein-Gordon equation is used to develop relativistic analogues of the quantum adiabatic approximation and the adiabatic dynamical and geometric phases. This method provides a precise definition of an adiabatic evolution of a Klein-Gordon field in a curved background spacetime. The purpose of this article is to employ the results of Ref.  in the investigation of geometric phases due to a spatially homogeneous background spacetime.

The phenomenon of geometric phase in gravitational systems has been previously considered by Cai and Papini , Brout and Venturi , Venturi , Casadio and Venturi , Datta , and Corichi and Pierri . The motivation for these studies extends from the investigation of the study of weak gravitational fields  to various problems in quantum cosmology and quantum gravity .

These studies were mostly plagued by the problem of constructing appropriate inner products on the space of solutions of the Klein-Gordon equation. This problem has so far been solved for stationary spacetimes. Therefore the existing results have a limited domain of applicability. The most important feature of the method developed in Ref.  is that it avoids the above problem by showing that indeed the adiabatic geometric phase is independent of the choice of an inner product on the space of solutions of the field equation. This enables one to investigate the phenomenon of the adiabatic geometric phase for spatially homogeneous spacetimes which are clearly non-stationary. One must, however, note that the results presented in this article are relevant to adiabatically evolving Klein-Gordon fields . This raises the question whether the adiabaticity of the evolution is compatible with the fact that the background spacetime is non-stationary. It turns out that the answer to this question is in the positive, i.e., there are Klein-Gordon fields in a non-stationary spacetime which do have adiabatic evolutions.

The organization of the paper is as follows. In section II, the results of Ref.  which will be used in this paper are briefly reviewed. The computation of the geometric phase for general spatially homogeneous (Bianchi) models are discussed in Section III. These are applied in the analysis of the geometric phase problem for Bianchi type I and type IX models in sections IV and V, respectively. A summary of the main results and the concluding remarks are given in section VI.

Throughout this paper the signature of the spacetime metric g is taken to be $`(,+,+,+)`$. Letters from the beginning and the middle of the Greek alphabet are associated with an arbitrary local basis and a local coordinate basis of the tangent spaces of the spacetime manifold, respectively. The letters from the beginning and the middle of the Latin alphabet label the corresponding spatial components and take values in $`\{1,2,3\}`$.

## II Two-Component Formalism and the Adiabatic Geometric Phase

Consider a complex scalar field $`\mathrm{\Phi }`$ defined on a globally hyperbolic spacetime $`(M,\mathrm{g})=(\mathrm{IR}\times \mathrm{\Sigma },\mathrm{g})`$ satisfying

$$\left(g^{\mu \nu }_\mu _\nu \mu ^2\right)\mathrm{\Phi }=0,$$
(1)
where $`g^{\mu \nu }`$ are components of the inverse of the metric $`\mathrm{g}`$, $`_\mu `$ is the covariant derivative along $`/x^\mu `$ defined by the Levi Civita connection, and $`\mu `$ is the mass.

Denoting a time derivative by a dot, one can express Eq. (1) in the form

$$\ddot{\mathrm{\Phi }}+\widehat{D}_1\dot{\mathrm{\Phi }}+\widehat{D}_2\mathrm{\Phi }=0,$$
(2)
where

$`\widehat{D}_1`$ $`:=`$ $`{\displaystyle \frac{1}{g^{00}}}\left[2g^{0i}_ig^{\mu \nu }\mathrm{\Gamma }_{\mu \nu }^0\right],`$ (3)
$`\widehat{D}_2`$ $`:=`$ $`{\displaystyle \frac{1}{g^{00}}}\left[g^{ij}_i_jg^{\mu \nu }\mathrm{\Gamma }_{\mu \nu }^i_i\mu ^2\right].`$ (4)
A two-component representation of the field equation (2) is

$$i\dot{\mathrm{\Psi }}^{(q)}=\widehat{H}^{(q)}\mathrm{\Psi }^{(q)},$$
(5)
where

$`\mathrm{\Psi }^{(q)}`$ $`:=`$ $`\left(\begin{array}{c}u^{(q)}\\ v^{(q)}\end{array}\right),`$ (8)
$`u^{(q)}`$ $`:=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\mathrm{\Phi }+q\dot{\mathrm{\Phi }}),v^{(q)}:={\displaystyle \frac{1}{\sqrt{2}}}(\mathrm{\Phi }q\dot{\mathrm{\Phi }}),`$ (9)
$`\widehat{H}^{(q)}`$ $`:=`$ $`{\displaystyle \frac{i}{2}}\left(\begin{array}{cc}\frac{\dot{q}}{q}+\frac{1}{q}\widehat{D}_1q\widehat{D}_2& \frac{\dot{q}}{q}\frac{1}{q}+\widehat{D}_1q\widehat{D}_2\\ & \\ \frac{\dot{q}}{q}+\frac{1}{q}+\widehat{D}_1+q\widehat{D}_2& \frac{\dot{q}}{q}\frac{1}{q}\widehat{D}_1+q\widehat{D}_2\end{array}\right),`$ (13)
and $`q`$ is an arbitrary, possibly time-dependent, nonzero complex parameter.

Next consider the eigenvalue problem for $`H^{(q)}`$. Denoting the eigenvalues and eigenvectors by $`E_n^{(q)}`$ and $`\mathrm{\Psi }_n^{(q)}`$, i.e.,

$$H^{(q)}\mathrm{\Psi }_n^{(q)}=E_n^{(q)}\mathrm{\Psi }_n^{(q)},$$
(14)
one has

$$\mathrm{\Psi }_n^{(q)}=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1iqE_n^{(q)}\\ 1+iqE_n^{(q)}\end{array}\right)\mathrm{\Phi }_n^{(q)},$$
(15)
where $`\mathrm{\Phi }_n^{(q)}`$ satisfies

$$\left[\widehat{D}_2iE_n^{(q)}(\widehat{D}_1\frac{\dot{q}}{q})\left(E_n^{(q)}\right)^2\right]\mathrm{\Phi }_n^{(q)}=0.$$
(16)
This equation defines both the vectors $`\mathrm{\Phi }_n^{(q)}`$ and the complex numbers $`E_n^{(q)}`$.

The following is a summary of some of the results obtained in Ref. .

* Eq. (16) reduces to the ordinary eigenvalue equation for $`\widehat{D}_2`$, if $`\widehat{D}_1=\dot{q}/q`$. In this case $`\mathrm{\Phi }_n^{(q)}`$ and $`E_n^{(q)}`$ are independent of the choice of $`q`$, and one can drop the labels $`(q)`$ on the right hand side of Eq. (15).
* $`\mathrm{\Phi }_n^{(q)}`$ belong to the Hilbert space $`_t=L^2(\mathrm{\Sigma }_t)`$ of square-integrable functions on the spacelike hypersurfaces $`\mathrm{\Sigma }_t`$ where the integration is defined by the measure $`[det(^{(3)}\mathrm{g})]^{1/2}`$ and $`{}_{}{}^{(3)}\mathrm{g}`$ is the Riemannian three-metric on $`\mathrm{\Sigma }_t`$ induced by the four-metric $`\mathrm{g}`$.
* Suppose that

+ $`\widehat{D}_1=\dot{q}/q`$,
    + $`\widehat{D}_2`$ is self-adjoint with respect to the inner product $`,`$ of $`_t`$,
    + $`\widehat{D}_2`$ has a discrete spectrum,
    + during the evolution of the system $`E_nE_m`$ iff $`mn`$, i.e., there is no level-crossings, and
    + $`E_n`$ is a non-degenerate eigenvalue of $`H^{(q)}`$, alternatively $`E_n^2`$ is a non-degenerate eigenvalue of $`\widehat{D}_2`$.

Then for all $`mn`$, $`\mathrm{\Phi }_m|\dot{\mathrm{\Phi }}_n=\mathrm{\Phi }_m|\dot{\widehat{D}}_2\mathrm{\Phi }_n/(E_n^2E_m^2)`$. The relativistic adiabatic approximation corresponds to the case where the latter may be neglected. In this case, an initial two-component Klein-Gordon field

$$\mathrm{\Psi }^{(q)}(0)=e^{i\alpha _n(0)}\mathrm{\Psi }_n^{(q)}(0)+e^{i\alpha _n(0)}\mathrm{\Psi }_n^{(q)}(0),$$
(17)
with $`n0`$ and $`\alpha _{\pm n}(0)\text{ }\mathrm{C}`$, evolves according to

$$\mathrm{\Psi }^{(q)}(t)e^{i\alpha _n(t)}\mathrm{\Psi }_n^{(q)}(t)+e^{i\alpha _n(t)}\mathrm{\Psi }_n^{(q)}(t),$$
(18)
where $`\alpha _{\pm n}(t)=[\alpha _n(0)+\alpha _n(0)]/2+\gamma _n(t)+\delta _{\pm n}(t)`$, $`\alpha _{\pm n}(0)`$ are arbitrary constants, $`\gamma _n(t)`$ is the geometric part of both $`\alpha _{\pm n}(t)`$ and $`\delta _{\pm n}(t)`$ is the dynamical part of $`\alpha _{\pm n}(t)`$. They are given by

$`\gamma _n(t)`$ $`=`$ $`{\displaystyle _{R(0)}^{R(t)}}𝒜_n[R],`$ (19)
$`\delta _\pm (t)`$ $`=`$ $`i\xi _n(t)\pm {\displaystyle \frac{\eta _n(t)}{2}},`$ (20)
where

$$𝒜_n[R]:=\frac{i\mathrm{\Phi }_n[R]|d|\mathrm{\Phi }_n[R]}{\mathrm{\Phi }_n[R]|\mathrm{\Phi }_n[R]}=\frac{i\mathrm{\Phi }_n[R]|\frac{}{R^a}|\mathrm{\Phi }_n[R]}{\mathrm{\Phi }_n[R]|\mathrm{\Phi }_n[R]}dR^a,$$
(21)
is the Berry’s connection one-form , $`R=(R^1,\mathrm{},R^n)`$ are the parameters of the system, i.e., the components of the metric, $`d`$ stands for the exterior derivative with respect to $`R^a`$,

$`\xi _n(t)`$ $`:=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^t}f_n(t^{})[1\mathrm{cos}\eta _n(t^{})]𝑑t^{},`$ (22)
$`f_n(t)`$ $`:=`$ $`{\displaystyle \frac{d}{dt}}\mathrm{ln}[q(t)E_n(t)],`$ (23)
and $`\eta _n`$ is the solution of

$$\dot{\eta }_n+f_n\mathrm{sin}\eta _n+2E_n=0,\eta _n(0)=\alpha _n(0)\alpha _n(0).$$
(24)
In view of Eqs. (6), (7), and (10), Eq. (18) can be written in the form

$$\mathrm{\Phi }=ce^{i\gamma _n(t)}(e^{i\delta _n(t)}+e^{i\delta _n(t)})\mathrm{\Phi }_n,$$
(25)
where $`\mathrm{\Phi }`$ is the one-component Klein-Gordon field, i.e., the solution of the original Klein-Gordon equation (1), and $`c:=\mathrm{exp}[\alpha _n(0)+\alpha _n(0)]`$.
* If in addition to the above adiabaticity condition one also has $`\dot{E}_n0`$, then an initial eigenvector $`\mathrm{\Psi }_n^{(q)}(0)`$ evolves according to

$$\mathrm{\Psi }^{(q)}(t)e^{i\alpha _n(t)}\mathrm{\Psi }_n^{(q)}(0),$$

where the total phase angle $`\alpha _n`$ again consists of a geometric and a dynamical part. This case corresponds to what is termed as ultra-adiabatic evolution in Ref. .

Suppose that conditions 3.1) – 3.4) are satisfied, but $`E_n`$ is $`𝒩`$ fold degenerate. Then the condition for the validity of the adiabatic approximation becomes $`\mathrm{\Phi }_m^I|\dot{\mathrm{\Phi }}_n^J0`$, for all $`mn`$ and $`I,J=1,2,\mathrm{}𝒩`$. Here $`\mathrm{\Phi }_n^1,\mathrm{}\mathrm{\Phi }_n^𝒩`$ are orthogonal eigenvectors spanning the degeneracy subspace of $`_t`$ associated with $`E_n^2`$. In this case, $`\mathrm{exp}[i\alpha _{\pm n}(t)]`$ become matrices of the form $`e^{i\delta _{\pm n}(t)}\mathrm{\Gamma }_n(t)`$ where

$$\mathrm{\Gamma }_n(t):=𝒫\mathrm{exp}[i_{R(0)}^{R(t)}𝒜_n],$$
(26)
$`𝒫`$ is the path-ordering operator, $`𝒜_n`$ is a matrix of one-forms with components

$$𝒜_n^{IJ}[R]:=\frac{i\mathrm{\Phi }_n^I[R]|d|\mathrm{\Phi }_n^J[R]}{\mathrm{\Phi }_n^I[R]|\mathrm{\Phi }_n^I[R]}.$$
(27)
In this case, a solution of the one-component Klein-Gordon equation (1) is given by

$$\mathrm{\Phi }=\underset{I,J=1}{\overset{N}{}}(c_n^Ie^{i\delta _n(t)}+c_n^Ie^{i\delta _n(t)})\mathrm{\Gamma }_n^{IJ}(t)\mathrm{\Phi }_n^J,$$
(28)
where $`c_{\pm n}^I`$ are constant coefficients determined by the initial conditions.

If the parameters $`R`$ undergo a periodic change, i.e., $`R(T)=R(0)`$ for some $`T\mathrm{IR}^+`$, then the path-ordered exponential $`\mathrm{\Gamma }_n(T)`$ which is called the cyclic adiabatic geometrical phase, cannot be removed by a gauge transformation. Eq. (27) shows that the formula for the relativistic adiabatic geometric phase has the same form as its non-relativistic analogue . In particular for $`𝒩=1`$ (the non-degenerate case), one recovers Berry’s connection one-form (21). In this case $`\mathrm{\Gamma }_n(T)`$ reduces to an ordinary exponential and yields the Berry phase $`e^{i\gamma (T)}`$ for the Klein-Gordon field.

In the remainder of this paper, I shall consider the problem of the adiabatic geometric phase for a Klein-Gordon field minimally coupled to a spatially homogeneous gravitational field.

It is important to note that the cyclic geometric phase has physical significance, if one has a cyclic evolution222For a discussion of the meaning of a cyclic evolution of a Klein-Gordon field see Ref. .. For an adiabatic evolution, the evolving state undergoes a cyclic evolution, if the parameters of the system, in this case the components of the metric tensor, change periodically in time. This corresponds to the spatially homogeneous (Bianchi) cosmological models which admit periodic (oscillatory) solutions333For an example see Ref. .. This observation does not however mean that the connection one-forms $`𝒜_n`$ and their path-ordered exponentials $`\mathrm{\Gamma }_n(t)`$ lack physical significance for general nonperiodic Bianchi models. The cyclic adiabatic geometric phases have noncyclic analogues which occur in the evolution of any quantum state, .

A noncyclic analogue of the non-Abelian cyclic geometric phase has recently been introduced by the present author . In view of the results of Ref. , the noncyclic adiabatic geometric phase for an adiabatically evolving Klein-Gordon field is given by

$$\stackrel{~}{\mathrm{\Gamma }}_n(t):=w_n(t)\mathrm{\Gamma }_n(t),$$
(29)
where $`w_n(t)`$ is an $`𝒩\times 𝒩`$ matrix with entries:

$$w_n^{IJ}(t):=\mathrm{\Phi }_n^I[R(0)]|\mathrm{\Phi }_n^J[R(t)].$$
(30)
The noncyclic geometric phase has the same gauge transformation properties as the cyclic geometric phase. In particular, its eigenvalues and its trace are gauge-invariant physical quantities, .

If it happens that the connection one-form $`𝒜_n`$ is exact, i.e., it is a pure gauge, then there are two possibilities:

* The curve $`C`$ traced by the parameters $`R`$ of the system in time has a part which is a non-contractible loop. In this case the cyclic or noncyclic geometric phase will be a topological quantity analogous to the Aharonov-Bohm phase . Such a geometric phase will be called a topological phase. Topological phases can occur only if the parameter space of the system has a nontrivial first homology group.
* The curve $`C`$ does not have a piece which is a non-contractible loop. In this case, the geometric phase is either unity (the cyclic case) or it depends only on the end points, $`R(0)`$ and $`R(t)`$, of $`C`$ (the noncyclic case). Such a geometric phase will be called a trivial geometric phase.

## III Spatially Homogeneous Cosmological Models

Consider Klein-Gordon fields in a spatially homogeneous (Bianchi) cosmological background associated with a Lie group $`G`$, i.e., $`M=\mathrm{IR}\times G`$. In a synchronous invariant basis the spacetime metric g is given by its spatial components $`g_{ab}`$:

$$ds^2=g_{\alpha \beta }\omega ^\alpha \omega ^\beta =dt^2+g_{ab}\omega ^a\omega ^b,$$
(31)
where $`\omega ^a`$ are the left invariant one-forms and $`g_{ab}=g_{ab}(t)`$. Throughout this article I use the conventions of Ref. .

The first step in the study of the phenomenon of the adiabatic geometric phase due to a spatially homogeneous cosmological background is to compute the operators $`\widehat{D}_1`$ and $`\widehat{D}_2`$ of Eqs. (3) and (4) in the invariant basis. It is not difficult to see that with some care these equations are valid in any basis. One must only replace the coordinate labels ($`\mu ,\nu ,\mathrm{},i,j,\mathrm{}`$) with the basis (in this case invariant basis) labels $`(\alpha ,\beta ,\mathrm{},a,b,\mathrm{})`$, and interpret $`_a`$ as the action of the operators $`\widehat{X}_a`$ associated with the dual vector fields to $`\omega ^a`$. This leads to

$`\widehat{D}_1`$ $`=`$ $`g^{ab}\mathrm{\Gamma }_{ab}^0,`$ (32)
$`\widehat{D}_2`$ $`=`$ $`\mathrm{\Delta }_t+\mu ^2,`$ (33)
$`\mathrm{\Delta }_t`$ $`:=`$ $`g^{ab}_a_b=g^{ab}\widehat{X}_a\widehat{X}_b\mathrm{\Gamma }_{ab}^c\widehat{X}_c,`$ (34)
where $`\mathrm{\Delta }_t`$ is the Laplacian on $`\mathrm{\Sigma }_t`$, $`_a`$ are the covariant derivatives corresponding to the Levi Civita connection,

$$\mathrm{\Gamma }_{\alpha \beta }^\gamma :=\frac{1}{2}g^{\gamma \delta }(g_{\delta \alpha ,\beta }+g_{\beta \delta ,\alpha }g_{\alpha \beta ,\delta }+g_{ϵ\alpha }C_{\delta \beta }^ϵ+g_{ϵ\beta }C_{\delta \alpha }^ϵ)\frac{1}{2}C_{\alpha \beta }^\gamma ,$$
(35)
as derived in Ref. , $`g_{\alpha \beta ,\gamma }:=\widehat{X}_\gamma g_{\alpha \beta }`$, and $`C_{\alpha \beta }^\gamma `$ are the structure constants:

$$[\widehat{X}_\alpha ,\widehat{X}_\beta ]=C_{\alpha \beta }^\gamma \widehat{X}_\gamma ,$$
(36)
with $`\widehat{X}_0:=/t`$. In view of the latter equality, the structure constants with a time label vanish. This simplifies the calculations of $`\mathrm{\Gamma }`$’s. The only nonvanishing ones are

$`\mathrm{\Gamma }_{ab}^0`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{g}_{ab},`$ (37)
$`\mathrm{\Gamma }_{ab}^c`$ $`=`$ $`{\displaystyle \frac{1}{2}}g^{cd}(g_{ea}C_{db}^e+g_{eb}C_{da}^e){\displaystyle \frac{1}{2}}C_{ab}^c.`$ (38)
In view of these relations, the expression for $`\widehat{D}_1`$ and $`\widehat{D}_2`$ may be further simplified:

$`\widehat{D}_1`$ $`=`$ $`{\displaystyle \frac{}{t}}\mathrm{ln}\sqrt{g},`$ (39)
$`\widehat{D}_2`$ $`=`$ $`\mathrm{\Delta }_t+\mu ^2=(g^{ab}\widehat{X}_a\widehat{X}_bC_{ab}^bg^{ac}\widehat{X}_c)+\mu ^2,`$ (40)
where $`g`$ is the determinant of $`(g_{ab})`$. Note that for a unimodular, in particular semisimple, group $`C_{ab}^b=0`$, and the second term on the right hand side of (40) vanishes. The corresponding Bianchi models are knows as Class A models.

As seen from Eq. (39), $`\widehat{D}_1`$ acts by multiplication by a time-dependent function. Therefore, choosing $`q=i\sqrt{g}`$, $`\widehat{D}_1=\dot{q}/q`$. This reduces Eq. (16) to the eigenvalue equation

$$(\widehat{D}_2E_n^2)\mathrm{\Phi }_n=(\mathrm{\Delta }_t+E_n^2\mu ^2)\mathrm{\Phi }_n=0,$$
(41)
for the operator $`\widehat{D}_2`$ which being essentially the Laplacian over $`\mathrm{\Sigma }_t`$, is self-adjoint. This guarantees the orthogonality of $`\mathrm{\Phi }_n`$ and the reality of $`E_n^2`$. 444One can also show that, since $`q`$ is imaginary, the Hamiltonian $`H^{(q)}`$ is self-adjoint with respect to the Klein-Gordon inner product, i.e., the inner product (11) of Ref. .

The analysis of the $`\mathrm{\Phi }_n`$ is equivalent to the study of the eigenvectors of the Laplacian over a three-dimensional group manifold $`\mathrm{\Sigma }_t`$. The general problem is the subject of the investigation in spectral geometry which is beyond the scope of the present article. However, let us recall some well-known facts about spectral properties of the Laplacian $`\mathrm{\Delta }`$ for an arbitrary finite-dimensional Riemannian manifold $`\mathrm{\Sigma }`$ without boundary.

The following results are valid for the case where $`\mathrm{\Sigma }`$ is compact or the eigenfunctions are required to have a compact support555This is equivalent with the case where $`\mathrm{\Sigma }`$ has a boundary $`\mathrm{\Sigma }`$, over which the eigenfunctions vanish. :

* The spectrum of $`\mathrm{\Delta }`$ is an infinite discrete subset of non-negative real numbers.
* The eigenvalues are either non-degenerate or finitely degenerate.
* There is an orthonormal set of eigenfunctions which form a basis for $`L^2(\mathrm{\Sigma })`$.
* If $`\mathrm{\Sigma }`$ is compact, then the first eigenvalue is zero which is non-degenerate with the eigenspace given by the set of constant functions, i.e., $`\text{ }\mathrm{C}`$. If $`\mathrm{\Sigma }`$ is not compact but the eigenfunctions are required to have a compact support, then the first eigenvalue is positive.

Another piece of useful information about the spatially homogeneous cosmological models is that (up to a multiple of $`i=\sqrt{1}`$) the invariant vector fields $`\widehat{X}_a`$ yield a representation of the generators $`L_a`$ of $`G`$, with $`L^2(\mathrm{\Sigma }_t)`$ being the representation space, one can view the Laplacian $`\mathrm{\Delta }_t`$ as (a representation of ) an element of the enveloping algebra of the Lie algebra of $`G`$. Therefore, $`\mathrm{\Delta }_t`$ commutes with any Casimir operator $`𝒞_\lambda `$ and consequently shares a set of simultaneous eigenvectors with $`𝒞_\lambda `$. This in turn suggests one to specialize to particular subrepresentations with definite $`𝒞_\lambda `$. In particular for compact groups, this leads to a reduction of the problem to a collection of finite-dimensional ones.666Here I mean a finite-dimensional Hilbert space.

In the remainder of this article I shall try to employ these considerations to investigate some specific models.

## IV Bianchi Type I

In this case $`G`$ is Abelian, therefore $`X_a`$ are themselves Casimir operators and the eigenfunctions of $`\mathrm{\Delta }_t`$, i.e., $`\mathrm{\Phi }_n`$, are independent of $`t`$. Hence the Berry connection one-form (21) and its non-Abelian generalization (27) vanish identically, and the geometric phase is trivial.

## V Bianchi Type IX

In this case $`G=SU(2)=S^3`$. The total angular momentum operator $`\widehat{J}^2=_a\widehat{J}_a^2`$ is a Casimir operator. Therefore, I shall consider the subspaces $`_j`$ of $`_t=L^2(S_t^3)`$ of definite angular momentum $`j`$. The left-invariant vector fields $`\widehat{X}_a`$ are given by $`\widehat{X}_a=i\widehat{J}_a`$, in terms of which Eq. (36), with $`C_{ab}^c=ϵ_{abc}`$, is written in the familiar form:

$$[\widehat{J}_a,\widehat{J}_b]=iϵ_{abc}\widehat{J}_c,$$
(42)
with $`ϵ_{abc}`$ denoting the totally antisymmetric Levi Civita symbol and $`ϵ_{123}=1`$.

Eq. (41) takes the form:

$$(\widehat{H^{}}+k_n^2)\mathrm{\Phi }_n=0,$$
(43)
where $`\widehat{H^{}}`$ is an induced Hamiltonian defined by

$$\widehat{H}^{}:=g^{ab}(t)\widehat{J}_a\widehat{J}_b,$$
(44)
and $`k_n^2:=E_n^2\mu ^2`$. Therefore, the problem of the computation of the geometric phase is identical with that of the non-relativistic quantum mechanical system whose Hamiltonian is given by (44). In particular, for the mixmaster spacetime, i.e., for $`g_{ab}`$ diagonal, the problem is identical with the quantum mechanical problem of a non-relativistic asymmetric rotor, .

Another well-known non-relativistic quantum mechanical effect which is described by a Hamiltonian of the form (44) is the quadratic interaction of a spin with a variable electric field $`(E_a)`$. The interaction potential is the Stark Hamiltonian: $`\widehat{H}_S=ϵ(_aE_a\widehat{J}_a)^2`$. The phenomenon of the geometric phase for the Stark Hamiltonian for spin $`j=3/2`$, which involves Kramers degeneracy , was first considered by Mead . Subsequently, Avron, et al  conducted a thorough investigation of the traceless quadrupole Hamiltonians of the form (44).

The condition on the trace of the Hamiltonian is physically irrelevant, since the addition of any multiple of the identity operator to the Hamiltonian does not have any physical consequences. In general, one can express the Hamiltonian (44) in the form $`\widehat{H^{}}=\widehat{\stackrel{~}{H}^{}}+\widehat{H}_0^{}`$, where $`\widehat{\stackrel{~}{H}^{}}:=Tr(g^{ab})\widehat{J}^2/3`$,

$$\widehat{H}_0^{}[R]:=\underset{A=1}{\overset{5}{}}R^A\widehat{e}_A,$$
(45)
is the traceless part of the Hamiltonian, and

$`\widehat{e}_1`$ $`:=`$ $`J_3^2{\displaystyle \frac{1}{3}}\widehat{J}^2,\widehat{e}_2:={\displaystyle \frac{1}{\sqrt{3}}}\{\widehat{J}_1,\widehat{J}_2\},`$
$`\widehat{e}_3`$ $`:=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}\{\widehat{J}_2,\widehat{J}_3\},\widehat{e}_4:={\displaystyle \frac{1}{\sqrt{3}}}(\widehat{J}_1^2\widehat{J}_2^2),`$ (46)
$`\widehat{e}_5`$ $`:=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}\{\widehat{J}_1,\widehat{J}_2\},`$
$`R^1`$ $`:=`$ $`g^{33}{\displaystyle \frac{1}{2}}(g^{11}+g^{22}),R^2:=\sqrt{3}g^{13},`$
$`R^3`$ $`:=`$ $`\sqrt{3}g^{23},R^4:={\displaystyle \frac{\sqrt{3}}{2}}(g^{11}g^{22}),`$ (47)
$`R^5`$ $`:=`$ $`\sqrt{3}g^{12}.`$
As shown in Refs.  the space $`^{}`$ of all traceless Hamiltonians of the form (45) is a five-dimensional real vector space. The traceless operators $`\widehat{e}_A`$ form an orthonormal 777Orthonormality is defined by the inner product $`\widehat{A},\widehat{B}:=3Tr(\widehat{A}\widehat{B})/2`$. basis for $`^{}`$. Removing the point $`(R^A=0)`$ from $`^{}`$, to avoid the collapse of all eigenvalues, one has the space $`\mathrm{IR}^5\{0\}`$ as the parameter space. The situation is analogous to Berry’s original example of a magnetic dipole in a changing magnetic field, . Again, a rescaling of the Hamiltonian by a non-zero function of $`R^A`$ does not change the geometric phase. Thus the relevant parameter space is $`=S^4`$. Incidentally, the point corresponding to $`0\mathrm{IR}^5`$ which is to be excluded, corresponds to the class of Friedmann-Robertson-Walker models.

Following Berry  and Ref. , let us identify $`t`$ with the affine parameter of a curve $`C:[0,\tau ]S^4`$ traced by the parameters $`R^A`$ in $`S^4`$. Such a curve may be defined by the action of the group $`SO(5)`$ which acts transitively on $`S^4`$. Therefore the time-dependence of the Hamiltonian may be realized by an action of the group $`SO(5)`$ on a fixed Hamiltonian. As it is discussed in Ref. , it is the unitary representations $`𝒰`$ of the double cover $`Spin(5)=Sp(2)`$ of $`SO(5)`$ (alternatively the projective representations of $`SO(5)`$) which define the time-dependent Hamiltonian:

$$\widehat{H}_0^{}[R(t)]=𝒰[g(t)]\widehat{H}_0^{}[R(0)]𝒰[g(t)]^{}.$$
(48)
Here, $`g:[0,\tau ]Sp(2)`$, is defined by $`R(t)=:\pi [g(t)]R(0)`$, where $`\pi :Sp(2)SO(5)`$ is the canonical two-to-one covering projection and $`R(t)`$ corresponds to the point $`C(t)S^4`$. The emergence of the group $`Sp(2)`$ is an indication of the existence of a quaternionic description of the system, .

Let us next examine the situation for irreducible representations $`j`$ of $`SU(2)`$. As I previously described, $`\widehat{J}^2`$ commutes with the Hamiltonian. Hence the Hamiltonian is block-diagonal in the basis with definite total angular momentum $`j`$. For each $`j`$, the representation space $`_j`$ is $`2j+1`$ dimensional. Therefore the restriction of the Hamiltonian $`\widehat{H^{}}_0`$ to $`_j`$ and $`𝒰[g(t)]`$ are $`(2j+1)\times (2j+1)`$ matrices.

Let $`\{\varphi _{j_n}^I\}`$ be a complete set of orthonormal eigenvectors of the initial Hamiltonian $`\widehat{H}_0^{}[R(0)]`$, where $`I`$ is a degeneracy label. Then the eigenvectors of $`\widehat{H}_0^{}[R]`$ are of the form

$$\mathrm{\Phi }_{j_n}^I[R]=𝒰[g]\varphi _{j_n}^I,$$
(49)
and the non-Abelian connection one-form (27) is given by

$$𝒜_{j_n}^{IJ}=i\varphi _{j_n}^I|𝒰^{}d𝒰|\varphi _{j_n}^J.$$
(50)
For the integer $`j`$, bosonic systems, it is known that the quadratic Hamiltonians of the form (45), describe time-reversal-invariant systems. In this case it can be shown that the curvature two-form associated with the Abelian Berry connection one-form (21) vanishes identically . The connection one-form is exact (gauge potential is pure gauge) and a nontrivial geometric phase can only be topological, namely it may still exist provided that the first homology group of the parameter space is nontrivial. For the problem under investigation $`=S^4`$, and the first homology group is trivial. Hence, in general, the Abelian geometric phase is trivial. The same conclusion cannot however be reached for the non-Abelian geometric phases.

In the remainder of this section, I shall examine the situation for some small values of $`j`$.

* $`j=0`$: The corresponding Hilbert subspace is one-dimensional. Geometric phases do not arise.
* $`j=1/2`$: In this case, $`\widehat{J}_a=\sigma _a/2`$, where $`\sigma _a`$ are Pauli matrices. Using the well-known anticommutation (Clifford algebra) relations $`\{\sigma _a,\sigma _b\}=2\delta _{ab}`$, one can easily show that in this case

$$\widehat{H}^{}=\frac{1}{2}\underset{a}{}g^{aa}(t)\widehat{I},$$
(51)
where $`\widehat{I}`$ is the identity matrix. Therefore, the eigenvectors $`\mathrm{\Phi }_n`$ are constant ($`g_{ab}`$-independent), the connection one-form (50) vanishes, and no nontrivial geometric phases occur.
* $`j=1`$: In this case the Hilbert subspace is three-dimensional. The Abelian geometrical phases are trivial. The nontrivial matrix-valued geometrical phases may be present, provided that the Hamiltonian has a degenerate eigenvalue. Using the ordinary $`j=1`$ matrix representations of $`\widehat{J}_a`$, one can easily express $`\widehat{H}_0^{}`$ as a $`3\times 3`$ matrix. It can then be checked that in the generic case the eigenvalues of $`\widehat{H}_0^{}`$ are not degenerate.888This is true for all integer $`j`$, . However, there are cases for which a degenerate eigenvalue is present. A simple example is the Taub metric, $`(g_{ab})=\mathrm{diag}(g_{11},g_{22},g_{22})`$, which is a particular example of the mixmaster metric, . For the general mixmaster metric, the eigenvalue problem can be easily solved. The eigenvalues of the total Hamiltonian $`\widehat{H}^{}`$ of Eq. (44) are

$$\frac{1}{g_{11}}+\frac{1}{g_{22}},\frac{1}{g_{11}}+\frac{1}{g_{33}},\frac{1}{g_{22}}+\frac{1}{g_{33}}.$$

Therefore the degenerate case corresponds to the coincidence of at least two of $`g_{aa}`$’s, i.e., a Taub metric. However, even in the general mixmaster case, one can find a constant ($`g_{aa}`$-independent) basis which diagonalizes the Hamiltonian. Hence the non-Abelian connection one-form vanishes and the geometric phase is again trivial. This is not however the case for general metrics. In the Appendix, it is shown that without actually solving the general eigenvalue problem for the general Hamiltonian, one can find the conditions on the metric which render one of the eigenvalues of the Hamiltonian degenerate. Here, I summarize the results. Using the well-known matrix representations of the angular momentum operators $`\widehat{J}_a`$ in the $`j=1`$ representation  one can write the Hamiltonian (44) in the form:

$$\widehat{H}^{}=\left(\begin{array}{ccc}t+2z& \xi ^{}& \zeta ^{}\\ \xi & 2(t+z)& \xi ^{}\\ \zeta & \xi & t+2z\end{array}\right),$$
(52)
where

$`t`$ $`:=`$ $`{\displaystyle \frac{1}{2}}(g^{11}+g^{22})g^{33},z:=g^{33},`$
$`\xi `$ $`:=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(g^{13}+ig^{23}),\zeta :={\displaystyle \frac{1}{2}}(g^{11}g^{22})+ig^{12}.`$
Then it can be shown (Appendix) that the necessary and sufficient conditions for $`\widehat{H}^{}`$ to have a degenerate eigenvalue are

+ for $`\zeta =0`$: $`\xi =0`$, in which case, $`\widehat{H}^{}`$ as given by Eq. (52) is already diagonal. The degenerate and non-degenerate eigenvalues are $`t+2z`$ and $`2(t+z)`$, respectively. In terms of the components of the metric, these conditions can be written as: $`g_{11}=g_{22}`$ and $`g_{ab}=0`$ if $`ab`$. This is a Taub metric which as discussed above does not lead to a nontrivial geometric phase.
    + for $`\zeta 0`$:

    $$\zeta =𝒵e^{2i\theta },t=𝒵|\xi |^2/𝒵,$$
    (53)
    where $`\mathrm{exp}[i\theta ]:=\xi /|\xi |`$ and $`𝒵\mathrm{IR}\{0\}`$. In this case the degenerate and non-degenerate eigenvalues are $`2(𝒵+z)|\xi |^2/𝒵`$ and $`2(z|\xi |^2/𝒵)`$, respectively.

For the latter case, an orthonormal set of eigenvectors is given by:

$$v_1=\frac{1}{\sqrt{2}}\left(\begin{array}{c}e^{2i\theta }\\ 0\\ 1\end{array}\right),v_2=\frac{1}{\sqrt{1+2𝒳^2}}\left(\begin{array}{c}𝒳e^{i\theta }\\ 1\\ 𝒳e^{i\theta }\end{array}\right),v_3=\frac{1}{\sqrt{2(1+2𝒳^2)}}\left(\begin{array}{c}e^{2i\theta }\\ 2𝒳e^{i\theta }\\ 1\end{array}\right),$$
(54)
where $`𝒳:=|\xi |/(2𝒵)`$. In view of the general argument valid for all non-degenerate eigenvalues, the geometric phase associated with $`v_3`$ is trivial. This can be directly checked by substituting $`v_3`$ in the formula (21) for the Berry connection one-form. This leads, after some algebra, to the surprisingly simple result $`𝒜_{33}:=iv_3|dv_3=d\theta `$. Therefore, $`𝒜_{33}`$ is exact as expected, and the corresponding geometric phase is trivial. Similarly one can compute the matrix elements $`𝒜_{rs}:=iv_r|dv_s`$, $`r,s=1,2`$, of the non-Abelian connection one-form (27). The result is

$`𝒜`$ $`=`$ $`\left(\begin{array}{cc}1& \\ ^{}& 0\end{array}\right)\omega ,`$ (57)
$``$ $`:=`$ $`{\displaystyle \frac{2𝒳e^{i\theta }}{\sqrt{2(1+2𝒳^2)}}}={\displaystyle \frac{2ϵ\xi }{\sqrt{2+\left|\frac{\xi }{\zeta }\right|^2}}}={\displaystyle \frac{ϵ(g^{13}+ig^{23})}{\sqrt{1+\frac{(g^{13})^2+(g^{23})^2}{(g^{11}g^{22})^2+(2g^{12})^2}}}},`$
$`\omega `$ $`:=`$ $`d\theta ={\displaystyle \frac{g^{13}dg^{23}g^{23}dg^{13}}{(g^{13})^2+(g^{23})^2}},`$
where $`ϵ:=𝒵/|\zeta |=\pm 1`$. As seen from Eq. (57), $`𝒜`$ is a $`u(2)`$-valued one-form, which vanishes if $`g^{23}/g^{13}`$ is kept constant during the evolution of the universe.

It is also worth mentioning that the requirement of the existence of degeneracy is equivalent to restricting the parameters of the system to a two-dimensional subset of $`S^4`$. Thus, the corresponding spectral bundle  is a $`U(2)`$ vector bundle over a two-dimensional parameter space $`\stackrel{~}{}`$. The manifold structure of $`\stackrel{~}{}`$ is determined by Eqs. (53). In terms of the parameters $`R^A`$ of (47), these equations are expressed by

$$R^5=f_1R^4,R^1=f_2R^4+\frac{f_3}{R^4},$$
(58)
where

$`f_1`$ $`:=`$ $`{\displaystyle \frac{2R^2R^3}{(R^2)^2(R^3)^2}},f_2:=\pm {\displaystyle \frac{(R^2)^2+(R^3)^2}{\sqrt{3}[(R^2)^2(R^3)^2]}},`$
$`f_3`$ $`:=`$ $`{\displaystyle \frac{(R^2)^2(R^3)^2}{2\sqrt{3}}},f_4:=(R^2)^2+(R^3)^2.`$
Here $`f_4`$ is also introduced for future use. In addition to (58), one also has the condition $`(R^A)S^4`$. If $`S^4`$ is identified with the round sphere, this condition takes the form $`_A(R^A)^2=1`$. Substituting (58) in this equation, one finds

$$(1+f_2^2+f_3^2)(R^4)^4(1f_42f_2f_3)(R^4)^2+f_3^2=0.$$
(59)
Eq. (59) may be easily solved for $`R^4`$. This yields

$$R^4=\pm \frac{3[(R^2)^2(R^3)^2]^2}{8[(R^2)^2+(R^3)^2]^2}\left[1\frac{2}{3}[(R^2)^2+(R^3)^2]\pm \sqrt{1\frac{4}{3}[(R^2)^2+(R^3)^2]}\right].$$
(60)
Note that the parameters $`R^A`$ are related to the components of the inverse of the three-metric through Eqs. (47). Thus the parameter space $`\stackrel{~}{}`$ is really a submanifold of the corresponding minisuperspace. Fig. 1 shows a three-dimensional plot of $`R^4`$ as a function of $`R^2`$ and $`R^3`$, i.e., a plot of the parameter space $`\stackrel{~}{}`$ as embedded in $`\mathrm{IR}^3`$. Note that $`R^2=\pm R^3`$ renders $`f_1`$ and $`f_2`$ singular. The corresponding points which are depicted as the curves along which the figure becomes non-differentiable must be handled with care. The smooth part of $`\stackrel{~}{}`$ consists of eight connected components, each of which is diffeomorphic to an open disk (alternatively $`\mathrm{IR}^2`$).
* $`j=3/2`$: This case has been studied in Refs.  in detail. Therefore I suffice to note that it involves nontrivial geometric phases. Note that because of Kramer’s degeneracy, one does not need to restrict the minisuperspace to obtain degenerate eigenvalues. Every solution of the Bianchi type IX model involves a non-Abelian geometric phase.

## VI Conclusion

In this article I applied the method developed in Ref.  to investigate the existence of cyclic and noncyclic adiabatic geometric phases induced by spatially homogeneous cosmological backgrounds on a complex Klein-Gordon field. Unlike the examples presented in Ref. , here the freedom in the choice of the decomposition parameter $`q`$ turned out to simplify the analysis.

I showed that for the Bianchi type I models Berry’s connection one-form vanished identically. This was not the case for the Bianchi type IX models. Hence, for these models nontrivial non-Abelian adiabatic geometric phases could occur in general. A rather interesting observation was the relationship between the induced Hamiltonians in the Bianchi type IX models and the quadrupole Hamiltonians of the molecular and nuclear physics. I also showed that even for the integer spin representations nontrivial geometric phases could exist. This should also be of interest for the molecular physicists and chemists who have apparently investigated only the fermionic systems (half-integer spin representations.) A rather thorough investigation of the non-Abelian adiabatic geometric phase for arbitrary spin 1 systems has been conducted in .

As described in Ref.  the arbitrariness in the choice of $`q`$ leads to a $`GL(1,\text{ }\mathrm{C})`$ symmetry of the two-component formulation of the Klein-Gordon equation. In the context of general relativity where the Poincaré invariance is replaced by the diffeomorphism invariance, one can use the time-reparameterization symmetry of the background gravitational field and the geometric phase to absorb the magnitude $`|q|`$ of the decomposition parameter $`q`$ into the definition of the lapse function $`N=(g^{00})^{1/2}`$. In this way only a $`U(1)`$ subgroup of the corresponding $`GL(1,\text{ }\mathrm{C})`$ symmetry group survives. The $`GL(1,\text{ }\mathrm{C})`$ or $`U(1)`$ symmetry associated with the freedom of choice of the decomposition parameter seems to have no physical basis or consequences. It is merely a mathematical feature of the two-component formalism which can occasionally be used to simplify the calculations.

The application of the two-component formulation for the Bianchi models manifestly shows that this method can be employed even for the cases where the background spacetime is non-stationary. One must however realize that the present analysis is only valid within the framework of the relativistic adiabatic approximation . Although, the (approximate) stationarity of the background metric is a sufficient condition for the validity of the adiabatic approximation, it is not necessary. This can be easily seen by noting that for example in the case of Bianchi IX model, for spin $`j=1/2`$ states, one has $`\dot{\mathrm{\Phi }}_n=0`$, so $`\mathrm{\Phi }_m|\dot{\mathrm{\Phi }}_n=0`$. Therefore, although the spacetime is not stationary, the adiabatic approximation yields the exact solution of the field equation. This shows that in general for arbitrary non-stationary spacetimes, there may exist adiabatically evolving states to which the above analysis applies.

The extension of our results to the non-adiabatic cases requires a generalization of the analysis of Ref.  to non-adiabatic evolutions.

## Acknowledgments

I would like to thank Bahman Darian for many fruitful discussions and helping me with the computer graphics, Teoman Turgut for bringing to my attention one of the references, and Kamran Saririan for mailing me copies of a couple of the references.

## Appendix

In this Appendix I show how one can obtain the conditions under which the Hamiltonian (52) has degenerate eigenvalues without actually solving the eigenvalue problem in the general case.

The analysis can be slightly simplified if one writes the Hamiltonian (52) in the form:

$`\widehat{H}^{}`$ $`=`$ $`(t+2z)\widehat{I}+\widehat{\stackrel{~}{H}},`$
$`\widehat{\stackrel{~}{H}}`$ $`:=`$ $`\left(\begin{array}{ccc}0& \xi ^{}& \zeta ^{}\\ \xi & t& \xi ^{}\\ \zeta & \xi & 0\end{array}\right),`$ (64)
where $`\widehat{I}`$ is the $`3\times 3`$ identity matrix. Clearly, the eigenvalue problems for $`\widehat{H}^{}`$ and $`\widehat{\stackrel{~}{H}}`$ are equivalent. Computing the characteristic polynomial for $`\widehat{\stackrel{~}{H}}`$, i.e., $`P(\lambda ):=det(\widehat{\stackrel{~}{H}}\lambda \widehat{I})`$, one finds:

$$P(\lambda )=\lambda ^3+t\lambda ^2+(|\zeta |^2+2|\xi |^2)\lambda (t|\zeta |^2+\zeta \xi ^2+\zeta ^{}\xi ^2).$$
(65)
If one of the eigenvalues (roots of $`P(\lambda )`$) is degenerate, then

$$P(\lambda )=(\lambda l_1)(\lambda l_2)^2.$$
(66)
Comparing Eqs. (65) and (66), one finds

$$t=l_1+2l_2,l_2^2+2l_1l_2=(|\zeta |^2+2|\xi |^2),l_1l_2^2=(t|\zeta |^2+\zeta \xi ^2+\zeta ^8\xi ^2).$$
(67)
Furthermore since $`l_2`$ is at least doubly degenerate, the rows of the matrix:

$$\widehat{\stackrel{~}{H}}l_2\widehat{I}=\left(\begin{array}{ccc}l_2& \xi ^{}& \zeta ^{}\\ \xi & tl_2& \xi ^{}\\ \zeta & \xi & l_2\end{array}\right),$$
(68)
must be mutually linearly dependent. In other words the cofactors of all the matrix elements must vanish. Enforcing this condition for the matrix elements and using Eqs. (67), one finally finds that either $`\xi =\zeta =l_2=0`$ and $`l_1=t`$, or

$`\zeta `$ $`=`$ $`𝒵e^{2i\theta },t=𝒵{\displaystyle \frac{|\xi |^2}{𝒵}},`$
$`l_2`$ $`=`$ $`𝒵,l_1=(𝒵+{\displaystyle \frac{|\xi |^2}{𝒵}}),`$
where $`\mathrm{exp}[i\theta ]:=\xi /|\xi |`$ and $`𝒵\mathrm{IR}\{0\}`$.

Figure Caption:

Figure 1: This is a plot of $`R^4=R^4(R^2,R^3)`$. The horizontal plane is the $`R^2`$-$`R^3`$-plane and the vertical axis is the $`R^4`$-axis. The parameter space $`\stackrel{~}{}`$ is obtained by removing the intersection of this figure with the planes defined by: $`R^4=0`$, $`R^2=R^3`$ and $`R^2=R^3`$. The intersection involves the curves along which the figure becomes non-differentiable.

.