text:
# Perturbation Spectra in the Warm Inflationary Scenario

## I Introduction

One of the central paradigms in cosmology is the idea of inflation, and its elegant solution to the problems of standard cosmology. In addition inflation provides an explanation for the origin of structure from quantum fluctuations frozen in by crossing the Hubble radius during the rapid expansion. However in the standard inflationary scenario (see e.g. ref ) the rapid expansion phase and subsequent reheating of the universe back to a radiation dominated model are treated as two separate periods. The motivation for this was both simplicity and the belief that interactions with the field driving inflation would be minimal . The switching on of field interaction dynamics at the end of inflation results in violent re- and pre-heating effects such as a resonant amplification  which can lead to generic amplification of perturbations .

From the point of view of particle physics, interactions are inevitable, and they lead one to consider the effects on inflation of a phenomenological “friction” term describing the decay of the inflaton field into other fields during inflation . One result of this is that density perturbations are no longer generated via quantum fluctuations, but are dominated by the larger thermal fluctuations .

In this paper we consider the warm inflationary scenario  from the phenomenological point of view. Several types of phenomenological warm inflation models exist in the literature . Warm inflation dynamics also has been studied in quantum field theory  and some of the phenomenological models can be derived from first principles . Our purpose here is to determine generic observational features of the warm inflation scenario, without referring to specific models.

The paper is organised as follows. In Section II results are derive for evolution of the inflationary and radiation fields during the slow-roll phase of warm inflation. Assuming the dissipation rate dominates over the expansion rate, a relationship is found between the radiation energy density and the energy density of the inflationary field. The spectrum of adiabatic perturbations is calculated in Section IIIA and its spectral index in Section IIIB. Isocurvature perturbations are also generated during warm inflation and we present a new derivation of them in Section IIIC and of their spectral index in Section IIID. Tensor perturbations are discussed in Section IIIE. We compare our results with standard inflation in Section IIIF. The conditions for ending warm inflation are discussed in Section IV. Our results are illustrated by the specific cases of inflation with polynomial and power-law potentials in Section V. Our analytic expressions are tested against exact numerical models in Section VI, where we find good agreement between the two. In Section VII we discuss observational constraints on warm inflation and how it may be distinguished from the standard inflation scenario by future CMB experiments. Finally, our conclusions are presented in Section VIII. Two appendices are included in the paper. The first deals with a derivation of the thermal fluctuation of the inflationary field, while the second presents new results for generic polynomial and exponential potentials. Throughout this paper, the conventions of  are followed.

## II Evolution of the inflationary and radiation fields

The evolution of the inflaton field, $`\varphi `$, is given by

$$\ddot{\varphi }+(3H+\mathrm{\Gamma })\dot{\varphi }+V^{}=0,$$
(1)
where $`H\dot{a}/a`$ is the Hubble parameter, $`a`$ is the cosmological expansion factor, and $`\mathrm{\Gamma }`$ is the dissipation coefficient. $`V(\varphi )`$ is the potential of the inflationary field. We assume a spatially flat universe. A dot denotes differentiation with time and a prime indicates $`d/d\varphi `$. The energy–density of the radiation field, $`\rho _\gamma `$, obeys the transfer equation

$$\dot{\rho _\gamma }+4H\rho _\gamma =\mathrm{\Gamma }\dot{\varphi }^2,$$
(2)
The equations are completed by the Friedmann equation for the evolution of the universe;

$$H^2=\frac{8\pi }{3m_p^2}(V+\rho _\gamma ),$$
(3)
where $`m_p`$ is the Planck mass. During an inflationary era the potential field dominates both the kinetic energy of the inflationary field and the energy density of the radiation, so the Friedmann equation can be reduced to

$$H^2=\frac{8\pi }{3m_p^2}V.$$
(4)
We assume the slow–roll condition, $`\dot{\varphi }^2V(\varphi )`$, throughout the warm inflationary era, yielding

$$\dot{\varphi }=\frac{V^{}}{3H(1+r)},$$
(5)
where the ratio of the production rate of radiation to the expansion rate, $`r`$, is defined as

$$r\frac{\mathrm{\Gamma }}{3H},$$
(6)
and parameterises the importance of dissipation in the warm inflationary model.

For the dissipation–dominated regime of warm inflation, the dissipation rate, $`\mathrm{\Gamma }`$, is much greater then the expansion rate; $`r1`$. The evolution of the inflaton field during this phase is governed by

$$\dot{\varphi }=\frac{V^{}}{3rH}.$$
(7)
Hence the slow-roll condition can be met even for relatively steep potentials, due to the damping effect of dissipation on the evolution of $`\varphi `$.

If we also assume that during warm inflation radiation production is quasi-stable ($`\dot{\rho _\gamma }\mathrm{\Gamma }\dot{\varphi }^2`$), equations (2) and (7) give

$$\rho _\gamma \alpha T^4=3r\left(\frac{1}{2}\dot{\varphi }\right)^2,$$
(8)
where $`\alpha g_{}\pi ^2/30`$ is the Stefan-Boltzmann constant and $`g_{}`$ is the number of degrees of freedom for the radiation field. In the standard cosmology $`g_{}100`$. Equation (8) shows that the energy density in the radiation field is controlled by the kinetic energy of the inflationary field as well as the dimensionless dissipation rate. As the kinetic energy of $`\varphi `$ is small during the slow-roll regime we require $`r`$ to be large to generate substantial radiation energy-density. Equation (8) is not strictly true at the start and end of the warm inflationary phase, as the stable radiation production condition will be violated.

In contrast the energy density of the inflaton during the slow-roll phase is dominated by its potential energy, giving

$$\rho _\varphi =V(\varphi ).$$
(9)
The radiation density can be re-written as

$$\rho _\gamma =\frac{1}{12r}\left(\frac{V^{}}{H}\right)^2,$$
(10)
during warm inflation, since in the slow-roll phase the kinetic energy of the inflaton is controlled by the slope of the inflationary potential.

Introducing the dimensionless inflationary slow-roll parameter

$$\epsilon \frac{m_p^2}{16\pi }\left(\frac{V^{}}{V}\right)^2$$
(11)
we can relate the energy density of the radiation field to the energy-density of the inflaton by

$$\rho _\gamma =\left(\frac{\epsilon }{2r}\right)\rho _\varphi .$$
(12)
We have verified by numerical modelling (Section VI) that equation (12) holds very well during most of the warm inflationary phase.

We assume that warm inflation begins from a radiation dominated era, and starts when the energy density in the inflationary field is equal to that in the radiation field. Initially the radiation energy is redshifted away by the quasi-exponential expansion but enters a quasi-static phase when radiation production stabilises during the dissipation–dominated phase. Warm inflation ends when the universe heats up to become radiation dominated again, at a time when

$$\rho _\gamma =\rho _\varphi .$$
(13)
At this epoch the universe stops inflating and smoothly enters a radiation dominated big-bang phase. This is one of the most attractive features of warm inflation, and suggests that its predictions are robust in the sense that they are not altered during the smooth transition to the radiation-dominated epoch. The material components of the universe are created by the decay of either the remaining inflationary field or the dominant radiation field.

The approximate relation, equation (12), and equation (13) suggest that an approximate condition for warm inflation to both begin and end is

$$\epsilon =2r,$$
(14)
although we do not expect this relation to be accurate since the condition of stable radiation production is violated at both epochs. In particular at the beginning of warm inflation the ambient radiation field is redshifted away during the sharp transition to warm inflation and the radiation density drops rapidly. Then at the end of the warm inflation phase, the production of radiation must end and we expect $`\dot{\rho }_\gamma 3H\rho _\gamma \mathrm{\Gamma }\dot{\varphi }^2`$ during the turnover regime going from inflation-domination to radiation-domination.

In addition at both the beginning and end of warm inflation as the energy density of radiation is comparable to that of the inflationary potential, equation (3) should be used to describe the evolution of the universe rather than equation (4), introducing another approximation.

In Section VI we test the accuracy of our assumptions over the whole course of a warm inflationary phase by numerically integrating equations (1), (2) and (3) and find that for many models it appears that these approximations are accurate enough to be used to estimate the conditions over a large period of warm inflation, including the final few e-folds at the end of warm inflation.

The condition for a warm inflationary phase can be summarised by

$$\epsilon <2r.$$
(15)
It is straightforward to show that the same result can be derived from the requirement of accelerated expansion in the warm inflationary epoch. This requires that $`\ddot{a}>0`$ implying $`\dot{H}/H^2<1`$ or $`\dot{\varphi }H^{}(\varphi )/H^2<1`$. Using equations (4) and (7) we again find the inequality of equation (15). The slow-roll condition and the additional strong dissipation condition

$$r>1$$
(16)
then implies warm inflation. One immediate consequence is that the usual inflationary condition, $`\epsilon <1`$, is no longer required as inflation (quasi-exponential expansion) can be maintained for steep potentials due to the damping of dissipation.

Finally, one can differentiate equation (7) with respect to time and, following from the slow-roll condition $`\ddot{\varphi }V^{}`$, we find a constraint on the second derivative of the potential,

$$|\eta |3r^2$$
(17)
where

$$\eta \frac{m_p^2}{8\pi }\frac{V^{\prime \prime }}{V}$$
(18)
is the curvature slow-roll parameter for inflation. This should be regarded as an extra condition for sustained warm inflation, and is a generalisation of the usual inflationary constraint, $`|\eta |1`$. Hence warm inflation can take place for a wider range of potentials than standard supercooled inflation.

## III The perturbation spectra

### A Adiabatic perturbations

Density perturbations are calculated in the same way as for standard inflation , via the equation

$$\delta _H=\frac{2}{5}\left(\frac{H}{\dot{\varphi }}\right)\delta \varphi ,$$
(19)
where $`\dot{\varphi }`$ is given by equation (7). This equation holds for a mixture of inflationary and radiation fields, in the quasi-static regime of radiation production. It also can be shown that this result holds for any scalar field interacting adiabatically with another field, in the quasi-static limit, with the equation of state $`p=\omega \rho `$. In the presence of non-adiabatic perturbations (see Section III C) this relation does not necessarily hold. Specifically, if there non-adiabatic pressure perturbations, these can induce curvature perturbations to evolve on super-horizon scales . However, we expect these effects to be small for the scenario discussed here and equation (19) can still be used for the order-of-magnitude estimates presented here.

The amplitude of perturbations, $`\delta _H`$, is related to the final dimensionless linear power spectrum of mass perturbations by

$$\mathrm{\Delta }^2(k)=\left(\frac{k}{aH}\right)^4T^2(k)\delta _H^2(k)$$
(20)
where $`T(k)`$ is the matter transfer function relating initial spectra to the final, post-recombination linear spectra.

The fluctuations in the inflationary field are no longer generated by quantum fluctuations, but by thermal interactions with the radiation field;

$$\delta \varphi ^2=\frac{k_FT}{2\pi ^2},$$
(21)
where $`k_F=\sqrt{\mathrm{\Gamma }H}=\sqrt{3r}HH`$ is the freeze-out scale at which dissipation damps out the thermally excited fluctuations. The origin of the freeze-out wavenumber $`k_F`$ is explained in Appendix A based on the original derivation in . Note that in the warm inflationary epoch, the freeze-out length scale is generally much smaller than the Hubble radius. As one final point, when $`\mathrm{\Gamma }H`$ and $`TH`$, the amplitude of thermal fluctuations is equal to the amplitude of quantum fluctuations generated by Hawking radiation.

