\section*{Acknowledgments}
\author[G.~Kaehler]{Goetz~Kaehler\affil{1}\corrauth and Alexander~Wagner\affil{1}}
\address{\affilnum{1}Department of Physics, North Dakota State University, Fargo, ND 58108, U.S.A.}
\emails{{goetz.kaehler@ndsu.edu} (G.~Kaehler)}

\begin{abstract}
We analyze the Lattice Boltzmann method for the simulation of fluctuating hydrodynamics by Adhikari {\em et al.} [Europhys. Lett. {\bf71}, 473 (2005)] and find that it shows excellent agreement with theory even for small wavelengths as long as a stationary system is considered. This is in contrast to other finite difference and older lattice Boltzmann implementations that show convergence only in the limit of large wavelengths. In particular cross correlators vanish to less than $0.5\%$. For larger mean velocities, however, Galilean invariance violations manifest themselves through errors of a magnitude similar to those of the earlier implementations.
\end{abstract}

\maketitle

\section{Introduction}
Fluctuations are important for many hydrodynamic phenomena, from colloid diffusion to phase-separation close to the critical point. Particle based methods such as Stochastic Rotation Dynamics~\cite{SRD}, Lattice Gas~\cite{LG} or Molecular Dynamics simulations~\cite{MD} naturally give rise to stochastic noise.
In contrast the lattice Boltzmann (LB), or finite difference discretization of the Navier Stokes equations require fluctuations that have to be included manually.
The guiding principle for doing this is the theory of the fluctuating Navier Stokes equations \cite{ll9}. Despite the success of applying the Navier Stokes equations to very small-scale flows formally the hydrodynamic limit requires large wavelengths. For fluctuating hydrodynamics the constraint of large wavelengths becomes important and standard discretization will give results that are not in agreement with statistical physics for shorter wavelengths. For a detailed analysis of simulating fluctuating hydrodynamics using finite difference methods and some remedies to improve this situation see the recent manuscript of A.~Donev~\cite{donev1}. Similar deficiencies are found for implementations of fluctuating Navier Stokes equations using the Lattice Boltzmann approach introduced by Ladd~\cite{ladd1}. It is, however possible to use a more fundamental approach to include fluctuations in the LB method. Adhikari {\em et al.}\cite{adhikari1} introduced noise on all nonconserved modes, not only the hydrodynamic ones, leading to a scheme which shows good agreement with theory even for large wavelengths. Duenweg {\em et al.} rederived this noise implementation from detailed balance considerations of lattice gases \cite{duenweg}. Both approaches are numerically identical. In this paper we study the degree of improvement achieved and show that many of the deficiencies that plague finite difference discretizations of fluctuating Navier Stokes equations are absent in this Lattice Boltzmann implementation as long as we consider a system with vanishing mean velocity. For large mean velocities Galilean invariance is violated and errors of a similar magnitude to the earlier implementations are observed.

\section{Fluctuating Lattice Boltzmann with Ghost Noise}
Following the derivation of Adhikari {\em et al.} \cite{adhikari1} we start with the Lattice Boltzmann equation (LBE)
\begin{equation}
\label{eqn:LBE1}
f_i(\mathbf{x}+\mathbf{v}_i, t+1) = f_i(\mathbf{x}, t) + \sum_j \Lambda_{ij}\left\lbrack f_j(\mathbf{x}, t) - f_j^0(\mathbf{x},t)\right\rbrack + \xi_i(\mathbf{x}, t).
\end{equation}
Here the the $f_i$ are the particle densities at position $x$, time $t$ associated with with velocity $\mathbf{v}_i$. $\Lambda_{ij}$ is the collision matrix and $\xi_i$ are the noise terms. We use the standard local equilibrium distribution given by
\begin{equation}
\label{eqn:fnull}
f_i^0 = \rho w_i \left\lbrack 1 + \frac{3}{c^2} \mathbf{u}.\mathbf{v}_i + \frac{9}{2c^4}\left(\mathbf{u}.\mathbf{v}_i\right)^2 - \frac{3}{2c^2}\mathbf{u}.\mathbf{u}\right\rbrack,
\end{equation}
which is the discretized version of a Maxwell distribution \cite{he1997, qian1992}. In equilibrium the $f_i$ will fluctuate around this distribution. The noise terms $\xi_i$ must be chosen such that, in the case of isothermal Lattice Boltzmann (LB), the density $\rho=\sum_i f_i$ and momentum $\rho\mathbf{u} = \sum_i f_i \mathbf{v}_i$ are conserved, i.e. $\sum_i \xi_i = 0$ and $\sum_i \xi_i v_i = 0$. Furthermore a proper fluctuation dissipation theorem (FDT) corresponding to the collision operator $\Lambda_{ij}$ is obeyed. This implies that the $\xi_i$ are correlated. We can find a representation in which the noise terms are uncorrelated by transforming the LBE into moment space. The moments are given by
\begin{equation}
\label{eqn:transform1}
M^a(\mathbf{x}, t) = \sum_i m_i^a f_i(\mathbf{x}, t).
\end{equation}
So far this is a standard Multi-Relaxation-Time (MRT) representation \cite{dhumieres, lallemand1, adhikari2}. The back transform is given by $f_i(\mathbf{x}, t) = \sum_a n_i^a M^a(\mathbf{x}, t)$. However, in order to construct a proper FDT these transforms cannot be orthogonal as in other MRT methods~\cite{dhumieres, lallemand1}, so here we have $n_i^a \ne m_i^a$. Instead the transforms are chosen such that
\begin{equation}
\label{eqn:hermitenorm}
\sum_i w_i m_i^a m_i^b = \sum_i m_i^a n_i^b = \delta^{ab}
\end{equation}
with $n_i^a = w_i m_i^a$ while maintaining a diagonal moment space representation of the collision operator $\Lambda_{ij} =- \sum_a\sum_b n_i^a \frac{1}{\tau^a} \delta^{ab} m_j^b$.
Now the moment transformation matrices are orthogonal with respect to the Hermite norm. Such transforms with weighted norms were proposed before \cite{vergassola1990, dellar2002, dellar2003} in different contexts. The necessity of the Hermite norm is briefly outlined after Eq.~(\ref{eqn:decouplemoments}) below and allows for a convenient definition of the moment space noise terms $\xi^a$ as independent random variables. We can now rewrite the collision term of the Lattice Boltzmann equation in terms of the moments $M^a$ as
\begin{equation}
\label{eqn:collision}
f_i(\mathbf{x}+\mathbf{v}_i, t+1) = \sum_a n_i^a \left\lbrace M^a(\mathbf{x}, t) - \frac{1}{\tau^a} \left\lbrack M^a(\mathbf{x}, t) - M^{a, 0}(\mathbf{x}, t) \right\rbrack + \xi^a \right\rbrace.
\end{equation}
Adhikari {\em et al.}\cite{adhikari1} then obtain the FDT by performing a Fourier transform of the fluctuations from the mean of the moments $\delta M^a = M^a - \langle M^a \rangle$. They then use the $k$-independence of these for an ideal gas to obtain
\begin{equation}
\label{eqn:FDT}
 \left\langle \xi^a \xi^c \right\rangle = \frac{\tau^a + \tau^c - 1}{\tau^a \tau^c} \left\langle \delta M^a \delta M^c \right\rangle.
