Mapping dust in the giant molecular cloud Orion A

The Sun is located close to the Galactic mid-plane, meaning that we observe the Galaxy through significant quantities of dust. Moreover, the vast majority of the Galaxy's stars also lie in the disc, meaning that dust has an enormous impact on the massive astrometric, photometric and spectroscopic surveys of the Galaxy that are currently underway. To exploit the data from these surveys we require good three-dimensional maps of the Galaxy's dust. We present a new method for making such maps in which we form the best linear unbiased predictor of the extinction at an arbitrary point based on the extinctions for a set of observed stars. This method allows us to avoid the artificial inhomogeneities (so-called 'fingers of God') and resolution limits that are characteristic of many published dust maps. Moreover, it requires minimal assumptions about the statistical properties of the interstellar medium. In fact, we require only a model of the first and second moments of the dust density field. The method is suitable for use with directly measured extinctions, such as those provided by the Rayleigh-Jeans colour excess method, and inferred extinctions, such as those provided by hierarchical Bayesian models like StarHorse. We test our method by mapping dust in the region of the giant molecular cloud Orion A. Our results indicate a foreground dust cloud at a distance of 350 pc, which has been identified in work by another author.


INTRODUCTION
The giant molecular clouds Orion A and B are sites of continuing star formation.Indeed, they are the regions closest to Earth in which high-mass (M > 8 M ⊙ ) star formation is occurring.Both Orion A and B are filamentary structures.Orion A, which is comet-like, consists of a dense head (containing the Orion Nebula, at galactic coordinates l = 206 deg, b = −16.4deg, and its largest cluster of young stars, the Orion Nebula Cluster) and a diffuse tail extending south-west on the sky for several degrees.Orion B is more uniform.It contains the Flame Nebula at its eastern extremity (at galactic coordinates l = 209 deg, b = −19.4deg) and extends north-west on the sky, again for several degrees.
Photometric surveys of young stars and young stellar objects, which trace the gas of the molecular cloud in which they have been concieved and born, have charted the gross on-sky distribution of Orion A. These surveys include those by: Meingast et al. (2016), Meingast et al. (2018), and Großschedl et al. (2019), as part of the VISION surveys, which used near-infrared observations made by the VISTA telescope augmented by visible and near-infrared observations by Pan-STARRS and mid-infrared observations by Spitzer; Megeath et al. (2012Megeath et al. ( , 2015)), which used mid-infrared observations made by Spitzer and near-infrared observations made by 2MASS; Carpenter (2000), which used near-infrared observations made by 2MASS; and Wright et al. (2010), which used mid-infrared observations by WISE.Astrometric surveys of the same young stars and young stellar objects have allowed us to chart the gross three-⋆ E-mail: a.gration@surrey.ac.uk dimensional distribution of Orion A. These surveys include those by Kounkel et al. (2018) and Großschedl et al. (2018), using Gaia Data Release 2 (DR2).They have established that the head of Orion A is a roughly spherical object, with a diameter of 15 to 20 pc, lying at a distance of s = 400 pc at l = 209 deg, b = −19.5 deg, and that its tail extends away from the Sun for a distance of 75 pc at an angle of 70 deg to the plane of the sky.
To chart the fine three-dimensional distribution of Orion A we can use dust, which is coupled to the gas component of the interstellar medium (ISM) and hence traces the densest gas of the molecular cloud in which it sits.The dust is accessible to us by its extincting effect on light passing through it.But the task of charting the dust of Orion A is very different from the task of charting its stellar population.It is a matter of inference rather than direct observation.Given extinctions for stars in the region of Orion A we can infer the dust density at an arbitrary point in the cloud, a process known as 'three-dimensional dust mapping'.A number of teams have already made three-dimensional dust maps of Orion A. These include: Schlafly et al. (2015), using only photometry from Pan-STARRS1; Rezaei Kh et al. (2018), using astrometry from Gaia Data Release 1 together with photometry from 2MASS and WISE, Rezaei Kh et al. (2020), using astrometry from Gaia DR2 together with photometry from 2MASS and WISE;and Dharmawardena et al. (2022), using astrometry from Gaia DR2 together with photometry from 2MASS, WISE, and Gaia DR2.All three identify a filament of gas consistent with the tail of young stellar objects and young stars.Rezaei Kh et al. (2020) report a cloud lying in the foreground of Orion A, at a distance of s = 350 pc, which they tentatively associate with a stellar cluster identified by Bouy et al. (2014).However, this cloud is not reported by either of the other two groups.
Although the distributions of the gas and dust are of interest in their own right, the dust is also something of a nuisance since it obscures and reddens observations of stars within and beyond it.Because the Sun is located close to the Galactic mid-plane we observe our galaxy through significant quantities of dust.Moreover, the vast majority of the Galaxy's stars also lie in the disc, so extinction has an enormous impact on stellar surveys and, in turn, on attempts to fit chemical and dynamical models of the Galaxy to stellar observations.To fully exploit the data from stellar surveys we require good three-dimensional dust maps of the whole Galaxy.A number of teams have made such three-dimensional dust maps.These include: Arenou et al. (1992); Drimmel et al. (2003); Marshall et al. (2006); Sale et al. (2014) (based on the method of Sale 2012); Chen et al. (2014) (based on the method of Berry et al. 2012); Lallement (2015), Lallement et al. (2014), Lallement et al. (2018), Lallement et al. (2019), and Vergely et al. (2010) (based on the method of Vergely et al. 2001); Green et al. (2015), Green et al. (2018), and Green et al. (2019) (based on the method of Green et al. 2014); and Leike & Enßlin (2019) and Leike et al. (2020) (based on the method of Enßlin & Weig (2010) and Enßlin & Frommert 2011).
To date, however, dust maps have suffered from either artificial inhomogenities known as 'fingers of God' (Marshall et al. 2006;Green et al. 2018;Lallement et al. 2019), or have been resolution-limited (Sale et al. 2014;Rezaei Kh et al. 2020;Leike et al. 2020).We propose a new method for creating three-dimensional dust maps and use it to map Orion A. Orion A is a good testing ground for dustmapping methods that might be used to map the whole Galaxy.
When making a three-dimensional dust map we wish to infer the density, ρ(r) or, equivalently, the extinction, A(r), at an arbitrary point, r, given knowledge of the positions, r 1 , . . ., r n , and extinctions, A 1 , . . ., A n for some set of stars.In this paper we treat the ISM as the realization of a random field, and compute the best linear unbiased predictor of the extinction at a point given extinctions for a set of stars in the region of that point.This is equivalent to forming the generalized least squares estimator of the mean of the extinction and then performing a nonparameteric fit to the residuals.We then differentiate this predictor to give a predictor of the density at the same point.Our method requires minimal assumptions about the statistical properties of the ISM.Indeed it requires us to know only the covariance of the densities of the ISM and to have a model of its mean.It does not require us to assume a distribution for the dust density, and in particular it does not require us to assume a log-normal distribution, as is common in three-dimensional dust mapping.We follow Sale & Magorrian (2014) in using a physically motivated model of the covariance function that captures the turbulent structure of the ISM.The parameters of this covariance model form the hyperparameters of our predictor of extinction, which we optimize using the method of leave-one out cross validation.We make two maps of dust in the region of Orion A out to 500 pc: one using extinctions computed using the Rayleigh-Jeans colour excess method, and one using extinctions from the StarHorse catalague.The resulting maps display no fingers of God, and may be constructed at arbitary resolution although the effective resolution is always limited by the density of observed stars.
Our method relies on knowledge of only the broadest statistical properties of the density field (i.e. its first and second moments) and we begin, in Section 2, by summarizing these.We then introduce the best linear unbiased predictor, in Section 3, testing it against simulated data in Section 4 before using it to map Orion A, in Section 5. Our maps broadly agree with those made by Rezaei Kh et al. (2020) and, like them, show a foreground dust cloud at a distance of 350 pc.In Section 6 we discuss the consequences of our assumptions about the statistical properties of the density field and show that our method is insensitive to them so long as they correctly account for the fractal nature of the ISM.

