text:
# Directed avalanche processes with underlying interface dynamics

## I Introduction

Avalanche phenomena are common in nature. Examples range from accumulating snow on mountain slopes, slow shearing between continental plates Carlson1989 , rerouting in river networks, to creeping magnetic flux lines in super conductors Bassler1998 . Following the work by Bak et al. Bak , physicists aim to capture the essential aspects of such dynamical systems with simple automaton processes, commonly referred to as sandpile models and self-organized criticality (SOC). Impressive successes have been achieved, like reproducing power-law distributions in avalanche events similar to those observed in nature, and the start of a classification scheme of such processes in terms of so-called universality classes Ben-Hur1996 . Unfortunately most of these are numerical in nature. Analytical exact results remain rare.

Directed avalanche phenomena form a subclass of these SOC processes. Dhar and Ramaswamy introduced the first directed sandpile model and solved it exactly Dhar1989 . This was possible because in their model the avalanche propagation is governed solely by its two edges, and those two follow independent random walk dynamics. Tadić and Dhar Tadic1997  introduced a directed model in which particles are allowed to pile up beyond the critical height, by replacing the automaton’s deterministic toppling rule by a stochastic one Tadic1997 . The density of critical sites tunes itself and at distances far from the driving edge the propagation of active sites approaches the directed percolation (DP) Kinzel1983  threshold. The scaling properties of the avalanche distributions are thus linked to the critical exponents characterizing the DP universality class. Another example of a stochastic directed avalanche process is the model introduced and studied numerically by Pastor-Satorras and Vespignani Pastor-Satorras2000 . Similar as in the above model by Dhar and Ramaswamy, the stable landscape configurations (between avalanche events) lack internal correlations in the stationary state. This allowed Paczuski and Bassler Paczuski2000  and also Kloster et al. Kloster2001  to link this dynamic process to so-called Edwards-Wilkinson Edwards1982  (EW) interface growth and to derive the exact scaling exponents of the avalanche distributions.

This novel world sheet type connection between avalanche dynamics and interface growth is particularly promising, because interface dynamic processes like EW and Kardar-Parisi-Zhang Kardar1986  (KPZ) growth are very well understood, in particular in 1+1 dimensions (1+1D) where the scaling properties are known exactly. However, the above models that are linked to EW type growth are rather poor examples, because EW growth is described by a simple linear stochastic (diffusion type) Langevin equation; correlations factorize, and important caveats in the relation to avalanche dynamics can be obscured by this simplicity.

We set out to generalize this approach to nonlinear interface dynamic processes, and recently introduced a directed unloading sandbox model Chen2002  in which the two dimensional (2D) avalanche dynamics relates to 1+1D KPZ type interface growth. We derived exponent relations between the avalanche and interface growth scaling properties, which are generic, and valid beyond our specific model. Our numerical results for the avalanche distributions (for length, width, depth, and mass) follow indeed these exponent relations. Moreover, the avalanche critical exponents obey the predicted KPZ values within a few percent, an accuracy typical to avalanche simulations. However, our numerical accuracy is better than that; mostly because of a careful finite-size scaling (FSS) analysis. The exponents seem to converge to values that are slightly different from the KPZ values.

This left us with a puzzle. What is the origin of these small deviations? Is this a fundamental effect; or do the exponents ultimately converge to the KPZ values, but with unusually large corrections to scaling. In this paper we address these issues. We also provide a more detailed discussion of these world-sheet-type relationships between avalanche and interface growth dynamics. Our first paper was short and did not include many of the details that are crucial for the analysis presented here.

The fundamental difference between conventional KPZ interface growth and avalanche dynamics arises from the averaging process over KPZ type space-time world sheets. In normal Monte Carlo (MC) simulations of interface growth the distribution functions are determined in terms of ensemble averages over a set of totally uncorrelated space-time MC runs. In contrast, the avalanche dynamics gives rise to KPZ world sheets that are strongly correlated. Two subsequent MC runs are identical except inside a single avalanche area. This difference in averaging, uncorrelated versus avalanche correlated MC runs, therefore emerges as a key issue for understanding the scaling properties of avalanche dynamics. This issue did not arise in the earlier EW type avalanche models due to the linear nature of the EW process. However for nonlinear dynamics, like KPZ, avalanche-correlated-type sampling could well lead to novel interface scaling exponents.

Speaking against a shift in the values of the exponents, are arguments like: the KPZ stationary state, i.e., the sand surface profile far way from the driving edge, can not be affected by the avalanche-correlated-type averaging, because large avalanches that span the entire width of the box occur periodically. These completely refresh the surface far way from the driving edge regularly, and thus wipe out all correlations between MC runs. This suggests that we are only dealing with much larger than usual corrections to scaling. The details are more complex than this simple argument, but we will establish that indeed the exponent values do not change.

The paper is organized as follows. In the next section we present the unloading sandbox model. In Sec. III we comment on how directed avalanche dynamics can be linked to interface growth in one lower dimension. Next, in Sec. IV, we show that in the interface growth interpretation our specific model belongs to the KPZ universality class. In Sec. V we derive the generic exponent relations between interface growth and directed avalanche dynamics, and in Sec. VI we test this numerically for our specific model.

In the second half of this paper we address the small deviations in the numerical values of the exponents from those of conventional KPZ growth. In Sec. VII we present numerical results detailing how the traces left on the surface profile by previous avalanches influence both the avalanche exponents and the interface growth ones. These scars in the rough surface enhance the surface roughness. We cast this enhanced interface roughness in terms of corrections to scaling, and determine what value the critical dimension of the corresponding irrelevant operator $`O_{\mathrm{sc}}`$ (in the sense of renormalization theory) should have. Next, we identify the geometric meaning of $`O_{\mathrm{sc}}`$, starting with a study of the one dimensional (1D) version of our model where a similar phenomenon takes place, Sec. VIII. In 1D the interface growth process is a simple random walk, and the avalanche correlated sampling relates to the scaling properties of merging random walkers. $`O_{\mathrm{sc}}`$ represents the distribution of avalanche endpoints in the 1D surface, and can be studied directly from the rounding of the surface profile near the driving edge. In Sec. IX we return to the full 2D case. The scars of previous avalanches form lines on the surface. We identify $`O_{\mathrm{sc}}`$ with the angle these lines make with respect to the direction perpendicular to the driving edge, and confirm with an analytic argument that the critical dimension of $`O_{\mathrm{sc}}`$ is equal to $`x_{\mathrm{sc}}=z`$ with $`z`$ the KPZ dynamic exponent. Finally, we summarize our results in Sec. X.

## II An Unloading Sandbox

Imagine a box filled with granular material, as illustrated in Fig. 1. One of its four retaining walls is slowly lowered, such that the sand spills out from that side, and thus slowly unloads the box and establishes a sloped surface. In the quasistatic limit, the wall moves slow enough that the unloading events can be described as distinct avalanches. The box can be three dimensional, leading to 2D avalanche dynamics on a 2D surface, or can be 2D (like in a very narrow box) giving rise to 1D avalanches on a 1D surface.

Inspired by this we consider a so-called solid-on-solid model defined on a 2D lattice. Height variables $`h(𝐫)`$ are defined on a square lattice. We will consider two versions of the model. In the continuous height version, the heights are real numbers. In the discrete model, the heights are integers, $`h(𝐫)=0,\pm 1,\pm 2,\mathrm{}`$. The former corresponds to a continuous material without internal structure, but strong cohesion up to a specific length scale $`s_c`$, while the latter corresponds to layered material where the surface height is quantized.

The 2D lattice is rotated diagonally such that the propagation direction of the avalanche is along the diagonal direction denoted by $`y`$. This is the direction in which the avalanche will run. Throughout this paper the coordinate perpendicular to $`y`$ will be denoted by $`x`$. Figure 2 illustrates this geometry.

The configurations are subject to the following stability condition. The column of particles on site $`𝐫=(x,y)`$ is supported by the two columns, $`𝐫_l=(x1,y1)`$ and $`𝐫_r=(x+1,y1)`$ directly below it and is stable when its height is less than the minimum of the heights at these two supporting sites increased by a fixed amount

$$h(𝐫)\mathrm{min}[h(𝐫_l),h(𝐫_r)]+s_c.$$
(1)
$`s_c`$ is a constant. In the version of our model where the heights are continuous variables $`s_c`$ represents the only length scale in the $`h`$ direction and can be set equal to 1 without loss of generality. Throughout this paper we will also set $`s_c=1`$ in the discrete $`h`$ model.

Consider a stable configuration, after $`\stackrel{~}{t}1`$ avalanches. The $`\stackrel{~}{t}`$-th avalanche is triggered at the highest site $`𝐫=(𝐱_{\stackrel{~}{t}},0)`$, on the $`y=0`$ driving boundary (or, in the discrete height model, by randomly choosing one of the highest sites) and reducing its height by a random amount $`0<\eta _{\stackrel{~}{t}}s_c`$. This likely creates unstable sites in the next $`y=1`$ row. Those are updated by replacing their height by an amount equal to the lowest of the two supporting columns in the previous row and then adding an uncorrelated random amount $`0\eta (𝐫)s_c`$ with uniform distribution, as

$$h(𝐫)\mathrm{min}[h(𝐫_l),h(𝐫_r)]+\eta (𝐫).$$
(2)
This updating continues row by row until all the sites are stable again. Only after that the next avalanche is started. The toppling of a site only effects the stability of the two sites immediately above it in the next $`y`$-row. Therefore we can update the system row-by-row in increasing order of $`y`$.

Direct experimental realizations of this unloading sandbox model are not our immediate concern (the focus is on establishing a generic theoretical relationship between avalanche dynamics and interface growth), but we expect that this model is applicable to actual experimental unloading sandboxes. One of the most important issues in this context is the row-by-row nature of the toppling rule. This is a crucial feature for our purposes, allowing the identification with KPZ interface growth (in the next section). In real unloading sandboxes the sand removed from row $`y`$ rolls down hill and likely disturbs the already stabilized lower surface levels. Experimental realizations can avoid this from happening, e.g., by choosing very light grains (compared to the cohesion forces). Note that our dynamic rule does not allow the build-up of any pockets (deeper than $`s_c`$) on the surface that might trap such downward rolling grains.

Conservation laws are crucial to avalanche dynamics. Unlike most avalanche processes, our model does not conserve mass while the avalanche propagates. That might raise the specter of our model not being (self-organized) critical. The connection to KPZ growth (an intrinsic critical process) dispels this phantom. Moreover, the global slope of the surface is preserved during each avalanche run, and conservation of steps in the profile plays the role analogous to conservation of mass.

The analysis of the dynamics involves distribution functions of various characteristic features of the avalanches, The common examples are: length, width, depth, and mass. The avalanche length $`l`$ will be defined throughout this paper as the maximum distance $`y`$ the avalanche travels from the driving edge; the width $`w`$ as the maximum departure of the $`x`$-coordinate (perpendicular to the propagation direction) from the trigger point $`x`$-coordinate; the depth $`\delta `$ as the maximum height change the avalanche creates at any of the affected sites; and the mass $`m`$ as the total amount of material removed by the avalanche.

## III Avalanches versus epitaxial interface growth

