Influence of inhomogeneous stochasticity on the falsifiability of mean-field theories and examples from accretion disc modeling

Despite spatial and temporal fluctuations in turbulent astrophysical systems, mean-field theories can be used to describe their secular evolution. However, observations taken over time scales much shorter than dynamical time scales capture a system in a single state of its turbulence ensemble. Comparing with mean-field theory can falsify the latter only if the theory is additionally supplied with a quantified precision. The central limit theorem provides appropriate estimates to the precision only when fluctuations contribute linearly to an observable and with constant coherent scales. Here we introduce an error propagation formula that relaxes both limitations, allowing for nonlinear functional forms of observables and inhomogeneous coherent scales and amplitudes of fluctuations. The method is exemplified in the context of accretion disc theories, where inhomogeneous fluctuations in the surface temperature are propagated to the disc emission spectrum--the latter being a nonlinear and non-local function of the former. The derived precision depends non-monotonically on emission frequency. Using the same method, we investigate how binned spectral fluctuations in telescope data change with the spectral resolving power. We discuss the broader implications for falsifiability of a mean-field theory.


INTRODUCTION
Mean-field approaches have been widely applied in theoretical astrophysics including for models of stellar structure, convection (Cox & Giuli 1968), accretion discs (Shakura & Sunyaev 1973), magnetic dynamos (Roberts & Soward 1975;Krause & Raedler 1980), and interpreting large-scale structure in the Universe. In such approaches, physical quantities are typically decomposed into mean and fluctuating parts, and the former is explicitly solved for while the latter is modeled using closure methods. To facilitate an analytically tractable theory, an infinite scale separation between the mean and fluctuating fields is often assumed, allowing for Reynolds averaging rules: linearity of averaging, interchangeability between averaging and differentiation, and invariance of mean quantities under subsequent averaging (see, e.g., Zhou et al. 2018).
In contrast to this idealized assumption however, realistic astrophysical turbulent flows often exhibit finite and moderate scale separations. The ratio between the system and the turbulence scales is ∼ O(10) to ∼ O(100) for geometrically precision error bars of the theory, and a better assessment of whether they are stochastic or systematic can be deduced.
The question of computing precision of mean-field predictions can be mathematically formulated as follows: Let A(τ ) be some field with (i) mean A to predicted by some theory, (ii) fluctuations δA that modeled statistically; and τ is a set of formal variables representing spatial coordinates, time, or any other parameters that A may depend on. If K is an observable related to A via where V is some appropriate interval and k is some function of A, what is the fluctuation in K due to that of A? Zhou et al. (2018) focused on the context where A is a galactic magnetic field, K is the Faraday rotation, and the integration is carried out over a line of sight.In this case, k(A) is the density-weighted magnetic field along the line of sight, and thus a linear function of the mean magnetic field. The coherence length and the variance of the fluctuations were taken to be independent of τ , the position along the line of sight. The integral of k above is then replaced by a discrete summation and the central limit theorem readily applies. In the ensemble of different realizations of turbulence, the variance of K is then simply N times smaller than that of k where N is the number of turbulent eddies along the line of sight.
In the present work, we generalize the method to include cases where k(A) is a nonlinear function of A, and both the variance and the coherent scales of the fluctuations depend on τ . We apply the formalism to accretion disc theories. For a wide range of accretion models, angular momentum transport is modeled by solving for the dynamics of the mean velocity field with angular momentum transport mediated by an imposed mean stress. These models are manifestly mean-field theories, and predictions from these models should therefore be presented with precision error bars.
Our method herein enables efficiently quantifying fluctuations about the mean-field predictions for the disc spectrum at all photon frequencies, thereby providing a quantitative way of accessing the precision of the models considered. This is accomplished by propagating local disc perturbations with given amplitudes and coherent scales to the emission spectrum. At a given photon frequency, this requires including contributions from the entire disc, superseding the local treatment of Blackman et al. (2010). Quantifying the theoretical precision becomes particularly important when comparing to snapshot observations of a source because then the observed fluctuations cannot be identified by their temporal characteristics.
In Section 2 we introduce the needed method for computing fluctuations of an integral from local contributions with inhomogeneous amplitudes and coherent scales. In Section 3 we apply the method to a generic thermal disc spectrum, and more specifically to protoplanetary discs in Section 4, and dwarf novae in 5. In Section 6 we consider how the derived imprecision should be binned when taking into account comparison to observations with finite telescope resolving power. We summarize in Section 7.