Combining equations (4), (7), (19) and (21) we find

$$\delta _H^2=\frac{32}{75m_p^4}\frac{V}{\epsilon }\left[r^{5/2}\left(\frac{\sqrt{12}T}{H}\right)\right].$$
(22)
The temperature of the radiation field is related to the potential by equations (8), (9) and (12) yielding

$$T^4=\frac{\epsilon V}{2r\alpha }.$$
(23)
Expressing the amplitude of density perturbations in terms of inflationary parameters we find

$$\delta _H^2=\frac{16}{25\sqrt{\pi }m_p^3}\left(\frac{2}{\alpha }\left(\frac{r^3V}{\epsilon }\right)^3\right)^{1/4}.$$
(24)
The ratio of thermal to quantum fluctuations is $`r^{5/2}T/H`$. During strong dissipation, $`r1`$, while towards the end of the warm inflationary era $`T/Hm_p/V^{1/4}1`$. Hence the thermal spectrum will dominate over the quantum spectrum of fluctuations during a warm inflationary phase.

### B Spectral index of adiabatic perturbations

The spectral index of the adiabatic scalar perturbations is given by

$$n_s1=\frac{d\mathrm{ln}\delta _H^2}{d\mathrm{ln}k}$$
(25)
where the logarithmic interval in wavenumber is related to the number of e-folds from the start of inflation, $`N(\varphi )`$, and the scalar field by

$$d\mathrm{ln}k(\varphi )=dN(\varphi )=\frac{8\pi }{m_p^2}\frac{rV}{V^{}}d\varphi .$$
(26)
Hence the spectral index of adiabatic density perturbations in warm inflation is given by

$$n_s=1\frac{3}{4r}(3\epsilon 2\eta ).$$
(27)
As we require that $`r1`$ we expect the deviation from scale invariance to be suppressed in the warm inflationary model, until $`\epsilon r`$, at the end of the warm inflation phase. If $`\epsilon >(2/3)\eta `$, or $`\eta `$ is negative, the spectrum of adiabatic fluctuations will be red, with the spectrum suppressed on small scales, as in standard supercooled inflation which generally also produces slopes $`n_s1`$. However, in the opposite case, when $`\epsilon <(2/3)\eta `$, the spectrum of adiabatic fluctuations will be blue, in contrast with predictions from generic standard inflation models. With a blue spectrum, the interesting possibility exists of producing primordial black holes on small scales .

### C Isocurvature perturbations

Isocurvature perturbations also are generated in the warm inflationary era due to thermal fluctuations in the radiation field. These perturbations can be characterised by fluctuations in the entropy, $`S_X`$, where the species $`X`$ is undergoing thermal fluctuations relative to the number density of photons

$$S_X=\frac{3}{4}(\delta _X\delta _\gamma )$$
(28)
and we have assumed the $`X`$-field is the ambient radiation field excluding photons.

Assuming that the isocurvature fluctuations are thermal in origin, their perturbation spectrum is

$`\delta _{\mathrm{iso}}^2`$ $`=`$ $`{\displaystyle \frac{1}{25m_p^2}}{\displaystyle \frac{k_FT}{2\pi ^2}}`$ (29)
$`=`$ $`{\displaystyle \frac{1}{25\sqrt{\pi ^3}}}\left({\displaystyle \frac{2r\epsilon }{\alpha }}\right)^{1/4}\left({\displaystyle \frac{V^{1/4}}{m_p}}\right)^3.`$ (30)
This relation is true only for non-interacting fields, but again should hold for order-of-magnitude estimates. The ratio of isocurvature to adiabatic perturbations is

$$R_{\mathrm{iso}}\frac{\delta _{\mathrm{iso}}^2}{\delta _H^2}=\frac{\epsilon }{16\pi r^2}$$
(31)
so that when $`r1`$ and $`r\epsilon `$ adiabatic perturbations will dominate. In practice the scalar modes will form a quadratic sum of adiabatic and isocurvature modes,

$$\delta _s^2=\delta _H^2+\delta _{\mathrm{iso}}^2,$$
(32)
and which can be separated observationally since isocurvature generated fluctuations of the CMB are out of phase with adiabatic perturbations. Another discriminator is if super-horizon perturbations can be observed, for instance in the CMB . Since isocurvature modes do not contribute to curvature perturbations (in the comoving gauge) the contribution from isocurvature modes should be absent from super-horizon perturbations, leading to a drop in amplitude of the scalar perturbations at the horizon scale. Isocurvature perturbations have previously been discussed in the context of warm inflation by , although no detailed calculation of their amplitude was given.

### D Spectral index of isocurvature perturbations

Following the derivation of the spectral index of adiabatic perturbation, the spectral index of the isocurvature perturbations, defined by

$$n_{\mathrm{iso}}\frac{d\mathrm{ln}\delta _{\mathrm{iso}}^2}{d\mathrm{ln}k},$$
(33)
is

$$n_{\mathrm{iso}}=\frac{1}{4r}(\epsilon +2\eta ).$$
(34)
This differs from isocurvature fluctuations generated by quantum fluctuations, where $`n_{\mathrm{iso}}=2\epsilon `$. If $`r`$ is large, the spectrum will again tend to be flat, until $`\epsilon r`$. If $`\epsilon >2\eta `$ the thermal isocurvature spectrum will be red as it is in standard supercooled inflation.

### E Tensor perturbations

Tensor perturbations do not couple strongly to the thermal background and so gravitational waves are only generated by quantum fluctuations, as in standard supercooled inflation. Their spectrum is given by

$$A_\mathrm{g}^2=\frac{32V}{75m_p^4},$$
(35)
with spectral index

$$n_\mathrm{g}\frac{d\mathrm{ln}A_g^2}{d\mathrm{ln}k}=\frac{2\epsilon }{r},$$
(36)
and the inequality $`4<n_g<0`$ arising from the warm inflation condition $`\epsilon r`$.

The ratio of tensor to scalar perturbations is

$$R_g\frac{A_\mathrm{g}^2}{\delta _H^2}=\frac{\sqrt{\pi }}{3}\left(\frac{2\epsilon }{r^3}\right)^{3/4}\left(\frac{\alpha V}{m_p^4}\right)^{1/4}.$$
(37)
Hence, as we are in the strongly dissipative regime, where $`r1`$ and $`\epsilon r`$, we expect the contribution from gravitational waves to be small unless the magnitude of the potential is large or $`r`$ is order unity.

The ratio of the amplitude of isocurvature to tensor perturbations scales as

$$\frac{\delta _{\mathrm{iso}}^2}{A_g^2}=\frac{3m_p}{32\sqrt{\pi ^3}}\left(\frac{2r\epsilon }{\alpha V}\right)^{1/4},$$
(38)
which will be greater than unity unless the scale of the potential is large. Hence we can expect that the inequalities

$$R_g<R_{\mathrm{iso}}<1$$
(39)
in general will hold. We find that this is true for all the of the numerical models we consider in Sections V and VI.

### F Comparison with Standard Supercooled Inflation

A nice feature of warm inflation is that these results are significantly different from standard supercooled inflation case, where

$`\delta _H^2`$ $`=`$ $`{\displaystyle \frac{32V}{75m_p^4\epsilon }},`$ (40)
$`n_s`$ $`=`$ $`16\epsilon +2\eta ,`$ (41)
$`n_g`$ $`=`$ $`2\epsilon ,`$ (42)
$`R_g`$ $`=`$ $`\epsilon .`$ (43)
Thus near future observations should be able to distinguish between the two. An immediate consequence of our analysis is that there need be no consistency relationship between the scalar spectral index and the ratio of tensor to scalar perturbations, as in standard inflation. We discuss the observational differences between the models in Section VII.

## IV Conditions at the end of inflation

At the end of the warm inflationary epoch, the dissipation parameter, $`r`$, is related to the slow-roll parameter, $`\epsilon `$, by equation (14). Substituting in equations (6), (4) and (11) we form a differential equation which can be solved for $`V(\varphi )`$ to find the magnitude of the potential field at the end of warm inflation,

$$V(\varphi _{\mathrm{end}})=\frac{\pi }{6}\left(\frac{\mathrm{\Gamma }}{m_p}\right)^2\varphi _{\mathrm{end}}^4,$$
(44)
up to a constant, but independent of the shape of the inflationary potential. This selects out the $`V\varphi ^4`$ potential as a special case where warm inflation never ends (see Appendix B). This is due to the energy density of the radiation field “tracking” the energy density in the inflationary field, $`\rho _\gamma \rho _\varphi `$, and so the condition for warm inflation to end is never met. In this case warm inflation continues until the slow-roll conditions are no longer met and the inflaton field decays away, which then ends inflation. For $`n>4`$ and exponential potentials, the radiation density will always fall faster than the energy-density of the inflaton field (see Appendix B). In this case warm inflation will occur, in the sense that a sizable radiation component can persist for a considerable duration during inflation, well beyond the requirements for solving the problems of standard cosmology. However to end such a warm inflation regime, this behavior of the potential must be local, and at some point the vacuum energy must be released at a sufficiently fast rate to allow the radiation to dominate.

## V Specific Cases

In order to understand how warm inflation differs from the standard supercooled inflationary scenario, it is useful to consider a few specific examples. Here we consider models of warm inflation with polynomial and exponential potentials.

### A Polynomial inflation

We consider the simplest inflationary model

$$V(\varphi )=m^2\varphi ^2/2$$
(45)
for a massive scalar field with mass $`m`$. With this potential the slow-roll parameters are

$$\epsilon =\eta =\frac{m_p^2}{4\pi \varphi ^2}.$$
(46)
The Hubble parameter is given by

$$H(\varphi )=\sqrt{\frac{4\pi }{3}}\left(\frac{m}{m_p}\right)\varphi $$
(47)
while the dissipation parameter scales like

$$r=\frac{1}{\sqrt{12\pi }}\left(\frac{\mathrm{\Gamma }}{m}\right)\frac{m_p}{\varphi }.$$
(48)
From equation (48) we see that the ratio $`\mathrm{\Gamma }/m`$ controls the evolution of the system.

Solving equation (7) we find that the inflaton field evolves as

$`\varphi (t)`$ $`=`$ $`\varphi _0e^{m^2t/\mathrm{\Gamma }}`$ (49)
$``$ $`\varphi _0\left(1{\displaystyle \frac{m^2}{\mathrm{\Gamma }}}t\right).`$ (50)
Compared with standard inflation, where

$$\varphi (t)\varphi _0\left(\frac{mm_p}{\sqrt{12\pi }}\right)t,$$
(51)
we see that during warm inflation the inflaton field decays due to dissipation into the radiation field.

The energy densities of the inflationary potential is given by equation (9) and (45) while the radiation field is given by

$$\rho _\gamma =\sqrt{\frac{3}{\pi }}\left(\frac{m^3m_p}{8\mathrm{\Gamma }}\right)\varphi =\sqrt{\frac{6}{\pi }}\left(\frac{m^2m_p}{8\mathrm{\Gamma }}\right)\rho _\varphi ^{1/2}.$$
(52)
The number of e-folds from the end of warm inflation is given by

