{'fraction_non_alphanumeric': 0.03416112245570121, 'fraction_numerical': 0.01474211185033924, 'mean_word_length': 4.36779577116187, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Avalanche Spatial Structure and Multivariable Scaling Functions; Sizes, Heights, Widths, and Views through Windows

We introduce a systematic method for extracting multivariable universal scaling functions and critical exponents from data. We exemplify our insights by analyzing simulations of avalanches in an interface using simulations from a driven quenched Kardar-Parisi-Zhang (qKPZ) equation. We fully characterize the spatial structure of these avalanches- we report universal scaling functions for size, height and width distributions, and also local front heights. Furthermore, we resolve a problem that arises in many imaging experiments of crackling noise and avalanche dynamics, where the observed distributions are strongly distorted by a limited field of view. Through artificially windowed data, we show these distributions and their multivariable scaling functions may be written in terms of two control parameters, the window size and the characteristic length scale of the dynamics. For the entire system and the windowed distributions we develop accurate parameterizations for the universal scaling functions, including corrections to scaling and systematic error bars, facilitated by a novel software environment SloppyScaling.


I. INTRODUCTION
Systems that have crackling noise and avalanches exhibit scale invariance and power laws, which point to the notion of underlying universality [1]. These systems include many of the best-studied examples of nonequilibrium critical phenomena, and much progress has been made in a renormalization group context [2,3]. The renormalization group implies that the long length and time behavior near critical points is governed by universal exponents and scaling functions. However, the predictive power of these theoretical studies has hiterto been underutilized; the primary focus of experiments and numerical simulations has been on precise estimates of critical exponents, rather than on the universal joint predictions of properties involving several control parameters and/or measured quantities.
A wide variety of materials and natural systems have been studied in the context of non-equilibrium critical phenomena. Many of these systems exhibit avalanches which have power law size distributions. These include Barkhausen noise in ferromagnets [4][5][6][7][8], fluid imbibition into porous media, flux-line depinning [3,[9][10][11][12][13], and martensitic transformations [14,15], to name a few. In the first three of these systems, avalanches are the result of the jerky motion of an interface (domain wall, fluid front, flux-line) in a disordered environment, and can be described by the same family of front-propagation models.
In this manuscript, we study the spatial structures of avalanches in a front-propagation model in two dimensions, developing tools and methods needed for systematic study and extraction of these multiparameter universal scaling functions. To illustrate the utility and importance of these functions, we apply them systematically to a practical experimental problem -the size distributions of avalanches seen through a viewing window. This problem illustrates (a) the complexity and sophistication of the different emergent size distributions, (b) the relationships amongst the probability distributions of heights, widths, and sizes and their utility in generating predictions for windowed avalanches, and (c) the use of functional forms and least-squares fits to analyze and report on these multiparameter scaling functions.
Imaging experiments have been used in recent years to study a wide variety of systems exhibiting crackling noise or similar dynamics. Barkhausen noise is measured making use of the magneto-optical Kerr effect [5,6,16], allowing one to examine the domain wall motion in 2D thin films in real-time. In experiments on superconducting vortices [17], a magneto-optical (MO) setup is also used. In experiments of fracture [18], fluid imbibition [19] and granular systems [20,21], the dynamics are also followed with a camera.
These visualization experiments provide an unusual opportunity: we now can study the universal properties of the spatial morphology -various distributions of heights, widths, angles, local heights, etc of either the arXiv:1112.0714v1 [cond-mat.stat-mech] 4 Dec 2011 avalanches or the fronts. However, the measurements of these properties are often distorted by a limited field of view. We hereby take this problem and develop the scaling theory for the universal functions needed to characterize the results of a generic imaging experiment -the distribution of avalanche sizes seen through windows.
The limited field of view in experiments distort the size distributions of avalanches, and cause difficulties in characterizing the critical exponents. Naturally, there is a bias towards small avalanches; large ones are cut off by the boundaries of the window. It can also distort the size distribution if pieces of large avalanches cut off by window boundaries are counted as small ones. Experiments have taken a variety of approaches to deal with such windowing effects: "Laser reflectometry" [16] on Barkhausen noise uses the magneto-optical Kerr effect, but lacking spatial resolution, lumps fragments and avalanches together; meanwhile, other optical Kerr experiments have shown [6] that the effective size exponent τ for this lumped distribution depends strongly on the window width. Work by Kim et al. [22] report quite striking distributions but do not specify whether their data includes avalanches that touch the boundary. In experiments on superconducting vortices, or magnetic flux avalanches [17] avalanches exceeding a certain size are discarded. In fluid imbibition [19], the edges of the system are purposely left out to avoid any distortion produced by side walls. In granular systems, where avalanche dynamics in piles of rice are studied with real-time reconstruction [20,21], and fracture experiments [18], where the dynamics are followed with a high-speed camera, boundary effects are not considered but may also be important.
We will show comprehensively how to analyze all of the size data lying within a window, and how to use the different classes of avalanches to get independent measures of various critical exponents. Indeed, window-width finitesize effects need not be avoided, but properly treated may provide additional measures of the critical exponents and the spatial structure.
Characterizing spatial structures of avalanches must go far beyond the traditional focus on critical exponents. Many experiments report power laws, however through this study we emphasize that one can make predictions about both power laws and scaling shapes from data, as has been demonstrated in a previous study of avalanche temporal shapes [23]. Indeed, traditional scaling collapse methods fail for functions of more than two variables. To optimally extract the behavior and estimate errors, we need to do simultaneous analysis of many different properties. We thus introduce a software environment, SloppyScaling [24], which facilitates the exploration and development of simultaneous fits of multiple data sets with parameterized forms of universal scaling functions. With this approach we are taking the first steps towards the use of scaling methods as a practical engineering tool.

II. SUMMARY OF KEY RESULTS
Since our theme is multifaceted, readers may be interested in focusing on different aspects of this work. In this section we present an overview of key results in this paper and a summary of their corresponding sections.
1. Universal spatial structures of avalanches in directed percolation depinning. We provide a substantive analysis of the universal spatial morphology of avalanches in the quenched KPZ (qKPZ) model in 1+1 dimensions (the model is discussed in section III). Figure 1(a) shows avalanches in a typical simulation of this model. Analogous to magnetic systems, we have added a "demagnetization factor" k that parameterizes the restoring force, which allows us to access many metastable configurations of the front near the depinning transition, and controls the typical width of an avalanche, L k ∼ k −ν k . Figure 2 shows examples of these resulting avalanches from simulations of various k. In section IV, we thoroughly examine the spatial structure of the avalanches, including sizes s (the total area covered between pinned fronts, which would correspond to the total magnetization change in magnets, or the avalanche size of a rice pile), and also widths w and heights h (which measure the length of an avalanche in directions perpendicular and parallel to the direction of the motion of the front; this corresponds to studying the shapes of the magnetic domains or flux lines). We examine and fit the distributions of these sizes, heights and widths in section IV.
2. Avalanches in Windows. As mentioned in the Introduction, in many imaging experiments the limited field of view distorts the avalanche size distribution. This is illustrated in Fig. 1(b) and (c). In Fig. 1(b) most of the avalanches are cut off by the left and right boundaries-if one were to count the area of these avalanches, we would count large ones as smaller ones. The resulting size distribution inside the window (dashed line of Figure 1(c)) has a very different power law and shape compared to the full system distribution (solid line of Figure 1(c)).
In section V, we show in detail how avalanches which cross boundaries exhibit distinctly different size distributions and critical exponents. For simplicity, we consider a strip geometry where only the left/right boundaries are relevant (since avalanches in our system are flat and wide, few touch the top or bottom). We can therefore separate avalanches into four different categories, 00, 01, 10, and 11, depending on whether avalanches touch (1) or do not touch (0) the left or right boundaries: internal avalanches (00) (Figure 1(d)), split avalanches (10 and 01) (Figure 1(e) and (f)), and spanning avalanches (11) (Figure 1(g)). The internal (00) avalanches share the same power law as that of the full distribution, with a cutoff controlled by both the system size L k and the window size W . The split (10 or 01) avalanches have a modified power law-a smaller exponent as larger avalanches are counted as smaller halves. (See Figure 1)(c)). The spanning avalanches also exhibit a smaller exponent, although this is not obvious in the shape of its scaling function. We can see that it has both an outer cutoff due to L k and an inner cutoff due to the window size W , since avalanches must be large enough to span the window (purple dash-dotted line in Figure 1(c)). The internal and spanning avalanches also have distinct universal scaling functions with a cutoff controlled predominately by the window size for windows comparable to or smaller than the size of L k . In Section V and Section VI we give a thorough analysis of these modified power laws, the different scaling shapes, and the results of fits to data.
Having established a sophisticated method of analyzing both experimental and simulation data, we can utilize this analysis to enhance the collection of data in visualization experiments. Section VII has some suggestions for how to collect data and simultaneously analyze the scaling behavior of different magnifications, and extract multiple exponents.
3. Functional forms A main emphasis of our work is that we fit an entire functional form instead of power laws [25], this includes the shape of the scaling function, and analytic and singular corrections to scaling. The benefit of approaching a scaling problem this way is that it allows us to account for both universal and non-universal effects in a consistent way. Writing down the functional form that is given by the data for a certain universality class will also be useful for identifying and characterizing other systems that are thought to belong to the same universality class.
We have found that, to analyze the windowed distributions, we need to first thoroughly examine the spatial structure of avalanches for the full system. In particular, in order to analyze the avalanche pieces left inside the window, we need to define height and width distributions and also joint distributions of sizes and widths. Section IV discusses these and also the results of fits for such distributions for the qKPZ model we are studying. To focus on the scaling region, and minimize lattice effects, we will discuss these distributions in terms of fractional area distributions, A(S) ∝ SP (S), the average fraction of the system that a given size takes up.
A remarkable result is that size, height, and width distributions can be fitted with a nearly identical functional form. For example the size distribution is: Here we have the shorthand S k = S/L 1+ζ k . The width A(w|L k ) and height A(h|L k ) distributions are similar, with the form of 1 is identical in form to the one quoted above in Eq 1.
Fitting the size, height and width distributions at once, we can extract multiple exponents-not only the commonly measured size distribution exponent τ , but also the exponent ν k which, as mentioned in point 1, relates the typical width of an avalanche to k, and also the roughness exponent ζ. The roughness exponent ζ, which measures the fluctuations of the interface, is typically quoted for front propagation models as a means of characterizing the universality class. (In our analysis, we have found a range for the roughness exponent ζ (from 0.62 to 0.72) that differs from the literature value 0.63 [13]; we discuss this in Appendix B.) One must note that the functional forms we choose are a practical tool to summarize existing information. While they may be inspired by analytical calculations, and chosen to be consistent with known asymptotics, they should be trusted only in the ranges over which they have been measured.
4. Multivariable scaling problems. The size distribution in Eq. 1 has a scaling form with one scaling variable. However, in this paper we will consider many scaling forms with more than one variable, such as a joint size and width distribution (IV), or the windowed distributions (section V and VI). In these cases, two or more scaling variables are important for describing the shape of the distribution (as seen in Figure 1(c)). For example, the general form for the 11 windowed distributions is: Here the scaling functions become distributions with two scaling variables, the rescaled size s k and rescaled window width W k . The traditional "scaling collapse" methods become problematic when multiple scaling variables are simultaneously important; this has hitherto retarded the effective study and use of these powerful universal joint distributions. We present a systematic method for analyzing scaling problems with multiple control variables. In our problem, the two control parameters are the demagnetization factor k, and the window width W . We will show that the interplay between k and W is important for determining the shape of the avalanche size distributions. In particular, we can write scaling functions for the distributions of the avalanche pieces in terms of these two scaling variables s k = s/L 1+ζ k and W k = W/L k (as seen in Figure 1(c)). For example, the scaling function for the 11 distribution is: Section VI discusses in detail the functions A 00 (s k , W k ), A 10 (s k , W k ), and A 11 (s k , W k ) that we consider. One may also find the results of fits (figures and tables of parameters) in section V.

SloppyScaling
The analysis in this paper is done in the software environment SloppyScaling and includes a Bayesian analysis of systematic error bars, which are explained in Appendix A 4. SloppyScaling allows us to fit data without collapses, which as mentioned above is problematic when there is more than one scaling variable involved. Included in the software setup are automatic fits of data to theory functions with nonlinear leastsquares and ease of visualizing results. This software may be applied to many different multivariable scaling problems, making full use of universality and the predictions of the renormalization group. All of the fits in this paper and their corresponding figures (including axis labels) were generated directly and automatically using SloppyScaling.

III. MODEL
We use a model for imbibition fronts to produce avalanches, it is extensively studied, with well-established critical exponents.
We simulate crackling noise using a quenched KPZ model in 1+1 dimensions (see Appendix A 1 for details on the implementation) [10], with dynamics given by: where h(x, t) is the height of the front, F the driving force increasing quasistatically, linear and non-linear terms for the KPZ model controlled by the parameters γ and λ respectively, and η Gaussian quenched disorder. In the spirit of magnetic avalanche models [26], we have added a term analogous to a demagnetization field −k h , which allows us to have pinned fronts at many metastable configurations. Models like these have been simulated [9] mostly near the depinning transition (and with k = 0). With k = 0, we define the area between each pinned front as an avalanche of size S. Avalanches produced by this model are thought to belong to the directed percolation depinning (DPD) universality class [13,27]. The avalanches are self-affine, long and wide, with ζ < 1. The roughness exponent ζ characterizes the ruggedness of the front, and also governs the scaling of avalanche heights h with widths w, h ∝ w ζ [9,28]. In our model, k controls (c) k = 10 −7 (k is what controls the typical avalanche sizes in a given system, giving a characteristic width L k .) Larger demagnetization fields stop avalanches more strongly, hence large k corresponds to smaller avalanches. The colored regions represent avalanches. The fronts are moving from bottom to top. Notice that the two simulations are statistically similar to one another apart from a rescaling of heights and widths. Note that most of the area is covered by the largest avalanches.
the typical size of the avalanche (Figure 2), the larger the k the smaller the typical size of an avalanche. We define the characteristic width of an avalanche in the full, unwindowed system to be L k = k −ν k . We will write all the scaling forms in terms of the length L k rather than the demagnetizing factor k directly to emphasize the analogies to finite-size scaling, as we are also studying finite window sizes compared to the size of L k .

IV. AVALANCHE SHAPES AND DISTRIBUTIONS OF THE FULL SYSTEM
In this section we will introduce various avalanche spatial distributions and their scaling forms. These forms will allow us to motivate the windowed distributions, and also serve as an example where traditional collapses may lead to questionable results.

A. Area-weighted size distributions
Traditionally, to describe an avalanche size distribution, we write the probability distribution as a power law times a universal scaling function. For example: L k is the characteristic width of an avalanche, and h ∼ L ζ k the typical largest height, and therefore the appropriate scaling variable to describe the area of an avalanche is S k = S/L 1+ζ k . When studying the spatial properties as in our case, the probability distribution P (S) (Eq. 5) is not the best choice, as its normalization is highly affected by nonuniversal effects at the lattice spacing ( Figure 3). In fact, for τ > 1, the normalization integral diverges at its lower (ultraviolet) limit a → 0, but not for L k → ∞. Although we could study scaling functions that include a lattice cutoff, it is more interesting to focus on the large (infrared) avalanche cutoff, which depends on L k in a universal way. To this end, it is more appropriate to make use of the first moment of P (S), and to consider the area-weighted size distribution A(S) A(S) ∼ S · P (S).
A(S)dS has a natural physical interpretation: it is the fraction of the full system area covered by avalanches with sizes between S and S + dS. Its scaling form is thus: We use the subscripts Sk to distinguish the scaling of the size distributions governed by k from those governed by other control variables. The power of L k we pull out of the scaling function is needed to normalize A(S) to one, since A(S) is sensitive to the long-distance cutoff.
In particular,  Table I.
where the last integral converges for τ < 2 [43] and becomes small as L k becomes large. Notice that the normalization of a power law must either diverge at the lower or at the upper limit. By studying avalanches weighted by their first moment, the normalization depends explicitly on L k , the infrared cutoff. This is the regime we are mainly interested in, since we would like to study the finite size effects imposed by both k and the window size W . Alternatively, as done in [29], one may also define a scale S m = S 2 2 S , redefine the sizes as S/S m and a corresponding p(S/S m ) which has normalization ∞ 0 dssp(s) = 1. Here p(s) is universal, but is also not a probability distribution in the conventional sense. This definition has an effect equivalent to what we do here: to make the function universal and insensitive to non-universal lattice effects on normalization. Namely, our definitions are related in the following way: A(S) = S S 2 m p(S/S m ) and A Sk (S/S m ) = (S/S m ) 2 p(S/S m ). Here we use S m ≡ S k = S/L 1+ζ k , which is consistent with their definition of S m up to a constant factor. We pre-fer to focus on the more directly interpretable fractional area distribution A(S).
Furthermore, notice that we have been unorthodox in writing the scaling form (8) for A(S|L k ) with a power of both L k and of S outside the scaling function. Normally one factors out a single variable from the scaling function. For example, one could in principle write In this form, A(S|L k ) dS is invariant under rescaling, and also clearly preserves normalization. However, B(S k ) = S 2−τ k A(S k ) vanishes as S k → 0, so this form of the scaling function disguises the power law behavior of the avalanche size distribution. By choosing the form A Sk which is defined to be finite and non-zero as S → 0, we make manifest both the avalanche size dependence and the system size dependence.

B. Maximum height and width distributions
In addition to the size distributions, we can also study the avalanche height and width distributions. We define height along the direction of front propagation, measuring the maximum height of an avalanche, and width perpendicular to heights, measuring the maximum width of an avalanche.
How do the height and width distributions scale? An avalanche of height h has size S ≈ h w ∼ h h 1/ζ = h (1+ζ)/ζ , so the system area A(h|L k )dh covered by avalanches with heights between h and h + dh scales as and since dS/dh ∼ h 1/ζ , where h k = h/L ζ k is the appropriate scaling variable to describe the avalanche height.
Similarly, we can write the scaling form for the width distributions as: where w k = w/L k is the scaling variable for widths. Notice the pattern in all these formulas: A(X|L) is a function with A(Y ) of scaling variables (combinations of X and L invariant under the renormalization group rescaling) multiplying a power of the scaling variable, divided by the independent variable X, making A(X|L)dX invariant under rescaling, allowing it to be universal.
where Table I). The first exponential term is the universal scaling function, while the second accounts for non-universal analytic corrections at small x due to lattice effects [30]. (See Appendix A 5).
One may ask why we choose this particular scaling form. When fitting data to a function, there are many parameterizations we could use to describe the data. This form is motivated from a functional renormalization group expansion by Le Doussal and Wiese for static avalanche size distributions in a linear model [31]. Our model differs in that there is a nonlinear KPZ term leading to anisotropy, so that our avalanches belong to a different universality class, the DPD universality class.
Le Doussal and Wiese find the avalanche size distribution for the linear model, for all static universality classes (random-bond, random-field, and random-periodic), to leading order in d = 4 − (where d is the dimension of the interface), to be: Here their scaling form includes the large scale cutoff S m . Le Doussal and Wiese claim that their results for both static and dynamic avalanches agree up to one loop for systems with ζ < 1. Static avalanches are separated by equilibrium configurations (or ground states), and dynamic avalanches are connected through a sequence of metastable states. Our avalanches result from a moving interface near the depinning transition, so they belong to a dynamic universality class. We thus use Eq. 14 as inspiration for the scaling form of Eq. 13. However, one may note that there is no theoretical basis that it should work, since it is from a distinctly different universality class.
One cannot determine the values of ν k and ζ independently, if we fit the size, height, or width distributions with Eqns. 8, 11, and 12 separately. For example in the size distribution we can only determine the combination ν k (1 + ζ). We determine the three critical exponents τ , ν k , and ζ, by jointly fitting the size, height and width distributions.
The results of our fits are reported in Table I and shown in Figures 3(b) and 4. In particular, we find ζ = 0.62 ± 0.02 which is close to the highly-precise value of ζ = 0.63 found in the literature [12,28,[32][33][34]. We note that the parameters δ x , which in principle control the aymptotic decay of the scaling function, are estimated here from a fit to the entire distribution. The quoted errors do not represent a confidence on the asymptotic decay -merely a confidence in the predictions over the range where data has been fit.  Table I. FIG. 5: Size distribution A(S|L k ) scaling collapse (color online) We collapse the size distribution data using the universal exponents of τ = 1.24, ν k = 0.45, and ζ = 0.62, the best fit values of the joint fit between A(S|L k ), A(h|L k ) and A(w|L k ).
Finally, the respective scaling collapses for the size, height, and width distributions are shown in Figures 5, 6, and 7. Although scaling collapses are very useful in verifying critical behavior, we argue that they may be problematic for the purpose of determining critical exponents, and one should fit and make use of functional forms. In Appendix A 5 we will show how scaling collapses are un- able to incorporate the effects of corrections to scaling, and how these corrections may cause a drift in the critical exponents.
parameter best fit standard errors systematic errors in linear approx.
Universal Exponents Shared Exponents Non-Universal Exponents The traditional standard error bars are typically ∼ 64 times smaller than the systematic error bars quoted; however, they are a gross underestimate of the actual errors expected since our theory is both highly nonlinear and sloppy [35]. We quote each parameter to the significant figure indicated by its standard error, since the parameters are strongly correlated, truncating each parameter to its significant figure would yield a poor fit.

D. Local height distributions
In our analysis of the next section, we will make use of another scaling function of the same form: not the maximum height of an avalanche, but the distribution of heights given by random cross sections of avalanches. Let A(h x |L k ) dh x be the fraction of the system area consisting of points (x, y) where the enclosing avalanche has vertical cross-sectional height at x in the range (h x , h x + dh x ), then, This distribution gives a different measure of the typical shape of an avalanche. The local height distribution is shown in figure 8. Here the curves show a fit we have generated with a scaling function of the fit to 11 spanning avalanches (Eq 32) which cross both window boundaries, taking the limit of W = 0. Our best measure of the local heights is equivalent to a 11 distribution with window size 1. Details for how this works are provided in the section V C.

E. Joint distributions and multivariable scaling functions
Once we have distributions for the measures of size S, width w and height h, we can also explore the forms of joint distributions. The area A(w, S|L k ) dS dw of avalanches in the range of size (S, S + dS) and widths (w, w + dw) will go to zero strongly if the size S becomes either much larger than or much smaller than the typical size w 1+ζ of an avalanche of width w -so we may factor out any combination of powers of w and S without changing the singularity. It still makes sense, though, to factor out the L k -dependence. If we choose to factor out powers of S, therefore, we find where again we have a power of a dimensionless scaling variable, times a scaling function, divided by Sw ∼ The last joint distribution that will be useful is related to right-most pieces of an avalanche. Consider the right-most piece of width x of an avalanche of total width w and total size S; let this segment have size s. (This will be the size measured by a window that cuts the avalanche at the left-hand window boundary at x.) Let A(s, w, S|x, L) be the fraction of the system covered by such avalanche pieces. Then, in the same logic as before, this five-variable distribution can be written as a power law times a universal function of four variables: One can clearly work out scaling forms for joint distributions of several variables and other combinations. The ones we have discussed here will be needed in our analysis of windowing effects.

V. WINDOW EFFECTS
Now that we have laid the groundwork for exploring the shapes of avalanches, we focus on analyzing avalanches inside a viewing window. In this section, we focus on how to define the right power laws and scaling; we also give results for fits, extracting critical exponents. In the next section, we go into more detail about the scaling shapes of these distributions -the universal scaling functions for avalanche sizes viewed through windows.
In imaging experiments one often runs into the problem of not being able to see the whole system, distorting the avalanche size distribution. In particular, for Barkhausen noise, typical magnetic avalanches span many decades in size, far beyond the spatial resolution of optical microscopes. The natural solution is to take measurements at a variety of magnifications. Even though at the weakest magnifications the window size W > L k and most avalanches avoid the window boundaries, the effects of the boundaries will always dominate at the highest magnifications. The analysis in this section not only provides a method to correct for finite-size-like window effects on exponents, but allows us to actively make use of all the data for a range of magnifications.
We show in detail in this section and the next how the avalanches which cross different boundaries exhibit distinctly different size distributions and critical exponents ( Figure 1). As described in section II, we consider the avalanches measured in an infinite strip of width W (Figure 1(b)), for a system with characteristic length L k . We separate avalanches into different categories: internal avalanches (00), split avalanches (10 and 01), and spanning avalanches (11). Let us call A 00 (s|W, L k ) ds the area fraction covered by such avalanches with sizes in the range s, s + ds. (For A 00 , the segment size equals the total size.) The split (10 or 01) avalanches will have area fraction A 10 (s|W, L k ) ds for each s. The distribution A 01 (Figure 1(f)) of avalanches touching the right boundary naturally equals A 10 on average. A 11 (s|W, L k ) ds is the fraction of the strip spanned by 11 spanning avalanches. We mentioned in section II that the 00 avalanches have a power law that matches the full system, whereas the 10, 01, 11 avalanches all have modified power laws with a smaller exponent τ .
Besides different power law scaling, the universal scal-ing functions for these different avalanche distributions are also distinct. In particular, the cutoff dependence on window size is different for internal avalanches and split avalanches, while the spanning avalanches have both an outer and an inner cutoff due to the window size (since avalanches must be large enough to span the window).
We present the fits of these universal scaling functions in this section and discuss their shapes in more detail in the next.
We know that all avalanches in the window are of one of the 00, 10, 01, 11 types, so As in the previous section, consider how each of these distributions A zz rescales under a coarse-graining by a factor b. The zz denote our indices for the various windowed distributions (00, 10, 01, 11). Each A zz ds, being a geometrical quantity (a fractional area), must be invariant under rescaling (with two invariant scaling variables): However, this is clearly not the form which makes the size and window-width dependence of the avalanche sizes manifest. We are allowed to factor powers of the invariant scaling variables W/L k and s/L 1+ζ k out of the scaling function B: The appropriate powers of s, L k , and W to pull outside depend upon which of the three distributions we are considering. For the distributions A 00 and A 10 , A 01 , we choose τ zz and υ zz (powers of the invariant scaling variables X = s/L 1+ζ k and Y = W/L k ) to make the resulting scaling function go to a constant at small X and/or Y . This way the power laws we pull out describe the behavior of the limit of small s, and the way in which the avalanches are cut off by the window size (as s approaches W 1+ζ ) are described by the scaling function. On the other hand for A 11 there are no small avalanches (they have to be large enough to span the window), and as L k → ∞, all avalanches span the window and become 11 avalanches, and in this limit the distribution will not go to zero, so for this distribution we instead pull out powers of W and s.

A. Internal Avalanches
First let us consider A 00 (s|W, L k )ds, the window area spanned by avalanches of sizes in (s, s + ds) that do not touch the boundaries. This can be computed explicitly from the function A(w, S|L k ) (Eq. 16) which gives the area covered by avalanches of width w and size S (note that for 00 avalanches, the segment pieces s = S): where (W − w)/W is the probability that an avalanche whose center lies within the window is entirely contained in the window (i.e., the avalanche center lies within (W − w)/2 of the center of the window). Changing variables from w to Ω = w/s 1/(1+ζ) , (22) with no explicit dependence on the window width W (so υ 00 = 0), and the same critical exponent τ 00 = τ that is given by the non-windowed distribution. Note that Eq. (22) is of the general form given by Eq. (19) and (20). This scaling equation is also consistent with our numerics: the normalization of the distribution for small avalanches is independent of W , and τ 00 = 1.26 ± 0.02 is consistent with the bulk τ = 1.24 ± 0.04.
Using this scaling form, we fit the 00 data (jointly with 11 and 10 data) with a scaling function given by a parameterized functional form: with s k = s/L ( k 1 + ζ). Our fit for the parameter n 00 is 1.62, which we believe to be 1+ζ (see next section). This makes the term s k /W n00 k = S/W 1+ζ which is another natural invariant scaling variable. Figure 9 shows the results of a nonlinear least squares fit, with shaded areas as estimations of fluctuations in the theory corresponding to systematic error bars on our parameters. In Eq. 23, τ , ζ, (and ν k which is hidden in L k ) are universal exponents shared amongst the three different distributions, A 00 1 and A 00 2 are (non-universal) analytic corrections to scaling reflecting lattice effects on small avalanches, and the other parameters encapsulate the shape of the universal scaling function A 00 . The fitted results for the other universal and non-universal parameters are quoted in Table II. We describe the scaling shapes and their motivation in more detail in section VI.
One may note the exponent ζ in our fits is fixed to the literature value of 0.63. If we allow for a free fit on all the parameters, it shifts to ζ = 0.68 ± 0.02. Although the free fit is 2.5σ away from the accepted value of ζ, the fit with fixed ζ = 0.63 has only a 50% higher cost than the free fit minimum, suggesting an average of 0.5σ drift on the parameters instead of 2.5σ as seen in Table II. This suggests three cautions (1) The estimate on our systematic error (0.02) is a lower bound estimate, and in fact the systematic error should be higher. (2) The fact that our systematic error should be higher also implies that the scaling functions are imperfect and may be improved upon. (3) There could be corrections due to a crossover that depends on both λ and k which we have not accounted for, which are distorting the fit. The subtleties and nuances in measuring the exponent ζ, and the possible origins of this drift are discussed in more detail in Appendix B.    Table II).