Propagating fluctuations from integrand to integral
To demonstrate how to propagate fluctuations to an integrated quantity, suppose K represents the total luminosity from a one dimensional object in a Cartesian geometry. If this is given by k(r) per unit length at position r, then where a and b specify the boundary. Now consider a local luminosity fluctuation δk(r) with a vanishing mean but allow for inhomogeneous variance σ 2 k (r), i.e., δk 2 (r) = σ 2 k (r), where the angle brackets denote an ensemble average over all possible realizations of the fluctuation. We take the fluctuations in the neighborhood of position r to be coherent over a scale l(r), and assume that l is much smaller than the local variation scales of k, σ 2 k , and l itself; that is l ≪ |k/(dk/dr)|, l ≪ |σ 2 k /(dσ 2 k /dr)|, and l ≪ |l/(dl/dr)|. In the simplest case, for which both σ 2 k and l are independent of position, the total luminosity becomes the sum of N = (b − a)/l number of independent and identically distributed random variables, and thus the variance of K is according to the central limit theorem. We now generalize to cases where both σ 2 k and l are smooth functions of r. We provide two methods which give identical results. We use the first method here and a second is given in Appendix A. We discretize the system into grids with the size of the ith cell being li ≡ l(ri), and the grid can be constructed in the following way: r1 = a; ri+1 = ri + li; i = 1, 2, · · · .
(4) K is given by the sum of contributions from fluctuating cells, where ki = k(ri). Because l(r) ≪ |k/(dk/dr)|, li can be the line element. The mean of K 2 is Since each cell covers a local coherent length and varies independently, we have Thus where the last step we have used the condition l(r) ≪ |σ 2 k /(dσ 2 k /dr)|. The variance of K is thus We test the analytical formula (9) against stochastically generated data using an example where we have σ 2 k (r) and l(r). Specifically we take Equation (38) for an accretion disc model, which is derived and explained in more detail in Section 5, but here this simply motivates use of the mathematical example σ 2 k (r) = r 37/8 e 0.2r 3/4 e 0.1r 3/4 − 1 4 ; l(r) = l0r 9/8 with two choices l0 = 0.03 or l0 = 0.05. Here length scales are in dimensionless units, normalized by the disc inner radius. The form of σ 2 k (r) is derived based on a black-body spectrum with a mean disc surface temperature ∝ r −3/4 , with r being the non-dimensional disc radius.
Using Equation (10), we compare the prediction of σ 2 K from Equation (9) with that derived from numerically generated data sets. For the latter, we discretize the interval [a, b] and assign each mesh point a random number (representing k) according to a multivariate Gaussian probability distribution function (PDF) with a zero mean and a covariance matrix where ri,j are the coordinates of a pair of cells on the grid. For each realization of k, its values on all mesh points are summed to give K, and different realizations of k construct an ensemble, over which the variance of K is then calculated. With a = 1 and different values of b, the results are plotted in Figure 1. Importantly, at a given b the analytical formula is on average ∼ 2000 times faster than averaging over an ensemble with 512 members 1 , and provides an efficient way of estimating σ 2 K . Comparing the l0 = 0.03 and l0 = 0.05 cases, it is evident that the larger deviation between data points and the theoretical curve at larger b is due to the breakdown of the condition l ≪ |σ 2 K /(dσ 2 K /dr)|. In fact the two length scales becomes comparable at b 600 with l0 = 0.03, and b 400 with l0 = 0.05. In general, when the correlation length becomes large enough to violate the aforementioned condition, Equation (9) delivers less accuate results. In the current case and later examples, the main contribution to σ 2 K comes from the region where l ≪ |σ 2 K /(dσ 2 K /dr)| still holds, and therefore yields the correct order of magnitude.
Having verified Equation (9) numerically, we can use it to formulate the following rules for propagating precision to an integrated quantity from its integrand: (i) Replace the integrand by the variance of the integrand multiplied by the correlation length.
(ii) Keep the line element of the integration unchanged.
The result can be readily extended to included coordinatedependent metric factors in the integration measure.