The focus of this paper is on how the above avalanche dynamics relates to interface growth in one lower dimension. Each stable sloped surface configuration of a directed sandpile can be reinterpreted as a world sheet (space-time configuration) of an interface in one lower spatial dimension. The direction in which the avalanches propagate plays the role of time and the perpendicular coordinates the role of space. Our 2D unloading sandbox is equivalent to a 1D growing interface. Such an interpretation makes sense only when the stability condition and the avalanche dynamic rule is directional and local in space-time, such that causality is not violated in the interface growth interpretation. The stability condition (1) and toppling rule (2) of our model are row-by-row in nature and therefore indeed Markovian in this sense.

Every stable configuration of the sandpile represents a possible interface growth life line (space-time-evolution interface world sheet). The conventional procedure for determining the scaling properties of growing interfaces is to average over a large set of completely independent MC runs. This would mean, in sandbox language, an ensemble average over completely refreshed surfaces, each totally uncorrelated from the previous one (except typically for the initial condition in row $`y=0`$). The toppling rule (2) is applied to all sites in every row, and repeated row-by-row, instead of only the unstable sites created by toppling only the highest site in the initial row.

In avalanche dynamics, however, two subsequent growing interface life lines in this ensemble differ only inside the avalanche area. From the interface growth perspective this represents a rather peculiar and dangerous correlated type MC run averaging procedure. The MC runs of KPZ space-time configuration are strongly correlated, and this raises the specter of a change in the interface roughness scaling properties. The numerical evidence, presented below is sufficiently ambiguous that this issue will preoccupy us in the second half of this paper.

## IV KPZ growth

In this section we demonstrate that the interface growth model conjugate to the unloading 2D sandbox belongs to the 1+1D KPZ universality class. The time evolution of the interface is governed by the the toppling rule of the sand model with $`y`$ in Eq. (2) representing time $`t`$,

$$h(x,t+1)=\mathrm{min}[h(x+1,t),h(x1,t)]+\eta (x,t).$$
(3)
In the conventional global type interface evolution (i.e., totally refreshing non-avalanche-type uncorrelated MC runs) every site in row $`t+1`$ is updated according to this rule.

Figure 3 illustrates the interface dynamics for one time step, $`tt+1`$. Conceptually the time step can be split into two parts; the deterministic $`\mathrm{min}[]`$ operator part and the stochastic random deposition $`\eta `$ part.

Note that because of the diagonal orientation of the square lattice (see Fig. 2), the lattice sites are not “stationary in time”. The conceptually easiest interpretation to resolve this flip-flopping is to first double the number of lattice sites and then to require them to be paired alternately with their right or left neighbors at even and odd times; at even times sites $`2n`$ and $`2n+1`$ are fused to be at equal heights and at odd times the $`2n1`$ and $`2n`$ sites.

The upper panel shows the deterministic first half of the update (from the drawn to the dash line). The partners switch and the $`\mathrm{min}[]`$ operation equalizes their heights by choosing the lowest of the two. so this step always removes material.

This can be interpreted also in terms of a movement of the steps in the interface. All up-steps move to the right and all down-steps to the left; while up and down steps merge when they meet at one site.

The lower panel illustrates the second half of the update. The height of each fused pair increases by a random amount $`0\eta s_c`$.

Deposition-type interface dynamics like this typically belongs to the KPZ universality class Kardar1986 . Indeed, Eq. (3) can be rewritten as

$`h(x,t+1)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[h(x+1,t)+h(x1,t)\right]`$ (4)
$``$ $`{\displaystyle \frac{1}{2}}\left|h(x+1,t)+h(x1,t)\right|+\eta (x,t),`$
and from this easily be identified to be a discrete form of the KPZ Langevin equation,

$$\frac{h}{t}=^2h\frac{\lambda }{2}\left(h\right)^2+\eta .$$
(5)
The crucial point is that the coefficient of the nonlinear term $`\lambda `$ is clearly present. There is no hidden special symmetry of some kind that makes it vanish by accident. At $`\lambda =0`$, the KPZ equation would reduce to EW growth.

To confirm the KPZ nature and make sure that the $`\lambda `$ is large enough that corrections to scaling from the EW point ($`\lambda =0`$) are not obscuring the KPZ scaling, we perform MC simulations on the interface dynamics as illustrated in Fig. 3. The MC runs are completely independent.

We measure the time evolution of the interface width $`W`$ defined as

$$W^2(L_x,t)\overline{(h\overline{h})^2}$$
(6)
with over bars (angle brackets) indicating average over $`x`$ (ensemble). Starting from, e.g., a flat initial condition it should scale as

$$Wt^\beta $$
(7)
at intermediate times $`0tL_x^z`$, and saturate at

$$WL_x^\alpha $$
(8)
for $`tL_x^z`$; with $`L_x`$ the length of the 1D interface. The exponents for the KPZ universality class in 1+1D are known exactly with $`\alpha =1/2`$, $`\beta =1/3`$, and $`z\alpha /\beta =3/2`$.

The numerical results are shown in Fig. 4. The values of $`\alpha (L_x)`$ are obtained from the saturated interface widths by imposing the scaling form (8) at adjacent values of the system size $`L_x`$. Similarly, the values of $`\beta (t)`$ are obtained from the transient interface widths by imposing the scaling form (7) at nearby times $`t`$. We like to remind the reader that simple log-log plots of $`W`$ versus $`L_x`$ and $`t`$ look typically impressively straight but are notoriously inaccurate. The construction of effective exponents, in the above manner might at first glance look less impressive (the data appears noisier), but this brings the analysis to a higher level where the leading corrections to finite-size and finite-time scaling become visible.

The approach to $`L_x\mathrm{}`$ in Fig. 4 is consistent with the leading correction to scaling exponent $`y_{\mathrm{ir}}=1/2`$ expected from the EW term $`^2h`$ in Eq. (5). The corrections to FSS are stronger when the height variables are discrete than when they are continuous. This is consistent with the smaller growth rate in the discrete height interface, and the fact the growth rate is typically proportional to the nonlinear term $`\lambda `$. On average, more material is removed during the first deterministic part of the update process when the surface heights are discrete.

## V Scaling Properties of 2D avalanches

In this section we derive the exact relations between the scaling properties of the avalanches and 1+1D KPZ interface growth. However, in the latter the world sheets are sampled in the correlated manner as outlined in Sec. III.

The characteristic feature of SOC is the lack of typical avalanche length, width, depth, or mass scales. The probability distributions follow power laws. For example, the distribution of avalanche widths scales as

$$P_ww^{\tau _w}$$
(9)
with scaling exponent $`\tau _w`$. Similarly, the avalanche length, depth, and mass distributions scale as power laws with exponents $`\tau _l`$, $`\tau _\delta `$ and $`\tau _m`$. We can summarize this in a metadistribution function $`P(l,w,\delta )`$; the probability to find an avalanche of a specific width $`w`$, length $`l`$, and depth $`\delta `$, obeys the scaling relation

$$P(l,w,\delta )=b^\sigma P(b^zl,b^1w,b^\alpha \delta )$$
(10)
with $`b`$ an arbitrary scale parameter. The exponents $`\sigma `$, $`z`$, and $`\alpha `$ are expected to be robust with respect to details of the dynamic rule, and thus are characteristic of the universality class to which this avalanche dynamics belongs. Single parameter distributions, such as $`P_w`$, follow by integrating out the other variables. This implies the following expressions for the $`\tau `$ exponents,

$$\tau _l=\frac{\sigma 1\alpha }{z},\tau _w=\sigma z\alpha ,\tau _\delta =\frac{\sigma 1z}{\alpha },$$
(11)
or inverted,

$$z=\frac{\tau _w1}{\tau _l1},\alpha =\frac{\tau _w1}{\tau _\delta 1},\sigma =\tau _w+z+\alpha .$$
(12)
Let’s presume that the avalanches are compact, i.e., that the inside and the boundaries of an avalanche are well defined and distinguishable (unlike in certain fractal structures), and that the sizes of the holes (unaffected regions) inside the avalanche do not scale with the avalanche size. This can be checked visually from typical simulation configurations, and both assumptions are indeed satisfied in our dynamics at least qualitatively. In that case, the mass of the avalanche must scale as $`mlw\delta `$, such that the critical exponent of the distribution of avalanche masses $`P_mm^{\tau _m}`$ obeys the identity

$$\tau _m=\frac{\sigma }{1+z+\alpha }.$$
(13)
There is one more relation between these critical exponents (leaving only two independent ones). The avalanche is initiated by lowering the bar at the driving edge of the box. In the stationary state the average surface profile is invariant, and therefore it shifts down at the same rate as the lowering bar. Thus we know how much mass drops out of the box on average.

To be more precise, during each avalanche event, the height of only one single boundary site at $`y=0`$ is lowered by, on average, an amount $`s_c/2`$. For a sandbox of width $`L_x`$ the boundary row is lowered by $`s_c/2`$ after $`L_x`$ avalanches. In the stationary state, the entire surface matches this lowering speed, such that the amount of removed sand is on average equal to $`L_xL_ys_c/2`$. Therefore, the average mass of each avalanche must be equal to

$$m=\frac{1}{2}s_cL_y.$$
(14)
The scaling properties of the mass distribution function tie into this because

$$m=m^{}P_m(m^{})𝑑m^{},$$
(15)
which can be evaluated using the metadistribution function as

$`m`$ $``$ $`{\displaystyle _0^{L_y}}𝑑l{\displaystyle _0^{\mathrm{}}}𝑑w{\displaystyle _0^{\mathrm{}}}𝑑\delta lw\delta P(l,w,\delta )`$ (16)
$`+`$ $`m_{L_y}{\displaystyle _{L_y}^{\mathrm{}}}𝑑l{\displaystyle _0^{\mathrm{}}}𝑑w{\displaystyle _0^{\mathrm{}}}𝑑\delta P(l,w,\delta ).`$
This equation incorporates finite-size effects. The box is presumed to be wide and deep enough, such that the the length $`L_y`$ of the box (in the direction perpendicular to the driving edge) is the only limiting finite-size factor. The first term in the above equation accounts for all avalanches that fit inside the box and the second term for the ones that reach the $`L_y`$ edge, and thus are prematurely terminated. The first integral scales as $`L_y^{(\sigma +2+2z+2\alpha )/z}`$ for large $`L_y`$. The second term scales with the same power, because the second integral scales as $`L_y^{(\sigma +1+z+\alpha )/z}`$ while the mass factor in front of it scales as $`mlw\delta L_y^{(1+z+\alpha )/z}`$. The result

$$mL_y^{(\sigma +2+2z+2\alpha )/z},$$
(17)
when compared to Eq. (14), yields the exponent identity

$$\sigma =2+z+2\alpha .$$
(18)
The validity of these exponent identities goes well beyond our KPZ type unloading sandbox. For example, the EW type directed avalanche models by Paczuski and Bassler Paczuski2000  and Kloster et al. Kloster2001  obey our Eq. (11) when we substitute for $`z`$ and $`\alpha `$ the EW values ($`z=2`$, $`\alpha =1/2`$). The scaling exponents of the original Dhar-Ramaswamy model can be described by the same equations with $`z=2`$, $`\alpha =0`$ as well.

## VI Numerical results for 2D sandbox avalanches

The discussion of the previous section leaves us with two independent avalanche critical exponents, $`\alpha `$ and $`z`$. The notation anticipates their identification with the scaling properties of a rough interface in interface growth. There, $`\alpha `$ is the scaling exponent of the interface width and $`z`$ the dynamic critical exponent. Indeed, the interface width relates to the depth of the avalanche, and time to the the length of the avalanche. We expect therefore that $`\alpha `$ and $`z`$ take same values as in 1+1D KPZ growth, $`\alpha +z=2`$ and $`\alpha =1/2`$.