STATISTICS OF THE INTERSTELLAR MEDIUM
The density of the ISM at a point, r, may be represented by the random variable, ρ(r).The set of all such random variables forms a random field, ρ := (ρ(r)) r∈R 3 , which represents the density of the ISM as a whole.Similarly, the extinction of light emitted by a source at a point, r, may be represented by the random variable where s := |r| is the line-of-sight distance of the source and r is the unit vector parallel to r.This gives another random field, A := (A(r)) r∈R 3 , which represents the extinction that would be undergone by light emitted at any point in space.The properties of these fields are usefully summarized by their means and variances.Since the expection and covariance functions are both linear we immediately find that the expectation of A(r) is and that the covariance of where s i := |r i | and s j := |r j | are the line-of-sight distances of the two sources and ri and rj are the unit vectors parallel to r i and r j .
Although the ISM is a complicated multiphase medium, the fluctuation in its density, can be approximated as being stationary and isotropic (Draine 2011).In these circumstances, its first and second moments are invariant under translations and rotations.Consequently, the covariance of the density fluctuations at two points, ∆ρ(r i ) and ∆ρ(r j ), is a function of their separation, |r j − r i |, and is given by the autocovariance function of the density fluctuations, k ∆ρ , such that By Bochner's Theorem (Adler 1981) the Fourier transform of k ∆ρ exists, and is related to the power spectrum of ∆ρ, denoted P ∆ρ , via the equation where k is the three-dimensional wavenumber.
It is common to assume that fluctuations in the density are described by Kolmogorov's theory of turbulence (Draine 2011), i.e. that its power spectrum is described by a power law, with index −11/3, over a wide range of wavenumbers.2This power law must fail at very small and very large wavenumbers.In the Kolmogorov theory of turbulence this happens at the boundaries of the inertial regime, i.e. at scales larger than the energy-injection scale, L 0 , or smaller than the dissipation scale.The energy-injection scale is that of stellar feedback, while the energy-dissipation scale is that of viscous forces.Although astronomical observations clearly probe the energy-injection scale, they do not probe the energy-dissipation scale, meaning that we need only worry about the failure of the power-law power spectrum for small wavenumbers.Motivated by this fact, Sale & Magorrian (2014) proposed a split-power law for the power spectrum of ∆ρ, with the location of the break being determined by the energy-injection scale, L 0 .It is given by for positive numbers Ω, γ, and R, a normalizing constant defined implicitly by the equation Above the wavelength associated with the energy-injection scale, L 0 , the spectrum is a power-law with index −γ, and below that wavelength it is a power law with index Ω.Sale & Magorrian call this power spectrum 'Kolmogorov-like', regardless of the value of γ.
The mean value of the density is in general unknown, but its covariance is equal to the covariance of the density fluctuations, since The density, ρ, is often taken to be lognormal, meaning that the extinction at a point, A(r), is the integral of a lognormal field.This has a probability density function (PDF) that cannot be expressed in closed form, but which is itself approximately lognormal (Ostriker et al. 2001).
If we assume that the density (as well as its fluctuations)3 is stationary and isotropic, with mean µ ρ and constant variance σ 2 ρ then the covariance of the densities at two points, ρ(r i ) and ρ(r j ), is also a function of their separation, |r j − r i |, and is given by the autocovariance function of the density, k ρ , such that The autocovariance of the density is then equal to the autocovariance of the density fluctuations, i.e. k ρ = k ∆ρ , and we have that where k ∆ρ is given by Equations 8 and 9.For each choice of autocovariance function, we need only a model of the mean of the extinction, which has a single free parameter, µ ρ .Despite this simplifying assumption, neither cov(A(r i ), A(r j )) nor cov(ρ(r i ), A(r j )) has closed form, and both must be computed numerically at some expense.For this reason we define the functions f and g such that where s 1 and s 2 are the line-of-sight distances to points r 1 and r 2 , and θ 12 is the angle subtended by them at the origin, which we precompute on a regular lattice of arguments and evaluate using linear interpolation.(We discuss these functions further in App.D).

