{'fraction_non_alphanumeric': 0.0355237770968108, 'fraction_numerical': 0.013430181969507812, 'mean_word_length': 4.440214035809837, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Disks as Inhomogeneous, Anisotropic Gaussian Random Fields

We model astrophysical disk surface brightness fluctuations as an inhomogeneous, anisotropic, time-dependent Gaussian random field. The field locally obeys the stochastic partial differential equation of a Mat\'ern field, which has a power spectrum that is flat at large scales and falls off as a power law at small scales. We provide a series of pedagogical examples and along the way provide a convenient parameterization for the local covariance. We then consider two applications to disks. In the first we generate a movie of a disk. In the second, by integrating over a movie of a disk, we generate synthetic light curves and show that the high frequency slope of the resulting power spectrum depends on the local covariance model. We finish with a summary and a brief discussion of other possible astrophysical applications.


INTRODUCTION
Astrophysical disks experience surface brightness fluctuations that can in principle be predicted by numerical solution of well-known governing equations. Three dimensional physical simulations are computationally expensive, however, and may still be missing important physical processes and be unable to resolve important physical lengthscales. What if the main features of the fluctuations could be captured in a simpler, easier-to-compute statistical model? In this paper we consider a model in which the surface brightness fluctuations are treated as a Gaussian random field (GRF).
GRFs are widely used in astrophysics to model correlated noise. In cosmology they are used to model initial conditions (e.g. Bardeen et al. 1986). In studies of quasars the light curve of the unresolved source is commonly modeled as a damped random walk (Kelly et al. 2009;MacLeod et al. 2010), which is a Gaussian process. In both cases the field is homogeneous and so a realization can be easily generated by drawing uncorrelated Fourier amplitudes from a Gaussian distribution with variance given by the power spectrum.
In astrophysical disks the correlation length and correlation time likely vary with local radius, possibly by orders of magnitude. It is not easy to see how Fourier techniques might be generalized to capture this inhomogeneity. In addition, the surface brightness correlation is likely anisotropic, as for example in the trailing spiral structures that populate flocculent spiral galaxies. How can one efficiently generate a realization of an inhomogeneous, anisotropic GRF to model fluctuations in astrophysical disks?
In practice we wish to sample the GRF on a finite mesh of points, which is equivalent to sampling a multivariate Gaussian with a prescribed covariance matrix. For small sets of sample points one might sample directly from the multivariate Gaussian. A simple but computationally suboptimal way of doing this is to transform to a basis that diagonalizes the covariance matrix (Karhunen-Loève transformation), draw independent amplitudes, and transform back to the original basis. A more efficient approach is to use a Cholesky decomposition of the covariance matrix. This approach is not practical for the large numbers of points considered here.
One might also generate an inhomogeneous, anisotropic GRF by distorting an initially homogeneous, isotropic GRF through a coordinate transformation (Sampson & Guttorp 1992). This is, effectively, how weak lensing acts on the microwave background. The technique is limited to certain geometries and boundary conditions. Suppose, for example, that one is modeling a flocculent spiral galaxy in which the number of spiral arms changes with radius. There is no coordinate transformation that can map the associated covariance function onto a rectangular, periodic domain with a homogeneous, isotropic GRF.
The technique we use here, which is borrowed from geostatistics, generates an inhomogeous, anisotropic field as the solution to a stochastic partial differential equation (SPDE; Whittle 1954Whittle , 1963Lindgren et al. 2011;Fuglstad et al. 2015). This method has a well-defined notion of a local covariance function, although the covariance function is limited to a particular -but useful -functional form. To our knowledge this is the first application of this technique in an astronomical context. The plan of this paper is as follows. In §2 we describe the SPDE technique. In §3 we provide simple examples of anisotropy, inhomogeneity, and time-dependent fields. In §4 we apply the technique to generation of synthetic movies of differentially rotating disks. In §5 we explore an application to generating broad band noise for disk light curves. §6 contains a summary and a guide to the main results.

GAUSSIAN RANDOM FIELDS AND SPDES
A random field f on a space X is a function such that for every x ∈ X, f (x) is a random variable. A Gaussian random field (GRF) is a random field such that the joint probability distribution on any set of N points {x 1 , . . . , x N } is an N -dimensional multivariate Gaussian distribution with mean vector µ(x n ) and covariance matrix C(x n , x m ), n, m ∈ {1, . . . , N }. Thus, a GRF is completely defined by a mean function µ(x) and a covariance function C(x, x ), x, x ∈ X. Since µ can always be subtracted off, we assume without loss of generality that all GRFs have zero mean.
GRFs have a simple definition and useful analytic properties and have therefore been used as a statistical model for a wide variety of phenomena. The central limit theorem implies that a superposition of independent, identically distributed random fields is a GRF. Thus, GRFs arise naturally in cosmology, since primordial fluctuations are thought to be the result of an interaction-free (or nearly interaction-free) scalar field during the inflationary period.
Homogeneous, isotropic GRFs such as those that arise in cosmology are readily generated in a Fourier basis. Homogeneity implies that the covariance function depends only on ∆x, so where f (x) is the expected value of f (x) andf (k) is the Fourier transformf (k) = f (x)e ikx dx. Thus, the Fourier transform of the covariance function C(∆x) is the power spectrum P (k). Thus one can generate a realization of a homogeneous GRF by drawing independent, normally distributed Fourier amplitudes with variance given by P (k). If the GRF is also isotropic then the power spectrum depends only on |k| (e.g. Bardeen et al. 1986). Inhomogeneity (or non-stationarity) 1 makes the generation of GRFs more difficult because the Fourier modes are no longer delta-correlated. Since the covariance matrix on a rectangular mesh with N points along each of d dimensions has N 2d components and Cholesky decomposition of an M xM covariance matrix is an O(M 3 ) operation, generating a GRF by directly sampling the resulting multivariate Gaussian would require O(N 3d ) operations. This is too costly for even modest grid size. Various approaches have been developed to make statistical modeling using inhomogeneous GRFs feasible, including low-rank approximations (e.g. Cressie & Johannesson 2008) or covariance tapering (e.g. Furrer et al. 2006).
The method we use here (Lindgren et al. 2011) takes advantage of a relationship between a particular covariance function known as the Matérn covariance, a stochastic partial differential equation (SPDE), and Gaussian Markov random fields (GMRFs). The Matérn covariance is stationary and isotropic and has the form x, y ∈ R d . Here λ is a scaling parameter, ν is a differentiability parameter, and K ν is the modified Bessel function of the second kind, order ν. For small r/λ ≡ |x − y|/λ, C ν ∼ σ 2 (1 − Ar 2ν ), where A is a constant. Notice that for ν = 1/3 this matches the covariance of a passive scalar in Kolmogorov turbulence. 2 For large r/λ, C ν → 0, that is, the field decorrelates on scales large compared to λ.
GMRFs are discrete Gaussian fields (e.g., random fields sampled on a Cartesian lattice) where the probability distribution of the field at a point depends only on its neighbors, that is, the field has the Markov property. Most numerical representations of Gaussian fields are discrete fields. Rue & Tjelmeland (2002) demonstrated that a Matérn field is well approximated by a GMRF, and Lindgren et al. (2011) showed that such a GMRF can be constructed efficiently using SPDEs for certain values of ν. In what follows we specialize to ν = 2 − d/2, where d is the number of dimensions; a generalization is discussed in the Appendix.
A field f (x) with a Matérn covariance for ν = 2 − d/2 can be expressed as a solution to the SPDE where d is the number of dimensions, N is a normalization constant, σ 2 is the variance of the field, and W is a standardized Gaussian white noise process (Whittle 1954(Whittle , 1963. The power spectrum for solutions to (3) is Since Gaussian white noise is easy to generate on a lattice, this changes the computational task from manipulating a large covariance matrix to solving a finite difference approximation to an elliptic partial differential equation. The numerically obtained solution of the SPDE on a finite grid or irregular lattice is a GMRF that represents the underlying GRF.
The key advantage of the SPDE method is that it can be readily generalized to inhomogeneous, anisotropic fields at little extra cost. By taking the SPDE as the definition of an inhomogeneous Matérn field, we can construct anisotropic and non-stationary variants by introducing position-dependence in the parameters, e.g. ( where the matrix Λ(x) introduces position-dependent anisotropy and correlation lengths, σ 2 is the local variance of the field, N is a normalization constant, and x ∈ R 2 or R 3 (since ν > 0). Because varying the correlation lengths also changes the variance of the field, the variance at each point is approximately normalized by the factor of (det(Λ)) 1/4 . This normalization holds when the correlation lengths varies on scales large compared to the correlation lengths.

EXAMPLES
Before applying this model to astrophysical problems we consider a sequence of models demonstrating the ability to model anisotropy, inhomogeneity (see also Fuglstad et al. 2015), and time dependence in two spatial dimensions.
The local covariance is controlled by Λ, which we parameterize as where λ 1 and λ 2 are the correlation lengths along the axes specified by the 2-d spatial unit vectors u 1 and u 2 and u 1 · u 2 = 0. Notice that (detΛ) 1/2 = λ 1 λ 2 . If the field is homogeneous, the resulting covariance function is where Evidently Λ −1 acts as a metric on the space. If λ 1 = λ 2 = λ (i.e. Λ = λ 2 I) the SPDE solution is a homogeneous, isotropic GRF. Anisotropy is introduced by choosing an anisotropy direction u 1 and setting λ 1 = λ 2 . Inhomogeneity is introduced by allowing Λ to vary across the domain.
The model can be extended to include time variation by building a GRF with two space and one time dimension, and a suitable specification for Λ (the Appendix describes an alternative approach to time dependence). Here, the vector x = (t, x, y). Introducing the velocity field v = (0, v x , v y ) and correlation time λ 0 , we write so that where ∆x = (∆t, ∆x, ∆y), u 0 = (1, v x , v y ), u 1 = (0, cos θ, sin θ), and u 2 = (0, −sin θ, cos θ). Notice that while u 1 and u 2 are still unit vectors, u 0 is not. Furthermore, u 0 is not orthogonal to u 1 or u 2 . Nevertheless, (detΛ) 1/2 = λ 0 λ 1 λ 2 .
As in the 2D case, Λ acts as a metric. In particular, along a spatial slice ∆x = (0, ∆x, ∆y), the form of r is identical to the 2D case, and thus λ 1 and λ 2 are the correlation lengths along the major and minor spatial axes at any point in time. On the other hand, setting ∆x = (∆t, v x ∆t, v y ∆t) = ∆tu 0 , we have s 2 = (∆t/λ 0 ) 2 . Thus λ 0 corresponds to the correlation time following the flow. Figure 2 shows an example that takes the homogeneous, anisotropic model from above and generates a time dependent model with a velocity field in the −x direction, advecting it to the left. Figure 2. A realization of an anisotropic, homogeneous field advected to the left by a velocity field. While the field is generated probabilistically, the model encodes a spatiotemporal correlation that is aligned with the velocity field. Thus, across a single correlation time, features appear to be moving across the domain before decorrelating. Notice that the model (9) is quite flexible and can be used in conjunction with the SPDE technique to introduce velocity fluctuations around any velocity field. The Appendix describes a generalization to three and more spatial dimensions.

APPLICATION TO RESOLVED DISKS
Next we consider a model for a resolved disk with surface brightness fluctuations. Possible applications including disk galaxies, disks around young stars (which appear to be notably lacking in turbulence), and planetary rings. Here we are motivated by Event Horizon Telescope observations of the disk around the black hole in M87 ((EHT) Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f). We will consider a simple image-plane model that assigns mean velocities to points on the image. The model does not include radiative transport, lensing, or other key physical processes. The local covariances are drawn from numerical simulations of disks, and we adopt the drastic simplification that the mean velocities follow a Keplerian profile.

Local Correlations in Disks
The local covariance of velocity fluctuations in disks can be derived from local model simulations of disks. Guan et al. (2009), for example, give an analytic fit for the covariance function on a plane lying at fixed altitude in a local numerical model of a disk. Their power spectrum is constant (white noise) at small k and forms a tilted ellipse at large k, asymptotically decaying as k −11/3 . Thus, the local power spectrum can be modeled by which corresponds to a covariance function where K 1/3 is a modified Bessel function of the second kind and s 2 = ∆x · Λ −1 ∆x. Here, x is a 3D spatial vector, and Λ parameterized as a spatially 3D extension of §3.
For simplicity we work in 2D, in the midplane of the disk, and approximate the Guan et al. (2009) covariance with a ν = 1/2 Matérn covariance. We then extend the spatial covariance to 2+1D as described in §2 and §3. This preserves key features of the field: an anisotropic covariance function with steep decay at small scales, decorrelation at large scales, and a power spectrum of the form here k = (ω, k x , k y ).

Global Model
We can now construct a time-dependent global disk model. First, we produce a realization of the fluctuation field f , which we interpret as fractional variation in disk surface brightness, using (9).
The velocity field is v = Ω Kẑ × x, with Ω K ∝ |x| −3/2 , corresponding to Keplerian rotation. The major axis of the correlation tensor u 1 is chosen to lie at a contant 20 • angle to a circle of constant radius. This is the opening angle of spiral features in the GRF, and is consistent with local model results (Guan et al. 2009). We let λ 0 , λ 1 , and λ 2 vary with radius but not azimuth, since the model is on average axisymmetric. It is natural to take λ 0 ∝ 1/Ω K . The correlation lengths are expected to be proportional to the disk scale height, which is in turn assumed proportional to the local radius (this assumption can be relaxed), so λ 1 ∝ r and λ1 λ2 = const. Finally, we set σ = 1. This completes the specification of f . We produced a realization by generating Gaussian white noise on a grid and solving the elliptic equation 3. We used the preconditioned conjugate gradient method with a semicoarsening multigrid preconditioner provided by hypre, a library of parallel solvers for linear systems.
To generate a movie of the disk we need to relate f to the surface brightness µ. The mean surface brightness of the disk is given by an "envelope function" g(r). In this example x ≡ r 0 /r and g(r) = x 4 exp(−x 2 ).
This has a "shadow" in the middle, like the EHT image of M87, and a surface brightness that drops off as r −4 at large radius. The surface brightness is then The control parameter (or function) n controls the fluctuation amplitude (this could also be done when the random field is generated). The envelope function is independent of the GMRF and can be chosen for convenience. Figure 3 shows the resulting non-stationary, anisotropic GRF with differential rotation and position-dependent correlation lengths. The model was run on a Cartesian mesh in x, y, t at a resolution of (N t , N x , N y ) = (1024, 256, 256), with periodic boundary conditions in time and Dirichlet boundary conditions in space. It was run on an Intel Xeon 6140, and required less than two minutes to converge to a solution. The fluctuation field parameters are λ 0 (r) = 2π Ω K , where Ω K = 1 r 3/2 , λ1 r = 5, and λ2 λ1 = 0.1. The envelope has a single parameter r 0 that was set to 1/10th of the width of the grid (≈ 26 mesh points).
The code used to produce this example is publicly available at https://github.com/AFD-Illinois/inoisy. Figure 3. A realization with differential rotation and position-dependent correlation lengths. The field, displayed in log scale, goes to zero in the center and near the edges due to the envelope function. As before, the appearance of advection is created by a spatiotemporal correlation that is applied along the velocity field, which for this realization is Keplerian. The position-dependent angle in the spatial correlation creates the spiral arm features.

APPLICATION TO UNRESOLVED DISK LIGHT CURVES
AGN light curves are commonly modeled as a damped random walk (Ornstein-Uhlenbeck process) (Kelly et al. 2009;MacLeod et al. 2010). Ornstein-Uhlenbeck processes are the only nontrivial stochastic processes that are Gaussian, Markov, and stationary. The power spectrum of an Ornstein-Uhlenbeck process is where ω is the angular frequency and τ is a characteristic timescale. The power spectrum is flat (white noise) at low frequency and scales as ω −2 at high frequency. Not all AGN light curves are well modeled by a damped random walk, however. Some have steeper spectra at high frequency (e.g. Smith et al. 2018). Can we use the GRF model of 4.2 to produce a more general model that connects the local structure of the disk, encoded in Λ, to the power spectrum of the light curve? A light curve is a time series constructed by integrating the source's fluctuating surface brightness over space at each instant. For a GMRF realization f (t, x), where here x is a 2D spatial vector, the light curve is where g(x) is the envelope function.
First, consider a simplified model calculation that shows that this integration can lead to interesting, nontrival results. If n is large then where L 0 is the luminosity when f = 0. To make the problem analytically tractable we use a stationary, isotropic GRF generated using the method of §3 (Equation 3 solved in 2 spatial dimensions plus time), with power spectrum Here N is a normalization constant. The power spectrum of the light curve P L (ω) is the k = 0 mode of the power spectrum P h (ω, k) of h(t, x) = f (t, x)g(x)/n. Then, using a hat to denote the Fourier transform, Using a Gaussian envelope the integrals can be done analytically and the light curve power spectrum is whereω = ω/ω 0 , q 2 = (σk 0 ) 2 , and Q 2 = q 2 (1 +ω 2 ). Whenω 1, P L approaches a constant, and when q 2 (1 +ω 2 ) 1, the power spectrum asymptotes to Thus, similar to the original power spectrum P f , the light curve's power is flat at low frequencies and falls off as ω −4 at high frequencies. However, in the intermediate range whereω 1 but Q 1, This intermediate regime occurs when q 2 = (σk 0 ) 2 = 2πσ λ 2 is small, i.e. when the width of the envelope function is much smaller than the correlation length. Figure 4 shows the shape of P L with varying q 2 . As q 2 becomes smaller, the intermediate regime becomes more prominent. The light curve power spectrum can also be calculated analytically for the more realistic envelope function applied to the same homogeneous, isotropic GRF. Like the envelope used in §4.2, this envelope is ring-like with a depression in the center, although it instead goes as r 3 at small r and falls off exponentially at large r.
The power spectrum of the light curve produced by (24) has the form where a i are integer coefficients andω, q, and Q are defined as before. As for the Gaussian envelope, the power spectrum is constant forω 1 and falls off as ω −4 at high frequencies, but in the intermediate regime whereω 1 but Q 1 is small, the power spectrum falls as ω −3 . Once again, this regime is determined by q, the ratio of the envelope width to the correlation length. When the envelope is much larger than the correlation length, we recover the slope of −4, but as the envelope becomes narrower, the power spectrum becomes shallower. Now consider the lightcurve of the global disk model generated in §4.2. Figure 5 shows the light curve and power spectrum of the realization shown in Figure 3. The slope of the power spectrum is shallower than ω −4 due to the combined effect of the envelope function and the inhomogeneity of the underlying GRF. Because the fluctuation field is anisotropic and inhomogeneous, the preceding simple models do not apply, but the qualitative behavior is similar.
It is interesting to ask whether the shape of the power spectrum, which is comparatively easy to observe, contains information about the underlying model. Figure 6 shows how the high frequency slope of the lightcurve power spectrum for a family of realizations similar to §4.2, with λ 1 /r varying from 5 to 195. The slope is averaged over 5 realizations. All other parameters are identical to §4.2. Evidently the slope of the power spectrum exhibits behavior that is consistent with the simple models described above: the power spectrum becomes shallower as the ratio of the correlation length to the envelope width increases.

SUMMARY
We have explored a model in which surface brightness fluctuations on an astrophysical disk are treated as a Gaussian random field. Realizations of the anisotropic, inhomogeneous fluctuation field can be generated by solving the stochastic partial differential equation 3. We provided pedagogical examples of anisotropic, inhomogeneous, and time-dependent anisotropic inhomogeneous fields in §3. The method requires that one specify an anisotropy tensor at every point in the domain. The key, useful result is a simple parameterization of the anisotropy matrix (9): one need only specify the anistropy and orientation of the covariance function in the two spatial dimensions, a velocity field, and a correlation time.
We applied the method to realize time-dependent, resolved images of a statistical disk model. The example shown in Figure 3 is difficult to distinguish from an animation of disk flow based on a physical simulation. The method presented here enables an inexpensive statistical simulation of a disk (even in three spatial dimensions, as described briefly in the Appendix); the realization shown here was generated in a few minutes.
[The physical inputs for the model are the correlation lengths and correlation times as a function of radius, which may be determined from simulations] Our statistical disk models provide a complement to physical simulations. Although physical simulations solve the governing equations and can therefore be predictive, they are subject to uncertainties related to physical and numerical approximations. Statistical models -if they can provide a good approximation to the simulations for particular parameter values -enable one to treat the physical simulations as a point in a larger parameter space of models, and therefore provide a universe of models to test physical simulations against.
We also used the statistical disk model to generate light curves for unresolved disks, and found that the models naturally produce an f −2 − f −3 power spectrum. This provides a means of connecting the space and time correlations of surface brightness fluctuations, which might be measured in a local model simulation, to the lightcurve.
[ Evidently it is interesting to examine how to extract model parameters from a dataset (inference). ] The technique used here seems likely to be useful elsewhere in astrophysics. For example, in modeling turbulent fluctuations in the interstellar medium (Sale & Magorrian 2018), providing realizations for turbulent inflow boundary conditions, or providing statistical models for any turbulent flow in which a mean flow field is known and for which a local covariance can be derived.

ACKNOWLEDGMENTS
This work was supported by the NSF grants AST-1716327, OISE-1743747, and a Romano Professorial Scholarship. We thank Gil Holder and Aviad Levis for their useful comments, and the referee for the helpful report that greatly improved the paper.
One motivation for considering a higher order model like this is for modeling a time-dependent process in three spatial dimensions (d = 4); for ν = 1/2 and d > 3 the variance of the GRF diverges, and the ν = 1/2 model is unsatisfactory , although this defect can be repaired by replacing the white noise process W by a red noise process, possibly realized as the solution to a separate SPDE. Notice that any finite difference operator with a compact stencil will have the Markov property, and the associated SPDE generates a GRF.

B. ALTERNATIVE APPROACH TO TIME DEPENDENCE
In §3 we introduced time dependence by treating the time coordinate on the same footing as the space coordinates, and orienting one axis of the correlation ellipse along a velocity vector. Here we consider a distinct procedure based on the SPDE where τ is a characteristic timescale and A is a constant. This SPDE was proposed by Lindgren et al. (2011) using slightly different notation. The power spectrum is from which one can see (because the integral over ω, k diverges in two spatial dimensions) that the pointwise variance of this model is formally infinite. The variance of the associated GMRF on a finite grid is not, however. We have implemented this form in the noisy code, available at https://github.com/AFD-Illinois/inoisyB. Another SPDE that produces a field that is equivalent to the example in §3 is where d is the number of spatial dimensions. The power spectrum is P ω,k ∝ τ 2 (ω − k · v) 2 + 1 + λ 2 k 2 ) 2 −1 .
Evidently the associated GRF has finite variance in two spatial dimensions, but not in three.