$`N(\varphi )`$ $`=`$ $`\mathrm{ln}k(\varphi )`$ (53)
$`=`$ $`{\displaystyle \frac{8\pi }{m_p^2}}{\displaystyle _\varphi ^{\varphi _{\mathrm{end}}}}𝑑\varphi ^{}{\displaystyle \frac{rV}{V^{}}}`$ (54)
$`=`$ $`\sqrt{{\displaystyle \frac{4\pi }{3}}}\left({\displaystyle \frac{\mathrm{\Gamma }}{m}}\right){\displaystyle \frac{1}{m_p}}(\varphi \varphi _{\mathrm{end}}).`$ (55)
Figure 1 shows the evolution of the energy densities of the radiation and inflationary fields as a function of number of e-folds. The upper line is the inflation energy density and the lower line the radiation energy density. The dashed lines are the results of our analytic model, from equations (45), (52) and (57), while the solid lines are the results of our numerical model (Section VI) of the evolution of the field through the warm inflationary era.

From the condition that $`\epsilon =2r`$ at the end of warm inflation we find that the magnitude of the inflationary field at this time is

$$\varphi _{\mathrm{end}}=\frac{1}{2}\sqrt{\frac{3}{4\pi }}m_p\left(\frac{m}{\mathrm{\Gamma }}\right)$$
(56)
where the initial field must have $`|\varphi _0|>\varphi _{\mathrm{end}}`$. In standard polynomial $`(n=2)`$ inflation this condition is $`|\varphi |_{\mathrm{end}}>m_p/\sqrt{4\pi }`$, so we find that for warm inflation this is weakened by the factor $`m/\mathrm{\Gamma }`$. If $`\mathrm{\Gamma }>m`$ then warm inflation can begin well below the Planck mass.

The total number of e-folds from the end of warm inflation is

$$N_{\mathrm{tot}}=\sqrt{\frac{4\pi }{3}}\left(\frac{\mathrm{\Gamma }}{m}\right)\left(\frac{\varphi _0}{m_p}\right)\frac{1}{2}.$$
(57)
As we require $`N70`$ e-folds to comfortably solve the problems of standard cosmology this implies

$$\varphi _035\left(\frac{m}{\mathrm{\Gamma }}\right)m_p.$$
(58)
The requirement of strong dissipation also leads to a constraint on $`\mathrm{\Gamma }`$. From equation (48) with equation (58) we find that at the beginning of inflation

$$r(\varphi _0)\frac{1}{215}\left(\frac{\mathrm{\Gamma }}{m}\right)^21$$
(59)
resulting in the requirement that $`\mathrm{\Gamma }15m`$ for warm inflation to commence in polynomial ($`n=2`$) inflation.

Calculating the spectrum of adiabatic perturbations yields

$$\delta _H^2=\frac{2}{75}\left(\frac{2}{3\alpha \pi ^5}\right)^{1/4}\left(\frac{\mathrm{\Gamma }}{m_p}\right)^{3/2}\left(\mathrm{ln}\frac{k}{k_{\mathrm{end}}}\right)^{3/4},$$
(60)
where $`\mathrm{ln}k_{\mathrm{end}}=\sqrt{4\pi /3}\mathrm{\Gamma }\varphi _{\mathrm{end}}/(mm_p)1/2`$ is the fiducial wavenumber at the end of warm inflation. Note that the amplitude of adiabatic perturbations is independent of the magnitude of the inflationary potential. In Appendix B we show that the amplitude of adiabatic perturbations in warm inflation is independent of the magnitude of the inflationary potential for all polynomial potentials.

Figure 2 shows the amplitude of adiabatic, isocurvature and tensor perturbations for the $`n=2`$ polynomial. In Appendix B we derive results for the more general polynomial and exponential cases. We find that the ratio of the amplitude of perturbations agrees with the inequality expressed by equation (39). For the $`n=2`$ case the isocurvature and tensor perturbations are strongly suppressed below the adiabatic perturbations. This difference can be reduced by increasing the magnitude of the inflationary potential.

The spectral indices are related by

$$n_s1=n_{\mathrm{iso}}=\frac{3}{8}n_g=\frac{3\epsilon }{4r},$$
(61)
and are all red. In Appendix B we derive more results for the general case of $`\varphi ^n`$ potentials. There we show that for polynomial potentials the spectral index is always red unless $`n>4`$, or for the case of an exponential potential (see Appendix B), in which case the spectrum is blue. The latter open up the interesting possibility of primordial black hole production in the early universe . A scale invariant spectrum, $`n_s=1`$, is found only for the $`n=4`$ model.

The dissipation ratio, $`r`$, and the slow-roll parameter $`\epsilon `$ are related by

$$r=\frac{1}{\sqrt{3}}\left(\frac{\mathrm{\Gamma }}{m}\right)\epsilon ^{1/2}$$
(62)
so that the ratio of isocurvature to adiabatic perturbations, given by $`\epsilon /r^2`$, is $`3(m/\mathrm{\Gamma })^21`$.

Figure 3 shows the spectral indices of the adiabatic, isocurvature and tensor spectra as a function of wavenumber. In all cases the spectral indices starts as a flat spectrum, when $`\epsilon r`$. In the case of $`n=4`$ spectral indices are exactly flat while for the exponential potential the spectral indices are almost flat.

## VI Numerical Models

While our analytic results simplify the whole analysis of warm inflationary models, it is important to verify the range of validity of this analysis and find where the assumptions used in their derivation are accurate. To test this we have numerically evolved the inflationary and radiation fields through a warm inflationary epoch, and calculated the adiabatic and tensor perturbation spectra. We numerically solve equations (1) and (2) for the evolution of the field and equation (3) for the evolution of the background universe using the methods of .

The results for the evolution of the energy-density of the fields for the case of a polynomial $`n=2`$ potential are shown in Figure 1. The magnitude of the inflationary and radiation energy-density fields are both normalised to unity in units of the potential at the start of inflation. We set $`\rho _\gamma =\rho _\varphi `$ initially, as it is reasonable to assume that inflation begins from a radiation-dominated epoch. We tested a number of other initial conditions including $`\rho _\gamma =0`$ and found little variation in our final results, as the models soon settled down to the same evolution after the first few e-folds.

To parameterise the models we used the more general expression for polynomial inflation from Appendix B, and choose parameters so that the amplitude of adiabatic perturbations $`\delta _H=2\times 10^5`$ at the $`50^{th}`$ e-fold, and allowed a total of 70 e-folds of warm inflation to elapse.

Figure 1 verifies that not only the shapes of the curves are well approximated, but also the relative magnitudes of the radiation and inflationary energy-densities are accurately reproduced. The main assumptions underlying our analytic results are that of slow-roll inflation and quasi-stable radiation production ($`\dot{\rho }_\gamma \mathrm{\Gamma }\dot{\varphi }^2`$). The latter assumption breaks down at the start of warm inflation as the radiation field initially is redshifted away. But radiation production by dissipation of the inflaton field soon catches up with the redshifting and stability is reached after a few e-folds.

We found accurate agreement between our analytic models and the full numerical models for polynomial potentials in the range $`0.5<n<4`$. In most cases the agreement was $`<1\%`$ over most of the range, except near the beginning of inflation. In the $`n=4`$ case the energy densities tracked each other as expected from our discussion in Section IV. We also tested our results for the exponential potential discussed in Appendix B, and again found similar accuracy, although here again the radiation field drops faster than the energy density in the inflationary field and warm inflation has to end by the natural decay of the inflaton, or by some other means.

We have also computed the amplitude and scale dependence of the perturbation spectra for adiabatic perturbations directly from equations (19) and (21). These are plotted for the case of $`n=2`$ in Figure 2. Again close agreement was found, with deviations less than $`1\%`$ over the range of wavenumbers. Again the poorest agreement was over the few e-folds at the beginning and end of warm inflation.

Finally we tested the predicted spectral indices for adiabatic perturbations (Fig. 3). Again we found good agreement between the analytic and numerical models, except for the first and last few e-folds. The sudden change in $`n_s1`$ at the beginning of warm inflation ($`N70`$) is real and due to the sharp transition to warm inflation when the change in the radiation density, $`\dot{\rho }_\gamma `$, is important.

## VII Observational Constraints

In standard inflation, there are four observable constraints ($`\delta _H`$, $`A_g`$, $`n_s`$, $`n_g`$) for the three parameters of slow-roll inflation, $`V_0`$, $`\epsilon `$, $`\eta `$, where $`V_0`$ is the magnitude of the inflationary potential. With four constraints there is a redundancy in the observations allowing a consistency relation . This is usually expressed as a relationship between the tensor-to-scalar ratio and the slope of the tensor spectrum,

$$R_g=\frac{n_g}{2}.$$
(63)
This appears to be generally true for slow-roll inflation, although recent work on pre-heating at the end of supercooled inflation may effect this result .

In warm inflation the dissipation factor, $`r`$, becomes an active parameter (since it is present at the end of supercooled inflation, but assumed inactive) giving four constraints for four parameters. Hence we do not expect the consistency relation of standard inflation to hold in warm inflation. This relation can only arise by coincidence. If the consistency relation is not observed it may indicate a warm inflationary phase. The dissipation factor can be estimated from the observables by the relation

$$r^2n_g\left(\frac{25\alpha \pi ^2}{32}R_g^4A_g^2\right)^{1/3}.$$
(64)
To discriminate between warm and standard inflation requires a measurement of all four observables. The MAP and Planck satellite missions will probe the CMB and provide strong constraints on the amplitude and shape of the scalar perturbation spectrum. Both MAP and Planck have polarisation detectors which are sensitive to tensor modes, but Planck has a much broader range of frequencies which will aid the subtraction of foreground contamination of the CMB signal. In practice the tensor modes may be so small as to be unobservable, hence it may not be easy to distinguish between the two models.

Current observations, from the COBE satellite, and subsequent airborne detectors (see  for a recent analysis), can be used to place constraints on the amplitude of scalar perturbations and the slope of their spectra. The low COBE multipoles indicate that $`\delta _H^210^{10}`$, giving the constraint

$`V^{1/4}`$ $``$ $`7.6\times 10^4(\alpha ^{1/3}\epsilon r^3)^{1/4}m_p,`$ (65)
$``$ $`7.6\times 10^{15}(\alpha ^{1/3}\epsilon r^3)^{1/4}\mathrm{GeV}.`$ (66)
If $`r^3\alpha ^{1/3}\epsilon `$, these energy scales are pretty much in keeping with those of Grand Unified Theories. However if $`r1`$ as required for warm inflation we open up the interesting possibility that the energy scale required can be much lower.

## VIII Conclusions

In this paper we have investigated the phenomenology of the warm inflationary scenario. In the slow-roll, warm inflation epoch we have found the general relationship between the energy-density in the inflaton and radiation fields. This has led us to the general criteria for warm inflation, $`\epsilon <r`$, $`|\eta |<3r^2`$ and $`r>1`$. The usual inflation condition, $`\epsilon <1`$, is generalised as dissipation can help sustain accelerated expansion of the universe for a wider range of inflationary potentials.