LINEAR PREDICTION AND MAP MAKING
When making a three-dimensional dust map we wish to infer the density or, equivalently, the extinction at an arbitrary point given knowledge of the extinctions and positions for some set of stars.Of course, neither the distance nor the extinction are observed directly.Instead, they must themselves be inferred from other stellar properties.For a limited number of stars (those for which we have observations in the infrared) we may estimate the extinction using the near-infrared colour excess (NICE) method of Lada et al. (1994) or its successors (the NICER method due to Lombardi & Alves 2001, the NICEST method due to Lombardi 2009, and the Rayleigh-Jeans colour excess method due to Majewski et al. 2011).The Rayleigh-Jeans colour excess (RJCE) method, for example, exploits the fact that all stars are well approximated as black bodies in the lowfrequency limit (where their spectral radiance is described by the Rayleigh-Jeans law), meaning that their low-frequency colours are approximately independent of spectral class.Majewski et al. (2011) showed that colours spanning the NIR and MIR are approximately constant for a range of spectral types, and that the colour H − [4.5 µ] exhibits especially little variation.They found, using the extinction law of Cardelli et al. (1989), that the K s -band extinction in magnitudes is with a scatter of 0.1 mag for F, G, and K stars, and 0.4 mag for B to M excluding late type dwarfs.However, due to observational error, both distances and extinctions estimated in this way may be unphysically negative.However, it is often necessary to infer the extinction and distance of a star using Bayesian methods of the kind pioneered by Burnett & Binney (2010), Bailer-Jones (2011), and Binney et al. (2014).These methods assume that a star's observed properties are parameterized by its instrinsic properties, and use Bayes's theorem to compute the posterior PDF of the intrinsic properties given the observed properties.We might wish to infer a single property or a tuple of properties, in which case we can find the marginal PDF of any element of that tuple.Several catalogues of stellar properties have been compiled in this way.They include the StarHorse catalogue (Santiago et al. 2016, Queiroz et al. 2018, Anders et al. 2019, Anders et al. 2021), and the catalogues of Bailer-Jones and his collaborators (Bailer-Jones 2015, Astraatmadja & Bailer-Jones 2016a,b, and Bailer-Jones et al. 2018, 2021).
The StarHorse catalogue takes the intrinsic properties of a star to be mass, age, distance, and V-band extinction.It takes the observable properties of a star to be effective temperature, surface gravity, overall metallicity, magnitude, and parallax, as catalogued by Gaia, Pan-STARRS1, 2MASS, and AllWISE.The likelihood of these observable properties is assumed to be Gaussian with mean given by a stel-lar model, and variance given by the observational errors.It assumes as a prior PDF for distance a physically motivated four-component model of the Galaxy (comprising the thin disc, thick disc, bulge/bar and halo).As a result, inferred distances are always positive.In the first release of the StarHorse catalogue, using Gaia DR2 (Anders et al. 2019), the prior PDF for extinctions is assumed to be uniform in the case of stars with low-noise parallax observations (signal-tonoise of five or greater) and a top hat, such that A V /mag ∈ [−0.3, 4], in the case of stars with high-noise parallax observations (signal-tonoise of less than five).In the second release, using Gaia Early Data Release 3 (EDR3) (Anders et al. 2021), the prior PDF for extinctions is assumed to be given by the maps of Green et al. (2019) or Drimmel et al. (2003) according to coverage.The full posterior PDF is not publicly available, but rather the 0.05, 0.16, 0.6, 0.84, and 0.95 quantiles.
The catalogues of Bailer-Jones provide only distances.They take the instrinsic property of a star to be distance, and its observable property to be parallax.The most recent of these catalogues (Bailer-Jones et al. 2018, 2021) assumes that the likelihood of the parallax is Gaussian with mean given by the reciprocal distance and variance given by the observational errors, and that the prior PDF for distance is that for a generalized gamma distribution.As a result, inferred distances are again always positive.Again, the full posterior PDF is not publicly available, but rather the 0.05, 0.16, 0.6, 0.84, and 0.95 quantiles.

The best linear unbiased predictor
For a list of directly computed or inferred extinctions A = (A 1 , . . ., A n ) at locations r 1 , . . ., r n we wish to predict the extinction, A(r), or density, ρ(r), at some arbitrary position, r, given the assumptions represented by Eqs 14 and 15.To predict the extinction we will use the best linear unbiased predictor (BLUP) of A(r), which we will denote Â(r).This is a linear combination of extinctions A 1 , . . ., A n with coefficients chosen so as to minimize the meansquare error, (where angle brackets again indicate the expectation) subject to the constraint that Â(r) be unbiased, i.e. subject to the constraint that ⟨ Â(r)⟩ = ⟨A(r)⟩.It is given by Goldberger (1962) as where the n × 1 vector-valued function γ is given by The expression for Â(r) is the sum of two terms.The first of these, γ(r) t A, is the generalized least squares (GLS) estimator of the mean of A(r), which we may rewrite as μρ s, where is itself the GLS estimator of the mean of the density, µ ρ , and is unbiased.Similarly, the expression Γ t A gives the GLS estimators of the means of A 1 , . . ., A n .The second term in the expression for Â(r) is then a weighted sum of the residuals of the observed values A and the GLS of their means. 4o predict the density we will use the derivative of the BLUP of A(r), which is also unbiased, since ⟨ρ(r) − ρ(r)⟩ = ⟨∂(A(r) − Â(r))/∂s⟩ = ∂ ⟨A(r) − Â(r)⟩ /∂s = 0.Because only γ and σ are functions of s it is easy to compute (see App.A for details).
Note that to find Â(r) and ρ(r) we need only to know the covariance of (A(r), A) and to assume that the expectation of A(r) is linear in s.We do not need to make any further assumptions about the distribution of the ISM.In particular, we do not need to assume that its distribution is lognormal.5

Computing the BLUP in practice
When using extinctions inferred with the RJCE method we can make the decomposition A i = A(r i ) + E i where A(r i ) is the intrinsic extinction at r i and E i is the uncertainty in its value.In this case the covariance matrix is assuming that A(r i ) and E j are independent for i j, and assuming that A(r) and E j are independent.When using extinctions inferred by Bayesian techniques that make stronger assumptions about he stars' neighbourhoods A i is the posterior prediction of a source's extinction, and both cov(A i , A j ) and cov(A(r), A i ) are properties of the Bayesian model used to construct it.Since neither is available as part of a catalogue we will approximate them by Equations 25 and 27, on the understanding that var(E i )δ i j no longer represents an observational error but instead quantifies the discrepancy between our approximation of the covariance and its true value.We will take var(E i ) to be equal to the variance of var(A i ).
To compute Â(r) we must know the position of each star.However, as we have already noted, the distance to a star is inherently uncertain.To circumvent this problem we may assume that the distance to a star is certain but that the uncertainty in the extinction is increased, using the fact that In this case, the variance of A i is increased by dA 2 i .Under this assumption the i j-th element covariance matrix, [Σ] i j , becomes [Σ] i j + dA 2 i δ i j , i.e. we add the square of the uncertainties, dA 1 , . . ., dA n , to the diagonal elements of Σ.If d s i ≪ L 0 then we may make the linear approximation, Of course, if the condition d s i ≪ L 0 is not met then this additional term in the variance will be too great.