B. Split Avalanches
Next, consider the avalanches that are split by one side of the window, say the left side, with the distribution A 10 (s|W, L k ). Physically, for small avalanches s and large L k /W this is clearly proportional to 1/W : the small avalanches extend only a small distance into the window (smaller than the window width), so the fractional area covered by them is proportional to one over the width of the window. This leaves us with a scaling form with τ to be determined. Note again that Eq. (24) is of form Eq. (19).
We can also write A 10 in terms of the distribution of right-most pieces A(s, w, S|x, L k ) from Eq. 17, integrating over all possible sizes S, all possible widths w, and all possible pieces x (from lattice size a to window width W ): dx/W is the relative probability that the avalanche intersects the left-hand boundary, and we have changed the integration limits at ∞ to the avalanche length scale L k . (For w < x < W , the original distribution is naturally zero.) After we rewrite the integration variables in terms of the invariant scaling variables, we can organize the form of the scaling function into the form of Eq. 24. This tells us that τ = τ . These results are consistent with our numerical fits: W has an exponent of minus one, and τ is equal to the system τ .
With the correct power laws pulled out, we can now write down a function to describe the data and cutoff: Again, our best fit n 10 is nearly 1 + ζ, so s k /W n10 k ∼ s/W 1+ζ . Also note that this distribution has the same functional form as A 00 in eqn 23 aside from a factor of s 1/(1+ζ) /W in front. Figure 10 shows the results of a joint nonlinear least squares fit with the 00 and 11 avalanche data, with shaded areas representing estimations of fluctuations in the theory corresponding to systematic error bars on our parameters. Here, as in the 00 distributions, τ , ζ (and ν k which is included in L k and the scaling variables s k and W k ) are universal exponents shared amongst the three different distributions, A 10 1 and A 10 2 are (non-universal) analytic corrections to scaling reflecting lattice effects on small avalanches, while the other parameters describe the shape of the universal scaling function A 10 (s k , W k ). Fitted results for the other universal and non-universal parameters are quoted in Table II. We describe the scaling shape and its motivation in more detail in section VI.    Table II).