Application to accretion discs
Applying the formalism to an an axisymmetric disc, we have for the total luminosity where k(r) is the luminosity per unit area from an annulus at radius r. Here the line element is dr, and the variance of the integrand is 4π 2 r 2 σ 2 k . Therefore For a two-dimensional non-axisymmetric disc, the luminosity is given by Assuming statistical axisymmetry, let the fluctuations of k be coherent over scales of lr(r) and l φ (r) in the radial and azimuthal directions, respectively. The luminosity per unit length of an annulus is and its variance is Since K = K1(r)dr, we obtain Equation (17) recovers the one-dimensional case (13) by taking l φ = 2πr. For quantities involving an integration in time, the method can be similarly extended. Equation (17) also reflects the intuitive expectation that the variance decreases with decreasing correlation lengths because neighboring fluctuations rapidly cancel out. For vanishing correlation lengths, our method is not valid and stochastic calculus must be invoked for a rigorous treatment. Below we explore two applications of Equation (17), to the spectra of geometrically thin, optically thick discs. We focus on finding spectral fluctuations at some given snapshot, where snapshot indicates a short time scale compared to eddy turnover times.

VARIANCE OF PREDICTED EMISSION SPECTRUM FROM A THERMAL DISC
Consider the thermal spectrum from a geometrically thin, optically thick disc. Given some surface temperature T (r, φ) that includes both a mean and a fluctuating part, the blackbody emission per unit frequency per unit area is where hP is the Planck constant, ν is the photon frequency, c is the speed of light, and kB is Boltzmann's constant. The total power per unit frequency from one side of the disc is where r * and rout are the inner and outer radii of the disc, respectively. We separate the total temperature T into a mean T and fluctuating part δT . The mean part T is assumed to be modeled by some mean-field theory which we assume to be axisymmetric and stationary so that T = T (r). The residual δT varies on turbulent time and length scales. All fluctuations are assumed to obey statistical axisymmetry although this can be generalized. In what follows, we give a unified formalism of calculating precision errors using the error propagation formula (17). A more formal and systematic description of three important different types of errors is presented in Appendix B.
For any local fluctuation in T at location (r, φ), whether from T or δT , with a variance σ 2 T (r) of its amplitude, the variance of the local fluctuation in f is As we are computing the expected variance of fluctuations around the mean, the term inside the angle brackets is evaluated at its mean value T = T . We assume that the fluctuation in f has the same coherent scales as those of fluctuations in T , as determined by the physical processes causig the latter. If the fluctuation of f is coherent over lengths lr(r) in the radial direction and l φ (r) in the azimuthal direction, the variance of the fluctuation in F will be determined by Equation (17) and given by Here, T ∂f /∂T as a function of T is determined by Planck's law (18), whereas the exact expressions of T itself, σ 2 T , lr, and l φ , will be model dependent.
We can also determine the fluctuation in the total luminosity for a frequency bandwidth over which fluctuations in the disc spectrum are coherent, in analogy to lr and l φ which are coherent scales in configuration space. Since the total luminosity is L = 2 ∞ 0 F dν, the variance of its fluctuation is where ∆ν is a frequency-dependent coherent bandwidth. One such example is given in Section 5. For convenience we define several dimensionless variables here for later use. The dimensionless disc radius isr = r/r * . The mean temperature can be written as T = T * T (r), where T * is the mean temperature at r * and includes all dependence on other disc parameters such as the central object mass, mass accretion rate, etc. The dimensionless frequency is then We denote the disc scale height by h(r), and the height-toradius ratio at the inner boundary by θ = h(r * )/r * .