Validation
Although our statistical model of the ISM is physically motivated, we would still like reassurance that Â(r) and ρ(r) are good predictors of A(r) and ρ(r).In the case of extinction we may test the performance of the BLUP using leave-one-out crossvalidation (LOOCV).We partition A into a set containing a single element, A i , and another containing n − 1 elements, We then find Â(r i ), namely the BLUP of A(r i ) based on A −i , and compute the residual, A i − Â(r i ).The behaviour of these n residuals should be consistent with the assumptions we have made about their distribution.Since we have only assumed that the mean is linear and that we know the autocovariance function, we need only check that the residuals obey Chebyshev's inequality, where std(A(r) − Â(r)) is the standard deviation of the residuals and λ is a positive number, although we would hope for them to do considerably better.They should be uniformly good across the mapped region, and display no trend in distance, latitude or longitude.However, we always expect our predictors to underperform near the boundaries of the mapped region, where they are constrained by fewer data.It is also useful to define the LOOCV score, which should be small compared to, say, the mean variance of the elements of A. Since we are unable to validate the BLUP of the density of the ISM in this way, we will take validation of Â(r) to be validation of ρ(r).

Choosing the model parameters
We have observed that for a given autocovariance function, our model of the mean of the extinction has a single free parameter, µ ρ .However, we do not in practice know the values of the parameter tuple that specifies the power spectrum, and hence the autocovariance function.We will assume that γ = 11/3 (to ensure Kolmogorov turbulence in the inertial regime), and that Ω = 1 (arbitrarily).However the variance of the density, σ 2 ∆ρ , and the length scale, L 0 , are less certain.We might wish to find their maximum-likelihood estimate (MLE).But this would require us to assume a model for the distribution of the extinction, contrary to the assumptions of the BLUP, so that by doing this we would lose one of the principal benefits of the method.Moreover, in computing the BLUP of the extinction we are necessarily computing the GLS estimator of the mean, and we do not want a second estimate of it.Instead, we may follow a procedure much used in machine learning, and choose the pair σ 2 ∆ρ and L 0 so as to minimize the LOOCV score, R LOO (eq.31).This does not allow us to formally estimate the parameter tuple in the way the method of maximum likelihood does.We should, instead, see it as a way of tuning the predictor so that it passes validation.Ultimately, our trust in the BLUP is justified only by its performance in validation, and we may choose the power spectrum's parameter tuple arbitrarily, so long as the resulting predictors pass this test.

Computational expense
The memory requirements for computing the BLUP are dominated by the storage of square matrices of size n × n and are therefore of order O(n 2 ).The computational complexity is dominated by the inversion of the covariance matrix, Σ.For small data sets, like those needed to map the dust in Orion A, we may do this using Cholesky decomposition, which has computational complexity of order O(n 3 ).However, for significantly larger data sets, like those needed to create global dust maps, this complexity is prohibitive.Such scaling problems are common to all three-dimensional dust mapping methods and a number of solutions are available.In our case it is best to use a low-rank or sparse approximation to the covariance matrix or to use the iterative inversion method described by Wang et al. (2019), which has complexity of order O(n 2 ).We will pursue this in future work.

METHODOLOGICAL TESTS USING SYNTHETIC DATA
We will shortly use our method to map dust within Orion A and its environment.In fact, we will follow Rezaei Kh et al. (2020)  , which contains approximately 10 000 stars.We will do this using two data sets: first, extinctions computed using the RJCE method; and second, extinctions reported by the StarHorse catalogue.But before doing so we will test our method using synthetic data that mimic those generated by the RJCE method.
To generate our synthetic data we realize the density field, ρ, on a 501 × 101 × 45 regular lattice in s, l, and b, under the assumption that ρ is a stationary and isotropic lognormal random process having a Kolmogorov-like power spectrum (Eq.9) with mean µ ρ = 1 × 10 −3 mag pc −1 , variance σ 2 ∆ρ = 1 × 10 −6 mag 2 pc −2 , energy-injection scale L 0 = 100 pc, and power-law indices Ω = 1 and γ = 11/3. 6 To do this, we note that the expectation and covariance of ln(ρ(r)) are given by the well-known formulae (see, for example, Coles & Barrow 1987) and generate a realization of a Gaussian random process, ln(ρ), which we then exponentiate, to give ρ.We then generate a realization of A on the same 501 × 101 × 51 lattice by approximating the line-ofsight integral using the trapezium rule, and choose a uniformly distributed random sample of this realization, size n = 10 000, to which we add normally distributed noise with zero mean and standard deviation of 0.1 mag.This forms a set of of synthetic observations, A = (A 1 , . . ., A n ) at locations r 1 , . . ., r n where we work in Galactic coordinates such that r i = (s i , l i , b i ).To each distance we add normally distributed noise with zero mean and standard deviation of 5 per cent.We then propagate errors from distance to extinction under the assumption that the linear approximation holds (Sec.3.1.1).Since the mean and maximum errors on distances are 19 pc and 25 pc respectively this approximation holds marginally.
We now map the entire area by evaluating Â and ρ for all points on the lattice.To do this, we assume that Ω = 1 and γ = 11/3, but optimize our choice of the remaining autocovariance parameters, σ 2 ∆ρ and L 0 , by minimizing the LOOCV score, R LOO (Sec.3.3).This is expensive, since it involves the inversion of the matrix Σ for every trial parameter pair.We therefore sample R LOO on a 25×25 logarithmically spaced rectangular lattice covering the region σ 2 ∆ρ /mag 2 pc −2 ∈ [10 −8 , 10 −4 ] and L 0 /pc ∈ [10,250].Doing this, we find that σ 2 ∆ρ = 3.0 × 10 −6 mag 2 pc −2 and L 0 = 98 pc (Fig. 1).These are close to the true values but we may not say they are consistent with the true values in any formal sense since we have no confidence region for them.We then perform LOOCV (Section 3.2) by inspecting the standardized leave-one-out residuals (Fig. 2).These have approximately unit standard deviation, and exhibit no trends in s, l, or b.Accordingly, we say that our predictors pass validation.Immediately, we find that the generalized least squares estimate of the density is μρ = 1.2 × 10 −3 mag pc −1 with root mean-square error (RMSE) 3.36 × 10 −4 mag pc −1 (App.B), consistent with the true value of 1 × 10 −3 mag pc −1 .
In Figures 3-6 we show our predictions for the plane b = 3.5 deg, and in Figure 7 our predictions for the line of sight l = 7 deg, b = 3.5 deg.Note that we cannot expect the BLUP to reproduce features on scales smaller than the mean stellar separation.However, it is large fluctuations in the density on this interstellar scale that result in the largest increases in extinction.These largest density fluctuations cannot be probed, and as a result the density map is smoothed out.Nonetheless, the predictions for extinction and density are everywhere consistent with their true values.Note that Â is increasing in s, even though we do not constrain it to be.Similarly, note that ρ is non-negative for all r.These facts are most clearly seen in Figure 7.The prediction interval for extinction increases with distance and is inhomogeneous, just as the extinction is.The prediction interval for density does not increase with distance and is homogeneous, just as the density itself is.Note that predictions are even good close to the boundary, and in particular within the region for which s < 200 pc, where the predictor is most poorly constrained.