Generalising the analysis of inflation in the slow-roll regime, we have derived expressions for the perturbation spectra of adiabatic and isocurvature fluctuations generated by thermal fluctuations during the warm inflationary epoch. The adiabatic spectrum produced by thermal fluctuations dominates that produced by quantum fluctuation. Its amplitude and spectral slope are both different from quantum generated perturbations and can lead to red or blue spectra. Isocurvature perturbations also are generated due to internal fluctuation between radiation fields. These are generally suppressed compared to the adiabatic spectrum.

Gravitational waves are still generated via quantum excitation, and have the same amplitude but different slope from standard inflation. The tensor-to-scalar ratio is suppressed by the relatively high amplitude of the adiabatic spectrum.

No consistency relation exists in warm inflation as there are now four parameters (the standard three of inflation and an active dissipation factor) and in general only four observables. This is probably the most important observational test of warm inflation, but may be difficult, even for Planck, if the tensor modes are strongly suppressed. In this case measurement of the isocurvature modes relative to adiabatic would provide an alternative test. However, verification of the consistency relation would rule out substantial interactions during inflation.

There are two important differences between standard and warm inflation. First, the slope of the inflationary potential does not need to be small, as dissipation will allow slow-roll even for steep potentials. Second, the energy scale of warm inflation may be much lower than required for standard inflation, due to the higher amplitude of thermal adiabatic perturbations and the constraints of Sachs-Wolfe effects observed in the CMB.

To investigate the warm inflationary model we need to study further the evolution of scalar perturbations during the transition phase from warm inflation to radiation domination when the energy densities of the radiation field and inflationary field are comparable. This will occur more smoothly than in standard inflation and should lead to observational differences between the models. In addition the model we have presented here for isocurvature perturbations is crude. A more complete understanding of the generic prediction of isocurvature perturbations for warm inflation is needed before detailed predictions can be made.

The model we have presented here is an extreme one, in that we have assumed that it is dissipation–dominated. A whole range of models lie between that of warm inflation and the standard supercooled inflation, and it would be interesting to explore their properties in the general framework we have outlined here. The development in  also would be beneficial for this.

Finally these results should be propagated to the post-recombination era using general-purpose transfer codes such as CMBfast . However to make definite, high-quality predictions for the CMB, such codes must be generalised to remove any of the relations of standard supercooled inflation, such as the consistency relation. This will allow us to predict with some confidence whether the MAP and Planck missions will be able to probe the high redshift universe with sufficient accuracy to distinguish between warm and supercooled inflation.

Acknowledgements: ANT and AB acknowledge the support of PPARC, and thank David Wands for useful discussion about the evolution of non-adiabatic super-horizon perturbations.

## Appendix A

The origin of the “freeze-out” wavenumber $`k_F`$ in equation (21) is reviewed here following the original and more detailed explanation in . To account for fluctuations, the inflaton now is written as $`\varphi (𝒙,t)=\varphi _0(t)+\delta \varphi (𝒙,t)`$. Assuming the universe remains near thermal equilibrium during warm inflation, the fluctuations of the inflaton field obey a fluctuation-dissipation relation. The evolution equation for the fluctuations, $`\delta \varphi (𝒙,t)`$, is obtained from the linearized deviation of the zero mode equation of motion equation (1), with inclusion of the spatial Laplacian term and the addition of a white–noise random force term

$$\mathrm{\Gamma }\frac{d\delta \varphi (𝒌,t)}{dt}=[𝐤^2+V^{\prime \prime }(\varphi _0)]\delta \varphi (𝒌,t)+\xi (𝒌,t),$$
(67)
where $`V^{\prime \prime }(\varphi _0)=d^2V(\varphi _0)/d\varphi _0^2`$ and we have Fourier transformed to momentum space. The white–noise represents the action of the thermal heat bath on the inflaton field, with the spectrum of the this noise term governed by the fluctuation-dissipation theorem,

$$\xi (𝒌,t)=0$$
(68)
and

$$\xi (𝒌,t)\xi (𝒌^{},t^{})_\xi \stackrel{T\mathrm{}}{=}2\mathrm{\Gamma }T(2\pi )^3\delta ^{(3)}(𝒌𝒌^{})\delta (tt^{}).$$
(69)
The solution of equation (67) is

$`\delta \varphi (𝒌,t)`$ $``$ $`{\displaystyle \frac{1}{\mathrm{\Gamma }}}e^{(tt_0)/\tau (\varphi )}{\displaystyle _{t_0}^t}e^{(t^{}t_0)/\tau (\varphi _0)}\xi (𝒌,t^{})𝑑t^{}`$ (70)
$`+`$ $`\delta \varphi (𝒌,t_0)e^{(tt_0)/\tau (\varphi _0)},`$ (71)
where $`\tau (\varphi )=\mathrm{\Gamma }/(k^2+V^{\prime \prime }(\varphi ))`$. The first term on the right hand side acts to thermalize $`\delta \varphi `$, whereas the latter is the memory term for the initial value of $`\delta \varphi `$, which over time becomes negligible.

In the cosmological case, equation (67) must be interpreted for physical wavenumbers, since the thermal effects as expressed in this equation act in accordance with the physical wavenumbers of the system. Thus for a given fluctuation mode, $`\delta \varphi (𝒌_c)`$, at comoving wavenumber $`𝒌_c`$, equation (67) expresses its dynamics within the time interval $`1/H`$ when the physical mode is $`k_{phys}k_c\mathrm{e}^{Ht}k`$. To be consistent with the thermalization conditions, during the time interval $`1/H`$ the mode $`\varphi (𝒌_c)`$ must thermalize at physical scale $`k`$. In terms of the solution equation (71), this requires the memory term to be negligible within a Hubble time, thus

$$\frac{k^2+V^{\prime \prime }(\varphi _0)}{H\mathrm{\Gamma }}>1.$$
(72)
The freeze-out wavenumber $`k_F`$ is at the point where this condition first holds, which for $`V^{\prime \prime }(\varphi _0)<\mathrm{\Gamma }H`$ is

$$k_F=\sqrt{\mathrm{\Gamma }H}.$$
(73)
If $`V^{\prime \prime }(\varphi _0)>H\mathrm{\Gamma }`$ then the freeze-out wavenumber is the same as for supercooled inflation, $`k_F=H`$. However, in general this regime never occurs during warm inflation and equation (73) is the governing condition. Note that $`H`$ varies during warm inflation, thus based on equation (73), so too will $`k_F`$. However, the variation of $`H`$ is typically very small, less than a factor 10, during the entire warm inflation period. Thus, up to order one factors, one can treat $`k_F`$ as a constant with $`H`$ evaluated at some appropriate time during warm inflation, such as the beginning. In our analytic approach we have included the $`\varphi `$-dependence of $`k_F`$.

In our treatment here we have taken the dissipative coefficient $`\mathrm{\Gamma }`$ to be independent of $`k`$, whereas in  it was mentioned that in general $`\mathrm{\Gamma }`$ will depend on $`k`$ and generally will decrease with increasing wavenumber. In this case, the condition equation (73) still is valid except it must now be solved for $`k_F`$ due to the $`k`$-dependence in $`\mathrm{\Gamma }`$.

## Appendix B

### Polynomial potentials

In this appendix we explore the properties of the warm inflationary scenario for the family of models with the functional form

$$V(\varphi )=\lambda m^4\left(\frac{\varphi }{m}\right)^n$$
(74)
where $`\lambda `$ is a dimensionless parameter, which allows more freedom in the normalisation of the potential.

The slow-roll parameters are given by

$$\epsilon =\frac{n^2m_p^2}{16\pi \varphi ^2},\eta =\frac{n(n1)m_p^2}{8\pi \varphi ^2}$$
(75)
with the relation

$$\epsilon =\frac{n}{2(n1)}\eta .$$
(76)
The expansion of the universe is governed by the Hubble parameter

$$H(\varphi )=\sqrt{\frac{8\pi \lambda }{3}}\frac{m^2}{m_p}\left(\frac{\varphi }{m}\right)^{n/2},$$
(77)
where the dimensionless dissipation parameter is

$$r=\frac{1}{\sqrt{24\pi \lambda }}\left(\frac{\mathrm{\Gamma }}{m}\right)\left(\frac{m_p}{m}\right)\left(\frac{\varphi }{m}\right)^{n/2}.$$
(78)
The energy density of the radiation field is given by

$$\rho _\gamma =\sqrt{\frac{3\lambda ^3}{2\pi }}\frac{n^2}{8}\frac{m_p}{\mathrm{\Gamma }}m^4\left(\frac{m}{\varphi }\right)^{(43n)/2},$$
(79)
and related to the energy-density of the inflation field by

$$\rho _\gamma =\sqrt{\frac{3}{2\pi }}\frac{n^2}{8}\frac{m_p}{\mathrm{\Gamma }}m_p^4\left[\lambda \left(\frac{m}{m_p}\right)^{4n}\right]^{2/n}\left(\frac{2\rho _\varphi }{m_p^4}\right)^{\frac{(3n4)}{2n}}.$$
(80)
This last relation tells us that for an $`n=4`$ potential the energy densities are proportional, $`\rho _\gamma \rho _\varphi `$, and that warm inflation cannot end by the condition given by equation (14). When the potential is steeper than this, $`n>4`$, or for the exponential potential, the radiation field drops more rapidly than the inflation field. Although warm inflation occurs, the effects of the radiation field will diminish and inflation will begin to evolve in the same way as for standard inflation. Hence for polynomial inflationary potentials we require $`n<4`$ for strong warm inflation to take hold.

The number of e-folds, $`N(\varphi )`$, the universe expands during inflation is

$$N(\varphi )=\sqrt{\frac{2\pi }{3\lambda }}\frac{8\mathrm{\Gamma }}{n(4n)m_p}\left[\left(\frac{\varphi _0}{m}\right)^{\frac{(4n)}{2}}\left(\frac{\varphi }{m}\right)^{\frac{(4n)}{2}}\right]$$
(81)
where we count the number of e-folds from the start of warm inflation (note in Section V we counted e-folds from the end of warm inflation). The number of e-folds can be related to $`\epsilon `$ and $`r`$ by

$$N(\varphi )=\left(\frac{2n}{4n}\right)\left(\frac{r_0}{\epsilon _0}\frac{r}{\epsilon }\right).$$
(82)
If $`n<4`$ then $`N`$ is dominated by the ratio of $`\epsilon /r`$ at the start of the inflationary phase, and to a good approximation the total number of e-folds is

$$N_{\mathrm{tot}}=\left(\frac{2n}{4n}\right)\frac{r_0}{\epsilon _0}.$$
(83)
The adiabatic spectrum of density perturbations is

$$\delta _H^2=\frac{8}{75}\left(\frac{2^7}{3^5\pi ^7\lambda ^3}\right)^{1/8}\left(\frac{\mathrm{\Gamma }^9}{\alpha n^6m_p^9}\right)^{1/4}\left(\frac{\varphi }{m}\right)^{3(4n)/8}.$$
(84)
Substituting in equation (81) we find that the amplitude of adiabatic perturbations in warm inflation is independent of the magnitude of the inflationary potential for polynomial potentials, when expressed as a function of e-folds, or scale.

The spectral index of adiabatic fluctuations is related to the slow-roll parameter, $`\epsilon `$, and the dissipation parameter, $`r`$, by