\end{equation}
One particular result of the derivation is that the moment fluctuations $\xi^a$ decouple because
\begin{eqnarray}
\left\langle\delta M^a \delta M^b\right\rangle & = & \sum_i \sum_j m_i^a m_j^b \langle \delta f_i \delta f_j \rangle \nonumber \\
& = & \sum_i \sum_j m_i^a m_j^b \bar{f}_i \delta_{ij}\nonumber \\
& = & \sum_i m_i^a m_i^b \bar{\rho} w_i \nonumber \\
\label{eqn:decouplemoments}
& = & \bar{\rho} \delta^{ab} .
\end{eqnarray}
Here we used $\langle \delta f_i \delta f_j \rangle = \bar{f}_i \delta_{ij}$ with $\delta f_i = f_i - \bar{f}_i$ where $\bar{f}_i$ is the spatially uniform global equilibrium distribution function \cite{ll10}. Adhikari also assumed that $\mathbf{u} \ll 1$ so that $\bar{f}_i = \bar{\rho} w_i$. This allows us to use the orthogonality relation of Eq.~(\ref{eqn:hermitenorm}) in the last step of the calculation above. For a different transformation we would obtain non-diagonal elements in the fluctuation matrix which will then require correlated noise terms which are more cumbersome to implement. For practical applications it is important to note that the $\mathbf{u} \ll 1$ condition for the noise introduces a non-Galilean invariant contribution. We comment on this in our validation section.
Inserting Eq.~(\ref{eqn:decouplemoments}) into Eq.~(\ref{eqn:FDT}) leads to a noise expression of
\begin{equation}
\label{eqn:noiseamp}
\xi^a = \frac{1}{\tau^a} \sqrt{\bar{\rho} \left(2 \tau^a - 1\right)}N,
\end{equation}
where $N$ is a random variable with zero mean and a variance of one.

Note that the moments $M^a$ are chosen to include the hydrodynamic moments. In the isothermal case discussed here they are comprised of the conserved quantities $\rho$ and $\mathbf{j}$, and the stress modes $\mathbf{\Pi}$.  The remaining degrees of freedom are often called ghost modes as they do not appear in the isothermal Navier Stokes equations. However, the key result of the Adhikari {\em et al.} paper \cite{adhikari1} was that they need to be taken into account when including noise. Thus we add noise on all non conserved quantities, i.e. stress and ghost modes, in Eq.~(\ref{eqn:collision}) according to Eq.~(\ref{eqn:noiseamp}).

In practice we implement this algorithm by calculating the moments by means of Eq.~(\ref{eqn:transform1}), performing the collision on the moments, adding the noise term and then transforming back into $f$-space as indicated in Eq.~(\ref{eqn:collision}). The streaming step is then done in $f$-space. This algorithm is almost as efficient as the standard LB implementation. The additional computational cost for calculating the ghost modes and the random numbers results in a computational overhead of less than $20\%$.

\section{Correlators in a D2Q9 Implementation}

\begin{figure}
\subfigure[Without ghost noise]{\includegraphics[scale = .4]{d2q9ngrr.eps}}
\subfigure[With ghost noise]{\includegraphics[scale = .4]{d2q9grr.eps}}
\caption{$S_{\mathbf{k}}(\rho)$ averaged over $2\times10^8$ iterations in a  $\tau^a=1$ for all $a$, $V=20^2$, fluctuating D2Q9 ideal gas without and with active ghost noise. Note that different scales are used to visualize the slight deviations seen in the ghost noise case.}
\label{fig:d2q9grr}
\end{figure}

\begin{figure}
\subfigure[Without ghost noise]{\includegraphics[scale = .4]{d2q9nguxux.eps}}
\subfigure[With ghost noise]{\includegraphics[scale = .4]{d2q9guxux.eps}}
\caption{$S_{\mathbf{k}}(u_x)$ averaged over $2\times10^7$ iterations in a  $\tau^a=1$ for all $a$, $V=20^2$, fluctuating D2Q9 ideal gas without and with active ghost noise. Note that different scales are used to visualize the slight deviations seen in the ghost noise case.}
\label{fig:d2q9guxux}
\end{figure}

To evaluate this method we present results for the D2Q9 (two dimensions, 9 base velocity vectors) LB model. The results are similar for other models, in particular we also tested D1Q3 and D3Q15. As D2Q9 base velocity set we use $\{v_i\} = \{(0,0), (1,0), (0,1),$ $(-1,0), (0,-1),$ $(1,1), (-1,1),$ $(-1,-1), (1,-1)\}$ and the $\{w_i\} = \{4/9,$ $1/9,$ $1/9,$ $1/9,$ $1/36,$ $1/36,$ $1/36,$ $1/36\}$
.
The matrix elements $m_i^a$ in transform~(\ref{eqn:transform1}) are then given by
\begin{equation}
\label{eqn:d2q93a}
\{m_i^a\} = \left(
\begin{array}{ccccccccc}
 1 &  1 &  1 &  1 &  1 &  1 &  1 &  1 &  1 \\
 0 &  \sqrt{3} &  0 &  -\sqrt{3} &  0 &  \sqrt{3} &  -\sqrt{3} &  -\sqrt{3} &  \sqrt{3} \\
 0 &  0 &  \sqrt{3} &  0 &  -\sqrt{3} &  \sqrt{3} &  \sqrt{3} &  -\sqrt{3} &  -\sqrt{3} \\
 0 &  \frac{3}{2} &  \frac{-3}{2} &  \frac{3}{2} &  \frac{-3}{2} &  0 &  0 &  0 &  0 \\
 0 &  0 &  0 &  0 &  0 &  3 &  -3 &  3 &  -3 \\
 -1 &  \frac{1}{2} &  \frac{1}{2} &  \frac{1}{2} &  \frac{1}{2} &  2 &  2 &  2 &  2 \\
 0 &  -\sqrt{\frac{3}{2}} &  0 &  \sqrt{\frac{3}{2}} &  0 &  \sqrt{6} &  -\sqrt{6} &  -\sqrt{6} &  \sqrt{6} \\
 0 &  0 &  -\sqrt{\frac{3}{2}} &  0 &  \sqrt{\frac{3}{2}} &  \sqrt{6} &  \sqrt{6} &  -\sqrt{6} &  -\sqrt{6} \\
 \frac{1}{2} &  -1 &  -1 &  -1 &  -1 &  2 &  2 &  2 &  2
\end{array}
\right).
\end{equation}
The corresponding elements $n_i^a$ of the back transform are defined by the requirement $n_i^a = m_i^a w_i$.
The zeroth moment then is the density $\rho$, the first and second are (up to a factor) the components of the momentum, the third and fourth the components of the shear stress and the fifth resembles the bulk stress \cite{kaehler2010}. The remaining three moments are the ghost modes. Thus the equilibrium moments $M^{a,0} = \sum_i m_i^a f_i^0$ are
\begin{equation}
\label{eqn:d2q94}
\begin{array}{ccc}
M^{0,0} & = & \rho\\
M^{1,0} & = &  \sqrt{3}\rho u_x \\
M^{2,0} & = &  \sqrt{3}\rho u_y \\
M^{3,0} & = &  \frac{3}{2} \rho (u_x^2 - u_y^2) \\
M^{4,0} & = &  3 \rho u_x u_y \\
M^{5,0} & = &  \frac{3}{2} \rho (u_x^2 + u_y^2) \\
M^{6,0} & = & M^{7,0} = M^{8,0} =  0 .
\end{array}
\end{equation}

We present here results for $k$-independence of the moment fluctuations predicted by Eq.~(\ref{eqn:decouplemoments}). In particular we consider the normalized static structure factor
\begin{equation}
\label{eqn:scorrelator}
S_{\mathbf{k}}(M^a) = N^a \left\langle \delta M^a(\mathbf{k}) \delta M^a(-\mathbf{k}) \right\rangle
\end{equation}
where $\delta M^a(\mathbf{k}) = \sum_{\mathbf{x}} \delta M^a(\mathbf{x}) e^{-i \mathbf{k} \cdot \mathbf{x}}$ is the discrete spatial Fourier transform of $\delta M^a$ and $\sum_{\mathbf{x}}$ is understood to be the summation over all discrete lattice sites. The normalization constant $N^a$ such that $S_{\mathbf{k}}(M^a) = 1$ is equivalent to $\bar{\rho}$. I. e. for the density $N^\rho = \frac{1}{\bar{\rho}^{3} V}$ and velocity components $N^{u_\alpha} = \frac{1}{\bar{\rho} V k_b T }$ where $k_b T = \frac{1}{3}$ for the isothermal D2Q9 model employed. A value of $1$ throughout $k$-space for the structure factor of any of the moments given in Eq.~({\ref{eqn:d2q94}}) thus indicates agreement with Eq.~(\ref{eqn:decouplemoments}). The volume $V$ is just the number of lattice points $V = \sum_{\mathbf{x}} 1$ and the division by it is just a normalization artifact of the Fourier transform.