MAPPING DUST IN ORION A
Having tested our map-making method in this way we can now use it to create maps using each of our two data sets.extinction is 0.17 mag, and the standard deviation of the error in extinction is 0.073 mag.
We combine these extinction data with the line-of-sight distances inferred by Bailer-Jones et al. ( 2021) using Gaia EDR3 (Gaia Collaboration 2018), again using the Gaia Archive cross-matching algorithm.We take for our distance prediction the median of the posterior distribution (which we expect to be close to the mean of the predicted distance since the probability density functions are unimodal and approximately symmetric at such small distances), and for our error half the difference of the 0.16 and 0.84 quantiles.The mean error in distance is 6.3 pc, and the standard deviation of the error in distance is 11 pc.We propagate errors from distance to extinction under the assumption that the linear approximation holds (Sec.3.1.1).The final sample contains 10 180 stars, having an average stellar separation of 3.9 pc.Of these sources, 1 434 have negative measured extinctions.We show all extinctions in Fig. 8.
We again map the entire region by evaluating Â and ρ at all points on a 501 × 501 × 45 regular lattice in s, l, and b, repeating the analysis we performed using synthetic data.By minimizing the LOOCV score, R LOO , we find that that σ 2 ∆ρ = 6.2 × 10 −5 mag 2 pc −2 and L 0 = 98 pc (Fig. 9).(Note that the precision of these values is limited by the size of the grid we use to search parameter space.)The leave-one-out residuals have better than unit standard deviation, and again exhibit no trends in s, l, or b, so that we may again say that our predictors pass validation (Fig. C1).The GLS estimate of the density is μρ = 3.8 × 10 −4 mag pc −1 with RMSE 2.75 × 10 −3 mag pc −1 (see App. B for a discussion of our quoted errors).In Figures 10 and 11 we show our predictions for the plane b = −19.5 deg (corresponding to Fig. 4 in the paper by Großschedl et al. 2019), and in Figure 12 we show our predictions for a sequence of on-sky regions at increasing line-of-sight distances (corresponding to Figs 2 and 5 in the 2020 paper by Rezaei Kh et al.).Note that Â is not increasing in s, and that ρ is negative for some r.Nevertheless, prediction intervals for Â are consistent with its being increasing in s, and prediction intervals for ρ are consistent with its being non-negative for all r.
Our maps clearly identify the head of Orion A (at s = 390 pc, l = −19.5 deg, b = 209 deg) and its tail, which extends to a distance of s = 470 pc albeit varying in density along its length.They also clearly identify the foreground cloud first noted by Rezaei Kh et al. (2020).In the plane b = −19.5 deg (Fig. 11) this appears as a chevron, the apex of which is located at s = 350 pc, l = 210 deg.In the plane s = 345 pc (Fig. 12) it is seen to have two lobes, one centred at l = −19.5 deg, b = 211 deg, and the other centred at l = 206.5 deg, b = 17 deg.The lower of these is located at the apex of the chevron, and is bifurcated at smaller distances.Figure 13 shows our predictions for the lines of sight centred on these lobes (corresponding to Fig. 3 in the 2020 paper by Rezaei Kh et al.).
Let us consider possible systematic errors in our predictions.These might arise from inadequate models of the mean and covariance of the extinction field (Eqs 14 and 15) or from systematic errors in the distance and extinction data that we use.Such errors might be associated with stars of a particular spectral class, causing them to appear too close or too distant, too red or too blue.The distances to our sources have been inferred using only parallax information so we do not expect there to be any systematic errors of this type in distance.However, the RJCE method for computing extinctions relies on a linear fit to observed colours.This assumption of linearity is good for a large range of spectral classes but could result in the extinctions of very cool stars being underestimated and the extinctions of very hot stars being be overestimated.Sytematic errors associated with spectral class should appear as overdensities in the joint distribution of (i) spectral class and distance or (ii) spec- tral class and extinction, and such overdensities should alert us to the presence of systematic errors in the computed values of either quantity.We can use the dereddened colour (G − K s ) 0 (computed using using the exinction law of Cardelli et al. 1989) as a proxy for spectral class and inspect the histogram for each distribution (not shown).In neither case do we see overdensities associated with any particular value of colour.We therefore claim that our data are free of significant systematic errors relating to spectral class.