C. Spanning Avalanches
Finally, consider the spanning avalanche distribution A 11 . First, remember that most of the area in general is spanned by the largest avalanches (since τ < 2). Therefore, as L k → ∞, 100% of the area is covered by avalanches of widths much larger than W , and hence A 11 must integrate to one in this limit. This makes it natural to pull out only powers of W and s outside the scaling form for A 11 .
Second, notice that the size of an 11 avalanche is basically given by its height. More specifically, as W → 0, the avalanches have size s = W h x , where h x is the height of the vertical cross section of the avalanche. Hence we can write A 11 in the limit W → 0 in terms of the distribution A(h x |L k ) of randomly chosen vertical cross sections of avalanches (eqn 15), choosing h x = s/W : where dh x /ds = 1/W . Remembering from Eq.
x , we substitute s/W for h x and take the limit of W k → 0 in Eq. 27 to give: we cancel the two W 's and add the dependence on the second scaling variable W k = W/L k to derive the scaling form for A 11 : Here lim Y →0 A 11 (X, Y ) = A hxk (X) and thus A 11 (X, 0) dX = 1 (implied by the fact that almost all avalanche area touches both boundaries as W/L k → 0). Also notice that since: we can rewrite A 11 (s|W, L k ) as: Figure 11 shows the results of a joint fit with the simulation data of the previous 00 and 10 distributions. For the 11 distributions we use the functional form: τ (and also ζ and ν k which are hidden in the scaling variables s k and W k ) are universal exponents shared    Table II). .
amongst the three different distributions, A 1 11 is the (nonuniversal) analytic correction to scaling reflecting lattice effects on small avalanches, while the other parameters describe the shape of the universal scaling function A 11 (s k , W k ). Note that we don't include the term A 2 11 , as we have done in the 00 and 10 distributions; this term turns out to be the same as another term in the universal scaling function in the limit of W → 0, and so it is redundant. (See Appendix A 5). The best fit universal and non-universal parameters are given in Table II. We discuss in more detail the motivation and form of the scaling function A 11 (X, Y ) in section VI.
We also test the limiting case of our scaling function with our data sets of W = 1 in Figure 8 of section IV E. The curves drawn in Figure 8 are with the function given in equation 32, using the best fit values of the joint fit of A 11 , A 10 and A 00 . Notice that the predictions of equation 32 matches the data for W k = 0, indicating our function satisfies lim Y →0 A 11 (X, Y ) = A hxk (X) as expected.