$$n_s1=\frac{(4n)}{n}\frac{3\epsilon }{4r}=\frac{(4n)}{n}\frac{3}{8}n_g=\frac{3(4n)}{(5n4)}n_{\mathrm{iso}}.$$
(85)
where the last two expressions relate the adiabatic spectral index to the isocurvature and tensor spectral index. From this expression we see that the warm inflationary adiabatic spectrum is only exactly scale invariant ($`n_s=1`$) for potentials with $`n=4`$, and is red unless $`n>4`$.

The isocurvature spectrum is given by

$$\delta _{\mathrm{iso}}^2=\frac{1}{50}\left(\frac{\lambda ^5}{6^3\pi ^{15}}\right)^{1/8}\left(\frac{m}{m_p}\right)^2\left[\frac{n^2}{\alpha }\frac{\mathrm{\Gamma }}{m_p}\left(\frac{\varphi }{m}\right)^{(5n4)/2}\right]^{1/4}.$$
(86)
The ratio of the amplitude of isocurvature to adiabatic perturbations is given by

$$R_I=\frac{\epsilon }{16\pi r^2}=\frac{3\lambda n^2}{32\pi }\left(\frac{m}{\mathrm{\Gamma }}\right)^2\left(\frac{\varphi }{m}\right)^{n2}.$$
(87)
In the special case of $`n=2`$ this is a constant, as found in Section V.

### Exponential potential

Consider potentials of the form

$$V(\varphi )=V_0\mathrm{exp}\frac{\varphi }{\kappa },$$
(88)
where $`\kappa `$ is a free parameter with units of mass.

The slow-roll parameters are

$$\eta =2\epsilon =\frac{m_p^2}{8\pi \kappa ^2}.$$
(89)
The energy densities of the radiation and inflationary fields are related by

$$\rho _\gamma =\sqrt{\frac{3}{2\pi }}\left(\frac{m_p}{8\mathrm{\Gamma }\kappa ^2}\right)\rho _\varphi ^{3/2}.$$
(90)
Hence the energy density of the radiation field falls more rapidly than that of the inflationary field, as predicted. In the exponential model, warm inflation will end in the same way as discussed earlier for $`n>4`$ polynomial potentials.

The number of e-folds from the start of warm inflation is

$`N(\varphi )`$ $`=`$ $`\sqrt{{\displaystyle \frac{2\pi }{3}}}\left({\displaystyle \frac{4\mathrm{\Gamma }\kappa ^2}{m_p}}\right)(V^{1/2}(\varphi _0)V^{1/2}(\varphi ))`$ (91)
$`=`$ $`\left(\left({\displaystyle \frac{r}{\epsilon }}\right)_0\left({\displaystyle \frac{r}{\epsilon }}\right)_\varphi \right).`$ (92)
The amplitude of adiabatic scalar perturbations is

$$\delta _H^2=\frac{8}{75}\left(\frac{2}{\pi }\right)^{7/8}\frac{1}{3^{5/8}\alpha ^{1/4}}\left(\frac{\kappa }{m_p}\right)^{3/2}\left(\frac{\mathrm{\Gamma }}{m_p}\right)^{9/4}\left(\frac{V}{m_p^4}\right)^{3/8}$$
(93)
while that of isocurvature perturbations is given by

$$\delta _{\mathrm{iso}}^2=\frac{1}{50}(6^3\pi ^{15})^{1/8}\left(\frac{\mathrm{\Gamma }}{m_p}\right)\left(\frac{m_p^2}{a\kappa ^2}\right)^{1/4}\left(\frac{V}{m_p^4}\right)^{5/8}$$
(94)
with ratio

$$R_{\mathrm{iso}}=\frac{9\pi ^{11/4}}{16\sqrt{2}}\left(\frac{\mathrm{\Gamma }}{m_p}\right)^2\left(\frac{m_p}{\kappa }\right)^2\left(\frac{V}{m_p}\right).$$
(95)
The adiabatic and isocurvature spectral indices are

$$n_s1=\frac{3\epsilon }{4r},n_{\mathrm{iso}}=\frac{5\epsilon }{4r}.$$
(96)
For warm inflation with an exponential potential, the adiabatic spectrum is always blue in contrast with standard inflation while the isocurvature is always red. The spectral indices are related by

$$n_s1=\frac{3}{5}n_{\mathrm{iso}}=\frac{3}{8}n_g.$$
(97)


format: text

text:
# Scaling of avalanche queues in directed dissipative sandpiles

## I Introduction

In the past decade the sandpile type of cellular automata played a special role in understanding of self-organized criticality in nonlinear dynamical systems (for a recent review see Ref. ). In sandpile automata the properties of the dynamics which are essential for the occurrence of self-organized critical states can be monitored in a direct manner. Apart from the relaxation rules, these are the following properties: type of driving and time-scale separation; conservation law of the dynamics; direction of mass flow; role of boundaries. In addition, the Abelian nature of the toppling rules in some sandpile automata enabled derivation of certain exact results , in contrast to other dynamical systems where such calculations are not available. Numerous sandpile models, both deterministic and stochastic  are shown to exhibit dynamic critical states in the limit of “infinitely slow” driving (i.e., at zero driving rate $`r=0`$). In this limit a new avalanche is initiated only after the previous one has stopped, thus the time scale separation is exactly observed. On the time scale of perturbations, the avalanche evolution is seen as occurring instantly. The existence of the critical states in the case of directed Abelian sandpile automaton at zero driving rate has been proved exactly by Dhar and Ramaswamy . At this point it is interesting to mention that the model studied in Ref.  is characterized by local driving and deterministic conservative dynamics. The automaton with conservative stochastic dynamics, on the other hand, has been shown to belong to another universality class . The presence of nonconservative defects in Dhar-Ramaswamy automaton leads to a subcritical behavior .

The behavior of driven dynamical systems at finite driving rates ($`r>0`$) represents an important subject both for theoretical and practical reasons. A finite driving rate may appear either as a control parameter set from outside, or as a probability distribution originating from another coupled stochastic process. In practice, the systems are driven by an external field, which oscillates with a finite frequency. Examples are Barkhausen noise , integrate and fire oscillators , granular material in rotating drums , etc. Queuing jobs at a server , e.g., in teletraffic, is an example where the frequency of arriving jobs is given by a random process.

Theoretically at finite driving rates $`r>0`$ the probability that a new avalanche starts before previous one has stopped increases with increasing $`r`$. This obviously leads to different statistics of avalanches, where an avalanche is understood to represent a non-interrupted activity of the system. For large driving rates a continuous flow (an avalanche which never stops) is expected in sandpiles. Similarly, a single spanning cluster may occur in driven disordered systems. Therefore, a time scale separation becomes less and less apparent with increasing $`r`$. In addition, by increasing driving rates, the local driving looses its strict sense. Thus fast driven sandpiles are placed between strictly local driving, where the system is driven at a single (random) site, and global driving, where the same perturbation applies to all sites in the system. The role of the conservation law (conservation of number of grains in the interior of the system) is also expected to be changed at finite driving rates. In $`r0`$ limit, locally driven nonconservative systems appear to be subcritical , whereas when the driving is global the critical states may appear even if the dynamics is dissipative . So far neither an exact theory nor a renormalization-group analysis of fast driven critical systems has been done. A general questions as the existence of critical states at finite driving rates and universal scaling properties of the system in the limit of large distances and long times remain yet to be understood.

Recently two numerical simulations elucidated certain important properties of sandpiles at finite driving rates. In a $`1`$-dimensional ricepile model Corral and Paczuski  have first introduced avalanches for a finite driving rate as non-interrupted periods of activity and have shown that these avalanches diverge at rates $`rr_c(L)L^{0.20}`$ for a given $`L`$. The rational behind this conclusion is that an ever-running avalanche occurs for driving rates $`r1/t_0`$, where $`t_0L^{z(2\tau _t)}`$ is the average duration of avalanches in zero driving rate . In another example Barrat et al.  have shown that in a $`2`$-dimensional symmetric Abelian sandpile model mixing of time scales at finite driving rates leads to correlations which appear to violate the fluctuation-dissipation theorem.

In this work we study a simple $`2`$-dimensional model with strictly directed flow of grains and deterministic toppling rules. We add particles only at the top row with driving rate $`r`$. The driving rate $`r`$ is defined as a number of added particles per time step of avalanche evolution. We consider both conservative and nonconservative dynamics. A fraction of sites $`c`$ are considered as annealed nonconservative defects. By toppling at a defect site two grains are lost, thus affecting the propagation of avalanche below that site . When $`c=0`$ the dynamics is strictly conservative. In the $`r=0`$ limit and $`c=0`$ the model has been exactly solved by Dhar and Ramaswamy . The critical states where shown to consist of heights $`h=0`$ and $`h=1`$ occurring with equal probability 1/2 at each lattice site. The duration of an avalanche is given by the probability distribution for large $`t`$

$$P_1(t)t^{3/2}.$$
(1)
In the presence of nonconservative defects it has been shown  that the screening of the power-law distribution in Eq. (1) occurs as

$$P_1(t)t^{3/2}\mathrm{exp}(t/\xi );\xi 1/c.$$
(2)
In this paper we will refer to the distributions in Eqs. (1) and (2) as probability distributions of elementary avalanches, to be distinguished from the combined avalanches, which occur at finite driving rates and which consist of series of elementary avalanches. Up to relatively high driving rates $`r=1`$ the model has the property that successive elementary avalanches are running one after the other, in contrast to cases studied in Refs. , where merging of avalanches may occur at any finite $`r`$. After an elementary avalanche is over the system is characterized by statistically unchanged distribution of heights, owing to a weak correlation between the avalanches in the recurrent states. A finite probability of avalanche collision, which accelerates flow of grains, occurs in this model only for $`r>1`$. Here we restrict the study to the case $`r1`$, where the formation of avalanche queues is a dominant phenomenon which determines the scaling properties of the system.

The problem of avalanche queue in our model can be regarded as an example of the server queue  of the class $`E,1,GI/\mathrm{},1`$, which is studied as a model in the analysis of various applied problems, e.g., in telecommunications, insurance, etc. These analogies are made clear by noticing that the avalanches of a directed sandpile model have the corresponding terms in the language of the queue theory as follows: (i) elementary avalanche — customer; (ii) duration of elementary avalanche — service time; (iii) driving rate — frequency of arrivals of customers; (iv) number of elementary avalanches coexisting at a given moment of time — number of working servers; (v) duration of a combined avalanche — busy period. The notation $`E,1,GI/\mathrm{},1`$ means that we deal with customers arriving by one and being served by one. The letter $`E`$ means that the arrival times are generated by a Bernoulli process with the distribution $`Prob(t=k)=p^k(1p);k=0,1,2,\mathrm{};p>0`$ . The symbol $`GI/\mathrm{}`$ means that the service times are identical independently distributed (i.i.d.) random variables and the infinite number of servers provides a non-restricted number of customers which can be served simultaneously. Despite a huge literature devoted to this subject , most of papers focus on the distribution of $`q_n`$, the number of working servers at a given moment of time. The importance of the tail behavior of the busy period distribution for fluid queues in telecommunications, generalized processor sharing and other applications was recently pointed out in Ref. . Here we concentrate on some other properties of the queue: the scaling behavior of the busy period and dissipated energy distribution.