According to the argument put forth in \cite{adhikari1} we expect the mean square fluctuations of all moments $M^a$ to be unity throughout $k$-space. For the density $\rho$ this is confirmed to three orders of magnitude in Fig.~\ref{fig:d2q9grr}(b) for $S_{\mathbf{k}}(\rho)$ and in Fig.~\ref{fig:d2q9guxux}(b) for $S_{\mathbf{k}}(u_x)$. We find similar agreement for all nine moments of the D2Q9 model. For comparison we set the noise on the non-hydrodynamic modes ($M^6, M^7, M^8$) to zero, recovering the original Ladd method~\cite{ladd1} and, as seen in Figures~\ref{fig:d2q9grr}(a) and~\ref{fig:d2q9guxux}(a), the lack of noise on the ghost terms leads to drastic deficiencies for large $k_x, k_y$ values. Note that there are no deficiencies in Fig.~\ref{fig:d2q9grr} for $k_x = 0$ and $k_y=0$. The reason is that the projection of the D2Q9 model onto one coordinate axis yields a D1Q3 model. The isothermal ideal gas D1Q3 model, however, only has one stress mode and no ghost modes and thus there is no difference between the Ladd and Adhikari implementations in these projections. This is again observed in Fig.~\ref{fig:d2q9guxux}(a) where $S_{\mathbf{k}}(u_x)$ exhibits white noise along the $k_x$ axis even in the absence of ghost noise.
\begin{figure}
\subfigure[Without ghost noise]{\includegraphics[scale = .40]{d2q9nguxuy.eps}}
\subfigure[With ghost noise]{\includegraphics[scale = .40]{d2q9guxuy.eps}}
\caption{$R_{\mathbf{k}}(u_x u_y)$ averaged over $8\times10^6$ iterations in a  $\tau^a=1$ for all $a$, $V=20^2$ fluctuating D2Q9 ideal gas simulation with and without active ghost noise. Again, take note of the different scales.}
\label{fig:d2q9gvxvy}
\end{figure}
Motivated by private communication with A.~Donev who is developing a general finite volume scheme to solve the fluctuating Navier Stokes Equations~\cite{donev1} based on a third order Runge-Kutta integrator we also measured the cross correlator
\begin{equation}
\label{eqn:rcorrelator}
R_{\mathbf{k}}(u_x, u_y) = N^{u_x} \left\langle u_x(\mathbf{k}) u_y^*(\mathbf{k}) \right\rangle.
\end{equation}
According to Eq.~(\ref{eqn:decouplemoments}) this quantity is expected to vanish. This is again confirmed nicely in Fig.~\ref{fig:d2q9gvxvy}(b) to three orders of magnitude. In contrast measurements of $R_{\mathbf{k}}(u_x, u_y)$ in an implementation without ghost noise exhibits significant correlations of up to $0.25 \bar{\rho}$ for intermediate $k_x$ and $k_y$ ranges as seen in Fig.~\ref{fig:d2q9gvxvy}(a).

The required condition in Eq.~(\ref{eqn:decouplemoments}), $\mathbf{u} \ll 1$, suggests that this noise implementation may suffer from a lack of Galilean invariance. To estimate the magnitude of this violation we consider an imposed mean velocity in the $x$-direction. We measured correlators for a fluctuating system with large superimposed velocity of $u_x = 0.1$. The results in Fig.~\ref{fig:d2q9uxlarge} indicate that indeed the moment fluctuations do not decouple and Eq.~(\ref{eqn:decouplemoments}) is no longer fulfilled. Compared to the Ladd implementation these errors are still smaller, but they approach the same order of magnitude for maximal accessible velocities. A more comprehensive investigation of these effects is subject of a forthcoming publication.

\begin{center}
\begin{figure}
\subfigure[$S_{\mathbf{k}}(\rho)$]{\includegraphics[scale = .258]{d2q9gurr.eps}}
\subfigure[$S_{\mathbf{k}}(u_x)$]{\includegraphics[scale = .258]{d2q9guuxux.eps}}
\subfigure[$S_{\mathbf{k}}(u_x, u_y)$]{\includegraphics[scale = .258]{d2q9guuxuy.eps}}
\caption{Correlators $S_{\mathbf{k}}(\rho)$, $S_{\mathbf{k}}(u_x)$, and $R_{\mathbf{k}}(u_x, u_y)$ averaged over $5\times10^6$ iterations or a $\tau^a=1$ for all $a$, $V=20^2$ fluctuating D2Q9 ideal gas simulation with a constant velocity of $u_x=0.1$.}
\label{fig:d2q9uxlarge}
\end{figure}
\end{center}
\section{Discussion and Outlook}
We have shown here that the Adhikari approach to use an improved LB method presents a promising scheme to simulate fluctuating hydrodynamics. The ability to interprete the ghost degrees of freedom as resulting from discrete particle distributions gives us the ability to systematically introduce fluctuations. This approach recovers fluctuations not only in the hydrodynamic limit but also for much shorter wavelengths. However, this is only true in the absence of flow. Since lattice Boltzmann methods are not generally used in this regime one may wonder if Galilean invariance violations may not erase much of the improvement achieved by including noise in the ghost modes. This is a subject to which we will return in a forthcoming paper.


\section{Introduction}
\label{sec:intro}

Accretion discs are ubiquitous in astrophysics, often forming around stars and compact objects where
angular momentum would  inhibit direct gravitational infall of the plasma onto the central object.
A long standing theme of research in accretion theory has been to understand how how the discs transport angular momentum and remain quasi-steady accretors. This transport  likely involves different combination of  local and global mechanisms depending on the circumstances. However, there is presently an intellectual gap between the study of the fundamental theory of angular momentum transport in accretion
discs via numerical simulations and practical accretion disc models
that can be used to compare with observed spectra.
For reviews that sample a range of perspectives
see \cite{2002apa..book.....F,2003LNP...614..329B,2009apsf.book.....H,2010arXiv1005.5279S,2010AN....331..101B}.

The most commonly used practical accretion disc model  is based on local viscous transport of angular momentum in axisymmetric discs and hides the unknown physics of transport into an effective local ``turbulent viscosity'' \citep{1973A&A....24..337S,1981ARA&A..19..137P,2002apa..book.....F}.
 Because microphysical  molecular viscosity is typically much too small to drive the observed accretion rates given constraints on surface densities, some kind of enhanced transport is needed.
While these mechanisms routinely involve some form of turbulent motion (often with a magnetic origin), the standard practical model equations describing accretion discs do not specify or accommodate
the nuanced physics of transport, lumping the processes into a single dimensionless transport coefficient.
The assumption that the discs are axisymmetric can then at best apply in a spatially or temporally averaged mean-field sense since turbulence necessarily violates axisymmetry locally and turbulent transport is a stochastic process.