VI. SCALING SHAPES AND RESULTS
In the previous section we wrote down scaling functions for each type of avalanche inside a window. In principle there are many possible parameterizations we can choose that would be able to capture the behavior of the data. In this section we explain how and why we chose each one, and also discuss the scaling function in the limit of small windows.

A. Scaling shapes and functional forms
We would like to capture the scaling behavior of both the finite size of the avalanches, and the effect of the window size on the distributions. We choose forms inspired by a functional renormalization group expansion for static avalanche size distributions for all universality classes [31], and further motivated by heuristic arguments for the cutoff dependence on the two scaling variables, s k = s/L 1+ζ k and W k = W/L k . Similar to in section IV C, we start with an avalanche size distribution of the form: We expect that as W → ∞, we will not see the effect of the window size, and the scaling forms will go to the limit of our proposed avalanche size distribution in Eq. 33. Keeping this in mind, we would like to write a function of A(s k , W k ), including the effect of the window size. For the 00 and 10 distributions, we expect the avalanche sizes to be cutoff by the window size W when W/L k < 1, hence a cutoff dependent on s k /W n k should be expected, where n > 0. While when W >> L k , their scaling forms should go to the limit of our proposed avalanche size distribution in Eq. 33. Therefore for the 00 and 10 avalanches, we propose a scaling form of: We have heuristic arguments for what this n zz value should take. Figure 12, a schematic of an avalanche cut by different windows, is meant to help illuminate our discussion. For the 00 internal avalanches, the largest avalanche contained fully within the window should have a width w that is roughly W . And since s max ∼ w max h max ∼ w 1+ζ max , then s max ∼ W 1+ζ and it follows that s k ∼ W 1+ζ k , giving us an expectation value   Here are the results of our joint fit for the windowed A00, A11, A10 distributions. The corresponding scaling forms which were fit are quoted in the table alongside the parameter results. Here systematic error bars which account for errors in the theory (see Appendix A 4 for explanation) are given. The traditional standard error bars are typically ∼ 30 times smaller than the systematic error bars quoted. The last column is the drift in parameters seen when allowing ζ to be a free parameter. Notice that these numbers are more or less similar or smaller than the estimated systematic error, except for ζ. (The problems in measuring ζ are discussed in Appendix B) As in Table I, we quote several digits more than the error bars warrant for individual parameters, because the errors are strongly correlated; truncating each parameter to its significant figures would yield a poor fit. Here we are depicting cases where this avalanche (or its segment) is the maximum avalanche size for the 00, 10, and 11 at different window widths. Boxes of different widths and colors are used to show the cases in which this may happen. The main avalanche may count as a 00 avalanche for a wide window, while part of it would count as a 11 avalanche for a smaller window; it could also count as a 10 avalanche if it happens to cross the window boundary. This figure illustrates our arguments for the shape of the cutoff (the exponent nzz) given a window size W for the 00, 11, 10 cases. Another small avalanche is drawn for the 11 case to show that the minimum size to cross the 11 window also introduces a separate cutoff.
of n 00 = 1 + ζ. Numerically, we find n 00 ∼ 1.62 when we fix ζ = 0.63. The fit plotted against one of the scaling variables s k and also the contour plot of the scaling shape is shown in Figure 13.
For the 10 or 01 spilt avalanches, since we are effectively measuring the ends of avalanches that spill into the window, n 10 depends on what the shape of the avalanche is at the edges. The largest portion of an avalanche to spill into the window will again be limited by the size of the window W . Here the shape follows the roughness of the two fronts preceding and following the avalanche, where h(x) ∼ x ζ for each, so plausibly h edge = h af ter − h bef ore ∼ W ζ . The size is then limited by s max ∼ w max h edge ∼ W 1+ζ , giving us an expectation value of n 10 = 1.63. Numerically, we find n 10 ∼ 1.62 (ζ is estimated from 0.62 − 0.72 in our various measures), matching our expectation. The fit plotted against one of the scaling variables s k and also the contour plot of the scaling shape is shown in Figure 14.
Now we move on to discuss the 11 spanning avalanches. Here the situation is slightly more complicated than the previous two cases, due to the distribution being strongly cut off at two length scales, as one may note from the Each contour reflects a drop of a factor of e in the scaling function. The black dots are at locations of the simulated data points used in the fit indicating where the fit should be a reliable prediction. The red solid line is log 10 W k = n00 log 10 s k which is the slope at the large avalanche cutoff, with n00 = 1.62, the best fit value. (n00 = 1 + ζ = 1.63 is the expected value from our heuristic arguments).
shapes of the distributions shown in Figure 11. First of all the avalanches need to be large enough to cross the window, implying an inner cutoff that depends on W 1+ζ k /s k (i.e. the cutoff is for s min /s 1 and so s min /s = W 1+ζ k /s k ); here the argument for the minimum size 11 avalanches follows from a similar argument for the maximum size 00 avalanches. The smallest avalanche that is able to span the window will have a width w min = W , whereas h min ∼ w ζ min , and s min ∼ h min w min ∼ W 1+ζ . However, in this case s k should be in the denominator of the scaling variable, since for s max >> s > s min the probability of the having a spanning avalanche grows as s increases. For the outer cutoff we expect that the maximum size is given by the window size W and the typical maximum height, i.e. s max ≈ W h max which implies that s k max ≈ W k max (h max /L ζ k ). This rescaled height (h max /L ζ k ) is constant. Therefore the cutoff for the large size avalanches should depend on s k /W k . Hence, we propose the scaling form below for the 11 spanning avalanches: where C 11 controls the strength of the inner cutoff and D 11 the outer cutoff. Figure 15(a) shows the shape of the scaling function plotted against one scaling variable s k , and Figure 15(b) gives the contour plot of this function against both variables. The best fit value of n 11 is n 11 ∼ 1.65, whereas the expected was n 11 = 1 + ζ = 1.63.
Finally, one may note that in our system all the n zz turn out to fit to our expected 1 + ζ within the error bars of ζ. One may be tempted to set n zz as 1 + ζ and have fewer fit parameters in one's form; however, we recognize that our geometrical arguments do not hold for front propagation that is super-rough with ζ > 1, or for models that allow overhangs, so we leave n zz as a free parameter to signify this geometrical constraint.