We perform MC simulations on the sandbox avalanche model and measure the avalanche metadistribution function $`P(l,w,\delta |L_y)`$, see Eq. (10). The sandbox is always taken wide and deep enough such that the box length $`L_y`$ acts as the only FSS type limiting factor. We average over $`2^{31}`$ avalanches. The reduced distributions, such as $`P_ll^{\tau _l}`$, follow from the metadistribution from, e.g., summation over $`w`$ and $`\delta `$.

Figure 5 shows FSS approximates for the $`\tau `$ exponents. They are constructed as follows. Power-law-decaying objects such as $`P_ll^{\tau _l}`$ are almost always subject to crossover-scaling-type effects, i.e., subdominant additional power-law terms. In the language of renormalization theory they originate from so-called irrelevant scaling fields and also from nonlinear scaling field effects. This is well documented in equilibrium critical phenomena, but most recent nonequilibrium scaling studies ignore this systematic effect, e.g., by simply making a log-log plot of $`P_l`$ as function of $`l`$ and drawing a least-square-fitting-type straight line through the data. Such results show very little statistical noise, but can give rise to significant systematic errors. An example of the importance of corrections to scaling, was the large spread in reported values of the stationary state roughness exponent $`\alpha `$ between various 2D KPZ-type-growth lattice models, which was resolved using a similar FSS analysis as presented here Chin1999 .

In the limit of large $`l`$ the subdominant additional power-law terms fade away. So, more weight must be put on the large $`l`$ part of the data than on the short $`l`$ section. However, it is a balancing act, because at large $`l`$ the results become noisier, since few avalanches reach that far.

The total number of avalanches that reach beyond $`y`$ scales as

$$Q_l(y)=_y^{\mathrm{}}P_l(l)𝑑l\frac{A}{\tau _l}y^{\tau _l+1},$$
(19)
if the fraction of avalanches of length $`y`$ scales as $`P_lAy^{\tau _l}`$ (these are only the leading terms). We construct a $`y`$ dependent approximate for the exponent $`\tau _l`$ from the ratio of these two quantities, as

$$\tau _l(y)=\frac{lP_l(y)}{Q_l(y)}$$
(20)
The results are shown in Fig. 5. (We do the same for the other distributions.) Plots such as this are intrinsically noisier than conventional simple log-log type of plots of the distributions, but they contain much more information. The variation with $`y`$ reflects the leading corrections to scaling. The statistical noise at large $`y`$ could be suppressed by running the MC simulation longer. The simulation time is the only limiting factor. We used $`2^{31}`$ avalanches and in that case, $`L_y=512`$ is the optimal box size.

In Fig. 6, we replot the same data in terms of $`\alpha `$, $`z`$, and $`\sigma `$, following Eq. (12) and using the same type of FSS analysis. From the trend of the curves, we conclude that $`\alpha =0.46\pm 0.01`$, $`z=1.52\pm 0.02`$, $`\sigma =4.43\pm 0.05`$, and $`\tau _m=1.48\pm 0.01`$. This means that the exponent relations (13) and (18) are satisfied well within the statistical noise limitations, i.e., within a few percent.

Surprisingly, the actual values for $`z`$ and $`\alpha `$, although close, differ significantly from the exactly known 1+1D KPZ values, $`\alpha =1/2`$ and $`z=3/2`$. They deviate more than warranted from statistical noise alone, and do not converge smoothly if the KPZ values are correct. The approximates for $`\alpha `$ actually undershoot the KPZ value $`\alpha =1/2`$, and those for $`z`$ overshoot $`z=3/2`$. This systematic effect needs to be explained. It could be that the exponents differ in a fundamental manner from the conventional KPZ values, or that we are looking at unusually large and slow corrections to FSS. The smallness of the deviations makes the latter more likely (except when this happens to be a continuously varying exponents scenario).

We will blame the correlated MC averaging feature for this, but it should be noted that avalanche distributions are intrinsically more sensitive to FSS effects than global interface features. Many avalanches in the ensemble are small compared to the global box size, and therefore sample and average the KPZ scaling properties over much smaller lengths and shorter time scales than in a conventional global interface roughness analysis at a comparable space-time box size.

One option is to push the run button on the computer and out perform all corrections to FSS. Unfortunately it would require extremely long MC times, to create large numbers of such large avalanches. It is doubtful we would be able to get far enough in a reasonable time span. Moreover this approach is intellectually unappealing. We prefer to search for the origin of the deviations in the exponents.

## VII Avalanche correlated MC runs

The basic premise of our exponent identities is that avalanches are like any other fluctuation on a 1+1D KPZ type world sheet. Initially flat KPZ interfaces (the sand surface next to the driving edge) roughen in time (moving away from the driving edge) in such a manner that at (KPZ) time $`y`$ the stationary state roughness is established within a length scale $`l_xy^{1/z}`$. This defines a so-called spreading cone. The avalanches are expected to follow the same pattern. However, the avalanche cone seems to spread slightly faster, since the above avalanche value for $`z`$ slightly exceeds the conventional KPZ value, and inside the avalanche the surface seems to be slightly less rough, since the avalanche value for $`\alpha `$ is slightly smaller.

In this and the following section we will establish that this is caused by correlations with previous avalanches. The new avalanche does not run its course on a pristine fresh KPZ interface world sheet but on an aged one scarred by previous avalanches.

There are two obvious tests to address the effects of these scars. The first one is to determine the avalanche distributions for only the first avalanche on a fresh KPZ world sheet (the initial condition), i.e., to refresh the entire surface completely after each avalanche. The results are shown in Fig. 7. The first-avalanches likely follow normal KPZ exponents: $`z`$ converges now smoothly towards $`z=3/2`$; while the FSS approximates for $`\alpha `$, although still too small, start to turn towards $`\alpha =1/2`$, and do not cross that value anymore. It should be noted that the FSS corrections are expected to be larger, and that the data is noisier than in Fig. 6, because although we ran the same number of avalanches ($`2^{31}`$), the fraction of large avalanches is smaller, leading to smaller and noisier amplitudes in the power-law tails of the distributions.

The second test of the role of the scars is to measure the global interface roughness for avalanche type correlated MC runs instead of completely refreshing MC runs. The upper panel of Fig. 8 shows the global interface width $`W^2`$ as function of time for several $`L_x`$s. The drawn lines correspond to avalanche correlated MC runs and the dashed line to conventional uncorrelated MC averaging. The drawn lines have bumps, i.e., the avalanche correlated runs lead to rougher interfaces at intermediate times.

This enhanced interface roughness is caused by the scars left by earlier avalanches. The scars vanish at very large $`y`$ because avalanches reaching that far span the entire system in the $`x`$-direction. Figure 9 shows a typical configuration of scars. The lines are the traces of previous avalanches, i.e., their edges. Latter avalanches wipe them out partially.

For finite system sizes, the stationary state interface width follows from the plateaus at large times. There the avalanche correlated and uncorrelated MC curves coincide. This is to be expected, because the large avalanches that span the entire system (in the $`x`$ direction at large $`y`$) occur at regular MC time intervals, such that the large $`y`$ part of the surface (i.e., the stationary state of the growth process) is completely refreshed periodically and therefore sampled effectively like in uncorrelated MC runs. As a result, the roughness exponent $`\alpha `$, defined by Eq. (8), is the same for the both cases.

Most avalanches do not extend into that large $`y`$ part of the surface. They terminate in the scarred part of the surface. Therefore, we define an alternative roughness exponent $`\alpha ^{}`$, associated with the scaling of the bumps, in terms of the maximized width

$$W^{}\underset{y}{\mathrm{max}}W(L_x,y)L_x^\alpha ^{}$$
(21)
more relevant for the avalanche scaling properties. Note that for uncorrelated MC runs, $`\alpha ^{}=\alpha `$, since the interface width increases monotonically in time.

The conventional method for measuring the exponent $`\beta `$, involves the slope at times $`y<L_x^z`$, and thus is sensitive to the bumps in $`W`$ as well. The results are shown in Fig. 10. Compared to those in Fig. 4, they clearly converge less smoothly, with larger corrections to scaling and we should wonder if they converge to the conventional exact KPZ values, $`\alpha =1/2`$ and $`\beta =1/3`$, at all.

In the lower panel of Fig. 8 we plot $`\mathrm{\Delta }W^2`$ as function of time, the difference between the squared widths of avalanched correlated MC runs (the drawn lines in the upper panel) and completely uncorrelated MC runs (the dashed lines in the upper panel). For infinite system size, $`\mathrm{\Delta }W^2`$ scales as $`\mathrm{\Delta }W^2y^s`$ with an exponent that numerically is very close to $`s1/3`$. Since the width itself scales as $`W^2y^{2/3}`$, it follows that the bumps in the width curves are a transient FSS effect.

This settles our basic issue at the numerical level; the avalanche correlated nature of the MC runs does not change the interface scaling exponents, but only gives rise to slow corrections to FSS. In the next two sections we will identify these corrections to scaling with the scars on the surface left behind by previous avalanches.

We start this analysis here by casting the deviations into the framework of corrections to scaling from a so-called irrelevant operator in the sense of renormalization theory. Let $`O_{\mathrm{sc}}(x)`$ be that irrelevant operator and $`u`$ be its scaling field. This mounts to presuming that the avalanche correlation between MC runs can be represented effectively by adding to the KPZ Langevin equation (5), a term $`uO_{\mathrm{sc}}(x)`$. We will have to determine below how $`O_{\mathrm{sc}}(x)`$ is related to the density of scars on the interface space-time world sheet left by previous avalanches. According to scaling theory, the presence of such a term to the Langevin equation leads to corrections to scaling in the interface width as

$$W^2(L_x,y,u)=b^{2\alpha }W(b^1L_x,b^zy,b^{y_{\mathrm{sc}}}u),$$
(22)
i.e., in the infinite-size limit, $`L_x\mathrm{}`$, to

$$W^2(y,u)=y^{2\alpha /z}S(y^{y_{\mathrm{sc}}/z}u),$$
(23)
and by expanding the scaling function $`S`$, while assuming that $`y_{\mathrm{sc}}<0`$, such that $`u=0`$ is a stable fixed point, and the argument $`y^{y_{\mathrm{sc}}/z}u`$ is a small parameter, to

$$W^2(y,u)=y^{2\alpha /z}\left[S(0)+y^{y_{\mathrm{sc}}/z}uS^{}(0)+\mathrm{}\right].$$
(24)
The critical exponent $`y_{\mathrm{sc}}`$ of this irrelevant scaling field must take the value $`y_{\mathrm{sc}}=\alpha `$ to account for the $`\delta W^2y^{1/3}`$ corrections in the interface width we found above. Moreover the operator must scale as

$$O_{\mathrm{sc}}(x)b^{x_{\mathrm{sc}}}$$
(25)
with critical dimension $`x_{\mathrm{sc}}=z`$, since the KPZ equation (5), implies that the terms $`uO_{\mathrm{sc}}(x)`$ and $`h/t`$ must scale alike. In the following two sections we will trace down the geometric identity of this mysterious operator $`O_{\mathrm{sc}}`$, starting with the 1D version of the model.

## VIII Surface rounding in the 1D unloading sandbox