Application to StarHorse data
The StarHorse catalogue includes quantiles for the inferred probability density functions of a source's extinction and distance, and therefore neatly provides all the information we need to construct a dust map.We use the first StarHorse catalogue (Anders et al. 2019) rather than the second (Anders et al. 2021) since this second catalogue begs the question by assuming for a prior on extinction the maps of Green et al. and Drimmel et al.Each source in the StarHorse catalogue is associated with the compound flags SH_GAIAFLAG, which quantifies the quality of the Gaia astrometry and photometry used, and SH_OUTFLAG, which quantifies the quality of StarHorse's inferences (Anders et al. 2019).The SH_GAIAFLAG flag consists of three digits flagging (1) high renormalized unit weight error, warning of a spurious astrometric solution, (2) high colour excess factor, warning of spurious photometry, and (3) source variability.The SH_OUTFLAG flag consists of five digits flagging (1) overall unreliability, warning of physically inconsistent inferred stellar properties, (2) large distance, warning of a spuriously large inferred distance, (3) unreliable extinction, warning of a spuriously small or spuriously large inferred extinction, (4) large extinction uncertainty, warning of an inferred extinction uncertainty of unit magnitude or greater, and (5) small uncertainty, warning of spuriously small uncertainties in any inferred stellar property.We use sources for which SH_GAIAFLAG is 000 and SH_OUTFLAG is 00000, ensuring that potentially spurious statistics of each kind are not included in our sample.
For each source in this sample, we realize a sample of the posterior extinction by assuming its PDF to be uniform within quantiles.In order to do this, we assume a lower limit on A V of −0.3 mag, and an upper limit of 4 mag, consistent with the StarHorse prior PDF for extinctions with errors in distance greater that 20 per cent.We take the standard deviation of the posterior extinction to be half the dif-ference of the 0.16 and 0.84 quantiles.The mean error in extinction is then 0.17 mag, and the standard deviation of the error in extinction is 0.011 mag.
We again take for our distance prediction the median of the posterior distribution, and for our error half the difference of the 0.16 and 0.84 quantiles.The mean error in distance is then 13 pc, and the standard deviation of the error in distance is 9.7 pc.Again, we propagate errors from distance to extinction under the assumption that the linear approximation holds (Sec.3.1.1).The final sample contains 12 719 stars, having an average stellar separation of approximately 3.7 pc.Of these sources, 1 434 have negative extinctions.
By minimizing the LOOCV, R LOO , we find that that σ 2 ∆ρ = 2.2 × 10 −5 mag 2 pc −2 , and L 0 = 128 pc (Fig. 14).However, the leaveone-out residuals have very high standard deviation of 4.5 mag 2 , although they exhibit no trends in s, l, or b (Fig. C2).As such, they fail validation.This may be due to the presence of young stellar objects in the StarHorse catalogue, which we have made no attempt to identify and remove, and which violate the assumptions of our model.Nevertheless, we find that GLS estimate of the density is μρ = 6.2 × 10 −4 mag pc −1 with RMSE 1.94 × 10 −3 mag pc −1 and, as before, we plot our predictions for the plane b = −19.5 deg (Figs 15  and 16), a sequence of on-sky regions at increasing line-of-sight distances (Fig. 17 18).Note that StarHorse extinctions are given for the V-band rather than K s -band, and that according to the Cardelli et al. (1989) extinction law, A V /A Ks = 8.8.Again, we note that Â is not increasing in s, and ρ is negative for some r.Nevertheless, prediction intervals for Â are again consistent with its being increasing in s, and prediction intervals for ρ are again consistent with its being non-negative for all r.
The head of Orion A (at s = 390 pc, l = 209 deg, b = −19.5 deg) appears only faintly in our maps, although the tail and foreground cloud both appear clearly.In the plane b = −19.5 deg the foreground cloud is seen as a linear structure rather than chevron, though it is again seen to have two lobes in the plane s = 345 pc.Since these maps do not pass validation we do not use them to make comment on the distribution of dust in Orion A. In particular we do not use them to comment on the existence of a foreground cloud.
Here we have used a single realization of the StarHorse posterior extinctions.Alternative realizations of the posterior extinctions produce similar results.We could instead have made multiple realizations and generated an empirical distribution for the predictor Â(r).In practice, however, we find that the mean and standard deviation of the empirical distribution converge quickly, and that they are in fact well approximated by a single sample.

DISCUSSION
To make our maps we have used a physically motivated Kolmogorov-like autocovariance function to describe the density fluctuations of the ISM.When using the BLUP (Eq.20) it is common to use one of many off-the-shelf library functions that do not have direct physical motivation.For example, it is very common to use a powered-exponential autocovariance function, k PE , given by where σ 2 PE is the variance of the field, L PE is a characteristic length scale, and where α ∈ (0, 2].Particularly popular is the case of α = 2, when k PE is known as the 'squared-exponential autocovariance function'.Autocovariance functions with compact support are also at- and k G (δ) = 0 otherwise, where σ 2 G is the variance of the field and L G is a characteristic length scale for the field, above which the autocovariance vanishes.
Any such standard autocovariance function may be used with our method by integrating it to give an approximation to cov(A(r i ), A(r j )) according to Equation 5.In turn, we can approximate the functions f and g (Eqs 16 and 17).(See App.B for plots of these approximations.)Although library functions of this kind have no direct physical motivation, we do require them to have properties consistent with the physics of the ISM.In particular we require the random processes they define to be continuous, undifferentiable and to have physically meaningful characteristic length scales.We should be mindful of these properties when we specify the autocovariance function and when we interpret our predictions.A stationary and isotropic random field is continuous (in the mean square) if its autocovariance function is continuous at zero, and is differ- The leave-one-out cross-validation score, R LOO , is found to have a minimum for σ 2 ∆ρ = 2.2×10 −5 mag 2 pc −2 and L 0 = 128 pc.The covariance matrix, Σ, is nearly singular for some parameter combinations (region shown in white), meaning that R LOO cannot be computed.Contours are shown at intervals of 1 × 10 −3 mag 2 .(See Sec. 5 for discussion.)entiable (in the mean square) if its autocovariance function is also twice differentiable at zero (i.e. a stationary random field is m-times differentiable everywhere if its autocovariance function is 2m-times differentiable at zero).The Kolmogorov-like autocovariance function is undifferentiable at zero, and hence defines an undifferentiable random field.This roughness is clear in our synthetic dust maps (bottom panel of Fig. 7).Gneiting's function is also undifferentiable at zero, and therefore also defines an undifferentiable random field.The powered-exponential is differentiable only for α = 2, when it is infinitely differentiable, and hence defines an infinitely differentiable random field.For α 2, however, it is undifferentiable, and hence again defines an undifferentiable random field.Indeed, realizations of the random field increase in roughness as α decreases.Whereas Gneiting's function is suitable for specifying the dust's autocovariance, the squared-exponential kernel is not.Using it will result in maps that lack sharpness and exhibit features that are ill-defined.
To find a physically meaningful interpretation of the characteristic length scale we may use the integral length scale of a random field, which is much used in the study of turbulence (Tennekes 1972).Suppose that k is a autocovariance function with variance σ 2 and characteristic length scale L. The integral length scale of the stationary random field associated with k is8 In order to make explicit its dependence on σ 2 and L let us write k as k(•; σ 2 , L).We introduce the physical length scale, Λ, such that noting that, whereas k is the autocovariance function of interest, k ∆ρ is the autocovariance of the density fluctuations (eq.8) with power spectrum defined by Sale & Magorrian (2014) and given in equation 9.This may be rewritten as or, equivalently, which allows us to find the physical length scale associated with a given autocovariance function.It approximates the energyinjection scale of the equivalent Kolmogorov-like density field.The Kolmogorov-like, squared-exponential, and Gneiting functions have physical length scales L 0 , 2.51L PE , and 0.466L G .
We can illustrate the consequences of our choice of autocovariance function by making maps using the squared-exponential and Gneiting functions with the data of Rezaei Kh. & Kainulainen (2022) and Bailer-Jones et al. (2021) (Sec.5).For the squaredexponential function we find no minimum in the LOOCV score: for large regions of parameter space inversion of Σ is numerically unstable (and hence the LOOCV may not be computed), whereas elsewhere the LOOCV score reduces monotonically as L PE approaches 10 pc (the lower limit of the search region) for all values of σ 2 PE (Fig. 19, top panel).This reflects the fact that the autocovariance function is misspecified, and that the LOOCV score of the BLUP can only be reduced by its overfitting the data.We may nonetheless choose the parameters of our function arbitrarily.Suppose, then, that we choose σ 2 PE = 1 × 10 −5 mag 2 and L PE = 40 pc, such that L PE is the characteristic length we would expect if Λ = 100 pc.In this case our predictors do pass validation, and we are able to make the maps shown in Figures 20 and 21.Note that the features in our map (in particular, the head and tail of Orion A, as well as the foreground cloud) are no longer resolved and that the map itself is much less sharp than before, reflecting the fact that the squared-exponential function results in smooth random fields that do not have the rough properties of the ISM, and that because of this its length scale cannot be given a physical interpretation.Moreover, our results are sensitive to our choice of length scale, meaning that a scaling of L PE will result in a scaling of the predicted features.With a characteristic length scale of L PE = 15 pc our predictors again pass validation, and the features in our map are resolved, but still lack sharpness.(We do not show this effect in our figures.)For Gneiting's function we find two minima in the LOOCV score, at σ 2 G = 1.78×10 −5 mag 2 , L G = 13.1 pc, and σ 2 G = 5.62×10 −5 mag 2 , L G = 191 pc (Fig. 19, bottom panel).Both pairs of parameters result in predictors that pass validation.However, it is the second pair that corresponds to the physical length scale we expect, since in this case we find that Λ = 89 pc.Maps made using this second pair of parameters are shown in Figures 22 and 23.They are all but indistinguishable from maps made using the first pair of parameters, indicating that our predictors are somewhat robust against misspecification of the length scale.Moreover, they are very similar to the maps made using the Kolmogorov-like autocovariance function (Figs 10 and 11), indicating that Gneiting's function provides a good approximation to the Kolmogorov-like function.