APPLICATION TO PROTOPLANETARY DISCS
Gas in protoplanetary discs is subject to turbulence, as indicated by molecular line observations (Hughes et al. 2011;Guilloteau et al. 2012) and dust distributions (Pinte et al. 2016). The underlying turbulence drivers include MRI (Velikhov 1959;Chandrasekhar 1960;Balbus & Hawley 1991), self-gravitation (Toomre 1964;Shlosman & Begelman 1987), and hydrodynamic instabilities. The driving mechanism may also vary between young and old objects, and between inner and outer regions, or midplane and surface layers in a single object. Regardless of the origin, if we assume that the eddy turnover timescale of dominant eddies is comparable to the local Keplerian timescale, the eddy turnover time can reach ∼ 30 yr just at 10 AU with a central solar-mass object. Thus for most parts of a protoplanetary disc, the turbulence timescale exceeds exposure timescales of telescopes, or even timescales of multi-epoch observations. For such a "snapshot" image, some turbulent eddies are brighter and some are dimmer than the average profile predicted by a mean-field theory, and they contribute to the observed thermal spectrum at all wavelengths. It is therefore necessary to ask, whether a deviation between observations and theory at a specific wavelength has a truly systematic or a merely stochastic origin. In this section, we quantitatively incorporate the effect of turbulent fluctuations on the mean-field prediction of the disc thermal spectrum. We isolate this turbulent effect by assuming all other parameters in the problem, e.g., mass accretion rate and the α parameter, remain time-independent.
Let a fluctuation in the surface temperature due to turbulence be δT , which also generates a fluctuation in the luminosity L ed of a turbulent eddy. Although a realistic probability distribution description for turbulence is still lacking, we capture the properties of δT by assuming that the luminosity L ed for each turbulent eddy, is drawn from an ensemble with a uniform PDF, so that PL ed (x)dx is the probability to find L ed between x and x + dx. Here L is the mean of L ed , and equals the luminosity of the mean-field disc model from an area equal to that of the eddy. In general, L depends on the disc radius in a mean-field model, and so does PL ed . This simple uniform PDF allows us to proceed simply, but more comprehensive statistical prescriptions may also be used (e.g., Lee & Gammie 2021). For optically thick discs, the PDF of T can be deduced from the relation L ed = σSBT 4 , so that where σSBT 4 max = L, and Tmax = 5T /2 9/4 , with T being the mean temperature field solved from a mean-field theory. Equation (25) then defines a PDF for the disc surface temperature whose mean is T . For a given mean-field disc model, P T (x) is known at all radii and the mean and the variance of f in the ensemble can be calculated by and respectively. Note that the ensemble mean of f is different from the flux obtained by using the mean of T , namely Consequently, there is a difference between the ensemble mean of F , and the value that would be obtained using first the mean of T . The latter is what is done for most mean-field disc models. The difference between these two approaches leads to what we call the "mismatch error" (ME): where F is given by Equation (19). The ME reflects the disagreement between the fluxes that arise from the following two theoretical approaches: (i) solve for mean fields such as the mean disc temperature T , and then calculate the disc spectrum using F | T =T ; (ii) solve for both the mean and the fluctuation in temperature, the latter statistically in a mean-field model, then combine to obtain total temperature T and the associated total spectrum, and take the average to get F . Approach (i) is commonly adopted, but approach (ii) is what should be used by theorists to more accurately compare to what observations measure.
The variance of F can be computed by propagating that of f . The correlation length of the latter is assumed to be isotropic and identified with the turbulence scale l in the model. In an α disc model we have and thus l ≃ α 1/2 h (Blackman 1998). The variance of F is then using Equation (17).
We now adopt a specific profile of T to qualitatively compute ME and FE. In mean-field models the mean temperature is typically related to the disc radius via a power-law relation, T = T * r p . In the context of protoplanetary discs, p varies from −1 to ∼ −1/2 depending on whether disc heating is dominated by star irradiation or viscous heating. The exponent also likely varies with the disc radius if a transition of heating source occurs. For simplicity, we consider a constant p here. A particular model is taken from Edgar et al. (2007) for which the heat of the central plane is solely due to viscous dissipation, and the mean surface temperature T and scale height h are solved to be T = T * r −21/40 , h = h * r 21/20 .
Combining Equations (29), (31), and (32), we compute the ME and the FE as shown in Figure 2. While FE originates from turbulence and measures the corresponding stochastic fluctuation around the mean spectrum, ME is systematic as is clear from Equation (29). Consequently, fitting averaged observational data (approximately F ) with typical theoretical predictions that amount to F | T =T , introduces a bias of inferred parameters. In the present protoplanetary disc example the bias in the maximum temperature at the disk inner radius T = T (r * ) is 6.5%, that in the exponent −∂ ln T /∂ ln r is 1.0%, and that for the disc outer radius is −5.4%. The numerical values are relatively small. They do depend on the choice of the PDF for T , but the result shows that the thin disk mean field model in this context are quite precise, inasmuch as observational uncertainties are larger.