Given the duration of elementary avalanches in Eq. (1), we may conclude that the directed sandpile model for $`r1`$ represents a special case of the queue theory with the power-law distribution of the service time with the exponent $`\nu =3/21=1/2`$. This implies that the average service time per customer diverges $`t\mathrm{}`$. In practice, service times are restricted to finite values, which corresponds to the power-law distribution with $`1<\nu <2`$  . This may explain why the queue with the distribution of the type given in Eq. (1) has not been studied so far. In our model a finite average duration of elementary avalanches $`t<\mathrm{}`$ is achieved in two cases: (a) In the case of distribution in Eq. (1) when the system size $`L`$ is finite, hence the distribution is truncated at $`t=L`$; (b) In the case of finite dissipation $`c>0`$, where the distribution in Eq. (2) has a characteristic duration $`\xi <\mathrm{}`$ for all finite $`c`$ values (see Refs.  and below).

For a finite $`L`$ we find a continuous flow phase (F) for low dissipation and large driving rates, and three regions with intermittent behavior of avalanche queues. These are regions with subcritical (S), nonuniversal (N), and multifractal (M) behavior, shown schematically on the phase diagram in Fig. 1. When the length separation $`L\xi `$ holds, we find a line in the $`(r,c)`$-plane where loss of particles on defects in the interior of the pile becomes “compensated” by fast adding of particles from outside. In the region above the compensation line the avalanche queues exhibit a scaling behavior with the scaling exponents depending on the driving rate: The slopes decrease whereas the fractal dimensions increase with driving rate. Cut-offs with a stretch-exponential behavior appear. In the limit $`c0`$ and when the system size $`L\xi `$ is varied multifractal scaling properties describe the scaled distributions, rather than a simple finite-size scaling.

In general, the distributions of combined avalanches are characterized by a scaling function of two arguments of the form

$$P(X,r,L_c)X^{\tau _X}𝒫(XL_c^{D_X},rL_c^{1/2}),$$
(3)
where $`X`$ represents either duration, $`t`$, or number of topplings , $`n`$, and $`L_c\mathrm{min}(\xi ,L)`$. The corresponding fractal dimensions $`D_X`$ are defined by

$$X_{\mathrm{}}\mathrm{}^{D_X},$$
(4)
where the average is taken over all combined avalanches of a selected length $`\mathrm{}`$ measured along the direction of transport.

Using analogy to the queue theory and the properties of the recurrent states we were able to derive an exact upper limit of the distributions of busy periods and to discuss the limit $`L\mathrm{}`$. We also derived the expression for the probability of continuous flow in the conservative limit.

The paper is organized as follows: In Sec. II we define the model and consider the case of conservative dynamics by numerical simulations on finite lattice. In Sec. III we present results of simulations in the case of finite concentration of nonconservative defects. In Sec. IV we present details of the analytical results. The paper contains a short summary of the results and discussion in Sec. V.

## II Multifractal queues of Dhar-Ramaswamy avalanches

The sandpile automaton model introduced by Dhar and Ramaswamy represents an example of a self-organized criticality with exact solution in the limit of zero driving rate . In this Section we consider the same model at finite driving rates $`0<r1`$. The dynamic rules of the model are summarized as follows : We consider a $`2`$-dimensional square lattice oriented downwards, with a dynamic variable—height $`h(i,j)`$—associated at each site. Grains are added at the top row only, and mass flow is only down. The toppling at a site $`(i,j)`$ occurs deterministically whenever $`h(i,j)h_c=2`$, and two grains are transferred downward, i.e.,

$$h(i,j)h(i,j)2;h(i+1,j_\pm )h(i+1,j_\pm )+1,$$
(5)
where $`(i+1,j_\pm )`$ represents two downward neighboring sites to the site $`(i,j)`$.

The probability distribution of duration of avalanches in zero driving rate $`P_1(t)t^{\tau _t^0}𝒫(tL^1)`$ with $`\tau _t^0=3/2`$ given in Eq. (1) becomes exact at large $`t`$ . In addition the dynamic exponent $`z^0=1`$. This implies that the average duration $`t_0L^{1/2}\mathrm{}`$. Similarly, the area $`s`$ enclosed by the boundary of an avalanche is given by the distribution $`D(s,L)s^{\tau _s^0}𝒟(sL^{D_s^0})`$, where $`\tau _s^0=4/3`$ and the fractal dimension $`D_s^0=3/2`$. Note that the number of toppled grains at each active site is two, then the distribution of the number of toppled grains within an avalanche, $`D(n)`$, is described by the same exponents, i.e., $`\tau _n^0=4/3`$ and $`D_n^0=3/2`$ at zero driving rate.

A finite driving rate $`r>0`$ is implemented as follows. An avalanche is started from the top and at each step of the avalanche progress a new particle is added with probability $`r`$ at a random site at first row. We also consider a deterministic addition of particles, i.e., we add a particle at regular intervals $`\mathrm{\Delta }t`$. Both approaches lead to the same results when the statistics is high enough. An added particle may trigger a new elementary avalanche before the previous one stops, thus making a pattern of active sites distributed over lattice. A snapshot of growth of a combined avalanche with marked active sites is shown in Fig. 2 (top). A combined avalanche—avalanche queue—is thus determined by a non-interrupted activity on the lattice and it stops when no more active sites occur. Then a new avalanche is started. It should be noted that when $`r>0`$ the number of added grains is higher than the number of combined avalanches. Another important remark is that the repeated toppling at a site may occur as soon as $`r>0`$. This leads to the inequality $`n>s`$, and thus $`D_n>D_s`$, and $`z>1`$. The numerical simulations confirm these conclusions (see below). Typically, we consider $`2\times 10^6`$ combined avalanches at each driving rate and lattice size. Periodic boundary conditions are applied in the perpendicular direction.

In this Section we perform numerical simulations for $`r>0`$ and finite lattice sizes $`L`$. The limit $`L\mathrm{}`$ will be discussed in Sec. IV. In main Figs. 3 and 4 are shown the integrated probability distributions of duration $`P(t^{}t,L)`$ and number of topplings $`D(n^{}n,L)`$ for fixed driving rate $`r=0.05`$ and various lattice sizes $`L`$. It should be noted that both slopes and cut-offs of these distributions appear to be different compared to ones of the elementary avalanches. In particular, slopes decrease with $`r`$ (see more detailed discussion in next Section). For instance at $`r=0.05`$ we find $`\tau _t=0.4`$ and $`\tau _n=0.3`$, in the steep part, and $`\tau _t=0.31`$ and $`\tau _n=0.2`$ in the flat part near the cut-off. A cut-off in the probability distribution of durations appear (cf. main Fig. 3). The characteristic jump at $`t=L`$ is related to the conditional probability: an activity lasts longer than $`L`$ steps only if the preceding avalanches last not shorter than $`t=L`$. The jump decreases with increasing lattice size.

It is interesting that these distributions can not be scaled according to a simple finite-size scaling (with new exponents), as one may naively expect. Instead, we find that a multifractal scaling applies according to the law

$$P(X,L,r)=\left(\frac{L}{L_0}\right)^{\mathrm{\Phi }_X(\alpha _X)};\alpha _X=\frac{\mathrm{log}(X/X_0)}{\mathrm{log}(L/L_0)},$$
(6)
where, as before, $`X`$ stands for $`t`$ or $`n`$. The corresponding spectral functions $`\mathrm{\Phi }_t(\alpha _t)`$ and $`\mathrm{\Phi }_n(\alpha _n)`$ are determined numerically for $`r=0.05`$ and shown in the insets to Figs. 3 and 4, respectively. The spectrum depends on the driving rate. For driving rates close to the line $`rL^{1/2}`$ extremely large avalanches may appear and the scaling fits fail. The origin of multifractality in the queues of Dhar-Ramaswamy avalanches can be found in the fact that an unrestricted multiple toppling at each site may occur, and that a toppling at a given site releases a local avalanche which propagates from that site downwards.

For comparison we show how a finite-size scaling fit of the data fails. In Fig. 5 we present an attempt of scaling collapse of the data shown in Figs. 3 and 4 above, according to the formula $`P(X,L)=L^\lambda 𝒫(XL^{D_X})`$. Note that the best collapse of the tails of distributions are obtained by fixing the fractal dimensions as $`D_t=1.20`$ and $`D_n=1.75`$, which correspond to $`\alpha _X`$ at the shoulder of the spectrum $`\mathrm{\Phi }_t(\alpha _t)`$ and $`\mathrm{\Phi }_n(\alpha _n)`$, respectively. Fixing a smaller (larger) value for the fractal dimension leads to systematic shifts of the distribution tails to the right (left) with increasing $`L`$. The best fit shown in Fig. 5 is obtained for $`\lambda =0.35`$, which satisfies (within numerical accuracy) the scaling relation $`\lambda =z(\tau _t1)=D_n(\tau _n1)`$ with $`\tau _X1`$ determined at the flat part of the distribution. However, as the Fig. 5 shows, fixing $`D_X`$ and $`\lambda `$ leads to the systematic shifts of the ’horizontal’ part of the distributions to the right with decreasing $`L`$. Fixing the exponents independently from the scaling relation results in crossing of the lines for different $`L`$ values.

For driving rates $`r>L^{1/2}`$ an ever-running avalanche may occur, representing a continuous activity on the lattice. The flow phase can be characterized by an average number of topplings per site, which is expected to have a nontrivial $`L`$-dependence. The probability of occurrence of the flow phase in the limit $`L\mathrm{}`$ will be discussed in Sect. IV.

## III Nonuniversality in dissipative dynamics

In the presence of dissipative defects $`c>0`$ the distribution of elementary avalanches, which is given in Eq. (2), appears to have a finite characteristic length $`\xi <\mathrm{}`$. Thus the average duration at zero driving rate is finite $`t_0<\mathrm{}`$. Precise value of the average duration is controlled by an external parameter—probability of dissipation $`c`$, and not by system size $`L`$, provided that $`L>\xi `$. The screening length $`\xi 1/c`$ was first estimated numerically in Ref. . A more precise analytical expression can be derived (see below and Ref. ) as $`\xi ^1\mathrm{ln}(1c)`$. Here we perform numerical simulations in the case $`L>\xi `$ at driving rates $`0<r1`$. In this range of driving rates we expect the role of lattice size in the analysis of Sec. II is to be replaced by the characteristic length $`\xi `$. In addition, competition between dissipation and driving rate leads to new phenomena.

In Fig. 2 (bottom) an example of a combined avalanche is shown for $`c=0.02`$ and driving rate $`r=0.5`$. It is remarkable that the number of topplings per site decrease with distance from the top. Intermittency of the dynamics as well as the occurrence of the long-range correlations can be seen by direct examination of the recorded activity of the system $`n(t)`$ at each time step. For the server queue the quantity $`n(t)`$ is interesting as the measure of the energy which is dissipated by the server at a given moment of time $`t`$. In Fig. 6 we show an example of the recorded signal for certain choice of parameters $`r`$, $`c`$, and $`L`$ corresponding to the region (N) of the phase diagram (cf. Fig. 1). A combined avalanche on this recording is represented by a set of peaks between two consecutive drops of the signal to the base line. The Fourier spectrum of the signal (shown on top panel in Fig. 6) exhibits a power-law behavior. The slope $`\phi 0.9`$ weakly increases with driving rate $`r`$.