CONCLUSION
We have made a three-dimensional dust map of the giant molecular cloud Orion A by computing the best linear unbiased predictor (BLUP, Eq. 20) of the extinction at every point on an arbitrarily dense lattice using extinctions for a set of observed stars in the region of this lattice.The BLUP requires us to make very few assumptions about the statistical properties of the ISM.It requires only that we know the covariance of the density at any two points, and that we have a model of the expected value of the density at any one point.Beyond this knowledge of the first and second moments of the density, we need know nothing of its distribution.
In practice, however, we do not know the covariance of the extinctions at two points, and instead have a model of it, just as we have a model for the expected value.Since we do not know the distribution of the dust density we are not able to use the method of maximum likelihood to optimize this model, but instead fit it to the data by choosing the parameter tuple that minimizes the LOOCV score.It is also convenient to assume that the density field is stationary and isotropic, and that it is described by Kolmogorov's theory of turbulence within the inertial range set by the energy-injection and eneryg-dissipation scales.Moreover, by assuming that the density is isotropic and homogeneous our model of the mean density has only one free parameter, µ ρ , which we can recover in the course of making our predictions.
This agnosticism is both a strength and a weakness of our method, since it does not allow us to encode the fact that density is al-  . 35).The leave-oneout cross-validation score, R LOO , is found to have a minumum for σ 2 ∆ρ = 1.78 × 10 −5 mag 2 , L 0 = 13.1 pc and σ 2 ∆ρ = 5.62 × 10 −5 mag 2 , L 0 = 191 pc.In both cases, the covariance matrix, Σ, is nearly singular for some parameter combinations (regions shown in white), meaning that R LOO cannot be computed.(Cp.Fig. 9.) Contours are shown at intervals of 5 × 10 −5 mag 2 .(See Sec.6 for discussion.)ways non-negative.Nonetheless, in experiments using synthetic data we find that predicted densities are, in practice, everywhere nonnegative.With observational data we find that predicted densities are, in practice, everywhere consistent with being non-negative.Negative predicted densities are the result of spuriously high predicted extinctions regressing to the mean along the line of sight, meaning that, along a line of sight, the extinction field is not monotone.Observational data also suffer from selection effects that cause the average observed extinction to decrease with distance (Sale 2015).The more heavily extinguished a star is the fainter it is and hence the less chance it has of appearing in any magnitude-limited catalogue.These negative densities are therefore due to either a lack of data or to model mispecification.In the latter case, either our assumptions about the statistics of the ISM fail to do it justice, or the data include sources with extinctions that are not due exclusively to the dust of the ISM (for example, the extinctions of young stellar objects is partially due to their dusty circumstellar discs and envelopes).Such sources are likely to be present in any catalogue, especially in catalogues containing nebulae like those of Orion A, even after cleaning.
The choice of autocovariance function is crucial to the quality of any three-dimensional dust map.We have used the physically motivated Kolmogorov-like function defined by Sale & Magorrian (2014).Other autocovariance functions may be used as substitutes so long as they have the appropriate properties, i.e. so long as they define a density field that is continuous and rough, in which case they seem to adequately describe the fractal behaviour of the ISM.These requirements exclude the popular squared-exponential autocovariance function, which defines continuous but smooth random fields.The use of a squared-exponential autocovariance function will result in smooth maps, the features of which will be washed out.Nevertheless, the Kolmogorov-like function has a physically meaningful parameter tuple that library functions like the powered exponential and the Gneiting function do not.However, by comparing the integral length scale of a substitute autocovariance function with that of the Kolmogorov-like autocovariance function we may find the physical length scale that corresponds to its characteristic length scale.This physical length scale approximates the energy-injection scale, and hence gives physical meaning to the characteristic length scale of the substitute function.
The maps we have made using RJCE extinction data combined with line-of-sight distances inferred by Bailer-Jones et al. ( 2021) broadly agree with those of Rezaei Kh et al. (2020) and Rezaei Kh. & Kainulainen (2022).In particular, they corroborate Rezaei Kh's claim that a foreground cloud exists at a distance of 350 kpc.However, the maps we have made using StarHorse data (Anders et al. 2019) fail even to identify the head of Orion A. Crucially, these maps fail validation, alerting us to the fact that they are untrustworthy.The unfiltered StarHorse catalogues, though wonderful, must be used with care.They contain many sources, such as young stellar objects, that violate the assumptions of its methodology, and which therefore have spurious extinction entries that are not flagged as such.
that we construct using RJCE extinction data combined with lineof-sight distances inferred by Bailer-Jones et al. ( 2021) pass validation.But those that we construct using StarHorse data (Anders et al. 2019) do not.We plot the standardized leave-one-out residuals for these predictors in Figures C1 and C2.We do not expect the residuals to have Gaussian distribution, but they must obey Chebyshev's inequality (see App. B).