APPLICATION TO DWARF NOVAE
Dwarf novae (DNe) are characterized by their regular outbursts, and thought to result from accretion disk instability (Osaki 1974;Hōshi 1979;Lasota 2001;Hameury 2020). The standard disk instability model relies on thermal instability and the fact that the disc opacity changes rapidly and nonlinearly with temperature at ∼ 10 4 K where hydrogen ionization takes place. The temperature is directly connected to the accretion rate and so where the temperature vs. surface density equilibrium curve is unstable, so is the accretion rate. Once the accretion rate increases in the disk to a value larger than the outer supply rate can accommodate the surface density and temperature drop until matter again builds up and the cycle repeats. How exactly the disc viscosity depends on heating is model dependent. Recently, Held & Latter (2021) demonstrated that if convection results from the strong opacity increase, its combined effect with the magneto-rotational instability (MRI) may lead to a significant increase in angular momentum transport, characterized by cyclic bursts of α, the stress-to-pressure parameter. The strengthened angular momentum transport in the simulations is speculated to result from convection generated magnetic fields reseeding the MRI, an effect most prevalent in cases with long cooling times and short resistive times. In these so-called strong MRI/convection cycles, α is enhanced in the MRI phase by approximately one order of magnitude. Several studies have also reported similar α bursts in stratified shearing boxes (Simon et al. 2011;Bodo et al. 2012; Hirose et al. 2014; Coleman et al. 2018), albeit some with different origins argued. In this section, we build our model based on the results of Held & Latter (2021), and explore how such fluctuations in α can affect the disc spectrum in the quiescent phase, during which the strong MRI/convection cycle is most likely because of the high resistivity. We assume that the representative temporal and spatial fluctuations in α are local, focusing on circumstances that apply before they lead to any global coherent structure over the entire disc. The typical cycle periods of DN bursts are observed to be O(10) orbital times in simulations, and therefore may occur multiple times during one hot or cold phase of the disc. The ratio between the Keplerian timescale at the outer radius and the quiescence time of the disc can be estimated from Equations (52) for an inside-out burst with large disc radii. Here, Tc is the disc mid-plane temperature, and δ is the difference between the logarithm of the maximum surface density in the lower equilibrium branch of temperature vs. surface density during the quiescent phase; typically δ 2 (c.f. Figure 11 in Lasota 2001). As such, at all radii the cycle period of the bursty α is much shorter than the disc quiescent time.
The relatively short cycle period leads to axisymmetric fluctuations, and accordingly the azimuthal coherence length is l φ = 2πr. The radial viscous diffusion time is much larger than the orbital time because which suggests that the α fluctuation can also be considered local in radius. In the simulations exhibiting the strong MRI/convection cycles Held & Latter (2021), α is defined by a volume average over a box of radial length 4h, based on which we assume a coherent scale lr = nh with a fiducial n = 5 in the radial direction. For the standard Shakura-Sunyaev model h/r = θr 1/8 (e.g., Frank et al. 2002), and thus lr = nθrr 1/8 . We now derive the temperature and flux fluctuations. For DNe in the quiescent state, the disc is not necessarily in the global viscous equilibrium (Lasota 2001), but we assume that local equilibrium between the mean surface black-body flux and the turbulent viscous dissipation holds, i.e., Since h, Σ, and r∂rΩ change on the viscous time scale, which is much longer than 10 Keplerian orbit times, an α burst event that spans a timescale of O(10) orbits causes a mean surface temperature fluctuation T ∝ α 1/4 . If the variance of fluctuations in α is σ 2 α , we have σ 2 We take σα/α = 0.5 as estimated from Held & Latter (2021).
Since σ 2 T = σ 2 T in this case, combining Equation (37), lr = nθrr 1/8 , and l φ = 2πr, we obtain from Equation (21) that For a standard Shakura-Sunyaev model we have T ≃ T * r −3/4 . The integral can be numerically carried out, and we show in Figure 3 the relative fluctuation σF /F and also the 1σ uncertainty around F by plotting F and F ± σF together.
To quantify the corresponding variation in the total luminosity, we need to determine the photon frequency range over which fluctuations in the disc spectrum are coherent. A reasonable estimate is obtained by assuming that photons at a frequency ν are emitted by a single annulus whose position is determined by Wien's displacement law: where β = hPν/kBT * is the dimensionless frequency defined. The coherence length lr of the temperature fluctuations is propagated to a coherent frequency width ∆β coh by ∆β coh = ∂β ∂r lr r * = 2.1nθr −5/8 = 0.89nθβ 5/6 .
One can verify that ∆β coh < F/|∂F/∂β| for all β ∈ [10 −3 , 10]. Using Equation (22), the relative fluctuation in the luminosity is found to be δL/L = 1%. If we identify the variation time tL of L as that of α for the most luminous annulus, we find tL/tKep,out = 0.0016 × O(10) where tKep,out is the Keplerian timescale at the disc outer radius.