The scaling properties of the distributions of avalanche queues depend on the mutual ratio of the driving and dissipation rates. In particular, the scaling function in Eq. (3) exhibits a nontrivial dependence on both arguments $`xt\xi ^1`$ and $`yr\xi ^{1/2}`$. (Another suitable choice of variables would be $`(tc,rt^{1/2})`$.) In Fig. 7 we show the distribution of the avalanche mass (dissipated energy) $`n`$ for fixed $`c=0.01`$ and various driving rates $`r`$. In general, the cut-offs of the distributions increase and slopes decrease with increased driving rates $`r`$. More detailed analysis of the slopes shows that for the range of values of driving rates the scaling behavior of the distributions can be described by the scaling exponents which depend on the driving rate. The slopes of various distributions and the corresponding fractal dimensions, which are defined in Eq. (4), are shown against driving rate in Fig. 8. The dissipation rate is fixed $`c=0.01`$. Note that $`\tau _{\mathrm{}}`$ in Fig. 8 represents slope of the distribution of the largest linear length reached in a combined avalanche. For a range of values of driving rates the scaling exponents decrease and the fractal dimensions increase with $`r`$, while the scaling relations between various exponents are found to be satisfied within numerical error bars.

The variation of the scaling exponents can be approximated by a linear dependence of $`r`$. It is interesting to note that a qualitatively similar behavior—linear variation of the scaling exponents with driving rate, has been measured experimentally in the case of Barkhausen noise in driven disordered ferromagnets . Our present analysis suggests that such behavior can be related to an interplay of driving rate and dissipation at defects, and that it applies more universally. Dependence of the fractal dimensions (and of the slope exponents via scaling relations) of the avalanche queues on driving rate $`r`$ can be linked to the $`r`$-dependence of the average length of the queue $`<N>=1+r<t>_0`$. It has been shown recently  that the fractal spectrum of the series of elementary signals in the case of transit times in the ricepile model varies as a power of the length of series. A precise $`r`$-dependence of the exponents in the case of avalanche queues requires additional work and will be given elsewhere.

The opposite effects on the exponents are obtained by increasing concentration of defects $`c`$ at fixed driving rate. In particular, the slopes of the distributions increase and fractal dimensions decrease with increasing concentration of defects in a limited range .

An exact expression for the scaling function $`𝒫(x,y)`$ in Eq. (3) can not be guessed. It appears that the cut of the surface $`𝒫(x,y)`$ at $`r=const`$ for well balanced values of $`c`$ and $`r`$ can be approximated by a stretch-exponential function, so that we have

$$P(X,c,r)=X^{\tau _X(r)}\mathrm{exp}(X^{\sigma _X}/\xi _X(r)).$$
(7)
Here $`X`$ stands for $`n`$ or $`t`$, and we find $`\sigma _n=1.14\pm 0.04`$ and $`\sigma _t=1.28\pm 0.06`$, for the distribution of energy and duration, respectively. The amplitudes $`\xi _X(r)`$ can be fitted (see top panel of Fig. 8) by the following function of driving rate $`r`$

$$\xi _X(r)=A_X(c)\mathrm{exp}(rB_X(c)),$$
(8)
for a fixed dissipation rate $`c`$.

The observed parameter-dependence of the probability distributions is also reflected in the behavior of the average duration and energy of combined avalanches. Notice that the average duration $`t`$ represents the average busy period of a server in the queue theory. Beside the average duration $`t`$, we also compute the average values of number of topplings (energy) $`n`$ and area $`s`$ (total number of distinct sites) affected by the processing of a combined avalanche. These average values are shown vs. $`r`$ in the inset to Fig. 7 for $`c=0.01`$. The values $`t`$ and $`n`$ increase with driving rate faster than an exponential function. Fitting the data in the inset to Fig. 7 we find

$$X=a_{0X}\mathrm{exp}(a_{1X}(c)r^\sigma ),$$
(9)
where $`Xn,t`$ and we estimate $`\sigma =1.2\pm 0.1`$. Note that, in contrast to durations and energies, the average area of an avalanche is bounded by the number of cells in the system $`sL^2`$.

As already mentioned above, for a given dissipation rate $`c`$ and $`L\xi `$ a driving rate $`r_0(c)`$ exists such that fast addition of grains compensates the losses in the bulk. In fact along an extremal line

$$r_0(c)\kappa \xi ^{1/2}\kappa \sqrt{c}$$
(10)
the coherence length remains constant. The existence of the compensation line Eq. (10) can be demonstrated by considering sets of data for average durations and energies against driving rate $`r`$, obtained for different characteristic length $`\xi 1/c`$. These data can be scaled according to the following scaling form

$$X\xi ^{D_X(2\tau _X)}=𝒢\left(r\xi ^{z(2\tau _t)}\kappa \right).$$
(11)
The corresponding scaling fits for the two cases $`Xt`$ and $`Xn`$ are shown in Fig. 9. Notice that the respective exponents $`D_n(2\tau _n)=1`$ and $`z(2\tau _t)=1/2`$ are exact values, thus leaving only one parameter, namely $`\kappa `$, to be determined by the fitting procedure. This is an advantage of having the exact solution for the elementary avalanches . ¿From the best fit we find $`\kappa =1.14\pm 0.1`$ in the given range of values of $`c`$ (see caption to Fig. 9). It is evident from Fig. 9 that the scaling function defined in Eq. (11) increases faster than an exponential.

In the simulations a continuous flow may occur in the case of dissipative dynamics at finite lattice size $`L`$ when the driving rate is increased. However, with increased system size $`L`$ the behavior is different from the one in the case of conservative dynamics discussed above. ¿From the numerical simulations alone we can not distinguish if the affected area of a continuous avalanche diverges with the system size $`L\mathrm{}`$, or it remains finite for the range of driving rates considered here. We will also discuss large $`L`$ limit in Sec. IV.

We have restricted our analysis to the case where the degree of dissipation is such that $`\xi L`$. For $`c0`$, however, we have that $`\xi \mathrm{}`$, thus the role of $`L`$ and $`\xi `$ is interchanged at some finite $`L`$. In the reverse limit when $`L\xi `$ the behavior is expected to be similar to the case of conservative dynamics at finite $`L`$, studied in Sec. II.

## IV Analytical results

We start the analytical description of the model with the assumption that individual avalanches, which form each combined avalanche, are statistically independent. By the definition of the model, each toppling in a single avalanche occurs later than those in the previous avalanches, so that the individual avalanches never intersect in the space-time points. Nevertheless, the next avalanche is sensitive to the configuration of occupation numbers left by the previous avalanches. In this way the individual avalanches are dependent on preceding avalanches. On the other hand, it is known from Abelian properties of the directed sandpile model , that the recurrent state is characterized by the independent distribution of occupation numbers zero and one at each site. Hence, one can expect that this property of recurrent state provides independent distribution of single avalanches at least for asymptotically large systems. This assumption allows us to consider the process of driving of the directed sandpile automaton as a sequence of independent identically distributed (i.i.d.) events. Another important consequence of the statistical independence is that stops of combined avalanches can be considered as recurrent events, i.e., the probability of two successive stops $`Prob(t_1,t_2)`$ at the moments $`t_1`$ and $`t_1+t_2`$ is given by the product $`Prob(t_1)Prob(t_2)`$.

We consider the probability distribution $`F(t)`$ that a combined avalanche starting at the moment $`t=0`$ stops for the first time at the discrete moment $`t`$. This means that $`F(t)`$ is the probability that the stop of all preceding elementary avalanches occurs till the moment $`t`$. Thus, $`F(t)`$ coincides with the probability distribution of duration of combined avalanches, when an ensemble of events is considered. Along with $`F(t)`$, it is useful to consider the function $`U(t)`$ which is defined as the probability that all preceding single avalanches stop to the moment $`t`$ regardless of how many stops of combined avalanches occurred before $`t`$. Noting that stops of combined avalanches are recurrent events, we can write for $`F(t)`$ and $`U(t)`$ the following identity

$$U(t)=F(1)U(t1)+F(2)U(t2)+\mathrm{}+F(t)U(0)$$
(12)
where it is convenient to put F(0)=0 and U(0)=1. For the generating functions defined by

$$u(s)=\underset{t=0}{\overset{\mathrm{}}{}}s^tU(t)$$
(13)
and

$$f(s)=\underset{t=1}{\overset{\mathrm{}}{}}s^tF(t)$$
(14)
one easily gets from Eq. (12) the known equation of the theory of recurrent events

$$f(s)=\frac{u(s)1}{u(s)}.$$
(15)
The total probability that a combined avalanche ever stops is given by $`f(1)`$. Therefore, the probability that a combined avalanche never stops, i.e., the probability of a continuous flow, which is given by

$$F(\mathrm{})=1f(1),$$
(16)
does not vanish if $`u(1)`$ in Eq. (15) is finite $`u(1)<\mathrm{}`$ .

### A Case $`c=0`$

In the case of conservative dynamics ($`c=0`$) the probability distribution of durations of elementary avalanches is given by Eq. (1). Then we can estimate the probability $`U(t)`$ as follows. Let $`\mathrm{\Delta }t1/r`$ be the average time interval between addition of successive particles to the first row. The probability $`Prob(xt)`$ that a single avalanche has duration less than $`t`$ is

$$Prob(xt)1\frac{b}{t^{1/2}}$$
(17)
for large $`t`$, where $`b`$ is a constant of the order unity when $`L`$ is large. Then for times $`t\mathrm{\Delta }t`$ we have

$$U(t)\left(1\frac{b}{t^{1/2}}\right)\left(1\frac{b}{(t\mathrm{\Delta }t)^{1/2}}\right)\mathrm{}\left(1\frac{b}{(\mathrm{\Delta }t)^{1/2}}\right).$$
(18)
Introducing $`k=t/\mathrm{\Delta }t`$ we can write Eq. (18) as the sum

$$\mathrm{ln}U(t)=\underset{n=1}{\overset{k}{}}\mathrm{ln}\left(1\frac{b}{(n\mathrm{\Delta }t)^{1/2}}\right)\frac{k^{1/2}}{(\mathrm{\Delta }t)^{1/2}}.$$
(19)
Approximating the sum by an integral leads to

$$U(t)\left(\frac{\sqrt{rt}b\sqrt{r}}{1b\sqrt{r}}\right)^{b^2r}\mathrm{exp}[2b(r\sqrt{t}\sqrt{r})]$$
(20)
for $`t1`$ and $`0<r1`$. For the infinite lattice there exists a constant $`c_2`$ such that

$$u(1)c_2\underset{t=0}{\overset{\mathrm{}}{}}\mathrm{exp}(r\sqrt{t})<\mathrm{}.$$
(21)
Therefore, we can find that for all finite driving rates $`r>0`$ there is a non-zero probability of continuous flow. In particular, the sum in Eq. (21) diverges at small $`r`$ as

$$u(1)\frac{1}{r^2}$$
(22)
leading to the probability of stop $`f(1)12r^2b^2`$, which decreases from unity as soon as a finite driving rate is applied. Then the probability of continuous flow Eq. (15) increases from zero by the same amount, i.e.,

$$F(\mathrm{})r^2,r0.$$
(23)
For large $`r`$ we expect that