The 1D version of the unloading sandbox shows the same type of differences between uncorrelated and avalanche-type correlated MC runs as the 2D version. We determined numerically the difference between the interface width for avalanche-correlated and uncorrelated MC runs, and found that it diverges as a power law $`\delta W^2y^{1/2}`$, with an exponent which is again (like in 2D) half the size of that for $`W^2y`$ itself. According the corrections to scaling formalism (24), the scaling dimension of $`O_{\mathrm{sc}}`$ must therefore be equal to $`x_{\mathrm{sc}}=z`$, just as in 2D.

The underlying interface dynamics becomes a zero dimensional growth model, i.e., a simple random walk in the $`h`$ direction with a nonzero drift velocity to account for the net tilt of the surface. The exponents of the various avalanche distribution functions must obey the same type of relations as in Sec. V

$$\tau _l=\frac{\sigma \alpha }{z},\tau _\delta =\frac{\sigma z}{\alpha },\tau _m=\frac{\sigma }{\alpha +z},$$
(26)
and

$$\sigma =z+2\alpha $$
(27)
Without loss of generality we can set $`\alpha =1`$ (measure all lengths in terms of $`\delta `$). These identities are satisfied exactly, and the exponents are the same for uncorrelated and avalanche correlated runs. From the interface dynamics perspective, a single directed random walker, the diffusion equation character of the dynamics implies that $`z=2\alpha =2`$. The values of all the other exponents follow from this, and are consistent with their values from the avalanche perspective. There, we are dealing with the statistics of merging random walkers. The number of walkers at a given “time” $`y`$ is equal to the number of avalanches of a length $`l`$ equal or larger than $`y`$ in the ensemble of MC runs. The density of the walkers decays as $`\rho (y)y^{1/2}`$ Hinrichsen1997 , such that the distribution of avalanche lengths obeys the form

$$P_l(l)=\left[\frac{}{y}\rho (y)\right]_{y=l}l^{3/2},$$
(28)
and therefore that $`\tau _l=3/2`$. The depth of the avalanche follows from the maximum separation between two subsequent walkers, and scales as $`\delta l^{1/2}`$, i.e., $`\alpha /z=1/2`$. The mass scales as $`ml\delta l^{3/2}`$, i.e., $`(\alpha +z)/z=3/2`$ and $`\tau _m=4/3`$.

This can be compared directly with the exponents of other 1D sandpile models, e.g., with results by Paczuski and Boettcher Paczuski1996  on the so-called Oslo sandpile model where $`\tau \tau _m1.55`$ and $`D(\alpha +z)/z2.23`$.

Let’s turn our attention now to the central issue, the difference between uncorrelated versus avalanche-correlated MC runs. Adding a term like $`uO_{\mathrm{sc}}`$ to the diffusion equation of motion creates a correction to the drift velocity of the random walk. This suggests we can identify the geometric meaning of $`O_{\mathrm{sc}}`$ directly by studying the deviations of the slope near the driving edge of the surface from its asymptotic value.

The average surface slope does not show any deviations (near the driving edge) from $`s_c/2`$ when we run the dynamics as a conventional random walk, which amounts to “completely refreshing” the surface after each MC run (uncorrelated MC runs). The avalanche-correlated runs do show a rounding of the surface near the driving edge

$$s(y)Ay^\kappa +\frac{1}{2}s_c$$
(29)
The numerical results for the exponent yield $`\kappa =0.98\pm 0.03`$, in accordance with $`x_{\mathrm{sc}}=z`$ from the interface width since $`\kappa =x_{\mathrm{sc}}/z`$ and $`z=2`$ for random walks.

This rounding originates from the distribution of termination points of the avalanches. A new random walk starts below the previous one and propagates until it meets the previous trajectory and terminates. The avalanche is the space between the trajectory of that new random walk and the already existent surface. The amount of rounding of the slope near the driving edge is proportional to the distribution $`\rho (y)`$ of merging points on the surface. Those are the scars from previous avalanches. Each random walker by itself does not contribute to the rounding, i.e., on average its walk has constant slope $`s_c/2`$ independent of time. However the merging process truncates each walk and does so in an upwards biased fashion. Each merging event causes the surface to drift upwards by a certain amount ($`s_c/2`$, in average, for the discrete $`h`$ version). Therefore the rounding of the surface is proportional to $`\rho (y)`$.

The entire process and the set of subsequent stable sand surfaces (Fig. 11) is therefore equivalent to a system of merging random walkers obeying the rule $`A+AA`$. That type of dynamics has received extensive attention recently and its various scaling properties are known exactly Hinrichsen1997 . There is little doubt that our 1D unloading sandbox is exactly soluble, using absorbing wall type random walk mathematics group . However, we will refrain from pursuing this path in this paper.

The critical dimension of $`O_{\mathrm{sc}}\rho (y)`$ can be estimated (for intuition building purposes) as follows. After adding a term $`uO_{\mathrm{sc}}`$ to the KPZ equation we should also write down an equation of motion for $`O_{\mathrm{sc}}`$ itself, to close the equations. The latter is not trivial, because the scars on the surface build up slowly in time, such that that the equation of motion for $`O_{\mathrm{sc}}`$ is highly nonlocal. On the other hand, the linear nature of the diffusion equation allows one to be somewhat frivolous with the order in which averages are taken, (without losing the essential physics, nor even the correct critical exponents).

Let $`\rho _{\stackrel{~}{t}}(y)`$ be the endpoint distribution after $`\stackrel{~}{t}`$ avalanches (MC time steps). During the last MC time step, one avalanche runs through the system. It refreshes the entire surface before its termination point $`y=l_{\stackrel{~}{t}}`$, such that $`\rho _{\stackrel{~}{t}}`$ at site $`y`$ does not change if the avalanche terminates before $`y`$; $`\rho _{\stackrel{~}{t}}(y)=1`$ if it terminates at $`y`$; and $`\rho _{\stackrel{~}{t}}(y)=0`$ if it extends beyond $`y`$:

$$\frac{\rho _{\stackrel{~}{t}}(y)}{\stackrel{~}{t}}=P_l(y)\rho _{\stackrel{~}{t}}(y)_y^{\mathrm{}}P_l(l)𝑑l$$
(30)
with $`P_l(l)`$ the probability that the avalanche terminates at distance $`l`$ from the driving edge. The stationary state endpoint profile therefore takes the form

$`\rho (y)={\displaystyle \frac{P_l(y)}{_y^{\mathrm{}}P_l(l)𝑑l}},`$ (31)
and $`P_l(l)l^{\tau _l}`$ yields

$$\rho (y)=\frac{1}{\tau _l1}y^1.$$
(32)
In other words, the surface curvature scales as $`\mathrm{\Delta }sy^{x_{\mathrm{sc}}/z}`$ with $`x_{\mathrm{sc}}=z`$, in agreement with the above results. Interestingly, this result is independent of the actual value of the scaling exponent $`\tau _l`$, provided that $`\tau _l>1`$, which has to be true for $`P_l`$ to be normalizable.

In conclusion, in 1D we identified the crossover scaling operator with the density of avalanche endpoints. These represent indeed the scars on the surface, the memory of previous avalanches.

## IX Avalanche rounding near the driving edge in 2D

As in the 1D model, the surface slope is modified by the iterated avalanche process. However, unlike in 1D, the average slope near the edge is not constant already in conventional interface dynamics (where the entire surface is being refreshed during each MC run). The surface slope is related to the growth rate of the underlying interface, and the rounding of the slope near the driving edge represents the transient growth rate of the KPZ interface from the initial configuration, e.g., a flat one:

$$s_\mathrm{f}(y)v_0+cy^{\kappa _\mathrm{f}}$$
(33)
with $`y`$ playing the role of time and the subscript, $`\mathrm{f}`$, denoting that the entire surface is refreshed. By direct numerical simulation of uncorrelated interface dynamics we find $`\kappa _\mathrm{f}0.7`$ (the left panel of Fig. 12). This is consistent with conventional KPZ scaling and power counting

$$sh/yy^{\alpha /z1}y^{2/3},$$
(34)
suggesting $`\kappa _\mathrm{f}=2/3`$.

We evaluate the surface slope profile $`s(y)`$ in avalanche correlated dynamics MC runs, in terms of the difference with respect to the uncorrelated case,

$$\mathrm{\Delta }s(y)=s(y)s_\mathrm{f}(y)y^\kappa $$
(35)
The FSS analysis for the exponent $`\kappa `$ (the right panel of Fig. 12) yields $`\kappa =1.05\pm 0.07`$. This is in agreement with $`x_{\mathrm{sc}}=z`$ and $`\kappa =x_{\mathrm{sc}}/z`$ implied by the corrections to scaling formalism (24).

Inside the bulk of an avalanche the interface is fully refreshed, and scales as in uncorrelated KPZ dynamics. At the avalanche boundaries, the slope of the surface is biased upwards, because of the merging with previous MC runs (which are on average shifted upwards by an amount $`s_c/2L_x`$ each time an avalanche is triggered). This means that the $`\mathrm{\Delta }s`$ is proportional to the density of scars in the surface. In 1D, the scars are point-like objects, the endpoints of the avalanches; but in 2D the avalanche boundaries are line objects. This nonscalar aspect makes that most line-segment contributions, when integrated along the boundaries of an avalanche, cancel out against each other.

To be more precise, $`s(y)`$ represents only the component of the slope in the $`y`$-direction, and the magnitude of those jumps depends on the local angle $`\theta `$ the boundary makes with the $`y`$-axis. This is an odd function, $`\mathrm{\Delta }(\theta )=\mathrm{\Delta }(\theta )`$, as illustrated in Fig. 13. The slope change is negative when the avalanche opens up and positive when it narrows down. The latter also implies that $`\mathrm{\Delta }(\theta )`$ has opposite sign for the left and right boundary of each avalanche. Notice that, while in the lattice model $`\theta `$ takes only two discrete values, it renormalizes to a continuous variable at larger length scales.

Let’s estimate the change in surface slope due to these scars in the same spirit as we did successfully in 1D. Consider one specific surface, and let $`s_{\stackrel{~}{t}}(y)`$ be the surface slope in a slice of the surface at distance $`y`$ from the driving edge, averaged over all $`x`$, after $`\stackrel{~}{t}`$ avalanches (MC time $`\stackrel{~}{t}`$). The last avalanche changes this as follows. Let $`w_{\stackrel{~}{t}}(y^{})`$ be the width of this avalanche, which terminates at $`y=l_{\stackrel{~}{t}}`$, in slice $`y^{}`$. The inside area of the avalanche is completely refreshed and therefore has the same average slope $`s_\mathrm{f}(y)`$ as in ordinary KPZ dynamics (totally refreshed subsequent world sheets). This leads to the following equation of motion,

$$\frac{s_{\stackrel{~}{t}}(y)}{\stackrel{~}{t}}=[\mathrm{\Delta }(\theta _L)\mathrm{\Delta }(\theta _R)]+w_{\stackrel{~}{t}}(y)[s_\mathrm{f}(y)s_{\stackrel{~}{t}}(y)]$$
(36)
The first term on the right hand side represents the creation of the two new avalanche edges, and the second term represents the refreshed surface inside the new avalanche. Note that $`s_{\stackrel{~}{t}}(y)/\stackrel{~}{t}=0`$ when this latest avalanche does not reach slice $`y`$, and that this is automatically taken care of because in that case $`\theta _L=\theta _R=0`$ and $`\mathrm{\Delta }(0)=0`$, while $`w_{\stackrel{~}{t}}(y)=0`$ for $`y>l_{\stackrel{~}{t}}`$. In the stationary state, after averaging over all possible avalanches, Eq. (36) leads to