ROLE OF TELESCOPE RESOLVING POWER IN DETERMINING FALSIFIABILITY
The imprecision we have computed above is independent of the finite resolving power and binning of data by given tele-scope that observational data is subjected to before it can be compared against theoretical predictions. We now show how to incorporate this. Let the spectral resolving power of a telescope be where ∆ν is the telescope resolving bandwidth at photon frequency ν. For snapshot measurements, defined by exposure time being smaller than flux fluctuation timescales, the observed flux F T is binned using the resolved bandwidth ∆ν = ν/RP, i.e., The flux spectral fluctuations are fully resolved when the resolved bandwidth is smaller than the coherent bandwidth of the fluctuations, i.e., when RP > ν/∆ν coh where ∆ν coh = ∆β coh kBT * /hP. The fluctuation in F T is then σ T = σF from Equation (38). In the opposite limit of RP < ν/∆ν coh , every neighboring number ∆ν/∆ν coh of fluctuating bandwidths are binned. Identifying K with F T and k with F in Equation (17) we have In Figure 4 we plot the relative precision σ T /F T for different constant RP values by evaluating Equation (43). The critical dimensionless photon frequencies above which RP > β/∆β coh are 10 −3.9 , 10 −0.3 , and 10 5.7 for the RP = 5, 20, 200 cases, respectively. Smaller RP reduces the effect of fluctuations because more fluctuation cells are averaged within a single resolving bandwidth.