B. The limit of small windows
Although there are noticeable imperfections in the theory function, the agreement is impressive between theory and simulation as seen in figures 13, 14, and 15 .
The scaling function for each of the three distributions is a competition between two types of terms: the rescaled size s k = S/L 1+ζ k and the rescaled window size W k = W/L k . Upon examining the fit parameters of the scaling function, all three distributions share the characteristic that at small W k the terms with pure powers of s 1/2 k and s δzz k become unimportant, leaving only s k /W nzz k for the 00 and 10 distributions, and for the 11 a s k /W k term. Notice that for the 00 and 10 distributions, since n zz is 1 + ζ, the L k dependence disappears for the universal scaling function at small W k . Therefore, the shape of the scaling function is cut off mainly by the window size W . In fact removing the s k terms for these functions does The upper orange solid line is the slope of the contour plot at the small avalanche cutoff, and has log 10 W k = 1/1.65 log 10 s k where 1.65 is the best fit n11 value. The orange dashed has log W k = 1/1.63 log 10 s k where 1.63 is the 1 + ζ value. The lower red solid line is log 10 W k = log 10 s k which is the slope at the large avalanche cutoff.
not affect the shape for W k < 5. For experiments that study systems in the same universality class as this one, this implies that data may be measured at large magnifications (small windows) and fit to extract exponents and scaling behavior without the extra, often unknown, scale of L k .