$$\overline{w_{\stackrel{~}{t}}(y)\left[s_\mathrm{f}(y)s_{\stackrel{~}{t}}(y)\right]}=\overline{\mathrm{\Delta }(\theta _L)\mathrm{\Delta }(\theta _R)}$$
(37)
Next, we perform an heuristic coarse-graining renormalization-type transformation. At large length scales, the average angle $`\theta `$ remains small, such that the right hand side can approximated as

$$\overline{\mathrm{\Delta }(\theta _L)\mathrm{\Delta }(\theta _R)}a\overline{\theta _L\theta _R}a\frac{\overline{w_{\stackrel{~}{t}}(y)}}{y}$$
(38)
Finally, we presume that in the stationary state it is not too bad to treat the KPZ height fluctuations deep inside the bulk of an avalanche and those near its edge as decoupled (at least in lowest order) such that

$$\mathrm{\Delta }s(y)=\overline{s_\mathrm{f}(y)s_{\stackrel{~}{t}}(y)}=a\frac{}{y}\mathrm{log}(\overline{w_{\stackrel{~}{t}}(y)}).$$
(39)
This yields $`\mathrm{\Delta }s(y)y^1`$, exactly the power-law decay we are looking for, and consistent with all the above numerical results.

The only requirement for the latter is that $`\overline{w_{\stackrel{~}{t}}(y)}y^\xi `$ decays as a power law. Again, like in Eq. (32) for 1D, the value the critical exponent $`\xi `$ does not matter. $`\overline{w_{\stackrel{~}{t}}(y)}`$ is equal to the average avalanche width in slice $`y`$ averaged over all avalanches. It is reasonable to expect, and we confirmed numerically, that this quantity scales with the same exponent as the average width of all avalanches longer than $`y`$, i.e., as

$$_y^{\mathrm{}}w(l)P(l)𝑑ly^{1/z\tau _l+1}$$
(40)
which yields $`\xi 1/3`$.

We are now ready to represent the crossover scaling operator $`O_{\mathrm{sc}}(x)`$ in terms of the scars on the surface. Consider time slice $`y`$. $`O_{\mathrm{sc}}(x)=0`$ when no scar line runs through site $`x`$, and otherwise is proportional to the angle the scar line makes with respect to the $`y`$-axis. However the sign also flips depending on whether this represents a left or right boundary of the original avalanche. The latter can be denoted by an arrow along the avalanche scar line. Alternatively, we can associate an age-field $`g(x,y)`$ to the entire surface, representing the age of the surface segments (how many MC time steps ago site $`x`$ was updated),

$$O_{\mathrm{sc}}\frac{\widehat{e}_yg}{|g|}$$
(41)
with $`\widehat{e}_y`$ a unit vector in the $`y`$-direction. The denominator arises because the magnitude of the age jump across the scar line $`|g|`$ does not play a role.

## X Summary

In this paper, we studied a directed avalanche model inspired by the unloading of a sandbox by means of a slowly lowering wall, and the wish to setup an avalanche dynamic rule belonging to the same universality class as KPZ type interface growth. The 2D sand surface represents the world sheet of the 1+1D growing interface.

The scaling exponents of the avalanche distributions are directly related to the dynamical and stationary state roughness exponents $`z`$ and $`\alpha `$ of KPZ growth in 1+1D, Eq. (11). However, we encounter one crucial difference. From the avalanche perspective the conventional uncorrelated MC runs correspond to completely refreshing the surface, i.e., an ensemble average over all possible initial conditions, without ever running an avalanche. From the KPZ perspective, the avalanche dynamics represents an unusual MC ensemble averaging procedure where subsequent interface world sheets only differ inside the single avalanche. This avalanche-correlated-type averaging enhances the interface roughness at time scales $`y<L_x^z`$, due to the scars of previous avalanches. It required a careful study, combining numerical and analytical tools, presented in the second half of this paper, to establish that these scars give rise only to larger than usual corrections to scaling and not to fundamentally different values of the global roughness scaling exponents $`z`$ and $`\alpha `$.

The effect of the scars can be represented by introducing an additional age field $`g(x,y)`$ to the height variables $`h(x,y)`$, that keeps track of how many MC runs ago site $`(x,y)`$ participated in an avalanche. This age-field couples into the KPZ equation (5) as an additional term of the form $`uO_{\mathrm{sc}}`$. The operator $`O_{\mathrm{sc}}`$ is proportional to the angle a scar makes with respect to the time-axis, and can be expressed in terms of the age field as shown in Eq. (41). We establish that the coupling of this age field to the KPZ equation is irrelevant in the sense of renormalization theory, both numerically and by writing down approximate equations of motion for $`uO_{\mathrm{sc}}`$. The scaling field $`u`$ renormalizes with exponent $`y_{\mathrm{sc}}=\alpha `$ and $`O_{\mathrm{sc}}`$ scales with critical dimension $`x_{\mathrm{sc}}=z`$.

We believe that the results of our work presented here can be generalized to most “Markovian” avalanche dynamic systems with local row-by-row type toppling rules, and that this is a promising route to improve our understanding of the scaling properties of avalanche dynamics in general.

This research is supported by the National Science Foundation under grant DMR-9985806.


format: text

text:
# DISCOVERY OF SOFT X-RAY EMISSION FROM IO, EUROPA AND THE IO PLASMA TORUS

## 1. Introduction

Imaging and spectral data from the infrared through the extreme ultraviolet provide important information on the makeup of the surfaces and atmospheres of the Galilean satellites Io, Europa, Ganymede, and Callisto (Carlson et al. 1999, Barth et al. 1997, Hall et al. 1998), discovered by Galileo Galilei in 1610, and of the Io plasma torus (IPT; Hall et al. 1994, Woodward, Scherb, & Roesler 1997, Gladstone & Hall 1998). Here we report the results of high-spatial resolution x-ray observations of the Jovian system with the Chandra X-ray Observatory showing that the Galilean satellites and the IPT are also x-ray emitters, and we discuss the physical processes that may contribute.

## 2. Observations

The Chandra X-ray Observatory (Weisskopf et al. 2000) observed the Jovian system on 25-26 Nov 1999 with the Advanced CCD Imaging Spectrometer (ACIS), in support of the Galileo flyby of Io, and on 18 Dec 2000 with the imaging array of the High Resolution Camera (HRC-I), in support of the Cassini flyby of Jupiter (Gladstone et al. 2001). During the ACIS observations, spanning 86.4 ks, the telescope focus and the planet were placed on the back-side illuminated CCD designated S3 in order to take advantage of this CCD’s sensitivity to low energy x-rays. The spacecraft was repointed four times to allow for the planet’s motion and was always oriented to allow the planet to move along the CCD’s second node. These data were corrected for charge-transfer-inefficiency effects (Townsley et al. 2000), and only the standard ASCA grades 0, 2, 3, 4, and 6 were retained in order to reduce the background induced by charged particles. Charge-transfer-inefficiency refers to the loss of charge during CCD readout for pixels furthest away from the readout node, leading to a low estimate for the event energy. Event grades depend on the distribution of charge induced by the event about the center pixel. X-ray and charge particle events produce different event grade distributions, allowing a powerful way to reduce unwanted background events. One ACIS pixel is 0.492 arcsec by 0.492 arcsec and the half-power diameter of the Chandra system level point-spread function is about 0.8 arcsec. No repointings were necessary during the HRC-I observations, spanning 36.0 ks (approximately one rotation of Jupiter), because of the instrument’s larger field of view and the observation’s shorter length. One HRC-I pixel is 0.1318 arcsec by 0.1318 arcsec.

Although the ACIS S3 optical filter is optically thick to optical light, there is an interference peak in its transmission near $``$9,000Å ($``$1.4 eV) which is likely to affect ACIS S3 observations of solar system objects and nearby bright late-type stars. Indeed, the ACIS S3 data for Jupiter were compromised by the throughput in low energy channels of large numbers of optical photons from the planet’s bright disk. To eliminate this problem for our analysis, all events within 34.4 arcsec ($``$1.4 times Jupiter’s radius) of the planet’s center were removed from the data. In any event, removal of this data is required for analysis of the Galilean satellites in order to avoid contamination by events associated with Jupiter. The planet’s optical emission also produced a bump in the ACIS bias frame, and all events within 34.4 arcsec of the center of this bump were removed from the data. For sufficiently bright x-ray sources, x-ray events occurring during the 41 ms necessary to read out each 3.2 s ACIS accumulation appear as a streak, or trailed image, in the data. For our ACIS Jupiter observations, optical photons from the planet’s disk do not contribute to a trailed image, as 41 ms is too short a time for sufficient accumulation to mimic an x-ray event. The trailed image from Jupiter’s true x-ray events is a small effect for E $`>`$ $``$440 eV, where we have reasonably accurate knowledge of Jupiter’s x-ray spectrum. From an extrapolation of Jupiter’s x-ray spectrum to lower energies, we infer that the trailed image is also not a dominant effect for our analysis of the Io Plasma Torus at energies down to 250 eV. The trailed image is unimportant for our analysis of x-ray emission from the Galilean satellites. Since we wished to remove x-ray events associated with Jupiter, we removed events within 34.4 arcsec of the planet’s center from the HRC-I data.

### 2.1. The Galilean satellites

In order to search for emission from each Galilean satellite, we removed Chandra’s orbital motion and transformed the data into a co-moving frame using an ephemeris obtained from the Jet Propulsion Laboratory. The size of the detect cells was set equal to the root-sum-square of the satellite’s radius and 2 arcsec. For the background, we determined an average number of events per detect cell in a larger region around the satellite in its own frame. This region was 98.9 by 98.9 arcsec for ACIS and 98.7 by 98.7 arcsec for HRC-I. The number of events within the detect cell centered on each satellite compared with this background value determines a probability of chance occurrence according to Poisson statistics. The results of this analysis appear in Table 1. Since these objects move across the field of view, Chandra’s spatial resolution at their position varies, so in each case we also compared the radial distribution of events around the satellite with that around 100 points chosen at random in the surrounding region. Io is detected in both the ACIS and HRC-I data, while Europa is detected in the ACIS data only. Ganymede appears at the 99% confidence level in the ACIS data. A simple calculation of joint probabilities confirms the detection of Io and Europa, with results that we find tantalizing for Ganymede and Callisto. The nominal energies of the x-ray events range from 300 eV to 1890 eV and seem to show a clustering between 500 and 700 eV. The mean energy, weighted appropriately by Chandra’s effective area, of the events are 515 eV and 664 eV, for Io and Europa, respectively. The corresponding estimated energy flux at the telescope and emitted power from the satellite are $`4.1\times 10^{16}`$ erg/s-cm2 and 2.0 MW for Io, and $`3.0\times 10^{16}`$ erg/s-cm2 and 1.5 MW for Europa. Figure 1 shows smoothed ACIS images of Io and Europa.

### 2.2. The Io Plasma Torus