CONCLUSION
The scale separations between mean and fluctuating fields in astrophysical flows are finite, in contrast to the idealized assumption of infinity, which affects both the accuracy and the precision of a given mean-field theory. In particular, spatially or temporally averaged observational data unavoidably includes contributions from small-scale fluctuations, and thus when fitting data, the inferred model may fail to match those from an accurate theory, or misleadingly appear to agree with an inaccurate theory, if fluctuations around mean-field results are not properly estimated. Falsifiability of mean-field theories by comparing to observations requires careful distinction between disagreements that result from accuracy with those that result from imprecision.
While improving accuracy means increasing the fidelity of the input physics that account for the finite scale separation, in this work we focus on calculations of precision of mean-field theories, defined as the variance of the fluctuations of meanfields propagated from small-scale fields. We have considered the general case of statistically inhomogeneous small-scale fluctuations, and the derived Equation (9) is shown to be consistent with data from a numerically realized ensemble, yet significantly preceding the latter in terms of efficiency.
We then exemplified the method by computing imprecision in the prediction of accretion disc thermal spectra, induced by (i) turbulent fluctuations, and (ii) meso-scale fluctuations in α, respectively. For both, we imposed fluctuations in temperature and its coherence length allowing both to be local functions of space. The consequent error propagation of these fluctuations to global observables, namely spectra and luminosity, were then derived. Although with small magnitudes, the derived spectra fluctuations indeed depend on photon frequency, and suggests that accurate falsification of a meanfield theory may require a likelihood function or data weighting that reflects the fidelity of the theoretical mean-field values in different regions of the parameter space, as the photon frequency in the current examples.
In contrast to previous dynamical models (Lyubarskii 1997;Dexter & Agol 2011;Cowperthwaite & Reynolds 2014;Turner & Reynolds 2021), each of which considers different models for the time evolution of imposed surface density or temperature fluctuations and pursues their observational signature, our work focuses on an efficient semi-analytical method of computing the propagation of inhomogeneous fluctuations to synthetic observables to compare with snapshot observations. We have assessed how precise the predictions of standard mean-field disc theories are when subjected to such fluctuations of a given amplitude.
In addition to the stochastic contributions, we have also shown that the nonlinear relation between basic physical quantities solved in theories (i.e., surface temperature) and observables (i.e., disc spectrum) can lead to a finite systematic mismatch between the meaning of the quantity that is predicted and that is observed and then averaged. This leads to a systematic bias when backing out disc parameters of a few per cent, and a more comprehensive treatment of turbulence would offer a more accurate determination of this bias.
that the theoretical mean-field model is fully accurate and instead focus on its precision.
For a member of D th , the theory predicts a mean value of the observable, F th . For a mean-field theory, F th is the same for all members of D th , equal to its ensemble average in D th : F th th = F th ; here, the superscript of the angle brackets means the average is drawn from D th . From each element of D th , we may construct synthetic observations which represent specific predictions processed from the mean-field theory to match what a given telescope would measure. We denote this by F obs,th . Note that F obs,th represents a predicted snapshot measurement of a turbulent object and so we can also construct the ensemble average F obs,th th . The two quantities F th and F obs,th th differ in general, and indeed differed for the example model in Section 4.
Having clarified the notation, the difference between an actual observed value and theoretical mean-field value for the observable is where we have decomposed the right side into three differences: Difference (B1) is how much a real observation deviates from its multi-epoch average; difference (B2) vanishes if the two ensembles D th and D obs are identical, quantifying the accuracy of the theory; difference (B3) measures whether the quantity predicted (F th ) means the same thing as the quantity observations actually measure ( F obs,th th ). We have assumed an accurate theory, so we set F obs = F obs,th .
Furthermore, we need not distinguish in which ensemble fields are averaged, and use · obs = · th = · in what follows. Consequently, Term (B2) vanishes. Since δF fluctuates in the ensemble, it will be meaningful to quantify its mean δF and variance σ 2 δF . Using Terms (B1) to (B3), we define the "filtering error" (FE) as where σ 2 X denotes the variance of the quantity X in the ensemble, and the last step comes from our accuracy assumption F obs = F obs,th . The "mismatch error" (ME) is ∆F = F obs,th − F th = F obs,th − F th .
In addition, the error associated with the perturbations in F th induced by varying model parameters (e.g., boundary conditions, transport coefficients) contributes an "intrinsic error" (IE), σ 2 IE . We then have δF = ∆F, The IE, FE and ME can all be determined theoretically because they only involve F th and F obs,th , but not F obs . Correspondingly, F th + δF ±σ δF will be a mean-field prediction with error bars, giving a finite range of where we expect the observed value F obs to locate. Thus to facilitate a more appropriate comparison between theory and observations than what is commonly done, we must quantify the precision of the former so that the error can be added to the mean and the result compared against an observation when that obsevation is a member of an ensemble rather than an ensemble average.
The non-vanishing ME in thermal disc spectra can be elucidated in the following way. Consider a member from D th . Let its surface temperature, as if we could measure it, be T obs,th , which is turbulent and varies with the disc radius r and azimuthal position φ. The observed total power emitted per unit frequency from one side of the disc is F obs,th = 2π 0 dφ rout r * 2hPν 3 /c 2 e h P ν/k B T obs,th − 1 rdr where I is the spectrum functional. The ensemble mean of the observation is F obs,th . On the other hand, in a mean-field disc theory, the constructed equations are solved to obtain a mean surface temperature T th . The predicted total emission is Even assuming an accurate theory T th = T obs,th , we still have F obs,th = I T obs,th = I T obs,th = F th (B11) because of the nonlinear dependence on the surface temperature. The difference between the left and right sides is defined as the ME.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.