VII. SUGGESTIONS FOR EXPERIMENTS AND CONCLUSIONS
What does our analysis imply for current experiments? How should one conduct the experiment and analyze the data? Here we discuss for the particular case of magnetic avalanches in Barkhausen noise, how to take into account window effects and further enhance the collection of data.
There are two optical methods for detecting avalanche distributions for Barkhausen noise in 2D thin films, and both make use of the magneto-optical Kerr effect (MOKE). When a polarized beam of light reflects off a magnetized sample, the reflected polarization is affected depending on the magnetization. A second polarizer can be used to filter this signal, and then using either a photodiode [16], or an optical microscope [5,6], we can collect data about the avalanches from the signal. For experiments using a photodiode (let us call this "laser reflectometry"), one can only measure the total magnetization change over time, and not individual avalanches. For experiments using an optical microscope (let us call this "avalanche visualization"), one can resolve individual avalanches and their shape.
In laser reflectometry experiments [16] we only have information for the magnitude of magnetization as a function of time, and cannot see which avalanches touch the boundary. Furthermore, in current techniques the laser spots are Gaussian in shape, and do not have sharp boundaries. However, there seems to be no fundamental reason why the illuminated region could not be optically generated with uniform illumination and sharp edges, up to some diffraction limit depending upon the geometry of the experiment. (A typical avalanche of interest is a few microns in size, large compared to the wavelength of optical light which is 400-700 nm). If one could make the edges of the laser spots sharper, one could adjust the laser spots to flicker between two sizes, one with a radius slightly larger than the other. Events that occur with the same magnitude in both the large size measurement and the smaller size measurement would be 00 internal avalanches. More elaborate sequences of spot shapes could be used to further distinguish 01 and 10 split avalanches from 11 spanning avalanches.
For avalanche visualization experiments [5,6,22], it is straightforward to separate the data into 00, 10, 01, and 11 avalanches for systems which have our strip geometry and flat fronts with ζ < 1. Our analysis will remain valid with minor corrections due to real-world experimental circumstances. For example, sometimes the propagation direction of the front is not parallel to the top and bottom boundaries. In this case there would be corrections depending on the angle of the tilt θ and the size of the window. The local heights h x would need to be adjusted with the factor of cos θ. The 10, and 01 avalanches would be cut off at an angle, but for self-affine (short, wide) avalanches these size corrections are likely irrelevant.
For materials with dipolar interactions and zig-zag shaped fronts [5,36], we are less confident that our methods can be applied without modification. The large vertical extent of the zig-zag front suggests that all four boundaries of the window will matter; therefore we will need to divide avalanches into more categories (0000, 1000, 0100, 0010, 1010...). The analysis for these types of avalanches will be more complicated. This would be an interesting problem to pursue by simultaneous analysis of simulations and the experimental data.
Knowing how to deal with window effects can be an important tool for these visualization experiments. By combining data at several magnifications (corresponding to different window sizes), we resolve a larger range of length scales. Higher magnifications will show the small avalanches, while lower magnifications will allow us to both capture larger avalanches and explore more fully the 00 internal avalanches. For example, if our CCD camera recording the images has a resolution of 1000 2 pixels, and we have a magnification of up to 5x-50x, we can simultaneously explore a range of window sizes and extend our effective spatial resolution from 1000 2 to 10000 2 .
More generally, this paper has provided the tools needed to extract from the experimental data for systems of similar interface dynamics the critical exponents τ , ζ and the universal scaling functions A Sk , A hk , etc. For these experiments, we can also measure widths and heights and the average shape, giving us an independent measure of ζ, and the universal scaling function A hxk . Our detailed analysis and comprehensive methods presented in this paper enable a more powerful interpretation of current experiments, and improved construction of future experiments.
Appendix A: Methods and software

Numerical simulation
The avalanche simulations in this paper were produced using a quenched KPZ model [10,11], with dynamics given by: where h(x, t) is the height of the front, F the driving force, k the "demagnetization field", linear and nonlinear terms for the KPZ model controlled by the parameters γ and λ respectively, and η gaussian random noise. This was run for system sizes of width L 4096, 8192, and 16384. The simulations have been run in a strip geometry (4096 × 8192, 8192 × 16384, and 16384 × 32768) and the bottom half of the simulations have been truncated to avoid transient effects due to the initially flat front. The left-right boundaries have periodic boundary conditions. The simulations are done using a discrete cellular automaton model, in which the displacement of the string h, the time t and the space x are all discretized and take integer values [13,37]. For a configuration {h i }, we compute the local force F i at each site i, leading to a discretized version of Eq. A1 where the sum runs over all the nearest neighbors nn for the site i, a is the discretization length that we set to 1, and η i (h i ) is a random force. The automaton dynamics are as follows: (1) increment the external field until one site is unstable (F i > 0); (2) determine for each site along the interface whether it is stable (F i < 0) or unstable (F i > 0); (3) advance all unstable sites by one step h i = h i + a = h i + 1 in parallel, generating a new value of the pinning force η i (h i ); (4) repeat until no sites are unstable (the end of the avalanche); (5) repeat (1-4) until the front passes the top of the simulated window.