Figure 2 shows the HRC-I image zoomed back from Jupiter, in the planet’s frame of reference. There is a region of diffuse emission on the dusk side of the planet, and a fainter region of diffuse emission on the dawn side of the planet. This morphology is similar to the dawn-dusk asymmetry seen in EUVE data (Hall et al. 1994, Gladstone & Hall 1998). The paths traced by Io, Europa, and Ganymede during the HRC-I observation are also shown. Io and Europa pass through the dusk side region of diffuse x-ray emission. We associate these regions of diffuse emission with the IPT. The background subtracted HRC-I count rates are $`3.2\times 10^3`$ counts/s-arcmin2 on the dawn side of Jupiter and $`5.67.1\times 10^3`$ counts/s-arcmin2 on the dusk side of Jupiter. The larger rate for the dusk side region applies to a circular source region of 1.36 arcmin2, and the smaller to the larger elliptical source region of 2.43 arcmin2 suggested by the ACIS observations. The rate for the dawn side region applies to both the smaller and larger source regions.

Figure 3 shows the ACIS image zoomed back from Jupiter, in the planet’s frame of reference over the energy ranges 250–1,000 eV (left panel) and 504–621 eV (center panel). Because this image was produced from four overlapping observations, and also because of the instrumental effects discussed above, we show an exposure map (right panel) for the image in Figure 3. The area of low exposure around Jupiter is due to the removal of events within 34.4 arcsec of the planet’s center or within 34.4 arcsec of the center of the bump in the bias frame. The ACIS image also shows regions of diffuse emission on the dawn side and dusk side of the planet. The paths traced by Io, Europa, Ganymede, and Callisto during the ACIS observations are also shown. In this case Io and Europa pass through the dawn side region of diffuse x-ray emission. The ACIS 250–1,000 eV background subtracted count rates are $`1.7\times 10^3`$ counts/s-arcmin2 both on the dawn side and dusk side of Jupiter. However, the 504–621 eV image shows significant asymmetry with stronger emission on the dawn side where Io and Europa were located during the ACIS observations. The HRC-I image also shows stronger emission on the side (dusk in that observation) where Io and Europa were located during that observation. The narrow band ACIS image also shows emission to the north of Jupiter, which we do not associate with the IPT.

Both a simple power-law and thermal bremsstrahlung give statistically acceptable fits to the background subtracted 250–1,000 eV ACIS IPT spectrum (see Table 2 and Figure 4). However, adding a single gaussian line does significantly improve these fits. The central energy of the gaussian is very close to a strong O VII line at 574 eV. Lines from other oxygen charge states may also contribute. The background subtraction procedure correctly removed instrumental lines such as Si K$`\alpha `$ at 1.74 keV. The 250–1,000 eV energy flux at the telescope from the IPT (within the ovals marked on Figure 3) is $`2.4\times 10^{14}`$ erg/s-cm2, corresponding to a luminosity of 0.12 GW, and is approximately evenly divided between the dawn side and dusk side. The apparent x-ray line emission around 570 eV originates predominantly on the dawn side of the planet. EUVE observations in 1993 (Hall et al. 1994, Gladstone & Hall 1998) determined an IPT energy flux of $`(7.2\pm 0.2)\times 10^{11}`$ erg/s-cm2 over the interval 370–735 Å (17–34 eV). This flux is $``$3,000 times that inferred from the Chandra data over the interval 250–1,000 eV.

## 3. Discussion

We have observed x-rays that originate from the IPT, the Galilean satellites Io and Europa, and probably Ganymede. Two known energy sources are available to drive these processes: (1) the solar x-ray flux; and (2) the energetic particle populations in the region of the IPT. Furthermore, we must consider the x-ray production processes that originate from the surfaces of the Galilean satellites as distinct and different from the x-ray processes that lead to emission from the IPT.

### 3.1. The Galilean satellites

We first consider the surface emissions from Io and Europa. Based on the results of model calculations (Peres et al. 2000), the energetic particles that constitute the outer radiation belts of Jupiter, and are found within the IPT, provide energy fluxes to the surfaces of Io and Europa $``$20 and $``$110, respectively, times larger than the solar x-ray flux during solar maximum conditions with active flares. We show below that the relevant incident particles and energies are those which deposit most of their energy in the top $``$10 micron of the surface. Irradiation models show that ion energy losses dominate that from electrons in this layer (Paranicas, Carlson, & Johnson 2001, Cooper et al. 2001). Even along Europa’s trailing hemisphere, we suspect the total energetic ion dose dominates the energetic electron dose down to about tens of microns in depth, below which point the situation is reversed (Paranicas, Carlson, & Johnson 2001). This, together with the clustering of x-ray event energies between 500 and 700 eV, suggests that the observed emission cannot be due to electron bremsstralung. Since the cross-sections for interaction and the efficiency for x-ray production from energetic ion impact and from fluorescence excited by solar x-rays are comparable, we conclude that the x-ray emission from Io and Europa is ultimately powered by the energetic ion population in the IPT. The surfaces of the Galilean satellites are exposed to bombardment by energetic ions diffusing radially inward through the region of the IPT. The incident flux is made up of energetic H, O and S ions (Paranicas et al. 2001, Cooper et al. 2001). Although Io and Europa have ionospheres which deflect low energy ions, more energetic ions with large gyro-radii will reach their surfaces and produce x-rays through particle-induced x-ray emission (Johansson, Campbell, & Malmqvist 1995). In fact the large gyro-radii for the x-ray producing ions means the effective cross-section of the Galilean satellites to IPT ions is larger than their physical size (Pospieszalska & Johnson 1989).

Io’s atmosphere arises from sublimation of surface SO2 frost, sputtering of the surface by Jovian magnetospheric particles, and volcanic activity, with the latter perhaps being the dominant contributor. This atmosphere is patchy rather than uniform, but can be expected to provide a lower limit to the energy of ions impacting the surface. Approximately 5 keV H+, 40 keV O+, and 60 keV S+ ions penetrate an atmosphere with a column density of $`10^{17}`$ SO2 molecules/cm2 (Lanzerotti et al. 1982); these limits fall between 1 and 10 keV for a column density of $`10^{16}`$ SO2 molecules/cm2. For comparison, electrons with energy greater than about 10 keV penetrate a column density of $`10^{17}`$ SO2 molecules/cm2 (Michael and Bhardwaj 2000). The integrated ion energy flux at Europa for H, O, and S above 10 keV (Paranicas et al. 2001, Cooper et al. 2001) is $`10^{10}`$ keV/s-cm2. The energetic ion flux at Io is $``$5 times smaller, but is less well known than for Europa. The depth of the region of interest is limited by the absorption of escaping x-rays in the surface. For Europa’s water ice surface, absorption in oxygen dominates, while for Io’s surface, absorption in oxygen is complemented by that in silicon, sulfur, and other elements. Oxygen K$`\alpha `$ x-rays at 525 eV in ice have an optical depth of unity at about 9.5 micron, implying an oxygen column density of 3.2 $`\times `$ 1019 cm-2. For incident protons with energy $`>`$ $``$300 keV, and an oxygen target, the cross-section for x-ray production is flat with energy and is $`10^{21}`$ cm2, implying an x-ray yield per proton of $``$0.01 for Europa. Similarly, the presence of sulfur dioxide (atomic weight 64) rather than water (atomic weight 18) suggests an an x-ray yield per proton $``$0.04 for Io. These yield estimates are uncertain; accurate estimates require more detailed work. The production cross-sections for energetic O and S ions are larger but the corresponding penetration depths are smaller suggesting similar yields. Dividing the ion energy flux by an assumed average ion energy of 100 keV and multiplying by the x-ray yield per ion gives an x-ray surface emission rate of $`1.5\times 10^6`$ photon/s-cm2 for Europa and $`1.2\times 10^6`$ photon/s-cm2 for Io. If 10% of these photons are O K-shell x-rays with energies around 525 eV, then the predicted power in O K-shell x-rays is $`3.9\times 10^{13}`$ erg/s and $`4.2\times 10^{13}`$ erg/s for Europa and Io, respectively. The corresponding fluxes at the telescope are $`8.1\times 10^{16}`$ erg/s and $`8.7\times 10^{16}`$ erg/s for Europa and Io, respectively. The fluxes predicted by these simple estimates are within a factor of $``$three of the observed fluxes for Europa and Io. Despite these low flux levels, further observations with the Chandra and XMM-Newton x-ray observatories may eventually permit the identification of lines from individual heavy elements.

### 3.2. The Io plasma torus

Turning now to the IPT, we estimate its x-ray flux (in O K-shell lines from O, O+, and O++) from fluorescence caused by the x-ray flux from the Sun. We use the results of model calculations for x-ray emission from the flaring active Sun (Peres et al. 2000), normalized to the model x-ray luminosity in the energy band 0.1-10 keV, to determine the solar spectrum at Jupiter’s orbit (the planet was 4.14 AU from the Earth and 4.96 AU from the Sun during the Nov 25-26, 1999, ACIS observations, and 4.12 AU from the Earth and 5.04 AU from the Sun during the Dec 18, 2000, HRC-I observations). Photo-ionization cross-sections are from analytic fits to the results of calculations using the Hartree-Dirac-Slater method (Verner et al. 1993). For photo-ionization of oxygen ions at energies just above the K-edge, the cross-sections are of order few $`\times `$ 10-19 cm2, decreasing with increasing energy roughly as 1/E3. For the fluorescent yield we use the value for neutral oxygen (0.0083) (Krause 1979). Galileo measured electron (Bagenal et al. 1997) and ion (Crary et al. 1998, Frank & Paterson 2000) densities above 1000 cm-3 in some regions of the IPT. Assuming the volume of a torus with semi-major axis equal to that of Io’s orbit and diameter equal to Jupiter’s equatorial diameter, and number densities in neutral oxygen and in ionized oxygen (mostly O+) of 1000 cm-3, we find a photon flux at the Earth of 2.6 $`\times `$ 10-7 photons/s–cm2 for the flaring active Sun. Our spectral fits require photon fluxes in the apparent line near 570 eV of $``$3 $`\times `$ 10-6 photons/s–cm2, about an order of magnitude larger. We conclude that the observed x-ray emission from the IPT is not due to fluorescence excited by solar x-rays.

Charge exchange processes have been invoked to explain Jupiter’s x-ray aurora (Cravens et al. 1995, Bhardwaj & Gladstone 2000) and x-ray emission from comets (Cravens 1997, Dennerl 1999). As energetic ions diffuse inward through the region of the Io plasma torus they are constantly subject to both electron stripping and charge exchange. Electron stripping generally has a higher probability than charge exchange above some cross-over energy. The cross-over energy is different for the different charge states with higher charge states having a generally higher cross-over energy. The cross-over energy for an OVII ion, above which electron stripping is favored, occurs above 11 MeV (Cravens 1997). Ions below a few tens of keV rapidly charge exchange in the outer torus, become neutralized, and are lost from the Jupiter system (Mauk et al. 1998). These energetic neutrals have been recently observed by the Cassini spacecraft during its flyby of Jupiter (Mitchell et al. 2001). However, as the higher energy ions slowly diffuse radially inward through the plasma torus region, increasing numbers are driven to higher and higher charge states. In fact we argue that all oxygen ions with energies over 16 MeV are fully stripped by the time they diffuse radially inward to 6 Jupiter radii near Io’s orbit (where they are observed). Upon colliding with a torus neutral and undergoing charge exchange to O VII, such ions can produce an x-ray photon at energies near 570 eV. From the Galileo Heavy Ion Counter experiment (Cohen et al. 2000), we find an integrated flux of oxygen ions above 16 MeV of 1.5 $`\times `$ 104 ions/s-cm2 . Combining this with an estimated charge exchange cross-section (Cravens et al. 1995) of 10-18 cm2, an estimated neutral density of 30 cm-3, and an estimated x-ray emission probability of 0.065, we obtain an x-ray volume production rate of 4.5 $`\times `$ 10-13 photons/s-cm3. Taking the same IPT volume as assumed for our fluorescence estimate, we find a photon flux at Earth of $``$4 $`\times `$ 10-10 photons/s-cm2, approximately 4 orders of magnitude smaller than the required flux near 570 eV. Similar considerations show that ion stripping by dust grains leading to x-ray emission also fails to explain the apparent line emission.