APPENDIX D: THE AUTOCOVARIANCE FUNCTIONS
We have introduced the functions f and g (16 and 17), which give the covariance of two points drawn from the density and extinction fields in terms of their line-of-sight distances and angular separation.We plot these functions in Figures D1, D2, and D5 for Ω = 1 and γ = 11/3 under the assumption of unit variance, σ 2 ∆ρ , and unit length scale, L 0 .For a given angular separation, θ 12 , the functions f and g are both symmetric in s 1 and s 2 .In the case that θ = 0 deg, the diagonals give the variance at a line-of-sight distance s 1 = s 2 , while in the case of θ 12 > 0 deg the origins give the variance, and the diagonals give the covariance at a common distance.
For all values of θ, the function f is peaked on the diagonal, and approximately zero far from it.For θ = 0 deg it is constant along the diagonal, and for θ > 0 deg it decreases along the diagonal from a maximum at the origin.The greater the angular separation, the faster it decreases.This behaviour is seen more clearly seen in Figure D3, which shows f for fixed s 1 and θ 12 and variable s 2 .
For θ = 0 deg, the function g is linear in s 1 and constant in s 2 above the diagonal.Below the diagonal it is constant in s 1 and linear in s 2 .This is to say that, along a line of sight, the extinction at a point is strongly correlated with the extinction at a second point more distant since that extinction must be at least as great, and decreasingly correlated with the extinction at a point less distant, since that extinction must be smaller.There is a smooth transition across the diagonal, more clearly seen in Figure D4, the width of which is approximately L 0 .This smooth transition is due to the fact that densities, and hence extinctions, are strongly correlated at distances of order L 0 or less, whether the second point is more or less distant than the first.For θ > 0 deg, the function g above the diagonal is linear in s 1 up to some critical value and constant in s 1 beyond that critical value, whilst being constant in s 2 .Below the diagonal the function g is constant in s 1 , whilst being linear in s 2 up to the critical value and constant in s 2 beyond the critical value.The greater the angular separation, the smaller the linear regime.Along distinct lines of sight, the extinction at a point is less strongly correlated with the extinction at a second point more distant since that extinction need not be as great and decreasingly correlated with the extinction at a point less distant since that extinction must again be zero at the origin.Again, there is a smooth transition across the diagonal, more clearly in Figure D4, due to the fact that extinctions are strongly correlated at distances of order L 0 or less.

D1 Alternative covariance functions
We have also discussed the use of autocovariance functions that have no direct physical motivation, and in particular the use of the squared-exponential and Gneiting functions (Sec.6).In Figure D6 we show the graphs of these, alongside the Kolmogorov-like autocovariance function, all for the case of unit variance and unit characteristic length scale.

Figure 1 .
Figure 1.Recovery of the parameters of the autocovariance function for the case of synthetic data.The parameters of the autocovariance function, σ 2 ∆ρ and L 0 (Eqs 9 and 5), may be chosen so as to minimize the leave-one-out cross-validation score, R LOO (Eq.31).Here, a minimum in R LOO is found for σ 2 ∆ρ = 3.0 × 10 −6 mag 2 pc −2 and L 0 = 98 pc.The true values of the parameters are σ 2 ∆ρ = 1 × 10 −6 mag 2 pc −2 and L 0 = 100 pc.(See Sec. 4 for discussion.)

Figure 2 .
Figure 2. Validation of the best linear predictor for the case of synthetic data.The standardized leave-one-out residuals of the predicted extinction may be used to validate our maps.Here, these standardized residuals are seen to have approximately unit standard deviation, and exhibit no trends in s, l, or b.Our maps therefore pass validation.(For the sake of clarity these plots show a sample of size 1000, uniformly distributed in volume.)(See Sec. 4 for discussion.)

Figure 3 .
Figure 3.Our synthetic extinction field, A (left), and the field recovered from a noisy sample of it using our method, Â (right), in the plane b = 3.5 deg.(See Sec. 4 for discussion.)

Figure 7 .Figure 8 .
Figure 7.Our synthetic fields A (top) and ρ (bottom), along with the fields recovered using our method, along the line of sight l = 7 deg, b = 3.5 deg.In each case the synthetic field is shown as a dashed line, the recovered field is shown as a solid line, and the three-RMSE prediction interval shown as a shaded band.(See Sec. 4 for discussion.)

Figure 9 .
Figure 9. Recovery of the parameters of the autocovariance function for the case of RJCE and Gaia data.The leave-one-out cross-validation score, R LOO , is found to have a minimum for σ 2 ∆ρ = 6.2 × 10 −5 mag 2 pc −2 and L 0 = 98 pc.Contours are shown at intervals of 5×10 −5 mag 2 .(See Sec. 5 for discussion.) ), and the lines of sight l = 206 deg, b = −17 deg, and l = 210 deg, b = −19.5 deg (Fig.

Figure 14 .
Figure 13.Colour excess and Gaia data: extinction, A Ks (top), density, ρ (bottom), and their predicted values, ÂKs and ρ, for the lines of sight l = 206.5 deg, b = −17 deg and l = 210 deg, b = −19.5 deg.These pass through the two lobes of the foreground cloud.The three-RMSE prediction intervals are shown as shaded bands.(See Sec. 5 for discussion.)

Figure D1 .Figure D2 .
Figure D1.Covariance of the density of the ISM at two points, r 1 and r 2 , as a function of the distances, s 1 and s 2 , and the angular separation, θ 12 (the function f , eq. 16) for Ω = 1 and γ = 11/3.In each panel the angular separation is fixed, and the distances allowed to vary.The covariance is shown on an arcsinh scale with linear width 0.001.

Figure D3 .Figure D4 .Figure D5 .Figure D6 .
Figure D3.Covariance of the density of the ISM at two points, r 1 and r 2 , as a function of the distances, s 1 and s 2 , and the angular separation, θ 12 (the function f , eq. 16) for Ω = 1 and γ = 11/3.In each panel the distance s 1 and the angular separation are fixed, and the distance s 2 allowed to vary.The covariance is shown on an arcsinh scale with linear width 0.001.