Nonlinear least squares fitting
We use nonlinear least squares methods for fitting data to theory functions, minimizing a cost defined as: Here θ are the parameters, y the function value, and σ i the error on the data points. The weight σ in our case is determined by fluctuations from run to run. Namely, we bin the data (in equispaced log bins) for each of N simulations, calculate the standard deviation of this value across runs, divide by √ N − 1 to get the fluctuation in the mean of that bin. [44] In our distributions, small avalanches occur more often (leading to small error bars), but large avalanches are more important, and the smallest avalanches suffer from non-universal lattice effects. So during the fitting process there is a tradeoff between fitting the region where there is good data and where the variations are most important. We use a number of methods to compensate for this imbalance. (1) We set a minimum error bar (1% of the data value) on the data points, making the error bars on smaller avalanches larger, and therefore decreasing their weight. [45] (2) Analytic and singular corrections to scaling can also account for non-universal effects. We include analytic corrections to scaling for the lattice effects in our scaling functions. These corrections appear in all distributions we discuss. A more detailed discussion is in section E of this appendix. (3) One may also skip points that have non-universal behavior when fitting. For our fits in this paper, all points are included.
Another issue arises in regions where one has sparse data; there may be bins that do not have any observations. For these zeros, the error bar should not be zero! We can use maximum likelihood methods to estimate theoretical errors. Say we have N experiments, and bin the data with L i sizes in each bin i. With the median size in the bin S i , the expected probability for each size is ρ i = A(S i )/S i , where A(S) is the theoretical distribution of sizes, and so getting a zero in one of the bins will have the probability p = (1 − ρ i ) N Li . Using maximum likelihood, this results in a cost of C = − log p = −N L i log(1 − ρ i ) ∼ N L i ρ given that ρ is small. The residual we add to the total cost is therefore: This calculation is generalizable to sparse data that is non-zero-say one has n events in a bin for N measurements, then the probability is p = i , and the appropriate error bar we get is approximately σ = 1 N Li . One could then use the larger of the error bar given by this argument or the statistical error bars from the simulation. In practice, for our simulations, we find it sufficient to use the statistical error bars but to compensate with minimum error bars given by Eq. A4.
In non-linear least squares fitting, we can also include priors in the cost if we have assumptions or information a priori about the parameters. For example, in our problem, we put priors on the exponents inside the scaling functions (n 00 , δ 00 , etc) to prevent them from going to large values and forcing their corresponding coefficients to zero. As a result our cost function now becomes: (A5) In Eq. A5, n represents all the arbitrary exponents that occur inside the scaling functions (n zz , m zz , δ zz etc...). log(δ 2 11 ) is included to prevent δ 11 from going to zero. A clear understanding of these techniques for fitting is important for increasing the reliability of our results, and acknowledging its limitations. In a following subsection, section A 4, we will further discuss how to estimate the reliability of results inferred by fitting data to a theory, generating systematic error bars for fitting results.

Software for fitting
To facilitate the exploration of this problem, we have developed a software environment, SloppyScaling, in Python.
This code is downloadable at http://www.lassp.cornell.edu/ sethna/Sloppy/SloppyScaling/SloppyScaling.htm. The main features of this code include various nonlinear-least squares fitting methods [38], automatic plotting for visualization, and methods for generating systematic error bars on the theory.

Systematic Error Bars
We have quoted in our results systematic error bars instead of the more commonly used standard error bars in the parameters. Standard errors given by the covariance matrix are expected to be erroneous for our problem, since our problem is highly nonlinear in the parameters, and also sloppy-parameter combinations in the sloppiest directions can vary an infinite amount without affecting the fit. We use a method due to Frederiksen et. al [35] for Bayesian estimation of errors. This method involves assuming that given a theory (M) which is imperfect, a spread of parameters (each corresponding to different models) may fit the data(D) in an equally acceptable matter. We can define a probability of a certain model with where C(θ) is the cost at a given set of parameters θ, and the effective temperature T sets a scale for the fluctuations away form the best fit. Since the cost at the best fit parameters C bf is a measure of how well the theory is doing, we choose T = 2C bf /N , where N is the number of parameters with "equipartition" allowing each degree of freedom a contribution of 1 2 T to the total cost. Ideally, after defining such a probability, one should sample parameter space to determine the systematic error bars on parameters. However, in our functions sampling is non-trivial due to parameter evaporation [39], the "entropy" for parameters drifting to infinity overwhelms their cost in degrading the resulting fits. Therefore, we make a quadratic estimate of the fluctuations in predictions, essentially using propagation of error to estimate the systematic error. The covariance matrix gives an error σ stat that assumes the temperature of 1, corresponding to P (θ|D, M ) = exp(−C(θ)). Using propagation of error, we calculate the systematic error according to our effective temperature T = 2C bf /N : The shaded plots are generated by sampling according to the Hessian at the best fit, weighting our steps in each eigendirection by the inverse square root of the eigenvalue, and scaling the steps with a low temperature (T L ). Then for our ensemble of parameters we calculate the fluctuations in the theory (residuals δr ens ) corresponding to the ensemble. We scale up these fluctuations according to the temperature defined by the best fit, or We have also estimated systematic errors by removing large k curves (which we believe have larger corrections), and looking at the corresponding drift in exponents and parameters. The estimate of the systematic error that this procedure gives is often similar (sometimes smaller) to the one using our temperature-scaled propagation of error.

Corrections to scaling
Corrections to scaling play an important role in our scaling functions, their inclusion accounts for nonuniversal effects, and helps increase the reliability of universal predictions. In each of the scaling functions in this paper we have included analytic corrections to scaling that capture the lattice effects of our automata simulations. They are of the form: This expansion for small S corrects for the lattice effects on small size avalanches which we believe to be present in the distributions. One can imagine that experiments may have other origins of non-universal effects, such as a nonlinear signal amplifier or distorting lens; one should always attempt to account for and include these [30]. One caveat is that, just as adding extra free parameters does not necessarily increase the quality of one's fits, adding corrections is not a guarantee for increasing the accuracy of one's scaling function. One should be careful in checking that the terms included in the corrections to scaling behave as expected in the region of interest, are subdominant when taking the appropriate limits, and do not confuse the main universal scaling function, either by canceling out terms or having the same effect. An example of this complication is seen in our studies of using the limit of the 11 distribution at W = 0 for the local height distributions.
In sections IV D and V C, we take the view point of using the parameterized from of A 11 (s|L k , W = 0) as the proper limit of A(h x |L k ). However, this data should be matched as well by A 11 (s|L k , W = 1), since in effect this is what we are fitting for the local height distributions. We can view the ratio between the two functions as a multiplicative correction to scaling from lattice effects (without considering the corrections to scaling for lattice effects in equation A9): Here we define h s = h x /L ζ k . Note however that for the term C 11 L (1+ζ−n11)m2 k h −m2 x , the value of the fit for n 11 equals 1.65, which is within the range of error for 1 + ζ = 1.63±0.02. Therefore this term is nearly C 11 h −m2 x , where m 2 = 1.61 ± 0.6. This term then has the same effect as the correction to scaling term A 11 2 /s 2 at small W (since s ∼ hW for 11 avalanches at small W ). Since A 11 2 /s 2 is only significant in the range of small s, which only occur in the 11 distributions at small W , these two parameters serve the same purpose, and it is redundant to include both for the fits. We therefore remove the term A 11 2 /s 2 for the 11 distributions for our fits.
One can also check the corrections and see all powers of L k are negative, and the multiplicative correction approaches unity as we get closer to the critical point. With our current fit, the term with L −1/2 k dominates the corrections. Notice that, originally, to account for lattice effects, we have added corrections in integer powers of h x (exp(A 11 1 /h x + A 11 2 /h 2 x )), which are subdominant compared to L −1/2 k . This implies that there are more dangerous corrections to scaling than originally inferred. By this method, one might check if the corrections originally included are sufficient-in this case they are close, but could be improved upon by systematically adding similar terms.
In a continuum case, the concept of local heights h x should describe a smooth shape tracing the depinning line, in our automata, since the smallest width is naturally 1 lattice spacing, the smallest possible window width is W = 1. In measuring the local heights, this discreteness limits the smoothness of the shape of h x , and gives rise to the corrections we see. Here we've seen that if we "know" from other measurements (in our case the 11 spanning avalanches) the right limit the universal scaling form should take, we may find the form of the corrections.