A large fraction of the IPT x-ray flux may result from bremsstrahlung radiation from non-thermal electrons. Above a few tens of eV, the electrons’ energy distribution is well-described by a power law rather than a Maxwellian (Sittler & Strobel 1987, Meyer-Vernet, Moncuquet, & Hoang 1995). These non-thermal electrons play a critical role in setting the charge state of the Io torus, and the energy balance of the IPT. Estimates of cooling from electron-impact excited UV emissions and heating from Coulomb collisions with ions (Smith et al. 1998) show the need for a large, additional source of energy. The nature of this energy source remains speculative, but it may result from wave-particle heating (Barbosa 1994) or inward diffusion of energetic ions and electrons from Jupiter’s middle and outer magnetosphere (Smith et al. 1998). Previous remote sensing observations have provided limited data on these non-thermal electrons, since UV and longer wavelength emissions are primarily sensitive to electrons with energies below a few eV.

Bremsstrahlung emission of soft X-rays are, on the other hand, produced by the non-thermal electrons in the few hundred to few thousand eV range. Based on the in-situ Ulysses observations, we assume the electrons are in a kappa, or generalized Lorentzian, distribution with a temperature of 10 eV and an index of $`\kappa `$ = 2.4 (Meyer-Vernet, Moncuquet, & Hoang 1995). Using this distribution and the Bethe-Heitler cross-section for bremsstrahlung emission (with Sommerfeld-Elwert correction), we numerically calculate the source rate and spectrum. At a density of 2000 cm-3, typical of the IPT, the source rate is 6 $`\times `$ 10-18 ergs/cm3/s. Integrated over the entire IPT, the radiated x-ray power between 250 and 1000 eV is approximately 3 $`\times `$ 107 W, or roughly one third of the observed signal. The predicted shape of the bremsstrahlung spectrum is harder than the measured 1/E6.8 spectrum, indicating that some of the flux near 250 eV is probably an extension of the FUV emission from the IPT (Hall et al. 1994, Gladstone & Hall 1998).

## 4. Conclusions

Chandra X-ray Observatory observations have shown that Io and Europa, and probably Ganymede, emit soft x-rays with 0.25-2.0 keV luminosities $``$1-2 MW. This emission is likely to result from bombardment of their surfaces by energetic ($`>`$ 10 keV) H, O, and S ions from the region of the IPT. The IPT itself emits soft x-rays with a 0.25-1.0 keV luminosity $``$0.1 GW. Most of this emission appears at the low end of the energy band, but an unresolved line or line complex is apparent in the spectrum at an energy consistent with oxygen. The origin of the IPT x-ray emission is uncertain, but bremsstrahlung from non-thermal electrons may account for a significant fraction of the continuum x-rays.

Further x-ray observations and more detailed modelling are needed to probe more deeply into the origin and properties of the x-ray emission from the Galilean satellites and from the Io plasma torus, and their implications for the Jovian magnetosphere.

The authors thank L. Townsley for making CTI corrections to the ACIS data, and F. Bagenal for helpful communications.


text:
# Heat Transport in a Strongly Overdoped Cuprate: Fermi Liquid and Pure 𝑑-wave BCS Superconductor

## Abstract

The transport of heat and charge in the overdoped cuprate superconductor Tl2Ba2CuO6+δ was measured down to low temperature. In the normal state, obtained by applying a magnetic field greater than the upper critical field, the Wiedemann-Franz law is verified to hold perfectly. In the superconducting state, a large residual linear term is observed in the thermal conductivity, in quantitative agreement with BCS theory for a $`d`$-wave superconductor. This is compelling evidence that the electrons in overdoped cuprates form a Fermi liquid, with no indication of spin-charge separation.

A fundamental question about the rich and baffling behavior of electrons in high-temperature superconductors is whether or not it can be understood in the framework of Fermi-liquid (FL) theory, the standard theory of electrons in solids. Several authors believe that when the concentration of electronic carriers in these cuprate materials is sufficiently low, as in the so-called underdoped region of the doping phase diagram, the basic excitations of the electron system are not the usual Landau quasiparticles characteristic of FL theory. In one class of proposals  for example, the electron is thought to fractionalize into a neutral spin-carrying excitation, called a “spinon”, and a spinless charge-carrying excitation, called a “holon” or “chargon”. However to this day, such “spin-charge separation” has not been confirmed experimentally. On the other hand, after 15 years of intensive research it is still not known with any certainty whether or not the ground state of cuprates is a Fermi-liquid in any region of the phase diagram. It is widely assumed that in the metallic-like overdoped regime at high carrier concentration FL theory does hold, but there is little solid evidence to support this lore.

In this Letter, we present the results of a study which show that strongly overdoped cuprates do not undergo spin-charge separation and their ground state is most likely a Fermi liquid. By measuring the transport of both heat and charge in the normal state at very low temperature, we were able to verify that one hole-doped cuprate in the overdoped regime obeys the Wiedemann-Franz (WF) law. This universal law is a robust signature of FL theory, stating simply that the electronic carriers of heat are fermionic excitations of charge $`e`$. In addition, the thermal conductivity in the superconducting state is found to be in good agreement with BCS theory for a superconductor with a pure $`d`$-wave order parameter.

The particular compound chosen for this study is Tl2Ba2CuO6+δ, because it can easily be overdoped. In many ways, it is the ideal cuprate material. Its crystal structure is tetragonal, without the CuO chains that complicate the properties of the orthorhombic compounds YBa2Cu3O7-δ (Y-123) and YBa2Cu4O8 (Y-124), or the buckling that alters the unit cell of Bi2Sr2CaCu2O8 (Bi-2212). It is made of a stack of single CuO2 planes, and is therefore not subject to possible bi-layer effects such as encountered in Bi-2212. It has a high maximum critical temperature $`T_c^{max}`$ of 90 K, at optimal doping, much as in Y-123 and Bi-2212. In this sense, it is free of the possible concerns about the low $`T_c`$ found in single-plane La2-xSrxCuO4 (LSCO). Finally, the $`d_{x^2y^2}`$ symmetry of its superconducting state has been confirmed by phase-sensitive measurements, at least at optimal doping .

Mackenzie et al. measured the resistivity of a strongly-overdoped crystal of Tl-2201, with $`T_c=15`$ K . In zero magnetic field, $`\rho (T)`$ was found to follow roughly a power law of $`T^{1.8}`$ from room temperature down to $`T_c`$ and extrapolate to $`\rho _07\mu \mathrm{\Omega }`$ cm at $`T=0`$. The resistive upper critical field at $`T0`$, $`H_{c2}(0)`$, is between 12 and 16 T (for fields perpendicular to the conducting CuO2 planes, i.e. $`Hc`$), depending on the precise criterion.

The overdoped samples of Tl-2201 used in this study were rectangular single crystals with typical dimensions of $`0.4`$ mm and $`0.2`$ mm in the tetragonal basal plane and $`10\mu `$m along the $`c`$-axis. The voltage pads had a width of $`25\mu `$m and the spacing between the electrodes was $`0.3`$ mm. They were grown by the same technique as used by Mackenzie and co-workers in previous studies . They have $`T_c15`$ K, in zero magnetic field. Using the empirical formula $`T_c/T_c^{max}=182.6(p0.16)^2`$, this translates into a carrier concentration of $`p=0.26`$ hole/Cu atom. To obtain such critical temperatures, the samples were annealed in 1 bar of flowing O2 at 350 oC for two days. The resistivity of our samples is essentially identical to that obtained previously , with $`\rho _0=5.6\mu \mathrm{\Omega }`$ cm. Both heat and charge transport were measured using the same contacts, made by diffusing silver epoxy. A typical value for the contact resistance was 0.1 $`\mathrm{\Omega }`$ at 4 K. The thermal conductivity was measured down to below 100 mK with a standard one-heater two-thermometer technique in a dilution refrigerator. The magnetic field was applied along the $`c`$-axis. The geometric factor used to convert from resistance (electrical or thermal) to electrical resistivity $`\rho `$ or thermal conductivity $`\kappa `$ was set by requiring that $`\rho (300\mathrm{K})=180\mu \mathrm{\Omega }`$ cm, the value obtained in previous studies of numerous crystals with the same doping level . The uncertainty on this value is estimated at $`\pm 10\mu \mathrm{\Omega }`$ cm.

The resistivity is shown in Fig. 1, for fields ranging from zero to above $`H_{c2}(0)`$. A slight positive magnetoresistance is observed, in agreement with previous work . The resistivity below 30 K (and above $`T_c`$) is best fit by the function $`\rho =\rho _0+bT+cT^2`$, with a substantial linear term (i.e. $`bT>cT^2`$ for $`T<15`$ K). The fitting parameters are $`\rho _0(H)=5.84`$, 5.99 and 6.15 $`\mu \mathrm{\Omega }`$ cm at $`H=7`$, 10 and 13 T, respectively, and $`b=0.064\mu \mathrm{\Omega }`$ cm K-1, $`c=0.0054\mu \mathrm{\Omega }`$ cm K-2 at 13 T. This unusual dependence was reported previously  and interpreted as “non-Fermi-liquid” behavior, in the sense that no linear term is expected in conventional FL theory. Deviations from the standard $`T^2`$ dependence have been observed in a number of heavy-fermion materials, for example, $`T^{1.2}`$ in CePd2Si2 below 20 K . In these materials, this is associated with the proximity to a quantum critical point (QCP), where antiferromagnetic order sets in as a function of pressure or chemical composition. In the case of cuprates, the obvious QCP would be the onset of superconductivity at a critical concentration $`p_c`$ close to 0.3 hole/Cu atom, but a QCP has also been postulated to exist inside the superconducting region.

The thermal conductivity $`\kappa `$ is shown in Fig. 2. The data is plotted as $`\kappa /T`$ vs $`T^2`$ to separate the contribution of electrons from that of phonons, given that the asymptotic dependence of the former as $`T0`$ is linear in $`T`$ while that of the latter is cubic. In other words, in Fig. 2, the electronic contribution is the residual linear term $`\kappa _0/T`$ given by the intercept of a linear fit with the $`T=0`$ axis. The value of $`\kappa _0/T`$ obtained in this way is: 1.41, 2.76, 3.47, 3.75, 3.87, 3.90, 3.95, and 3.95 mW K-2 cm-1, at $`H=0,1,2.5,4,5.5,7,10,`$ and $`13`$ T, respectively. As explained above, the uncertainty on the overall absolute value is approximately $`\pm 5\%`$. However, the relative uncertainty, e.g. between different fields, is much lower, around 1 %. This high degree of reliability is due to the fact that in these samples electrons conduct much better than phonons, and hence the slope of $`\kappa (T)/T`$ in Fig. 2 is weak relative to the intercept. Note that at high fields, electrons scatter phonons very effectively and $`\kappa (T)`$ is entirely electronic below 1 K.