In order to properly compare standard axisymmetric accretion theory with observations, the limitations of the  theory must be quantified. The theory has a limited precision which implies that agreement between models and observation can only be approximate. Disagreement within the precision error may be consistent with a stochastic nature of the discs.  A key quantity, and the focus of this paper, is the relative precision error (RPE)\footnote{Compare the statements `babies sleep between one and twenty three hours per day' and `babies sleep \unit{\ensuremath{8\pm1}}{\hour} per day.' The first statement has a large RPE, but is almost certain to be accurate. The second statement has a low RPE, but may well be incorrect. However, its low RPE makes the second statement testable} in the luminosity at a given frequency.
\citet{1998MNRAS.299L..48B} discussed the RPE for accretion discs in the context of Advection Dominated Accretion Flows (ADAFs) and thin discs in active galactic nuclei (AGN) and  X-ray binarys (XRBs) in which the observation times are often longer than the orbit times at a given radius.
In contrast,  protoplanetary discs often involve small observational exposure times compared
to disc orbit times.   In this context, we revisit quantifying the precision of standard accretion disc theory and make predictions for multi-epoch observations. The predictions can help to diagnose whether
 a systematic variability  may also be present.

 In section 2, we compute the azimuthal and radial contributions to the luminosity RPE
  and combine them for a spatially
 unresolved source. In section 3 we discuss how to minimize the RPE and the role of telescope
 spectral resolving power. In section 4 we consider how to apply the results
 to the specific case of LLRL 31, one of the few
 protoplanetary discs in which there are multi-epoch observations.
 Although we find that the observed variability in that source is systematic
 and not stochastic, the application exemplifies how to use the results herein.
 We also discuss further  implications for  future multi-epoch observations.
of protoplanetary discs, and conclude in section 5.

\section{Calculation of  RPE  in Luminosity from Turbulent Dissipation }

We assume a standard $\alpha$-disc prescription \citep{1973A&A....24..337S} for the turbulent viscosity.The characteristic turbulent eddy scales have characteristic velocities \ved{} and sizes \led{} so that the effective viscosity satisfies
\begin{equation}
\nu \approx \alpha \cs{} h \approx \ved{} \led{} \sim {v_{ed}^2\over \Omega}
\label{1}
\end{equation}
where \cs{} is the disc sound speed, $h$ is the scale height,
and $\Omega$ is the orbital speed.
We take $0<\alpha<1$ to be a constant over the whole disc.
The last similarity in (\ref{1}) follows from the assumption that the
eddy turnover time $\tau_{ed}\sim {1\over \Omega}$  consistent with simulations of the magnetorotational instability
(MRI) (e.g. \cite{2003LNP...614..329B}). For a thin disc in hydrostatic equilibrium, Eq. (\ref{1}) also implies that
$v_{ed} \sim \alpha^{1/2} c_s$ and $l_{ed} \sim \alpha^{1/2} h$ \citep{1998MNRAS.299L..48B}.

In what follows, we focus on cases for which the
required observational exposure time is much less than the eddy
turnover time at all radii. We can then regard a measured spectrum as
 an instantaneous snapshot.
This applies  to \emph{Spitzer} observations of protoplanetary discs.
The discs range in size from a few tenths of an AU{} to several tens of AU
so that \ted{} ranges from days (in the inner disc) to years (in the outer disc).
In contrast, a typical \emph{Spitzer} exposure takes only a few minutes.

Since the sources of interest are spatially unresolved,
we calculate the contributions to the RPE by
averaging in the azimuthal and in the radial directions for the axisymmetric theory.
If  the disc is sufficiently optically thick then  the photosphere
occurs within the first layer of energy containing eddies.
An optically thin disc requires additionally averaging over the third dimension of eddies.

\subsection{Azimuthal Contribution}

\label{sec:simple}

Consider an annulus of width equal to one eddy scale centered at a radius $r$ within the disc.
The number of eddies producing observable emission in this  annulus is given by
\begin{equation}
\Ned{} \sim  \frac{2 \pi r}{\led{}}  \left[ 1+\left({h\over l_{ed} }-1\right)e^{-\tau_{\nu} }\right],
\label{2}
\end{equation}
where $\tau_\nu$ is the optical depth. In the limit that $\tau_\nu \rightarrow 0$, eddies over the entire thickness of the disc contribute to the emission. For $\tau_\nu \rightarrow \infty$ only the top layer of eddies contribute
to the emission.

Over its lifetime, the brightness of each dominant  eddy will rise to some maximum and then
fall as it cascades to small scales. Let us assume that the rate of brightening and dimming are constant so
that the eddy brightens to its peak luminosity $L_0$ at  $t_{ed}/2$ and then falls to zero
at $t=t_0$.
The probably distribution function for the luminosity during the rise and fall is then
\begin{equation}
p(L_e) d{L_e} = \frac{1}{L_0}dL_e.
\end{equation}
Then the mean luminosity for each eddy is $\langle L_e \rangle =L_0/2$
and the
mean squared luminosity is
\begin{equation}
 {\langle L_e^2\rangle} = {1\over 2}\left[
  {\int_0^{L_0}  {L_e^2\over L_0} dL_e
\over \int_0^{L_0} {1\over L_0} dL_e } +
 {\int_{L_0}^0  {L_e^2\over L_0} dL_e
\over \int_{L_0}^0 {1\over L_0} dL_e }\right]
= L_0^2/3.\end{equation}
The variance is then
\begin{equation}
\sigma^2(L_e) = \mean{L_e^2} - \mean{L_e}^2 = \frac{L_0^2}{12}
\end{equation}

If the value of the mean luminosity of each eddy is normally distributed about
$L_0/2$ then the standard deviation of the luminosity per eddy when averaged over
all $N_{ed}$ eddies in a given annulus satisfies $ \frac{\sigma(L)}{\sqrt{\Ned{}}}$.
Using $L_\nu$ to indicate the luminosity at frequency $\nu$ measured from the
annulus whose  peak temperature corresponds to that  frequency, the fractional luminosity variation of the entire annulus of $N_{ed}$ eddies is then
\begin{equation}
\left[\frac{\Delta L_\nu}{{L_\nu}}\right]_{az}
=\frac{N_{ed}\Delta L_e}{N_{ed}\mean{L_e}}
= {\sigma (L_e) / {\sqrt N_{ed}} \over   L_0/2}
= \frac{1}{\sqrt{3 \Ned{}}}
\label{5}
\end{equation}
where we have used Eq. (\ref{2}),  the relations below Eq. (\ref{1}), and the subscript
$az$ to indicate the azimuthal contribution to the variation.
In the optically thick limit Eq. (\ref{5}) gives
$\left[\frac{\Delta L_\nu}{{L_\nu}}\right]_{az}=\frac{\alpha^{\frac{1}{4}}}{\sqrt{6 \pi}} \parenfrac{h}{r}^{\frac{1}{2}}$.

\subsection{Radial Contribution}
\label{sec:radial}

In determining the radial contribution to the luminosity RPE we
consider two independent  contributions. The first comes from
a radial smoothing implicit to the mean field theory. The second comes
from turbulent velocity fluctuations that add macroscopic fluctuations about the mean orbital
velocity at each radius. We address each in turn.

In formulating a mean field theory for accretion discs,
one can in  principle simply take vertical and azimuthal averages to obtain a theory that depends only
on radius without any radial smoothing.  However in the  commonly employed
 framework which invokes a turbulent viscosity
 (Shakura \& Sunyaev 1973),  the scale of the turbulence is assumed to be
smaller than the scale of the mean quantities.  The theory is not really resolved on spatial
scales below an eddy scale so we therefore include a
radial smoothing scale, $\xi$, chosen such that $\led{} < \xi < r_{\text{disc}}$.
In addition, telescopes have a finite spectral resolution and a given frequency  in practice corresponds to a radius range in the disc, not a precise radius for a spatially unresolved source.
We therefore invoke a radial smoothing  of mean quantities of the form
\begin{equation}
{\overline X}(r)
       \approx \mean{X(r,t)}_{\xi}
       = \int_{-\xi}^{\xi} X(r+\lambda,t) d \lambda,
\label{6}
\end{equation}
where $X$ is an arbitrary quantity to be averaged and $\overline X$ is its mean.
The similarity follows  from the assumption that eddy statistics are steady in time.
Blackman (1998) discussed an additional temporal average but
as mentioned above, our present interest is for cases in which the required snapshot observation time is much less than an eddy turnover time at all radii.

The radial contribution to the luminosity variation can be written as
\begin{equation}
\left[\frac{\Delta L_{\nu}}{L_{\nu}}\right]_{rad}
\approx
\left | \frac{r}{L_{\nu}} \pderiv{L_{\nu}}{r} \right|
\frac{\Delta r}{r}.
\label{eq:VariabilityEquation}
\end{equation}
We must now calculate the contributions to the right hand side.

We write the total imprecision in radius  as a quadratic sum of
 the smoothing and velocity fluctuation contributions. That is,
\begin{equation}
\frac{(\Delta r)^{2}}{r^{2}} \approx
\frac{\resolution{\Delta r}^{2}}{r^{2}}+
\frac{\fluctuate{\Delta r}^{2}}{r^{2}}
\label{eq:radiusRPEdefine}
\end{equation}
where first term comes from the blurring induced by the radial smoothing scale $\xi$  and the second term comes from the macroscopic turbulent velocity fluctuations.

For the radial smoothing term of  Eq.  (\ref{eq:radiusRPEdefine}),
we use the fact that
\begin{equation}
r - \frac{\xi}{2 \Nres{\frac{1}{2}}} < r <
r +\frac{\xi}{2 \Nres{\frac{1}{2}}}
\end{equation}
are indistinguishable once $\xi$ is chosen, so that
\begin{equation}
{(\Delta r)_{rs}\over r}= \frac{\xi}{2r \Nres{\frac{1}{2}}},
\label{12}
\end{equation}
where
\begin{equation}
\Nres{} = 1 +\left({h\over l_{ed} }-1\right)e^{-\tau_{\nu}}
+ \frac{t_{\text{obs}}}{\ted{}}
\label{12a}
\end{equation}
and measures the effective number of averaging scales per $\xi$. The role of $\tau_\nu$  the same as that for Eq.
(\ref{2}). In the limit $t_{obs}<<t_{ed}$ and $\tau_\nu >> 1$, $N_{fl} =1$.
For the velocity fluctuation contribution to Eq. (\ref{eq:radiusRPEdefine}),
we note that  the dominant contribution to the local velocity at a given radius is  $\propto r^{-1/2}$ from
Keplerian rotation. Thus
\begin{equation}
\frac{\fluctuate{\Delta r}}{r} \approx 2 \frac{\fluctuate{\Delta v_0}}{v_0} \approx 2 \frac{\ved{}}{v_0 \Nfl{\frac{1}{2}}},
\label{9}
\end{equation}
where we have again used the standard result for error on the mean as applied to $N_{fl}$ eddies each of whose mean velocity is normally distributed with fluctuations of order $v_{ed}$.
  The number of eddies over which the radial averaging is performed  is given by
\begin{equation}
\Nfl{} = \frac{\xi}{\led{}}\left[1+\left({h\over l_{ed} }-1\right)e^{-\tau_{\nu}} \right]
+ \frac{t_{\text{obs}}}{\ted{}}
\label{eq:NflDefine}
\end{equation}
For $\tau_\nu >> 1$ and

  $t_{\text{obs}}/t_{ed}\sim 0$, this gives $N_{fl}= \xi/l_{ed}$.
In this limit, adding (\ref{12}) and (\ref{9}) in quadrature, we obtain
\begin{equation}
\frac{(\Delta r)^{2}}{r^{2}} \approx
4 \frac{v_{\text{ed}}^{2} l_{\text{ed}}}{v_{0}^{2} \xi}
+ \frac{\xi^{2}}{4 r^{2}}.
\label{eq:deltarexpand}
\end{equation}

If each ring of material radiates around a specific blackbody temperature $T$,  the luminosity of an annulus of
width  $l_{ed}$ is
\begin{equation}
L_{\nu}  = 2 \pi r \ensuremath{B_{\nu}^{\text{peak}}}{}\paren{T(r)} l_{ed},
\label{eq:LnuDefine}
\end{equation}
where
\begin{equation}
B_{\nu}^{\text{peak}}(T)=A T^3
\label{16a}
\end{equation}
 is the peak value of the Planck function $B_\nu$ and
$A$ is a constant that comes from evaluating
\begin{equation}
B_{\nu}
(T) = T^3 \frac{2 k^3}{h^2 c^2} \frac{x^3}{e^x -1}
\end{equation}
at  $x \equiv h \nu / k T
\sim 2.82$,  corresponding to the Wien displacement law.
Then
\begin{equation}
\pderiv{L_{\nu}}{r}
=  2 \pi A T^2 \bracket{l_{ed}T + 3 r l_{ed}\pderiv{T}{r}+ r T {dl_{ed}\over dr}}.
\label{eq:pderivLnu}
\end{equation}
If we assume that
 \begin{equation}
T = C r^{-p},
\label{eq:TradialPowerLaw}
\end{equation}
then  combining with the scalings below Eq. (1) gives  $l_{ed}=\alpha^{1/2}h\propto c_s/\Omega \propto r^{{1\over 2}(3-p)}$.
In combination  with Eqs~(\ref{eq:VariabilityEquation}), (\ref{eq:LnuDefine}) and (\ref{eq:pderivLnu}), we then obtain
\begin{equation}
\left[\frac{\Delta L_{\nu}}{L_{\nu}}\right]_{rad}
\approx
\left| 1 - 3p +{3\over 2}-{p\over 2} \right|
\frac{\Delta r}{r}=\left| {5\over 2}-{7p\over 2} \right|
\frac{\Delta r}{r}
\label{eq:DiscVariability}
\end{equation}
\citet{2007MNRAS.381.1280E} calculated $p$ for a variety of disc models, and  their values are given in Table~\ref{tbl:TemperaturePowerLaws}.
Eq. (\ref{eq:DiscVariability}) is strictly applicable  for the $\tau_\nu >>1$ regime so
here $p=9/10$.

\begin{table}
\begin{tabular}{r|l}
Disc Model  & $p$ \\
\hline
Viscously Heated, Optically Thick  & $9/10$  \\
Viscously Heated, Optically Thin   & $3/5$ \\
Irradiated by central star   &  $1/3$
\end{tabular}
\caption{Values for $p$ in Eq. (\ref{eq:TradialPowerLaw}) for various disc models \citep{2007MNRAS.381.1280E}}.
\label{tbl:TemperaturePowerLaws}
\end{table}
\section{Total RPE and Role of  Resolving Power}

Eqs. (\ref{5}) and (\ref{eq:DiscVariability})
give the respective  azimuthal and radial contributions to the
RPE of the predicted luminosity in terms of disc model properties. They must be added in quadrature
to get the total luminosity RPE. The choice of $\xi$ can be taken to minimize the RPE but whether this minimum
can be compared with data depends on the resolving power of the instruments.
We explain these points further below.

\subsection{Minimizing the Intrinsic RPE}
First we determine the minimum  of Eq. (\ref{eq:deltarexpand}),
which applies in the optically thick $\tau_\nu >>1$ and $t_{obs}=0$ limit.
Differentiating and setting equal to zero gives
\begin{equation}
\frac{\xi_{\text{opt}}}{r} = 2 \parenfrac{\ved{}}{v_0}^{\frac{2}{3}} \parenfrac{\led{}}{r}^{\frac{1}{3}}.
\label{eq:xioptdefine}
\end{equation}
The optically thin case gives the same result because the optically thin corrections to $N_{rs}$ and $N_{fl}$
just produce an overall multiplicative factor of $l_{ed}/h$ to both terms in
Eq. (\ref{eq:deltarexpand}) and thus do not change its minimum.

Choosing $\xi=\xi_{opt}$ will minimize the RPE
of the theory.  Increasing $\xi$ would average over more turbulent fluctuations, but at the price of poorer spatial precision. Reducing $\xi$ would increase the random error due to a smaller number of eddies.
If  $\xi_{\text{opt}} < \led{}$,  we must use $\xi = \led{}$ instead.

Substituting Eq. (\ref{eq:xioptdefine}) into Eq. (\ref{eq:deltarexpand}), gives
\begin{equation}
\left . \frac{\Delta r}{r} \right |_{\text{min}}
=
\sqrt{3}
\parenfrac{\ved{}}{v_0}^{\frac{2}{3}}
\parenfrac{\led{}}{r}^{\frac{1}{3}},
\label{eq:deltarOpt}
\end{equation}
where we have used  $\ved{} \led{} = \alpha \cs{} h$.
Note that  $\xi_{\text{opt}}$ of Eq. (\ref{eq:xioptdefine})
can only be used to compare with observations if the observations are resolved on
scales $\xi\le \xi_{opt}$ since only then can the data be binned accordingly.
We will come back to this point. First we calculate the RPE for the luminosity using
$\xi =\xi_{opt}$.

Substituting $l_{ed}= \alpha^{1/2}h$ and $v_{ed}=\alpha^{1/2}c_s$ into
Eq. (\ref{eq:deltarOpt}), and in turn substituting the result into Eq. (\ref{eq:DiscVariability})
gives
\begin{equation}
\left . \frac{\Delta L_{\nu}}{L_{\nu}} \right |_{\text{rad,opt}} \approx \sqrt{3} \left| {5\over 2}-{7p\over 2} \right|\alpha^{\frac{1}{2}}
\parenfrac{h}{r}.
\label{eq:LvaryOpt}
\end{equation}
Combining results from Eqs. (\ref{5}) and (\ref{eq:LvaryOpt}) by adding the radial and azimuthal
contributions in quadrature,we have for the net minimum RPE in the $\tau_\nu >>1$ limit
\begin{equation}
\left . \frac{\Delta L_{\nu}}{L_{\nu}} \right |_{\text{opt}}\approx\left[{{3}\over 4}( 5-7p)^2\alpha \parenfrac{h}{r}^2 + {\alpha^{1/2}\over 6\pi} \parenfrac{h}{r}\right]^{1/2},
\label{eq:LvaryOptCom}
\end{equation}
where in the optically thick limit \citep{2007MNRAS.381.1280E} $p=9/10$ and
 $h/r$ is
 given by
 \begin{equation}
\frac{h}{r}  =
1.8 \times 10^{-2} \left(\frac{\dot{M}_{\ast}}{10^{-8} M_{\odot} yr^{-1}} \right)^{1\over 5} \left(\frac{M_{\ast}}{M_{\odot}} \right)^{-7\over 20} \left(\frac{r}{AU} \right)^{1\over 20}.
\label{depth}
\end{equation}

\subsection{Expressing Total RPE as function of Resolving Power}

For a telescope, we can only choose $\xi$ to be as small as the  spatial or spectral
resolution allows. If  the smallest
$\xi$ that the telescope can resolve is smaller than $\xi_{opt}$, the data can be binned
to create an effective $\xi=\xi_{opt}$ that minimizes the RPE.

We can estimate the radial averaging scale from the spectral resolving power of the telescope.
The resolving power is given by
\begin{equation}
R \equiv \frac{\nu}{\Delta \nu}.
\end{equation}
From the Wien displacement law, we know that a blackbody spectrum has $\nu_{\text{peak}} \propto T_{\text{peak}}$, so the effective resolving power is then  $R_{ef} \approx T / \Delta T$.
From Eq. (\ref{eq:TradialPowerLaw}), we then obtain
\begin{equation}
|\xi_{\text{ef}}| \approx \frac{r}{pR_{\text{ef}}}
\label{eq:xirealdefine}
\end{equation}
Substituting into Eq. (\ref{eq:deltarexpand}),  we find
\begin{equation}
\left . \frac{\Delta r}{r} \right |_{\text{ef}}\approx\left( 4 \alpha^{3/2} \parenfrac{h}{r}^{3} pR_{\text{ef}} + \frac{1}{4p^{2}R^{2}_{\text{ef}}} \right)^{\frac{1}{2}}
\label{eq:deltarReal}
\end{equation}
Using this in (\ref{eq:DiscVariability}) we obtain
\begin{equation}
\left . \frac{\Delta L_{\nu}}{L_{\nu}} \right |_{\text{rad,ef}}\approx \left |{ 5\over 2}-{7p\over 2 }\right |\left( 4 \alpha^{3/2} \parenfrac{h}{r}^{3} pR_{\text{ef}} + \frac{1}{4p^{2}R^{2}_{\text{ef}}} \right)^{\frac{1}{2}}.
\label{eq:LvaryReal}
\end{equation}
The azimuthal contribution Eq. (\ref{5}), does not depend on the frequency explicitly
so combining it in quadrature with Eq. (\ref{eq:LvaryReal})
gives  the total effective RPE
\begin{equation}
\begin{array}{r}
\left . \frac{\Delta L_{\nu}}{L_{\nu}} \right |_{\text{ef}}\approx\left[\left |{5\over 2}-{7p\over 2 }\right |^2\left(4 \alpha^{3/2} \parenfrac{h}{r}^{3} pR_{\text{ef}} + \frac{1}{4p^{2} R^{2}_{ef}} \right) + {\alpha^{1/2} \over 6 \pi}\left({h \over r}\right) \right]^{1\over 2}.
\label{eq:LvaryRealCom}
\end{array}
\end{equation}
It might seem surprising that increasing $R$ in Eq. (\ref{eq:LvaryReal})
 can give less precise results but recall that for a fixed $p$ and $r$ there is a one-to-one correspondence between $R_{ef}$ and $\xi$.  But as discussed above, there is a value $\xi=\xi_{opt}$  that minimizes the RPE.
There is therefore a corresponding value of $R_{ef}$  that minimizes the RPE and the two turn out to be
related by  $\xi_{opt} =  {r\over p R_{opt}}$.
Figure 1 illustrates these basic principles in  plots of  Eq. (\ref{eq:LvaryRealCom})
for different values of $\alpha$. Note that the minima are less pronounced as $\alpha$ is decreased.

Note that  if $R$ of the telescope is sufficiently large and the data are not more coarsely binned
by hand, the data could represent a radial resolution less than the eddy scale of the theory. The  $\alpha$-disc model based on turbulence becomes ill-defined below such scales.
Coarse graining of the data has to be considered carefully  before sensibly comparing
with a mean field theory.

\subsection{Optically thin case}
Finally, note that
for $\tau_\nu<<1 $  Eqs. (\ref{2})
(\ref{12a}) and (\ref{eq:NflDefine}) take a different values compared to the $\tau_\nu >>1$ case.
The overall calculation of the RPE follows that which led to (\ref{eq:LvaryRealCom}) with the analogous
result being
\begin{equation}
\begin{array}{r}
\left . \frac{\Delta L_{\nu}}{L_{\nu}} \right |_{\text{ef}}\approx\left[\left |{5\over 2}-{7p\over 2 }\right |^2\left(4 \alpha^{2} \parenfrac{h}{r}^{3} pR_{\text{ef}} + \frac{\alpha^{1/2}}{4p^{2} R^{2}_{ef}} \right) + {\alpha \over 6 \pi}\left({h \over r}\right) \right]^{1\over 2}.
\label{30thin}
\end{array}
\end{equation}
where in this $\tau_\nu <<1$ limit \citep{2007MNRAS.381.1280E}, $p=3/5$ and
\begin{equation}
\frac{h}{r}  =
2.8 \times 10^{-2} \left(\frac{\dot{M}_{\ast}}{10^{-8} M_{\odot} yr^{-1}} \right)^{1\over 10} \left(\frac{M_{\ast}}{M_{\odot}} \right)^{-2\over 5} \left(\frac{r}{AU} \right)^{1\over 5}.
\label{depththin}
\end{equation}
To obtain minimum (optimal) value of (\ref{30thin}) (the analogue of (\ref{eq:LvaryOptCom})),
we substiute in (\ref{30thin}) the value
$R_{ef}= R_{opt}= {r/p \xi_{opt}}$
and use $\xi_{opt}/r $ from Eq. (\ref{eq:xioptdefine})
which applies for both the optically thin and thick cases,

\begin{figure}
\centering
\includegraphics{TVThinRadius10auRange2000.pdf}
\caption{RPE in luminosity from turbulent dissipation  for optically thin case at $r=10$AU for the parameters
of LLRL 31 as a function of resolving power $R_{ef}$ for $\alpha = 0.1 (\text{Thick lined Blue}), 0.01 (\text{thin lined Red}), 0.001 (\text{Dashed})$. The higher $\alpha$ has the most prominent minimum, which does not correspond to the maximum resolving power but an intermediate value.}
\label{resolvingpower}
\end{figure}

\section{Implications for Multi-epoch Observations}

\begin{figure}
\includegraphics{RadialVariabilityThin.pdf}
\caption{RPE in luminosity from turbulent dissipation vs. radius $r$ for the optically thin case, using  $\alpha = 0.1 (\text{thick lined Blue}), 0.01 (\text{thin lined Red}), 0.001 (\text{Dashed})$. From left to right the three vertical dashed lines
correspond  to three specific radii whose orbit period equals different inter-epoch observation periods for  LRLL 31
\citep{2009ApJ...704L..15M}: 1 week ($r = 0.087$ AU) and 4 months($r = 0.58$ AU), and
 $\sim$ 25 months ($r=2$ AU).  For the respective inter-epoch periods,  these vertical lines mark the upper radial limit of observational validity of the curves, and above which the expected variability will be sharply reduced.
For example, if the inter-epoch interval is 25 months, then the RPE would predict the regime of validity for the curves to be $r< 2\AU$. Above $r> 2$ AU the prediction would be that $\Delta L_\nu/L_\nu$ falls near zero because the epochs would sample data from the same statistical ensemble.}
\label{fig:variabilityvsr}
\end{figure}

A basic application of the RPE  described above arises for multi-epoch observations
of a single object when the time scale for a  snapshot observation is shorter than
the eddy time at any radius. Emission from  radii
at which the time scale between epochs of observation
 exceeds the eddy turnover time  would then be expected to vary stochastically between epochs.
In contrast,  at radii larger than the that at which the eddy turnover time matches the epoch
interval time, successive observations sample the same member of the statistical
ensemble and no stochastic variability corresponding to the the local energy-dominating eddies would be expected.

\subsection{Lessons from the specific case  of LRLL 31}

Presently, multi-epoch data are sparse. Programs such as YSOVAR \citep{2010AAS...21542907S} will
eventually produce more data.  Nevertheless, to focus the practical application on a
on a specific object, we consider  IR Spitzer multi epoch observations
of the source LRLL 31 in the star-forming region IC 348
\citep{2009ApJ...704L..15M}
.  For reasons discussed below, this object
is most likely not showing its primary variability from  stochastic effects  we discuss but our analysis
to reach this conclusion below  illustrates how to use the concepts of the present paper to arrive at that conclusion.

 The stellar mass of LRLL 31 is $M = 1.8 M_{\odot}$ and the estimated accretion rate is
 $\dot{M} = 1.5 \times 10^{-8} M_{\odot} {\rm  yr}^{-1} ({\rm from\ }\text{Pa} \beta)$.

 For the low spectral resolution mode used to observe this source
 in the 5-38 micron range \citep{2009ApJ...704L..15M},
 the Spitzer spectral resolving power is $R\sim 200$. The observation times are minutes, well below
 the smallest relevant orbit time and consistent with setting $t_{obs}\sim 0$ in
 (\ref{12a}) and (\ref{eq:NflDefine}) and out assumption throughout the previous sections.

Figure 2 shows Eq. (\ref{eq:LvaryReal}) plotted as a function of  radius for three different values of
$\alpha$ and a value of $R_{ef}= 200$ for the optically thin case.

  The vertical lines show, for three separate choices of the inter-epoch time,
   the radii above which there would be no variability expected because the inter-epoch time is less than the eddy turnover time above these radii.

Tables 2  and  3  exemplify  how use of $R_{ef}=R$,
the resolving power for the telescope, compares to
the use of $R_{ef}=R_{opt}$, the value that  minimizes the luminosity RPE,
 for different choices of $r$ and $\alpha$.
For these tables, the epoch time is taken to exceed
the eddy turnover time at all radii. In reality, as emphasized by the vertical lines in Figure 2,
there would be negligible expected variability at radii above which the orbit time ($\sim$ eddy turnover time) exceeds a chosen inter-epoch observation interval.
For cases  in which  the telescope resolution (column 6) exceeds the optimal value (column 4)
the data can  in principle be binned such that  the binned data can be compared to the optimal $\alpha$-disc model prediction.  The tables highlight two trends in the Tables:  (i) the larger the value of $\alpha$ the larger the predicted
variability between epochs at a given radius and (ii) the larger the radius for a fixed $\alpha$ the larger the predicted variability.

Muzerolle (2009) argue  that the disc of LLRL 31,  at  $ r < 10$ AU,
 is optically thin and substantially cleared out (perhaps via the action of planet) below $10$AU but
optically thick for a least some range of larger radii. Thus for this object,
the values in the Table 2 would apply for $r<10$AU and the values of table 3 may apply for
some range of larger radii above which the disc  again becomes optically thin.
The value at which the latter occurs requires
more detailed modeling.

 There is also  a transition radius
above which viscous dissipation would be dominated by irradiation from the central star \citep{2006ApJ...638..897S,2007MNRAS.381.1280E},
namely
\begin{equation}
r_{c} = 0.3 AU \left(\frac{\dot{M}_{\ast}}{10^{-8} M_{\odot} yr^{-1}} \right)^{\frac{3}{4}} \left(\frac{M_{\ast}}{M_{\odot}} \right)^{\frac{7}{6}} \left(\frac{L_{\ast}}{L_{\odot}} \right)^{-\frac{5}{6}}.
\end{equation}
For $r>r_c$  the photospheric temperature induced by viscous dissipation falls
below that associated from surface dust illumination by the central star.
The values of $p$ would then flatten to $p=1/3$ (see Table 1) and irradiation would dominate
at all larger radii.
Disc emission resulting from  stellar irradiation
 is less influenced by the stochasticity of turbulence than emission via turbulent dissipation.
Since the local blackbody emission varies with $T^3$, we would expect a much reduced stochastic variability in regions of the disc dominated by emission from reprocessed starlight.
This implies that even for  epoch
intervals longer than the orbit periods at $r>r_c$ that the stochastic variability
discussed herein would be reduced compared to regions of $r<r_c$
It is important to emphasize that the RPEs calculated in  Tables 2 and 3  is a stochastic variability
corresponding to that of  mean field $\alpha$-disc
resulting from the disc luminosity contribution by stochastic process
such as  turbulent viscous dissipation. As applied to the specific object LLRL 31 therefore,
these would only amount to a variability in a subdominant contribution to the total emission above $r_c$.

 An absence of observed stochastic fluctuations
 for $r>r_c$ would support  the expectation that the primary emission in those regions is not the result of turbulent dissipation. For $r<r_c$ the absence of stochastic variability would be expected only
if  the inter-epoch time scale is shorter than the characteristic eddy turnover time at the radius
dominating the radiation at the particular frequency measured.
Note also that  the amplitude of the predicted stochastic variability
would  not depend on the epoch interval as long as the epoch interval is larger than the orbit time
at the radius producing the corresponding emission.
Other systematic non-axisymmetry or variable accretion rates would  add to the variability and could
be distinguished form stochastic effects given enough multi-epoch observations

For the specific case of LRLL 31, the largest time scale between epochs  is of order 4 months
which means that according to Fig. 2, the stochastic variability would be expected
only for $r< 0.58$ AU. This is well within the dust depleted inner region of this transitional object where
the disc is optically thin and where  the emission would come from dissipation internal to the disc rather than
reprocessed starlight.  Thus there would be some stochastic contribution to the measured variability between epochs predicted by the RPE of the mean field models.
{ However,  the double-digit percentages of  the observed variability reported by
\citet{2009ApJ...704L..15M}
 in LRLL 31 are larger than the values predicted in our Tables that arise  from stochastic variability alone }

In addition, the presence  of a single pivot wavelength  of 8.5 microns (see Fig 1 of
\citet{2009ApJ...704L..15M})
  at which there is no variability in this source
and above which the sign of the variability between epochs changes sign from that at low wavelengths
strongly suggests the influence of a global geometric feature such as a thick wall or a warp.
 In short
the dominant variability in LRLL 31 object likely comes from systematic rather than stochastic effects.

\section{Conclusions}

For $\alpha$  models of turbulent accretion discs we have derived a predicted
RPE in the luminosity as  function of frequency, focusing on the case
in which the observational exposure times  are small compared ti other dynamical time scales.
We  have discussed the implications for interpreting multi-epoch spectral observations of spatially unresolved protoplanetary discs.

The RPE depends on the radial scale of averaging and there exists an
optimal scale that minimizes the error which also corresponds to an optimal
spectral resolution when the correspondence between peak frequency of emission and radius
is made.  The data can be binned to compare with the theory of minimum
RPE  only when the spectral resolution of the telescope exceeds the value $R_{opt}$
 which minimizes the RPE.
 An optimal resolution exists because if the spectral resolution is too high then the instrument ends up
 sampling noise and if the spectral resolution is too low, then the instrument
 samples overly coarsely binned regions of the disk.

 The stochastic variability expected between multi-epoch observations of discs is an implicit prediction of
alpha disc theory and  can be directly compared to observed variabilities.
Complementarily,  the nature of the observed variability can  be used to constrain whether or not the behavior
of the disc is consistent with luminosity produced by stochastic or systematic processes.

The predicted stochastic  variability from $\alpha$ discs has several distinct characteristics that would signature its prevalence:

$\bullet$   The stochastic variability is predicted primarily for regions where the luminosity is dominated by internal
turbulent dissipation in the disc. For regions dominated by surface dust-reprocessed starlight, the stochastic variability could be expected to be strongly reduced.

$\bullet$.  For regions of the discs dominated  by from turbulent dissipation,
the strength of the stochastic variabilty would increase gradually with decreasing emission
frequency (and thus increasing radius) down to a critical frequency
below which the associated disc radii have orbital times exceeding the inter-epoch observation time. There the expected variability from the RPE would drop sharply.

$\bullet$  For inner regions of the discs in transitional objects, where emission may be dominated by
turbulent dissipation but the discs are optically thin, the stochastic variability would be expected
to be lower than that expected from the same radii for optically thick turbulent dissipation dominated discs.

In applying the basic ideas herein to the specific case of LRLL 31, we find that systematic variability NOT stochastic variability dominates.  This is indeed consistent with the conclusion of \citet{2009ApJ...704L..15M}.
As more multi-epoch observations of broader samples of discs at different stages in their lifetimes
are obtained,  the present work  may help provide a tool to distinguish stochastic processes
in discs from systematic dynamical changes and the time scales on which these occur.

Two overall lessons from this analysis are: (1) the precision of axisymmetric accretion theories can be quantified and deviations the alpha-accretion disc theory
and observations at a given epoch cannot rule out the theory if the deviations between
epochs exhibit stochastic behavior and fall within the expected RPE of the theory.
The theory may be incomplete but because a mean field theory for a turbulent system is intrinsically imprecise,
but the imprecision must be quantified so that the user realizes its limitations. (2)
Arbitrarily high spectral resolution can lead to  misleading comparisons between
theory and observation and binning the data may be necessary to compare the data with the theory of
minimum RPE.

\begin{table}
\begin{tabular}{c|c|c|c|c|c}
r(au) & $\alpha$ & $R_{\text{opt}}$ & ${\Delta L_\nu \over L_\nu}$
(\%) & $R_{\text{tel}}$ & ${\Delta L_\nu \over L_\nu}$(\%) \\
\hline
0.1 & 0.001 & 1832 & 0.08754 & 200 & 0.09226 \\
& 0.01 & 579.3 & 0.2781 & 200 & 0.2816 \\
& 0.1 & 183.2 & 0.8914 & 200 & 0.8915 \\
1 & 0.001 & 1156 & 0.1103 & 200 & 0.1139 \\
& 0.01 & 365.5 & 0.3514 & 200 & 0.3531 \\
& 0.1 & 115.6 & 1.135 & 200 & 1.144 \\
10 & 0.001 & 729.3 & 0.1392 & 200 & 0.1417 \\
& 0.01 & 230.6 & 0.4449 & 200 & 0.4451 \\
& 0.1 & 72.93 & 1.454 & 200 & 1.513 \\
20 & 0.001 & 634.9 & 0.1493 & 200 & 0.1515 \\
& 0.01 & 200.8 & 0.4780 & 200 & 0.4780 \\
& 0.1 & 63.49 & 1.570 & 200 & 1.661 \\
\end{tabular}
\caption{Optically Thin Case.
In column 4, we list the RPE varibilities calculated using $R_{ef}=R_{opt}$ and
and in column 6  the values corresponding to the Spitzer telescope low mode   $R_{ef}=R_{tel}\sim  200$.
The RPEs  are computed from Eq. (\ref{30thin}) for different $r$ and $\alpha$ and  $h/r$ from Eq. (\ref{depththin}) using $p=0.6$.
The optimal resolution is calculated from
Eq. (\ref{eq:xioptdefine})  as described below
Eq. (\ref{depththin}).
Note that if $R_{opt}< R_{tel}$ data can be binned to compare with the $\alpha$-disc model of minimal RPE.}
\label{tbl:variabilityalpha}
\end{table}
\begin{table}
\begin{tabular}{c|c|c|c|c|c}
r(au) & $\alpha$ & $R_{\text{opt}}$ & ${\Delta L_\nu \over L_\nu}$(\%) & $R_{\text{tel}}$ & ${\Delta L_\nu \over L_\nu}$(\%) \\
\hline
0.1 & 0.001 & 625.8 & 0.6935 & 200 & 0.7111 \\
& 0.01 & 249.2 & 1.116 & 200 & 1.118 \\
& 0.1 & 99.19 & 1.835 & 200 & 1.881 \\
1 & 0.001 & 557.8 & 0.7355 & 200 & 0.7510 \\
& 0.01 & 222.1 & 1.186 & 200 & 1.186 \\
& 0.1 & 88.40 & 1.958 & 200 & 2.030 \\
10 & 0.001 & 497.1 & 0.7802 & 200 & 0.7935 \\
& 0.01 & 197.9 & 1.260 & 200 & 1.260 \\
& 0.1 & 78.79 & 2.091 & 200 & 2.200 \\
20 & 0.001 & 480.2 & 0.7942 & 200 & 0.8069 \\
& 0.01 & 191.2 & 1.284 & 200 & 1.284 \\
& 0.1 & 76.10 & 2.133 & 200 & 2.256 \\
\end{tabular}
\caption{same as Table 2 but for optically thick case. The RPEs  are computed from Eqs. (\ref{eq:LvaryOptCom}) and  (\ref{eq:LvaryRealCom}) for different $r$ and $\alpha$ and for the optically thick version of $h/r$ in Eq. (\ref{depth}) and $p=0.9$.}
\label{tbl:variabilityalphathick}
\end{table}