Scaling collapses and their limitations
We have argued at various points in the paper that scaling collapses are limited and may lead to questionable results. Here we will illustrate an example of this. Using the critical exponents given by a free fit of the windowed distributions (ζ allowed to vary), we collapse the sizes, heights and widths. Comparing the figures included here and the ones in Section IV ( Figures 5, 6, and 7.), one can see that the collapses are of similar quality.
The fact that such distinct values of ζ can yield similar quality collapses may imply (1) our "systematic error" bars on ζ, estimated to be ±0.02, are in reality much larger, (2) collapses do not incorporate non-universal corrections to scaling, and these may have an important effect, (3) collapses are not a reliable way of verifying the values of critical exponents. In particular, we expect corrections to scaling due to large k to be responsible for the drift in exponents. Here only the combination of ν k ζ affect the scaling collapse, the large shifts in ν k and ζ mostly cancel in the product, yielding similar collapses. Comparing this with the collapse shown in Figure 6, we see that the two collapses are comparable in quality, where in Figure 6 the large avalanche cutoff is collapsed nicely, and here the smaller avalanches are collapsed better.
With these software tools and analytical methods, data at critical points may be analyzed while including multiple scaling variables, allowing for the treatment of a broad range of experiments, and also allowing for a far more rigorous estimation of statistical and systematic er-  Figure 7. Here only the combination of ν k affects the scaling collapse and (τ − 2)(1 + ζ) affects the shape of the scaling collapse. Since ν k does not change significantly, the quality of the collapses are similar.
rors. By using automatic fits to entire scaling functions, instead of traditional collapses, and by estimating systematic error bars, we facilitate the interpretation of data with multiple scaling variables and analytic corrections to scaling. This advance will allow for better characterization not only of noise in magnetic thin films and similar avalanche dynamics, but should be broadly applicable to all applications of critical phenomena and scaling theories to experiments and simulations.

Appendix B: The roughness exponent ζ
In the investigations presented in this manuscript, we have found that the estimates of the critical exponent ζ have been problematic. In this appendix, we will discuss various means of measuring this exponent, the significance of the range of values we find from various measurement methods, possible origins of this range, and implications for future research. We emphasize that any value of ζ in the range we observe (0.62 ± 0.02 to 0.72 ± 0.02) can describe all of our data essentially equivalently well.
The shape of the front has been studied as an identifying feature for front propagation models, which is usually characterized by defining a roughness exponent ζ, which is measured through a height-height correlation function: The quenched KPZ model we use is conjectured to belong to the Directed Percolation Depinning (DPD) universality class [9], which is conjectured in turn to belong to the Directed Percolation(DP) universality class. For a pinned interface in DPD, the roughness exponent ζ = 0.63 ± 0.01 [10,11] matches that of the ratio of correlation length exponents in DP ν ⊥ /ν = 0.6326±0.0002 [40]. (Here ξ ⊥ ∼ |p − p c | ν ⊥ and ξ ∼ |p − p c | ν , where p is the branching probability, and p c is the percolation threshold.). One may note that for a moving interface, the picture is less well known; there has been numerical study that shows ζ = 0.70 [12,41], but also arguments that the interface under this condition is not self-affine [13], that the moving regions have ζ = 1 and the pinned regions have ζ = 0.63. One could imagine that our demagnetizing force, like the velocity in DPD, could lead to a heterogeneous mixture of different scaling regions, converging to ζ = 0.63 as k = 0. Figure 19 shows measurements of the height-height correlations in our model. Using finite-size scaling for a numerical fit, we see that as k is tuned away from zero, ζ falls between 0.64 and 0.68, increasing with k. Measuring the local log-slope, one can see clearly a drift in the measured exponent in Figure 20.
The range of our estimates (varying from ζ = 0.63 ± 0.02 to ζ = 0.72 ± 0.02) is large compared to our error estimates; however, best fits with ζ fixed within this range had costs within 1.5 times that of the best fit value, indicating both that our quadratic estimates for the systematic errors are too small and that it may be challenging to definitively measure ζ in either simulation or experiment.
In our fits we find that ζ = 0.62 ± 0.02 for the size, width and height distributions joint fit, and ζ = 0.68 ± 0.02 for the 00, 10, 11 joint distributions. Notice that although the direction of front propagation is in general along the y-axis in our problem, portions of the front will be at various angles to the y-axis. Since the local direction of the front propagation is not fixed, we can also choose a rotationally invariant definition of height and width: defining the height and width of an avalanche along the axes of the moment of inertia tensor. We define root-mean-square heights and widths as the square root of its eigenvalues. So if we fit rms heights and widths jointly, we get ζ = 0.72 ± 0.02 much higher than 0.62 measured along the global axes. The difference in these two exponents seems to indicate that the local avalanche shape has a different geometry than the global avalanche front. Many front propagation models spontaneously break rotational symmetry through the orientation of the front. (Envision a circular front growing from a point, with differing front orientations.) Note, however, that the qKPZ dynamics is anisotropic, breaking rotational invariance. We have also looked at fits with the windows scaling functions involving subsets of simulations with different k. Using smaller values of k generally lead to fits of ζ closer to 0.63. This may point to imperfections in our theory function (are we missing some of the scaling behavior dealing with k?), or corrections to scaling (analytic or singular). One possibility is that there is a residual crossover effect having to do with another relevant variable. In the anisotropic form of the qKPZ model, the nonlinear term (Eq. 4) λ is a relevant variable, and is non-zero under renormalization [3,9,42], and although we simulate the model at fixed λ, we have observed there is a crossover effect following a direction having to do with both k and λ. We note that without the nonlinear KPZ term the qKPZ model becomes the quenched Edwards-Wilkinson (qEW) model, and that for the qEW ζ > 1. Literature suggests there also may be a crossover effect due to a runaway fixed point [42]. We think these last two possibilities can be explored with more simulations done on different λ and k, examining a crossover to the linear version of the qKPZ model (the quenched Edwards-Wilkinson), to make a more complete picture of the phase space.  (2) of the height-height correlation function; this is a measure for the roughness exponents for various simulation sizes L and k. The results of the local-log slope ln[C(r)/C(r/2)]/ ln (2) is consistent with what is seen in Figure 19. The lower red-dashed line shows ζ = 0.63, corresponding to directed-percolation depinning (DPD) and suggested by literature to be the correct value for our model. Whereas the upper black-dashed line shows ζ = 0.68 which is the result of our fits of windowed avalanche distributions. In the blowup of the region of r/2 = 10 − 100, we can see there is a trend of larger ζ corresponding to larger k simulations.