Fundamentally, the linear term in $`\kappa `$ at $`T=0`$ reveals the presence of fermionic excitations in the electron system. We can then ask whether these excitations carry charge. This question can only be addressed in the absence of any superfluid that can also carry charge, which amounts to testing the WF law in the normal state. This law is one of the most fundamental properties of a Fermi liquid, reflecting the fact that the ability of a quasiparticle to transport energy is the same as its ability to transport charge, provided it cannot lose energy through collisions. It states that the heat conductivity $`\kappa `$ and the electrical conductivity $`\sigma `$ of a metal are related by a universal constant:

$$\frac{\kappa }{\sigma T}=\frac{\pi ^2}{3}(\frac{k_B}{e})^2L_0$$
(1)
where $`T`$ is the absolute temperature, $`k_B`$ is Boltzmann’s constant and $`L_0=2.44\times 10^8`$ W $`\mathrm{\Omega }`$ K-2 is Sommerfeld’s value for the Lorenz ratio $`L\kappa /\sigma T`$ . Theoretically, electrons are predicted to obey the WF law at $`T0`$ in a wide range of environments: in both three or two dimensions (but not strictly in one dimension), for any strength of disorder and interaction , scattering and magnetic field . Experimentally, the WF law does appear to be universal at $`T0`$: until recently, no material had been reported to violate it. The first exception was found in optimally-doped Pr2-xCexCuO4 (PCCO), an electron-doped cuprate .

It is in general difficult to test the WF law in cuprate superconductors because of their high upper critical fields. In our crystals, the superconductivity has completely vanished by 13 T, at which field we find $`\kappa _0/T=3.95\pm 0.04`$ mW K-2 cm-1 and $`\rho _0=6.15\pm 0.03\mu \mathrm{\Omega }`$ cm, so that $`L=\rho _0\kappa _0/T=0.99\pm 0.01L_0`$, in perfect agreement with the WF law. Note that the Lorenz ratio does not suffer from the 5% uncertainty associated with the geometric factor, as both transport measurements are performed using the same sample with the same contacts. The error bars are therefore on the order of 1 %. In Fig. 2, the transport of heat and charge are compared directly by reproducing the charge conductivity at $`13`$ T from Fig. 1. This is done by plotting $`L_0/\rho (T)`$ vs $`T`$ using the fit to the $`13`$ T data for $`\rho (T)`$ (inset of Fig. 1). The charge conductivity $`L_0\sigma (T)`$ is seen to be equal to the heat conductivity $`\kappa (T)/T`$ at $`13`$ T.

The basic implication of this result is that the fermions which carry heat also carry charge $`e`$ and are therefore indistinguishable from standard Landau quasiparticles. In particular, there is no evidence of any spin-charge separation. Indeed, if electrons were to fractionalize into neutral spin-carrying fermions (spinons) and charged bosons (chargons) , there would be no reason to expect the WF law to hold, as the heat-carrying fermions would not take part in the transport of charge. This result therefore imposes a constraint on theories of spin-charge separation (SCS): the critical hole concentration $`p_{SCS}`$ at which electron fractionalization starts to occur is not the QCP where superconductivity starts to occur (on the overdoped side of the phase diagram), but can only be lower. In other words, any hypothetical onset of SCS must obey $`p_{SCS}<0.26<p_c`$. It therefore appears that the mechanism for superconductivity in this overdoped region of the phase diagram is not the condensation of charge-$`e`$ bosons, but most likely Cooper pairing. Note that (barring any profound electron-hole asymmetry) this conventional picture is expected to break down with underdoping, as suggested by the violation of the WF law in PCCO near optimal doping .

Although the standard FL description fails, as revealed by the non-quadratic $`T`$ dependence of $`\rho (T)`$, the basic nature of the electronic excitations in the limit of zero energy is that of Landau FL quasiparticles. (A similar situation is seen in heavy-fermion materials .)

In the absence of a magnetic field, there is a large residual linear term in the thermal conductivity of Tl-2201, namely $`\kappa _0/T`$ = 1.41 mW K-2 cm-1. A similar term has also been observed in other hole-doped cuprates, albeit at optimal doping, where it is much smaller: $`\kappa _0/T`$ = 0.14, 0.15 and 0.11 mW K-2 cm-1, in Y-123  , Bi-2212  and LSCO , respectively. Within BCS theory applied to a $`d`$-wave superconductor, this residual heat conduction is expected, arising from zero-energy quasiparticles induced by impurity scattering near the nodes in the $`d_{x^2y^2}`$ gap function. In the clean limit, where the scattering rate $`\mathrm{\Gamma }k_BT_c/\mathrm{}`$, it is universal (in the sense that it is independent of $`\mathrm{\Gamma }`$) and it depends only on the ratio of the two quasiparticle velocities ($`v_F`$ and $`v_2`$) which govern the Dirac-like spectrum of nodal quasiparticles, $`E=\mathrm{}\sqrt{v_F^2k_1^2+v_2^2k_2^2}`$ :

$$\frac{\kappa _0}{T}=\frac{k_B^2}{3\mathrm{}}\frac{n}{c}(\frac{v_F}{v_2}+\frac{v_2}{v_F})$$
(2)
where $`n`$ is the number of CuO2 planes per unit cell of height $`c`$ (along the $`c`$-axis), and $`\stackrel{}{k_1}`$ and $`\stackrel{}{k_2}`$ are unit vectors pointing in directions normal and tangential to the Fermi surface at the node, respectively. In other words, $`v_F`$ is the Fermi velocity in the nodal direction and $`v_2`$ is proportional to the slope of the gap at the node, $`d\mathrm{\Delta }/d\varphi =\mathrm{}k_Fv_2`$, with $`k_F`$ the Fermi wavevector.

Applying Eq. 2 to Tl-2201, for which $`n=2`$ and $`c=23.2\AA `$, we get $`v_F/v_2=270`$. A rough estimate using Fermi surface parameters typical of cuprates, namely $`v_F=2.5\times 10^7`$ cm/s and $`k_F=0.7\AA ^1`$ (the values measured in Bi-2212 ), and the simplest $`d`$-wave gap function, $`\mathrm{\Delta }=\mathrm{\Delta }_0\mathrm{cos2}\varphi `$, with the weak-coupling relation for a $`d`$-wave superconductor, $`\mathrm{\Delta }_0=2.14k_BT_c`$, gives $`v_F/v_2=210`$. This shows that the magnitude of $`\kappa _0/T`$ is in good agreement with the simplest BCS analysis.

It should be recognized that even though the mean free path in these samples is rather long (in the range $`5001000\AA `$ ), the scattering rate $`\mathrm{\Gamma }`$ is not small compared to $`T_c`$. It may be estimated using the standard expression for the normal state conductivity: $`\kappa _N/T=\frac{1}{3}\gamma _Nv_F^2\tau `$, where $`\gamma _N`$ is the specific heat coefficient and $`\tau =1/(2\mathrm{\Gamma })`$. With $`\gamma _N3`$ mJ K-2 mole-1  and $`v_F=2.5\times 10^7`$ cm/s, one gets $`\mathrm{}\mathrm{\Gamma }0.4k_BT_c`$. At finite $`\mathrm{\Gamma }`$, corrections to Eq. 2 give an increase in $`\kappa _0/T`$ . Assuming $`\mathrm{\Delta }_0=2.14k_BT_c`$, the correction for $`\mathrm{}\mathrm{\Gamma }/k_BT_c=0.4`$ is by a factor of approximately $`1.5`$ . Thus the correct value of $`v_F/v_2`$ is probably closer to $`270/1.5=180`$ .

It will be interesting to investigate the doping dependence of $`\kappa _0/T`$ as a way of measuring the dependence of the gap function on carrier concentration, via $`v_2`$. In the absence of further data on Tl-2201, we may compare with optimally-doped Y-123 ($`T_c=93`$ K) or Bi-2212 ($`T_c90`$ K), for which $`v_F/v_2=14`$ and 19, respectively . (The value of 19 for Bi-2212 agrees very well with the value of 20 obtained from ARPES measurements of $`v_F`$ and $`v_2`$ separately .) Under the assumption of a doping independent $`v_F`$, verified in both Bi-2212  and LSCO , one immediately sees that $`v_2`$, or the magnitude of the gap (near the nodes), scales roughly with $`T_c`$. This strongly suggests that the standard BCS relation between gap magnitude and transition temperature, $`\mathrm{\Delta }_0T_c`$, holds in the overdoped regime. This is in striking contrast with what is found in the underdoped region of the phase diagram. Indeed, our measurements on underdoped Y-123 and LSCO  reveal that $`v_F/v_2`$ decreases as $`T_c`$ is reduced by underdoping (see also ).

Several authors have proposed the existence of a QCP within the superconducting dome in the phase diagram of cuprates, either as a theoretical prediction to explain the diagram itself or as suggested in various experiments. Its location is usually taken to be near optimal doping, in the neighbourhood of $`p=0.2`$. If it is associated with a change in the symmetry of the superconducting order parameter, Vojta et al. have argued that the most likely scenario is a transition from a pure $`d_{x^2y^2}`$ state to a complex order parameter of the form $`d_{x^2y^2}+ix`$, where $`x`$ can have either $`s`$ or $`d_{xy}`$ symmetry . Dagan and Deutscher have recently reported a split zero-bias anomaly in their tunneling on Y-123 thin films as soon as the material is doped beyond optimal doping, a feature which they attribute to the appearance of a complex component to the order parameter in the bulk . The presence of a subdominant component $`ix`$ in the order parameter causes the nodes to be removed, as the gap can no longer go to zero in any direction. The observation of a residual linear term in the thermal conductivity, a direct consequence of nodes in the gap, therefore excludes the possibility of any such subdominant order parameter. (More precisely, since our measurement goes down to 100 mK, it puts an upper bound on the magnitude of $`|x|`$ relative to $`|d_{x^2y^2}|`$ at about 0.5 %.) Moreover, there is no subdominant order parameter in Tl-2201 at optimal doping . In other words, if there truly is a QCP between optimal doping at $`p0.16`$ and the critical point $`p_c0.3`$, it does not appear to be associated with the onset of a complex component in the order parameter.

In summary, the low-temperature transport properties of Tl-2201 with $`T_c=15`$ K show that spin-charge separation does not occur in strongly overdoped cuprates. The normal state at $`T0`$ satisfies the Wiedemann-Franz law perfectly, demonstrating that the only electronic excitations carrying heat and charge are Landau quasiparticles. The superconducting state obeys BCS theory in that the residual heat conduction is of the expected magnitude for a $`d`$-wave gap and the dependence of the low-energy spectrum on doping strongly suggests that the gap scales with $`T_c`$ in the conventional way. Finally, the possibility of a sub-dominant order parameter ($`ix`$) is ruled out.

We are grateful to D. Hawthorn, C. Lupien, J. Paglione, F. Ronning, and M. Sutherland for help in various aspects of the measurements, and to S. Kivelson, C. Lannert and K. Behnia for useful discussions. This work was supported by the Canadian Institute for Advanced Reasearch and funded by NSERC of Canada.
$``$ Present address: Laboratoire National des Champs Magnétiques Pulsés, 31432 Toulouse, France.