$$1F(\mathrm{})\mathrm{exp}(r).$$
(24)
If the size of the system is finite ($`L<\mathrm{}`$) the probability of stops $`U(t)`$ is bounded from above by a finite value

$$U(t)(rL)^{b^2r/2}\mathrm{exp}(2brL^{1/2}),$$
(25)
which follows from Eq. (20) taking only the dominant $`t`$-behavior for $`t=L1`$.

### B Case $`c>0`$

In the case of finite dissipation rate $`c>0`$ the dissipation leads to a finite characteristic length $`\xi `$ in the distribution of elementary avalanches in Eq. (2). This can be easily demonstrated using mean-field arguments  in reaction-diffusion systems. Let us suppose that at each site of the lattice one of species $`A`$ or $`B`$ is living. These species represent two possible states of the original model: $`A`$ corresponds to the empty site, $`B`$ to the occupied site. Due to the external driving force, new particles $`\varphi `$ are added to the first row of the lattice at rate $`r`$. The propagation of particles can be described by the following rules

$$A+\varphi B,B+\varphi A+2\varphi .$$
(26)
The kinetic equations corresponding to this scheme of “chemical” reactions are

$$\dot{n}_A(\mathrm{})=n_\varphi (\mathrm{})[n_B(\mathrm{})n_A(\mathrm{})]$$
(27)
$$\dot{n}_B(\mathrm{})=n_\varphi (\mathrm{})[n_A(\mathrm{})n_B(\mathrm{})]$$
(28)
$$\dot{n}_\varphi (\mathrm{})=n_\varphi (\mathrm{})+2n_\varphi (\mathrm{}1)n_B(\mathrm{})(1c)$$
(29)
where $`n_A(\mathrm{})`$, $`n_B(\mathrm{})`$, and $`n_\varphi (\mathrm{})`$ are concentrations of species $`A`$, $`B`$ and $`\varphi `$, respectively, at the $`\mathrm{}`$-th row. In the steady state we have $`\dot{n}_A=\dot{n}_B=\dot{n}_\varphi =0`$ and Eqs. (27)-(29) lead to the simple conditions for the concentrations $`n_A=n_B=1/2`$ and

$$n_\varphi (\mathrm{})=r(1c)^{\mathrm{}}.$$
(30)
For $`c>0`$, the density of particles $`\varphi `$ and, hence, the number of topplings in an avalanche decay exponentially with the distance $`\mathrm{}`$ from the top as $`n_\varphi (\mathrm{})=r\mathrm{exp}(\mathrm{}/\xi )`$. This implies that the characteristic length of the avalanche is $`\xi ^1\mathrm{ln}(1c)c`$. Therefore the above results, in particular Eq. (25), obtained for the case of finite lattices $`L`$ and $`c=0`$ apply also for the case $`c>0`$ by the substitution $`LL_c`$ with

$$L_c=\mathrm{min}\{\xi ,L\}.$$
(31)
### C Bounds for the busy time

If the combined avalanches are finite (i.e., there is no continuous flow), we can estimate their average duration using the known theorem from the theory of recurrent events (Ref. , Ch. XIII, Theorem 3 ). According to this theorem, the inverse average time of combined avalanches $`t^1`$ coincides with the limit of the sequence $`U(t)`$ when $`t\mathrm{}`$. Using the bound for $`U(t)`$ given by Eq. (25), we get

$$t(rL)^{b^2r/2}\mathrm{exp}(2brL^{1/2})$$
(32)
The true asymptotics of $`t`$ possibly contains an additional prefactor $`L^{1/2}`$ like in the case for $`r=0`$. Notice that we get numerically that the average duration as a function of $`r`$ (cf. inset to Fig. 7) increases faster than the exponential, which agrees with Eq. (32).

The combination $`rL^{1/2}`$ appears as a characteristic parameter determining the duration of combined avalanches. Thus, for $`\xi <L`$ it follows from Eqs. (32) and (30) that the coherence length remains constant if $`r`$ varies with $`c`$ as $`r\sqrt{c}`$, i.e. the increasing driving rate compensates the dissipation.

Another interesting feature of the probability distributions at finite driving rates is the occurrence of a stretch-exponential cutoffs both in dissipative and non-dissipative case (cf. Figs. 3,4,7). Indeed, we can see from Eq. (12) that $`U(t)>F(t)`$ for all finite $`t`$. Therefore, for non-dissipative case we have an exponential decay of combined avalanches

$$P(t)F(t)<(rt)^{b^2r/2}\mathrm{exp}(2br\sqrt{t})$$
(33)
in the thermodynamic limit $`L\mathrm{}`$, which follows directly from Eq. (20) for large $`t`$.

For finite lattice sizes $`L`$ or finite dissipation $`c>0`$ the function $`U(t)`$ is bounded from above by a constant given in Eq. (25) with $`L`$ replaced by $`L_c`$. In this case we can find the origin of an exponential cut-off in the following way. Consider an enveloping process which corresponds to propagation of the front of combined avalanches. Duration of an elementary avalanche starting at $`t_{i0}`$ is a simple linear function of time, $`t_i=tt_{i0}`$. The enveloping process consists of those topplings which occur at the maximal distance from the time axis at each moment of time $`t`$. The position $`x`$ of the front is an one-dimensional random walk confined to the interval $`[0,\xi ]`$. Starting from $`x=0`$, the walk performs a step ahead with some probability, and a step back, the length of which is a random variable. The combined avalanche stops if the random walk returns to the origin where it is trapped. The probability of the return to the origin from an arbitrary position $`x`$ is not smaller than $`U(\xi )`$, where

$$U(\xi )=(r\xi )^{b^2r/2}\mathrm{exp}(2br\xi ^{1/2}),$$
(34)
and $`\xi 1/c`$ as above. The survival time $`t`$ of the random walk under consideration does not exceed the period of successful tests in the Bernoulli scheme with the probability of “success” $`1U(\xi )`$. The period of tests in the Bernoulli process has the geometrical distribution

$$Prob(x=t)=(1U(\xi ))^tU(\xi ).$$
(35)
Hence, the time distribution for combined avalanches is bounded from above by the exponential function in Eq. (35). Using Eq. (34), we obtain

$$P(t)F(t)<A(1U(\xi ))^tA\mathrm{exp}(tU(\xi )),$$
(36)
where $`A`$ is a constant and $`U(\xi )`$ is given by Eq. (34). Note that this expression represents an upper bound of the distribution of busy periods in Eq. (7). Therefore $`1/U(\xi )`$ plays the role of an effective correlation length at finite $`r`$, in a qualitative agreement with the numerical data and Eq. (8) of previous Section.

## V Conclusions and discussion

We have shown that a finite driving rate $`r`$ is a relevant perturbation, which alters self-organized critical states in the directed sandpile automata. In the case of conservative dynamics $`r`$ couples to $`t_0L^{1/2}`$, thus leading to enhanced effects when the length scale is increased $`L\mathrm{}`$. A continuous flow eventually occurs, in which critical long-range correlations are destroyed. On a finite length scale, $`L_c=\mathrm{min}\{\xi ,L\}<\mathrm{}`$ either due to finite screening length $`\xi `$ or finite system size $`L`$, the critical states occur with qualitatively new correlation properties, which is manifested in (i) a multifractal scaling of combined avalanches when $`L\xi `$, and (ii) occurrence of compensation between driving and dissipation along a line $`r_0(c)\xi ^{1/2}\sqrt{c}`$, when $`L\xi `$. How precisely the effective coherence length $`\xi _{eff}(r)`$ of combined avalanches increases with driving rate depends on details of the relaxation process and grain addition. In the case of a finite input current at each site of the system, we find a finite toppling rate at all scales, $`n_\varphi (\mathrm{})r/c`$, compatible with $`\xi _{eff}(r)\mathrm{}`$. However, if grains are added only at the top, the correlation length increases exponentially with $`r`$ in the range $`0<r1`$, but remains finite presumably up to large driving rates. Here we restricted the discussion to the case $`r1`$, where queues of Dhar-Ramaswamy avalanches occur. Avalanche queuing for this range of driving rates is peculiar to our model, due to strictly local critical height rule and the directed transport. In the rice-pile and in the symmetric Abelian models  a perfect queuing is prevented by the collision of avalanches, which occurs at any finite $`r>0`$. Owing to the exact solution for behavior of elementary avalanches , we were able to study properties of the queue in detail. In particular, we find that average busy period is bounded from below by an exponential function $`t(rL)^{b^2r/2}\mathrm{exp}(2brL^{1/2})`$. The avalanche queue can be regarded as a multifractal set, in which the average length is regulated by the driving rate as $`N=1+rL^{1/2}`$. There are no waiting times for elementary avalanches, therefore our model corresponds to the realization known in the queue theory as “infinite number of servers”. Hence, the average number of jobs $`q_n`$ that can be served in parallel is unlimited. It should be stressed that our cellular automaton represents a new example in the queue theory in which queuing jobs are distributed according to a power-law distribution with the exponent $`\nu <2`$ and average duration of jobs is limited by a control parameter $`L_c=\mathrm{min}\{\xi ,L\}`$. We hope that more practical examples of this class can be found. We also believe that the study of the scaling properties of the queues, as we have done in this work, adds a new aspect which has not been considered so far in the queue theory. The observed nonuniversal scaling properties of avalanche queues can be related to variation of the average length of the queue with driving rate. The scaling exponents are found to vary approximately linearly with the driving rate $`r`$. A similar $`r`$-dependence was observed experimentally in other driven disordered systems, and seems to apply more generally.

A continuous activity on the lattice, corresponding to a flow phase occurs for $`rL^{1/2}`$ in the case of conservative dynamics. On an infinite lattice $`L\mathrm{}`$ probability of continuous flow increases from zero as $`F(\mathrm{})2b^2r^2`$, whereas probability of an intermittent avalanche-like flow decreases from unity with the same rate $`f(1)12b^2r^2`$. For the case of dissipative dynamics $`c>0`$ an effective coherence length increases with $`r`$ for $`1/c<L\mathrm{}`$. Our results suggest that the probability $`f(1)`$ remains finite even at high driving rates. In the limit $`L\mathrm{}`$ the compensation line extends to the point $`c0`$, $`r0`$. For the range of driving rates studied in this work we expect that the transport properties of grains in this model remain unchanged (looked at the time scale of avalanche propagation), compared to the transport at zero driving rate . Collision of avalanches, which occurs first at rates $`r>1`$ may accelerate the grain transport, possibly resulting in a new scaling behavior of the distribution of transit times. Notice that due to local critical height rules and deterministic topplings the depth of the active zone (defined in Ref. ) does not change in the flow phase of our model.

The analytical results in Sec. IV are derived assuming that elementary avalanches may be considered as independent events. We checked by computing numerically correlation function between events in a queue for finite $`L`$ that rather weak correlations occur. The correlations increase with the “distance” $`\tau `$ between avalanches as $`\tau ^\eta `$, where $`\eta =0.05\pm 0.01`$ with statistical error bars.

###### Acknowledgements.

Work of B.T. was supported by the Ministry of Science and Technology of the Republic of Slovenia. Work of V.P. was supported in part by the International Slovenian-Russian project. B.T. wishes to thank A. Corral for discussions.