The complex multi-scale structure in simulated and observed emission maps of the proto-cluster cloud G0.253+0.016 (`the Brick')

The Central Molecular Zone (CMZ; the central ~ 500 pc of the Milky Way) hosts molecular clouds in an extreme environment of strong shear, high gas pressure and density, and complex chemistry. G0.253+0.016, also known as `the Brick', is the densest, most compact and quiescent of these clouds. High-resolution observations with the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed its complex, hierarchical structure. In this paper we compare the properties of recent hydrodynamical simulations of the Brick to those of the ALMA observations. To facilitate the comparison, we post-process the simulations and create synthetic ALMA maps of molecular line emission from eight molecules. We correlate the line emission maps to each other and to the mass column density, and find that HNCO is the best mass tracer of the eight emission lines within the simulations. Additionally, we characterise the spatial structure of the observed and simulated cloud using the density probability distribution function (PDF), spatial power spectrum, fractal dimension, and moments of inertia. While we find good agreement between the observed and simulated data in terms of power spectra and fractal dimensions, there are key differences in the density PDFs and moments of inertia, which we attribute to the omission of magnetic fields in the simulations. This demonstrates that the presence of the Galactic potential can reproduce many cloud properties, but additional physical processes are needed to fully explain the gas structure.

However, despite its high density, the Brick exhibits a very low star formation rate (SFR), which is a long-standing puzzle (e.g. Longmore et al. 2012;Rathborne et al. 2014a;Henshaw et al. 2019). The same applies to the CMZ as a whole (Longmore et al. 2013a; Barnes et al. 2017), even though some individual clouds, such as Sgr B2, are actively forming stars ) and stellar clusters (Ginsburg & Kruĳssen 2018). One explanation for the presently low SFR in the region is the idea that the CMZ undergoes bursts of star formation, followed by quiescent phases (Kruĳssen et al. 2014;Krumholz et al. 2017;Armillotta et al. 2020;Sormani et al. 2020;Tress et al. 2020). Longmore et al. (2013b) and Barnes et al. (2019Barnes et al. ( , 2020 reported an increase of star formation activity as a function of orbital position angle among the subset of CMZ clouds known as the 'dust ridge' (which includes the Brick). This idea implies that the 'dust ridge' clouds may follow an evolutionary sequence, where the Brick is positioned on the cusp of starting to form stars and potentially a young massive cluster. Indeed, recent work by Walker et al. (2021) and Henshaw et al. (2022) found signatures of early star formation within the cloud. Therefore, the Brick represents a unique opportunity for studying the initial developmental stages of a likely progenitor of a massive stellar cluster, similar to the Arches (Rathborne et al. 2015). Furthermore, the Brick can help us understand star formation in extreme environments, such as present-day galactic centres and high-redshift galaxies (which share similar pressures and temperatures), without the resolution restrictions present in studies of objects at so much greater distances (Kruĳssen & Longmore 2013).
The structure of the Brick has been studied with the Atacama Large Millimeter/submillimeter Array (ALMA) (Rathborne et al. 2015). These high-resolution observations revealed complex, hierarchical substructure, in what was previously seen as a smooth, dark cloud (also see Henshaw et al. 2019). Rathborne et al. (2015) detected the line emission from 17 different molecules, which trace different chemistry and physical processes within the Brick. They report high fractal dimensions and large velocity dispersions for all of their observations, which is consistent with a highly supersonically turbulent environment.
Additionally, there have been numerical efforts towards modelling the gas dynamics and star formation in the CMZ. Some authors have performed simulations that follow the gas flow towards the galactic centre and connect kiloparsec scales to the inner few hundred parsecs (e.g. Sormani et al. 2015;Armillotta et al. 2019Armillotta et al. , 2020Tress et al. 2020;Sormani et al. 2020;Moon et al. 2021). These simulations allow us to study the formation and large-scale properties of molecular clouds in the CMZ. However, they do not yet have enough resolution to model the small-scale structures within the individual clouds. For this reason, it is important to also consider simulations of individual clouds within the CMZ environment, as done recently by Dale et al. (2019) and Kruĳssen et al. (2019). These works have focused on individual CMZ clouds moving in an external gravitational potential chosen to closely represent the potential in the CMZ.
To determine the effect of the Galactic potential on the star formation process in the clouds, Dale et al. (2019) performed a set of isolated cloud simulations and compared them to several sets of simulations evolved on orbits through the external potential. They found that the isolated simulations had very different SFR, sizes and velocity dispersions in comparison to the simulations evolved in the potential. This indicates that the Galactic potential has an important role in shaping the CMZ clouds.
In addition, Dale et al. (2019) used two sets of initial conditions for the simulations evolved with an external gravitational potential. Previous cloud-scale simulations (e.g. Bertram et al. 2015;Federrath et al. 2016) created initial conditions by selecting a virial parameter (ratio of kinetic to gravitational potential energy) of the cloud. However, the gravitational potential adopted by Dale et al. (2019) generates a compressive tidal field at the clouds' orbits , implying that their overall behaviour and evolution cannot be predicted using the virial parameter of a cloud in isolation. Dale et al. (2019) constructed initial conditions in which the clouds have turbulent support against the compressive tidal field, and not only against their own self-gravity. In their work, they demonstrated that this 'tidally-virialised' set of clouds is a better match for the observed SFR in the CMZ than classical 'self-virialised' clouds, as the latter underwent rapid gravitational collapse, inconsistent with the observations of the Brick. Additionally, Kruĳssen et al. (2019) found that the tidally-virialised simulations reproduce known CMZ properties, such as the cloud inclinations in the plane of the sky, aspect ratios, velocity dispersion, column densities and retrograde line of sight velocity gradients.
The aim of this work is to investigate to what extent the influence of the Galactic potential may be responsible for shaping the CMZ gas on sub-cloud scales. We achieve this by comparing the cloud-scale simulations of the Brick (Dale et al. 2019;Kruĳssen et al. 2019) and the ALMA observations of Rathborne et al. (2015). For completeness, we use both the tidally-virialised and self-virialised simulations, and characterise their substructure. In order to facilitate the comparison, we post-process a snapshot of each simulation and create synthetic ALMA emission maps. The exact steps for generating the synthetic observations are discussed in Section 2. We then study the structure of the emission maps in Section 3, and perform a quantitative comparison to the observations in Section 4. Finally, our main conclusions are summarised in Section 5.

SPH simulations
In this work we have used two of the hydrodynamics simulations presented in Dale et al. (2019) and have created synthetic emission maps based on these. We have selected one snapshot per simulation, corresponding to the present day position of the real Brick cloud. The simulations were performed with the SPH (smoothed particle hydrodynamics) code gandalf (Hubber et al. 2018). The chosen simulations are the high-density, tidally-virialised model (TVir) and the high-density, self-virialised model (SVir) from Dale et al. (2019). The former was used as an analogue to the Brick in the analysis of Kruĳssen et al. (2019), due to its approximately appropriate size, surface density, and velocity dispersion. The simulated clouds both have the same initial gas mass of approximately 4.5 × 10 5 M ⊙ , equally divided between 10 6 SPH particles. Both clouds were selected from a set of randomly generated velocity fields to have negative spin angular momentum consistent with the shearing motions in the clouds observed upstream from the Brick, and have been evolved along the eccentric orbit of Kruĳssen et al. (2015) and Henshaw et al. (2016). Along their orbit, the clouds have been experiencing their own self-gravity, as well as the static external gravitational potential, derived from the observed stellar distribution in the CMZ (Launhardt et al. 2002). The difference between the two simulations is in their initial velocity dispersion. While the SVir model has been assigned velocities based on an assumed internal virial parameter, the TVir model's velocities have been chosen to balance out the compressive tidal field, caused by the external potential. As a result, the SVir model is under-supported against the tidal forces and it undergoes more rapid collapse (Dale et al. 2019). While the clouds have been allowed to collapse and form sink particles, corresponding to sites of star formation, stellar feedback has been omitted from the simulations in order to isolate the effects of the external potential and the eccentric orbit.
The density threshold used for the sink creation is sink = 10 −17 g cm −3 , which leads to sink accretion radii sink = 0.035pc. The chosen threshold density is lower than the maximum resolvable density (1.9×10 −17 g cm −3 ), as given by the Bate & Burkert (1997) criterion, which ensures numerical accuracy of the behaviour of gas near the Jeans limit. By the time the simulated clouds reach the present-day position of the Brick along its orbit (after 0.74 Myr of evolution), ∼ 55% (TVir) and ∼ 65% (SVir) of their gas mass is transformed into sink particles. Note that at the particle resolution of these simulations ( part = 0.4463 M ⊙ ) the sinks do not represent individual stars, but only dense clumps that have collapsed 1 . For the subsequent analysis in this paper we ignore the sinks and only consider the remaining gas particles. The gas has been evolved using an isothermal equation of state, with a temperature of 65 K, which is a typical value observed in the CMZ (Ao et al. 2013;Ginsburg et al. 2016;Krieger et al. 2017). The assumed mean molecular weight is = 2.35. In their initial state, the clouds have been assigned a velocity field, consistent with turbulent motion, with a power spectrum of the form ( ) ∝ −4 , and resulting in virial parameters = 3.2 (TVir) and = 2.6 (SVir). These values later drop to 1.24 and 1.16 for TVir and SVir, respectively, by the time the clouds reach the present-day position of the Brick (which corresponds to the analysed snapshots in this paper). Even though there is no explicit subsequent driving of turbulence, the presence of the strong shear in the external gravitational potential maintains a large velocity dispersion ).

Post-processing work flow
We have created two sets of synthetic ALMA emission maps, one per snapshot of each simulation. The snapshots correspond to the present day orbital position of the Brick. The emission maps are produced as a result of post-processing the snapshots with the radiative transfer code polaris (Reissl et al. 2016), and then further with the casa (Common Astronomy Software Applications) package (McMullin et al. 2007). To facilitate these processing steps, we have made several intermediate steps, as shown in Figure 1. First, the SPH particle positions of each snapshot have been used in order to generate a Voronoi tesselation grid. The grid has been regularised using Lloyd's algorithm (Lloyd 1982) to ensure that the gas structure is well resolved. After the grid construction, SPH properties such as density and velocity have been mapped onto the grid cells, using the mapping method described in Petkova et al. (2018). All of the above data have been combined to create an input file for polaris, which is used for performing the desired radiative transfer. More detailed descriptions of each of the individual steps of the process are given below.

Choice of Voronoi grid
In order to post-process the simulation snapshots with polaris, it is necessary to use a grid representation of the gas distribution, as the radiative transfer code cannot be executed directly on the SPH particles.
1 Each sink particle contains tens of gas particles. Therefore, with the particle mass being part = 0.4463 M ⊙ , the lightest sink has a mass much greater than 1 M ⊙  Figure 1. A schematic of the post-processing work flow. The SPH data are used for the construction of a Voronoi grid, which is regularised with Lloyd's algorithm (Lloyd 1982). The SPH densities and velocities are then mapped onto the grid and, together with the grid's neighbours lists, are combined into a polaris grid file. We additionally create a polaris command file, which contains our assumed abundances of molecular tracers (see Table 1). polaris is executed using the grid file, the command file, and information on molecular transitions from the LAMDA database (Schöier et al. 2005). Finally, the output emission maps from polaris are further processed with the casa package (McMullin et al. 2007) to create synthetic ALMA observations.
A Voronoi grid is an unstructured mesh built around a set of generating sites. The grid is constructed such that each cell encloses all points in space closer to its own generating site than to any other generating site. This results in all grid cells in 3D having the shape of convex polyhedra (Dirichlet 1850;Voronoi 1908).
The benefit of adopting a Voronoi grid to represent SPH data is that a large span of length scales can be covered by a small number of grid cells. A Voronoi grid can follow the SPH particle distribution and hence provide better resolution in regions that require it. The most intuitive way of constructing the grid is, therefore, by using the particle positions as generating sites (which we will refer to as a simple Voronoi grid).
This approach does have some disadvantages, as discussed by Koepferl et al. (2017). Due to the geometric property of a Voronoi grid to bisect the distance between two neighbouring generating sites, grid cells on the interface between a high and a low density region end up being artificially elongated and make the dataset noisy. The left-hand panels of Figure 2 illustrate this issue. We show three density slices through the Voronoi grid of the TVir snapshot (at = 55, 62, 65 pc), which all exhibit sharp edges and low-density "noise" around their periphery. To counter that problem Koepferl et al. (2017) proposed that additional generating sites are inserted at chosen locations. Here we propose an alternative approach, consisting of two steps. The first step is to insert a small number (TVir: 770; SVir: 796) of additional generating sites in the empty areas sur-  , with darker blue being higher density. The plots on the left use the SPH particle positions as generating sites, while the plots on the right have grid cells regularised with Lloyd's algorithm with 5 iterations, and represent the cloud structure more cleanly.

Mass-weighted PDFs
Simple Voronoi grid Regularised Voronoi grid SPH snapshot Figure 3. Column density PDFs of the simple Voronoi grid, the Voronoi grid regularised with Lloyd's iterations, and the SPH snapshot of the TVir cloud. The left-hand panel presents the volume-weighted PDFs, and the right-hand panel contains the mass-weighted ones. The three underlying column density maps have been projected onto the plane of the sky, with a resolution of 3008 × 3008 pixels, and a pixel size of 0.33 ′′ (0.013 pc). The vertical dotted lines indicate the approximate equivalent of the volume density threshold for sink formation in the simulations (see the text). rounding the cloud. We have arranged these on a coarse hexagonal grid. This is necessary because polaris requires a cubic simulation volume as an input (see Section 2.5), and the simulated clouds have one spatial dimension with much smaller range than the other two. The second step is to use Lloyd's algorithm to regularise the shapes of the grid cells (Lloyd 1982). Lloyd's algorithm follows an iterative procedure. As an initial step, we construct a grid which uses the SPH particle positions as generating sites. We then move each generating site to the centroid of its cell and reconstruct the grid.
The repositioning of generating sites and grid reconstruction are repeated for several iterations (here we have chosen to use five, ensuring visual convergence of the density distribution). This process regularises the cell shapes and it produces a clean final look, as can be seen in the right-hand panels of Figure 2. Note that with each iteration the grid loses some information about the SPH particle distribution. Therefore, in the limit of infinite iterations, the grid cells would approach a uniform distribution.
In addition to the visual comparison of the two Voronoi grids of the TVir snapshot in Figure 2, we show the volume-weighted and mass-weighted probability distribution functions (PDFs) of the column density maps of the two grid representations and the original SPH snapshot in Figure 3. We only show the PDFs of the TVir snapshot, since those of the SVir snapshot are almost identical. The first thing to notice about the PDFs shown in Figure 3 is that there is a good agreement between the three density distribution representations for column densities Σ ≥ 16 M ⊙ pc −2 . For lower column densities, however, there are some discrepancies. When considering the volume-weighted PDFs, we see that the shape of the SPH snapshot curve vaguely resembles a lognormal with a flattened tail towards low densities. The PDF peaks at log(Σ /(M ⊙ pc −2 )) ≈ 3, with a rapid decline towards higher densities. Below the peak, there is an approximately flat, but decreasing region, until log(Σ /(M ⊙ pc −2 )) ≈ −3, below which the density PDF decreases more rapidly. A similar trend is followed by the PDF of the regularised Voronoi grid, but the approximately flat region extends to lower column densities (log(Σ /(M ⊙ pc −2 )) ≈ −5). This excess of low-density pixels in the regularised Voronoi grid with respect to the SPH snapshot comes from the fact that some of the zero column density pixels in the SPH column density map have very low but non-zero densities in the column density map of the regularised Voronoi grid. This is the result of a smoothing effect, where a small fraction of a particle's smoothing kernel volume overlaps with a grid cell at the periphery of the cloud, and hence produces a non-zero mass (and density) in the cell.
By contrast, the simple Voronoi grid does not have any pixels with zero column density by construction (since each grid cell contains a particle). As a result, its density PDF has a second, higher peak at log(Σ /(M ⊙ pc −2 )) ≈ −1, and no pixels with log(Σ /(M ⊙ pc −2 )) < −2. When considering the mass-weighted column density PDFs, which highlight the densities where the mass is primarily contained, the three representations of the snapshots are considerably more similar. The density PDFs of the SPH distribution and regularised Voronoi grid are both single-peak curves with close to identical shapes. The simple Voronoi grid matches the position of this peak in its PDF, but it contains an additional feature at lower densities, corresponding to the mass contained in its elongated cells at the cloud periphery.
In addition, due to the finite resolution of the simulations, the high-mass end of the column density PDFs decreases more steeply than it should in reality. Gas above the volume density threshold of sink formation is trapped in sink particles and it is not included in the PDFs. We convert this threshold to column density using Σ sink = 2 sink sink ≈ 1000 M ⊙ pc −2 , and we include it in Figure 3 as vertical dotted lines.
Together, Figure 2 and Figure 3 demonstrate that the regularised Voronoi grid is a better match of the SPH snapshot than the simple Voronoi grid. Therefore, for the rest of this work we adopt the regularised Voronoi grid. We perform the grid construction using the C++ library Voro++ (Rycroft 2009).

Mapping of SPH parameters onto grid cells
In order to perform the line radiative transfer, we must assign density and three-dimensional velocity to each grid cell of the Voronoi grid. We do this self-consistently by calculating the average cell density and velocity directly from the SPH dataset. This method is an extension of the mapping presented in Petkova et al. (2018) and it can be derived directly from the SPH formalism.
Let us consider a fluid property , which is defined for each particle . We can express at an arbitrary position r as: where is the total number of particles, and are the mass and density of , r is the position vector of , and is the kernel function. The kernel function is approximately Gaussian-shaped and it equals zero for |r − r | greater than a certain distance, given as a multiple of the smoothing length ℎ . The most commonly used SPH kernel function, and the one adopted by this work, is the cubic spline, defined as (Monaghan & Lattanzio 1985): (2) If we now consider a grid cell, the volume-averaged value of in cell is given as: where is the volume of cell . Combining equations 1 and 3 we can write that: Using a volume-averaged value is appropriate for performing the density mapping, because that ensures mass conservation between the SPH snapshot and the Voronoi grid. Hence we have the following expression for the cell density, : Note that in the above expression, the sum gives the cell mass, which is comprised of contributions from different particles. In contrast to the density mapping, the velocity mapping is typically done as a mass-average (since that conserves momentum). Therefore, we will use the following expression for each velocity component, , ( = {1, 2, 3}): By using the equations 6 and 7, we reduce the density and velocity mapping to simple sums, in which the only remaining unknown is the integral of the kernel function. Calculating this integral numerically is computationally expensive, hence we use the analytic form of the integral derived in Petkova et al. (2018), and we include a summary of it in Appendix A.

Processing with polaris
We perform radiative transfer post-processing of the simulated Brick with polaris (Reissl et al. 2016;Brauer et al. 2017;Reissl et al. 2019). The code uses 3D Monte Carlo radiative transfer for dust scattering and temperature calculations, and employs a ray tracing method for imaging and line emission. The calculations can be performed on a variety of grid structures, including a Voronoi tessellation. For its line radiative transfer, polaris makes use of the LAMDA database (Leiden Atomic and Molecular Database; Schöier et al. 2005), which provides pre-computed molecular parameters, such as quantum numbers, energy levels, Einstein coefficients, and collision rates with H 2 .
In this work, we perform line radiative transfer, using the large velocity gradient (LVG) approximation to calculate the level populations (Sobolev 1957). This assumes that the velocity difference between two neighbouring grid cells exceeds the thermal broadening of the line, and therefore the locally emitted photons are able to easily escape the grid. Note that polaris imposes an upper limit of 1 pc to the Sobolev length.
In order to test the validity of the LVG approximation, we use the photon interaction length, given as (Ossenkopf 2002 where is the thermal broadening of the emission line, and d /d is the local velocity gradient. The LVG approximation should be used when the interaction length for a given cell is smaller than the cell size. In order to check that the condition holds, we consider each Voronoi cell wall, and the two cells that share it. For each of these cells, we compute the velocity component perpendicular to the cell wall. We denote the absolute difference between these two velocity components as . By computing | / |, we obtain in units of the distance between the two cell generating sites. Figure 4 shows a histogram of , calculated in this way, as a function of cell density for the TVir snapshot. In this example, we use HNCO and = 65 K, which gives ≈ 158 m/s. We can see that the majority of the data points (77%) lie below the dashed line, which marks = 1. Note that the value of depends on the molecular mass; for the species that we are about to introduce we have = 146 − 200 m/s. This range in has a very minor effect on the number of datapoints below = 1 (71 − 78%), which indicates that the LVG approximation remains appropriate for our simulation snapshots.
We generate line emission maps of CS, HCN, HCO + , HNCO, HNC, H 2 CS, HC 3 N and N 2 H + , which are commonly detected molecules in CMZ clouds (e.g. Jones et al. 2012;Rathborne et al. 2015;Callanan et al. 2021). The gas temperature is assumed to be 65 K for consistency with the hydrodynamics simulation. This value also matches existing observations of the average gas temperature of the Brick (Ao et al. 2013;Ginsburg et al. 2016;Krieger et al. 2017). In addition, we adopt a constant abundance for each molecular tracer throughout the cloud. This simplified assumption will likely affect some of the properties of the emission maps, such as the correlation/decorrelation of subregions, the shape of the brightness PDFs and the power spectrum slopes (introduced later in the paper), but it will provide a reasonable starting point of our analysis. Indeed, observational studies of the CMZ region, such as that of Tanaka et al. (2018, see their fig. 18), show only modest spatial variation of abundances in the plane of the sky. The variation is typically even smaller when we only consider the extend of the Brick cloud. We account for these variations in our modelling by considering three different abundances for each molecular tracer (see below). In addition to abundance variation in the plane of the sky, there can also be variation along the line of sight, which is harder to deduce from observations. However, Rathborne et   We have assumed that = 158 m/s (HNCO at = 65 K). The dashed line denotes an interaction length equal to the separation between two neighbouring cells. Approximately 77% of the data points lie below the dashed line, which indicates that the LVG approximation is appropriate.
sity ratio for a few tracers throughout the Brick. For an optically thin tracer (such as HNCO) this constant ratio is consistent with having approximately constant abundance along the line of sight, and it further confirms that our assumptions are sensible. We defer a more in-depth chemical modelling to a future study. The molecular abundances that we assume are obtained primarily from the CMZ observations of Tanaka Table 1. The abundances are chosen as follows: • HCN, HNC: M11 report abundances of a few times 10 −9 for both species. T18 find HCN ∼ 10 −7 and HNC ∼ 10 −8 , but R18 report HCN ∼ 10 −8 for most of their tabulated sources, with HNC typically a factor of a few smaller. Therefore, we choose an abundance of 10 −8 for HCN, and 3 × 10 −9 for HNC in order to reproduce the ratio of about 3:1 between the integrated intensity of HCN and that of HNC.
• HCO + : T18 report abundances of ∼ 10 −8 , while R18 report a value a factor of a few smaller. M11 find that the abundance strongly depends on the choice of cosmic ray ionisation rate. When using cosmic ray ionisation rates comparable to the CMZ, they obtain HCO + abundances of a few times 10 −8 . Therefore, we choose an abundance of 10 −8 .
• N 2 H + : T18 report an abundance of a few times 10 −9 for most of the CMZ, and slightly higher in Sgr B2, but R18 give a value of around 10 −9 for most sources. We assume 10 −9 in this work.
• CS: T18 find abundances between 10 −8 and a few times 10 −8 , while M11 report values of around 10 −8 for cosmic ray ionisation rates consistent with the CMZ. R18 also find values of ∼ 10 −8 , and hence we assume an abundance of 10 −8 . Table 1. Summary of molecular tracer data, used as an input for polaris. The abundances have been primarily obtained from the work of Tanaka et al. (2018), Riquelme et al. (2018) and Meĳerink et al. (2011), and are discussed in more detail in the text. The frequencies of the transitions are taken from the LAMDA database (Schöier et al. 2005).

Molecule Fiducial Abundance
Abundance Range Transition Frequency (GHz) • HC 3 N: T18 report an abundance of ∼ 3 × 10 −9 for most of the CMZ, increasing to ∼ 10 −8 in Sgr B2. R18 also report values closer to 10 −8 in most sources. In consideration of these studies, we assume a fiducial abundance of 10 −8 .
• HNCO: The abundance of this molecule was not included in T18 or M11, but R18 report values between 10 −8 and a few times 10 −8 . Churchwell et al. (1986) find a value of 3 × 10 −9 in Sgr B2, and hence we assume an abundance of 10 −8 .
• H 2 CS: The abundance of this molecule was not reported in T18, M11 or R18. Therefore, we assume a value of 10 −9 , to match the observations of local clouds (e.g. Minh et al. 1991).
In addition to the abundances listed above, we consider a plausible abundance range for each molecular species (listed in Table  1). We use the extremes of these ranges to create sets of minimum abundance and maximum abundance emission maps alongside our main ones, in order to check the robustness of our results.
Using the above setup, we create a set of emission maps in different velocity channels. The line of sight of the emission maps coincides with that of the real Brick as viewed from Earth. We use 101 velocity channels, covering the radial velocity range between −50km s −1 and 50km s −1 , which includes all kinematic parts of the clouds. We have also performed tests with narrower velocity channels, which do not produce significantly different results. Note that the symmetry of the velocity range about 0 is internally imposed by polaris, whereas the velocity distributions of the simulated clouds are biased towards positive values (with a mean of ∼ 20 km/s, which is consistent with the real Brick; see Appendix B). Additionally, we do not include any sub-grid turbulent broadening of the lines. The emission maps are then integrated over the velocity range to create the final polaris output images. These integrated emission maps have dimensions of 3008 × 3008 pixels, and a pixel size of ∼ 0.013pc, corresponding to 0.33 ′′ at the distance of the CMZ. We justify the choice of pixel size in Section 2.7.

Processing with casa
As the final step towards the creation of synthetic ALMA observations from our simulation snapshots, we process the integrated emission maps from polaris with the casa 5.4.1-31 package (Common Astronomy Software Applications McMullin et al. 2007). Since we are insterested in performing a direct comparison with the data set of Rathborne et al. (2015), we closely follow their observational input parameters. Therefore, for each molecular tracer we generate six observations using the same array configurations (and hence the same -coverage) as Rathborne et al. (2015). Each observation is performed with the simobserve task and it consists of a 72 ′′ ×162 ′′ mosaic, centred around the densest part of the cloud (see the righthand panels of Figure 5). Each observation is run for 40 minutes on-source, with a precipitable water vapour of PWV = 1.5 mm. We then combine the six observations per emission line and clean the resulting image with the tclean task. The synthesised beam of the final images has an angular size of ∼ 1.9 ′′ (∼ 0.08 pc).

Resolution
It is crucial to verify whether the finite resolution of the SPH snapshots or of the intermediate post-processing steps can affect the reliability of the final synthetic emission maps. In order to do this, we compare the spatial resolutions of the SPH simulations, the Voronoi grids, the polaris output, and the synthesised beam of the ALMA setup (listed in Table 2). While the latter two have constant resolution element sizes throughout the cloud, the former two consist of a wide range of particle and cell sizes, respectively. SPH is a Lagrangian method, in which the resolution is defined by a particle mass and not by a spatial scale. Nonetheless, each SPH particle has an associated smoothing length ℎ, which is dynamically adjusted based on the local gas density. The densest regions of the clouds consist of the most compact particles, where we measure the minimum smoothing lengths of ℎ min,Tvir = 0.014 pc (tidally-virialised cloud), and ℎ min,Svir = 0.017 pc (self-virialised cloud). These high-density regions are likely to have strong contributions to the synthetic emission maps, and the fact that their local resolution (of a few times ℎ) is comparable to the synthesised beam size ( maj = 0.080 pc) suggests that the SPH snapshots resolve the density peaks sufficiently in order to produce robust synthetic emission maps.
However, since the SPH snapshots are not directly postprocessed with radiative transfer, but are first mapped onto a grid, we also need to consider the cell sizes. Similarly to the SPH simulations, the Voronoi grids have smaller cell sizes in regions of high density. Here we define the cell size as the cube root of the cell volume, due to the irregular shapes of the grid cells. We find the minimum cell sizes of the two clouds to be min,Tvir = 0.021 pc and min,Svir = 0.022 pc (tidally-virialised and self-virialised cloud, respectively). These values are between their corresponding ℎ min and 2ℎ min , which means that each high-density SPH particle is mapped onto multiple Voronoi cells. This is a desirable property and it indicates some level of resolution preservation.
Finally, we also need to consider the pixel size of the polaris output, which is pix = 0.013 pc. This value is smaller than the minimum Voronoi cell sizes, ensuring that at least one ray passes through each grid cell. However, even if this were not the case, polaris can perform ray splitting to ensure that at least one ray passes through each cell, and hence no cells remain "invisible" to  the ray-tracing algorithm. Additionally, the polaris pixels need to be smaller than the ALMA beam size. This is true by construction, since we have pix < min < 2ℎ min < min < maj . The above considerations suggest that the final emission maps should be robust and reliable. We have also tested this by altering the sizes of some of the of the underlying resolution elements, such as the SPH smoothing lengths and the polaris pixel size. We have created two smoothed (lower resolution) versions of the SPH snapshot by multiplying the smoothing lengths by a factor of 1.5 and 2 respectively. Independently from these tests, we have also created synthetic emission maps using an increased polaris pixel size of 0.04 pc. All of the above tests caused only minor changes of the final results.

SYNTHETIC EMISSION MAPS
By applying the procedure described in Section 2, we produce integrated emission maps, as well as synthetic ALMA observations of the TVir and SVir simulation snapshots. The integrated emission maps produced with polaris are shown in Figure 6. Comparing each of the emission lines, we see that their corresponding maps are visually distinct from one another, with some of them combining diffuse emission with bright small-scale emission peaks (e.g. HCN, HCO + , HNC), whereas others trace larger structures (e.g. HNCO). The differences in the size scale of the structures result from the fact that each molecule traces different densities within the clouds. In addition to the size scale variations, the emission maps span a range of peak intensities (0.15 − 10 mJy/arcsec 2 km/s). If we compare the TVir and the SVir emission maps for a given molecule, we see that they agree with each other in terms of peak intensity and visual appearance. This is not surprising, since the two simulation snapshots contain a similar range of densities (see Figure 3 and 5). scale, diffuse emission that is well traced by HNC and HCO + in Figure 6, cannot be seen in Figure 7 and 8 due to the limited uvcoverage of the synthesised beam. Through visual inspection, we find that the relative brightness of peaks in the synthetic ALMA observations vary between the different molecular tracers. Furthermore, as in Figure 6, the maximum peak intensity also varies with the choice of molecular tracer. Table 3 provides a summary of the peak intensities from Figures 7 and 8, and also provides reference values from the observations of Rathborne et al. (2015). We see that the TVir and SVir snapshots produce ALMA emission maps with very similar maximum peak brightness. Additionally, there is good agreement between the values for HCN, HCO + , HNCO and H 2 CS obtained from observations, compared to the synthetic ALMA observations of the two simulation snapshots. Note that the emission maps for HCN, HCO + and HNCO presented in Rathborne et al. (2015) combine the ALMA observations with additional single dish data, which slightly boosts the intensity peaks. The synthetic HCO + emission maps contain brighter peaks than the real Brick observations, which may be due to overestimation of the molecular abundance in our modelling.
In addition to a simple visual inspection, we can quantify the differences in the synthetic ALMA emission maps by computing   Table 3. Maximum peak intensities of the synthetic ALMA observations of the TVir and SVir simulation snapshots, and comparison to the peak intensities reported by Rathborne et al. (2015). For each simulation snapshot, we give the peak intensities of the emission maps with the fiducial abundances, as well as those with the minimum and maximum plausible abundances (see Table 1). Additionally, we provide the peak intensities both in units of [  the cross-correlation coefficients between them. Figure 9 shows the cross-correlation coefficients of the TVir emission maps. The results for the SVir emission maps are omitted due to their similarity to the TVir results. We find that all of the molecular lines are well correlated, with coefficients in the range 0.47 − 0.95. We estimate the uncertainties of these coefficients to be < 0.1 using a bootstrap technique in which we calculate the coefficients for many sub-regions of the full emission maps. In Figure 9 we can identify two groups of molecular tracers, which have strong correlations within the group, and have only moderate correlations outside the group. The first group consist of CS, HCN, HCO + and HNC, while the second group consists of HNCO, H 2 CS, HC 3 N and N 2 H + . Upon visual inspection, we find that the latter group traces some of the larger structures, while the former group mainly traces emission from compact peaks. It is worth noting that while fitting in with the latter group, N 2 H + is somewhat of an outlier. This is because the N 2 H + emission maps are largely dominated by noise, due to their overall low brightness.
If we consider only the molecular species which were also observed by Rathborne et al. (2015) (i.e. HCN, HCO + , HNCO and H 2 CS), the range of cross-correlation coefficients that we calculate becomes 0.51 − 0.95. The observations of Rathborne et al. (2015) have a broader range of cross-correlation coefficients (and somewhat lower, with 0.28 − 0.81). Encouragingly, the coefficients are distributed in the same way between the molecular couples, i.e. the lower/higher cross-correlation coefficients in the simulated emission maps correspond to the lower/higher cross-correlation coefficients in the observations. Specifically, HCN and HCO + are most strongly correlated with each other, as are HNCO and H 2 CS. This reflects the difference in excitation energy between both pairs of lines (see Table 1 of Rathborne et al. 2015). Rathborne et al. (2015) also compute the cross-correlation coefficients between the molecular line emission maps and the dust continuum emission map, the latter of which is expected to trace the gas density. They find that none of their molecular emission maps are strongly correlated with the continuum and hence the gas density, with cross-correlation coefficients in the range 0.09 − 0.39. In order to compare to their results, we compute the cross-correlation coefficients between the gas column density and the synthetic emission maps (see Figure 10). In reality, the continuum emission does  not perfectly trace the column density, since there are temperature fluctuations within the Brick. However, at 3 mm we probe the Rayleigh-Jeans tail of the dust emission, which is only weakly sensitive to temperature. Under these circumstances, using the column density and the dust continuum emission interchangeably is a reasonable approximation. Figure 10 shows the cross-correlation coefficients between the column density map of the TVir snapshot and the molecular emission maps from Figure 7, which we find to be in the range 0.41 − 0.66 for the assumed abundances. In addition, we have reconstructed the maps in Figure 7 for the minimum and maximum plausible abundances (see Table 1). With the exception of N 2 H + , we do not see a strong variation of the cross-correlation coefficients with the molecular abundance. We find that the best tracer of the column density is HNCO, which is in agreement with the findings of Rathborne et al. (2015), even though they report a significantly weaker correlation. The fact that our synthetic observations correlate a lot more strongly with the density structure than the real Brick observations is likely due to one or more of the simplifying assumptions that we have made. In particular, we have assumed constant temperature and abundances throughout the simulated clouds, while the real Brick has variations in temperature and abundances (e.g. Mills et al. 2015). Our models provide a first attempt at reproducing the different emission lines of the Brick, and can be improved in future work by including cooling, heating, chemistry, and more detailed modelling of the molecular species.

CLOUD STRUCTURE
Similarly to the emission maps of Rathborne et al. (2015), our synthetic maps show complex, hierarchical morphology. In this section, we aim to quantify the structure of the simulated Brick, using a variety of quantitative metrics, such as density probability distribution functions (PDFs), the fractal dimension, the spatial power spectrum, and the moments of inertia (through the J-plots method; Jaffa et al. 2018). Each of these metrics probes a different substructural property of the synthetic emission maps. The density PDFs show the relative occurrence rate of pixels of different brightness, while the power spectrum is sensitive to the spatial separation of the brightness peaks. The fractal dimension is a way of quantifying complexity in substructures of the cloud. A high fractal dimension indicates a turbulent cloud with winding contours. The opposite of that would be a low fractal dimension cloud, where each contour is close to circular, and it likely encloses a gravity-dominated re-   Table 1). The uncertainties of all coefficients are estimated to be < 0.1 by using a bootstrap technique. gion. Finally, the moments of inertia measured through the J-plots method characterises the shapes of the cloud substructures, and it divides them into centrally condensed, centrally rarefied, and elongated/filamentary structures.

Brightness and density probability distribution functions
We compare the relative occurrence of pixels of different brightness (see Figure 11). Since it is more physically meaningful to consider column density or volume density PDFs, and continuum emission traces the density structure better than line emission, we also include the PDFs of the 3 mm continuum emission from Kruĳssen et al.
(2019) (using the TVir snapshot only) and Rathborne et al. (2015) in Figure 12. Note that the -axis of Figure 12 is given in terms of ln( ), where is the ratio of a given pixel value and the mean pixel value included in the PDF. In the case of the continuum PDFs, we assume that the gas column density is proportional to the dust column density, which in turn is proportional the pixel intensity, and hence also equals the ratio between the column density and the mean column density. Additionally, we have excluded pixels under a certain threshold (which differs between the images) in order to separate the background noise from the emission, and to ensure that the PDFs only contain contributions from closed contours within the images (Alves et al. 2017). Figure 11 shows good agreement between the PDFs of the TVir and SVir simulations (solid and dashed lines respectively), which can be explained by similarities in the underlying density distributions of the two snapshots (see Figure 3 and 5). The shapes of the PDFs are approximately log-normal with a flattening at lower intensities. We also include the brightness PDFs of the molecular species that are observed by Rathborne et al. (2015) for comparison (i.e. HCN, HCO + , HNCO and H 2 CS; dotted lines). We find good agreement between the simulated and observed brightness distributions for H 2 CS. The other molecules show some discrepancies in terms of the peak brightness (as already shown in Table 3) and the slope of the higher-brightness end of the curves. Many of these discrepancies can be overcome by tuning the assumed abundances, and by more generally improving the chemical modelling. In particular, by slightly increasing the abundance of HNCO and slightly decreasing that of HCO + , we can mostly compensate for the discrepancies between the simulated and observed curves. In contrast, the peak brightness of HCN in the synthetic images is close to that of the real Brick, but there is a mismatch of the high brightness slope that likely requires better chemical modelling. Alternatively, some of these discrepancies may be arising from differences in the gas structure between the simulations and the Brick.
We find similar log-normal shapes in the 3 mm continuum emission (Figure 12), with widths of ln ,TVir = 0.70, and ln ,R15 = 0.26. Note that Rathborne et al. (2014b) report a PDF width of 0.34 ± 0.02 for the same observed 3 mm continuum emission map. The difference here likely arises from a slightly different choice of fitting range for the log-normal curve. In order to consolidate the two results, we will use ln ,R15 = 0.30 ± 0.05 (and also ln ,TVir = 0.70 ± 0.05) in all calculations below. Finally, the observed 3 mm continuum emission map has a high-density tail in its PDF, which deviates from the log-normal shape. This feature is caused by a few bright peaks, which have been attributed to regions of ongoing star formation (Rathborne et al. 2014b;Walker et al. 2021). The feature is not present in the continuum PDF of the TVir snapshot, despite the more prominent gravitational collapse in the simulation. This is due to the fact that the smallest (sub-sink radius; ≥ 10 −17 g cm −3 ) scales of star formation remain unresolved in our models, and any density enhancements are rapidly converted into sink particles, which are not included in the emission maps.
There is an enormous body of literature showing that the width of a log-normal volume (not column) density PDF can be expressed in terms of a small number of parameters (e.g. Vazquez These are the turbulence driving parameter, 3 , the 3D sonic Mach number, M, and the turbulent plasma beta, . Therefore, in order to link ln ,TVir and ln ,R15 to the above physical parameters, we first need to obtain the volume density dispersions. To do so, we follow the procedure used in Federrath et al. (2016), who analyse the PDF of the 3 mm continuum emission map of the Brick, presented in Rathborne et al. (2014b). First, we find the column density dispersion = [exp( 2 ln ) − 1] 1/2 (e.g. Price et al. 2011), which gives ,TVir = 0.80 ± 0.07 and ,R15 = 0.31 ± 0.05. Next, Federrath et al. (2016) find the conversion factor between the column density and the volume density dispersion, R 1/2 = 0.28 ± 0.11, using the method of Brunt et al. (2010). Note that the large uncertainty on the value of R 1/2 comes from the fact that the method assumes an isotropic mass distribution in the cloud, whereas Federrath et al. (2016) find a moderate level of anisotropy in the Brick. They attribute the anisotropy to the cloud's strong ordered magnetic field (Pillai et al. 2015), which is not present in the simulations, and hence we do not expect a larger uncertainty when applying this method to the TVir snapshot. By dividing by R 1/2 , we get the volume density dispersions / 0 ,TVir = 2.8 ± 1.1 and / 0 ,R15 = 1.1 ± 0.5, where 0 is the average volume density. Our estimate for / 0 ,R15 is consistent with that of Federrath et al. (2016), who find the value 1.3 ± 0.5.
For the real Brick, Federrath et al. (2016) find that M = 11 ± 3 and = 0.34 ± 0.35, which let them deduce that = 0.22 ± 0.12, using the relation Using our value for / 0 ,R15 , we find = 0.20 ± 0.13, which is consistent with the above result and with solenoidal turbulence driving.
Within the simulations, the parameters on the right-hand side of eq. 9 are either known or can be measured. First, the simulations do not include a magnetic field, which means that → ∞, and hence (1 + −1 ) −1/2 = 1. Second, by using a Helmholtz decomposition, we can represent the velocity field of the TVir simulation as the sum of a solenoidal component and a compressive component. Using the kinetic energy of the two components, we deduce that the turbulence driving is solenoidal (caused by by the strong shear arising from the Galactic potential), and hence ≈ 0.3 (Petkova et al. in prep.). Finally, we calculate M = 7.7±0.6, using the size-linewidth relation of the simulated cloud (Petkova et al. in prep.). Combining these parameters, we find that the right-hand side of eq. 9 equals 2.2±0.2. This is in agreement with the measurement of / 0 ,TVir = 2.8±1.1, as expected.
The direct measurements of the above parameters in the simulations allow us to investigate what causes the discrepancy between / 0 ,TVir and / 0 ,R15 , by exploring possible variations in , M, and . First, the turbulence driving parameter is in agreement between the simulated and observed cloud, and it is unlikely to change unless we alter the gravitational potential in the simulation or add other sources of turbulence that cause cloud-scale compression. Second, the Mach number that we measure in the simulated cloud is lower by about 1 compared to the observed Brick. If we used the observed value of M, the right-hand side of eq. 9 would become even larger and shift away from the / 0 ,R15 value. Instead, the simulation would need to have an even smaller Mach number (by a factor of ∼ 2) in order to match the observed density PDF. By necessity, this would lower the velocity dispersion of the simulated cloud, and it would cause it to undergo rapid collapse and star formation. Finally, if we were to modify the simulation to include the same turbulent magnetic field strength as the real Brick, the density dispersion would be modified by a factor of (1 + −1 ) −1/2 ≈ 0.5. The resulting right-hand side of eq. 9 would then predict / 0 = 1.1, which is in excellent agreement with / 0 ,R15 = 1.1 ± 0.5. The caveat here is that the measured value of in the Brick has a large error bar, and also that by introducing a magnetic field in the sim- ulations, one might also end up altering the Mach number. Despite these concerns, our analysis indicates that the discrepancy in the widths of the density PDFs may be attributed to the lack of a magnetic field in the simulations.

Spatial power spectrum
There is no explicit turbulence driving in the simulations of Dale et al. (2019). Despite that, Kruĳssen et al. (2019) find that the clouds maintain large velocity dispersion along their orbit, which is attributed to the presence of strong shear. To study this behaviour we compute the spatial power spectra of the synthetic emission maps, using TurbuStat (Koch et al. 2019). The power spectra are constructed by performing a 2D Fourier transform of the emission maps, and plotting the 1D radial profile of the resulting images as a function of length scale. Figure 13 shows the power spectrum of the TVir column density (top panel), alongside the power spectrum of the HNCO emission map of the TVir snapshot (bottom panel). The power spectrum of the column density follows a broken power-law, with a clearly defined break at ∼ 0.32 pc. Note that the position of the break does not change if we vary the pixel size of the column density map that we use. In contrast, the HNCO emission map power spectrum has a more complex shape, which does not overall resemble a single power law. As a result of the interferometric nature of the synthetic emission maps, the large spatial scales are filtered out, whereas the small spatial scales are correlated due to the beam profile. Therefore, there is a narrow window of spatial scales where we can trust the power spectrum values. We exclude spatial scales below three beam sizes (using the geometric mean of the major and minor beam axes; see Koch et al. 2020), which corresponds to values above 5.16 pc −1 in the power spectrum. We also exclude spatial scales above the maximum resolvable scale of 10.5 ′′ (see ALMA Cycle 0, Extended Configuration, Band 3), corresponding to values below 2.38 pc −1 . This range is marked with vertical dashed lines in Figure 13, and we find that within it the power spectrum can be fitted with a single power law. The slopes of the remaining TVir and SVir emission maps are summarised in Table 4.
The power spectrum of the TVir column density has two linear slopes: 1 = −2.25 ± 0.01 for large spatial scales, and 2 = −4.09 ± 0.01 for small spatial scales. This is similar to the results of Rathborne et al. (2015), who find a break in the 3 mm continuum emission power spectrum of the real Brick at 0.12 pc, with 1 ∼ −2.5 and 2 ∼ −8.0. However, note that the break in the power spectrum of the real Brick occurs very close to the beam scale, which means that it is possible for it to be an observational artefact (see Koch et al. 2020). The differences in the power spectrum parameters between the real Brick and the simulation reflect differences in the underlying physical parameters. For example, one explanation for the slight discrepancy in 1 may be the fact that the simulations do not include the effects of magnetic fields. Indeed, Padoan et al. (2004) find that their simulation with equipartition of kinetic and magnetic energy produces a slope of −2.25±0.01, while their super-Alfvénic one generates a steeper slope of −2.71 ± 0.01. However, there are many factors that influence the slope of the power spectrum, and the responsible physical processes cannot be uniquely identified. Typically, a break in the power spectrum occurs at a size scale where energy is being injected into the system, or dissipated more efficiently. Since the effect of the shear (and hence turbulent energy injection) is negligible at this size scale, we consider energy dissipation due to thermal processes or gravitational collapse. A potentially relevant length scale for these effects is the Jeans length. The average volume density within the simulations is difficult to define due to their prevalent substructure, but sensible values are bracketed by the initial volume-averaged volume density and the evolved mass-weighted average volume density for them is in the range 0.4−6 × 10 4 cm −3 . This corresponds to a Jeans length of 0.13−0.50 pc. This result is similar to the findings of Rathborne et al. (2015), who calculated the Jeans length of the Brick to be 0.10 pc, by assuming an average volume density of ∼ 10 5 cm −3 . Rathborne et al. (2015) deduced that the break in the power spectrum of the Brick is indeed related to the Jeans length. The break that we see in the simulations is also consistent with a Jeans-like scale, however it could be caused by other factors, such as the sonic scale or numerical dissipation. Because the current paper focuses on structural properties and a full kinematic analysis is needed to distinguish these interpretations, we will examine this in-depth in a follow-up study that focuses on the kinematic properties of the simulations.
Finally, Table 4 shows that the synthetic emission maps have power spectrum slopes between −4.2 and −2.3 (if we exclude N 2 H + with its high levels of noise). The group of molecules with compact emission (i.e. CS, HCN, HCO+ and HCN) have shallower slopes (above -3.5), while those that have more extended emission have steeper slopes (below -3.5). These slopes are broadly consistent with the results of Rathborne et al. (2015), who find values about −3 for HCN, HCO + and HNCO (see also Henshaw et al. (2020)). Since the emission maps of Rathborne et al. (2015) consist of ALMA observations combined with single dish data, their power spectra were measured over a greater range of spatial scales than ours, and hence their slopes are more robust. This difference in available power law fitting range can account for some of the observed slope discrepancies.

Fractal dimension
Gas clouds exhibit self-similarity on a range of size scales. However, the self-similar density profile of a star-forming core appears very different from the self-similar structure of a supernova remnant. The former has smooth, spherical density contours, while the latter has a complex, turbulent structure. This difference can be captured by the concept of fractal dimension, which corresponds to the level of complexity within a structure.
We have calculated the fractal dimension of the synthetic emission maps of the Brick, using two different methods -the perimeter-area method and the box-counting method. Both of these methods consider contours at different pixel thresholds within an emission map, and extract the fractal dimension from their properties. Note that here we use the 2D fractal dimension, since we measure it using 2D emission maps. Therefore, the fractal dimension is found within the range between 1 and 2, where 1 corresponds to smooth, circular contours, and 2 to high structural complexity.

Perimeter-area method
One way of calculating the fractal dimension of an observed cloud is by considering the area, , and the perimeter, , of different contours in the image. These two quantities are linked by the fractal dimension, , in the following way (Vogelaar & Wakker 1994): This means that an object with non-complex structure (such as a 2D-Gaussian) has ∝ 1/2 and = 1, while an object filling up the space with the length of its winding contours has = 2. In practice, we have cloud images consisting of pixels, which create limitations on our ability to accurately calculate and . To obtain we select all pixels with values above a chosen level (that is varied) and we search for connected regions within that subset. A region is connected when each of its pixels shares a wall with at least one other pixel in the region. Counting the number of pixels in a connected region gives , and is obtained by counting the number of pixels externally adjacent to the connected region (i.e. sharing a wall). Regions containing fewer than a chosen number of pixels are considered to be poorly resolved and hence we have not included them in the calculation of . Vogelaar & Wakker (1994) demonstrate that this threshold should be at around 20 pixels. However, for an interferometric observation, the synthesised beam is the relevant resolution element, which is typically represented by multiple pixels. We limit the area of the smallest structures included in this measurement to be at least 4 times the area of the synthesised beam. This limit was empirically determined, based on selecting ( , ) pairs that follow a single power law. In the case of the TVir and SVir emission maps, this size is 200 pixels, while in the real observations of the Brick, it is 65 pixels, due to differences in pixel sizes. Finally, the way in which we chose to determine and is not unique. Connected regions can include pixels sharing a diagonal border and not only a wall. The perimeter can include the sum of two pixel walls instead of just counting the pixel itself in parts where Table 5. Fractal dimension, , computed with the perimeter-area method for the simulated and observed emission maps of the Brick. Note that we have only included the observational results for those molecules from Rathborne et al. (2015), which were also used with the simulation. The results reported by Rathborne et al. (2015) follow the addition interpretation.  Vogelaar & Wakker (1994) and Federrath et al. (2009), tweaking the parameters of this method does not significantly change the value of , but only the proportionality constant in the power law.
We encounter two types of interpretations of the perimeter-area method found in the literature. One calculates the area and perimeter of each individual connected region and uses it as a data point to fit a line through (we refer to it as scatter approach), while the other sums up the areas and perimeters of all connected regions for a given pixel level (we refer to it as addition approach). Examples of the first interpretation can be seen in Vogelaar & Wakker (1994) and Sánchez et al. (2005), and the second one was used by Rathborne et al. (2015) 4 . More importantly, the two approaches are not equivalent and hence the fractal dimensions that they yield should not be compared to each other directly. This point can be illustrated if we consider a scenario of the scatter approach where all points lie on the line log( ) = 2 log( ).
By adding the areas and perimeters of two of the data points, we get a new data point which does not, in general, lie on the same line: On a more fundamental level, the difference between the interpretations resembles the distinction between intensive and extensive properties of matter. An intensive property is one which is scale free (such as temperature or pressure), while an extensive one depends on the length scale over which it is measured (such as mass or volume). The scatter approach yields an intensive property of the cloud, as its output would remain the same if we were to include an additional cloud in the emission map, with the same intrinsic fractal dimension as the original data. On the other hand, the addition approach can, under certain circumstances, produce a different result by including an additional structure (see Appendix C). What separates the outcome of the addition method from a true extensive property, is that it is not additive, i.e. we cannot compute the fractal dimension in two sub-regions of the cloud and add them to obtain the total fractal dimension. Due to this fundamental discrepancy between the two interpretations, we recommend using the scatter 4 Most authors do not explicitly state their interpretation of the perimeterarea method. Since the addition and scatter interpretations typically have visually distinct data sets (the former consist of monotonically increasing ( , ) pairs, while the latter appears as scattered data points about a line; see Figure 14), we use that as a diagnostic tool.
approach when computing and comparing the fractal dimension of gas clouds. In order to compare our emission maps to the existing literature, we compute the fractal dimension using both interpretations of the perimeter-area method. Figure 14 shows the addition and scatter interpretations, applied to the HNCO emission map of the TVir snapshot. We see that the addition interpretation produces a much steeper slope and hence a higher fractal dimension than the scatter interpretation. This behaviour continues for the remaining synthetic and real emission maps of the Brick (see summary of the results in Table 5). We will consider the implications of the measured fractal dimensions further below.

Box-counting method
An alternative way of calculating the fractal dimension can be achieved by using the box-counting method. Similarly to the perimeter-area method, we select the pixels with values above a chosen threshold (specified below). We then cover the image with a grid of boxes and count the number of boxes, ( ), containing the selected pixels as a function of box size, . By varying , we measure different values of ( ), and for some range of , log ( ) increases roughly linearly with − log . We can use that property to obtain the fractal dimension , given by (Federrath et al. 2009): Unlike , the fractal dimension obtained with the boxcounting method depends on the pixel threshold. Figure 15 shows the fit for the HNCO emission map of the TVir snapshot. In the top panel of the figure we see fitted for a chosen pixel threshold, while in the bottom panel we see the change of as a function of pixel threshold. In order to compare between the different emission maps, we have fixed the threshold to the mean pixel value in the image (excluding the background pixels). Typically this threshold is around ∼ 10 − 20% of the maximum pixel value. The results of these calculations for the emission maps of all molecular species are shown in Table 6.

Fractal structure
We find three different ranges of fractal dimension values for our synthetic emission maps, depending on the method used. The two different interpretations of the perimeter-area method yield = 1.6 − 1.9 for the addition approach and  Figure 15. Fractal dimension, , obtained with the box-counting method (see eq. 13) for the HNCO emission map of the TVir snapshot. Top: calculated at a pixel threshold given by the mean pixel value of the emission map (≈ 15% of the maximum pixel value). Bottom: calculated for a range of pixel thresholds, which are given as a fraction of the maximum pixel count in the image. The shaded area gives a 1 error bar to the slope of the linear fit for . for the scatter approach. There is one outlier of the scatter approach, which is the HNC emission map of the TVir snapshot, with = 1.11. Upon visual inspection of the image, we see that it consists of individual peaks, with close to circular contours, which explains the results. It is interesting to note that this image yields a very high fractal dimension of 1.84 with the addition interpretation, which is an example of why we do not recommend using this approach. Additionally, we find the fractal dimension range from the box-counting method to be = 1.3 − 1.6 (excluding the emission maps of N 2 H + , which are noisy and hence the level at which is calculated is closer to the peak brightness). There is no significant variation in fractal dimension between the TVir and SVir emission maps. In comparison, Rathborne et al. (2015) report a slightly lower fractal dimension = 1.4 − 1.7 using the addition interpretation, with values between 1.4 and 1.5 for the four molecular tracers that overlap with our data set. However, we find inconsistent fractal dimension values for HCO + and H 2 CS when we perform the measurement using the addition interpretation. This discrepancy is due to the fact that the addition interpretation is sensitive to the selection of structures within the image, and we exclude structures with areas that are smaller than four beam sizes. When using the scatter interpretation on the Rathborne et al. (2015) data, we find values within the range of = 1.2 − 1.3, which are consistent with the synthetic maps.
The two distinct ranges of tell different stories about the fractal nature of the Brick. While the addition interpretation sees structures with relatively high level of complexity, the scatter interpretation sees a smoother cloud. From a methods perspective, we trust the results of the scatter interpretation, as discussed in Appendix C. The values for interstellar gas clouds, obtained by other authors using this method, fall in the range of = 1.2−1.6 (Beech 1987;Bazell & Desert 1988;Dickman et al. 1990;Falgarone et al. 1991;Vogelaar et al. 1991;Hetem & Lepine 1993;Vogelaar & Wakker 1994;Sánchez et al. 2005). This places the fractal dimension of both the observed and simulated emission maps of the Brick in the lower half of that range.
The fractal dimension, given by the box-counting method, is harder to compare to literature values, since it is a function of the pixel threshold. Figure 15 shows that as we increase the pixel threshold, we gradually lower from about 1.6, down to about 1.4. This trend arises from the fact that by selecting only the brightest peaks, we consider small spatial scales, where we lack the sufficient resolution to see fractal substructure. Therefore, the highest value of serves as an upper estimate of the fractal nature of the synthetic emission maps.
All of the above results indicate that neither the observed, nor the simulated Brick cloud have unusual fractal structure, compared to typical interstellar gas clouds.

Moments of inertia characterised with the J-plots method
We now quantify the cloud substructure by using the principal moments of inertia, characterised through the J-plots method (Jaffa et al. 2018). With this method, we can distinguish between centrallycondensed (e.g. cores), centrally-rarefied (e.g. rings), and elongated structures (e.g. filaments).
In order to apply the J-plots method to the emission maps generated with casa, we first identify 2D substructures in the maps by constructing dendrograms with the astrodendro Python package 5 .
Each dendrogram structure is a contiguous group of pixels with values above a given threshold (similarly to the structures used for measuring the fractal dimension). The dendrogram is constructed starting from the brightest pixel and gradually decreasing the threshold to identify more structures and their hierarchical connections. The dendrogram construction is controlled by three parameters, which are the minimum pixel value (lowest possible threshold), the minimum delta (minimum threshold distance between two hierarchically connected structures), and the minimum number of pixels in a structure. For our analysis, we set the minimum value to three times the root-mean-squared (RMS) noise of the image background. Note that the RMS is computed for an empty region at the periphery of each image. The minimum number of pixels must exceed the beam size, as the latter defines our resolution. For the synthetic observations we have used 200 as the minimum number of pixels (with the beam area covering ≈ 25 pixels), and for the real observations we have used 65 pixels. This accounts for differences in the pixel size between the two sets of images, and it ensures that the minimum number of pixels cover approximately the same physical area. The minimum delta is set to 10% of the brightest peak.
After constructing the dendrogram, we compute the J moments of each structure in the following way, as described by Jaffa et al. (2018). Let a structure consist of pixels, with values Σ , where ∈ [1, ]. The J plot method treats Σ as surface density, and hence uses the terms "mass" and "centre of mass". We will adhere to this language for consistency with Jaffa et al. (2018), but note that in our calculations Σ is the intensity.
The area, , and mass, , of the structure described above are given respectively as: where Δ is the area of a pixel. We can then compute the moments 0 = , , , , , , using the position of each pixel center, ( , ), 6 which leads to the center of mass: The moments about the Cartesian axes are The principal moments are then given by where the indices {1, 2} refer to the different signs (∓). From this definition, the J moments are defined for = 1, 2: where 0 = /4 . Note that we always have 1 ≥ 2 . After computing the J moments, we plot 2 against 1 and we interpret the J plot in the following way. The position of a datapoint projected along the 1 = 2 diagonal determines the degree of central condensation/rarefaction of the structure, with negative values of 1 and 2 corresponding to a centrally-rarefied, and positive values of 1 and 2 corresponding to a centrally-condensed structure. Additionally, the projected position along the 2 = − 1 diagonal indicates the aspect ratio of the structure, with 2 = 1 = 0 being a circle, and 1 > 0, 2 < 0 being an elongated structure. In this context, a curved filament would count as a (somewhat) centrally-rarefied structure. Figure 16 shows the J plots of the synthetic (TVir: top panel; SVir: middle panel) and real (Rathborne et al. 2015: bottom panel) integrated HNCO emission maps of the Brick. The size of each data point is proportional to the physical area of the corresponding structure, and the colour indicates the pixel threshold at which the structure has been identified within the dendrogram. The structures in Figure 16 occupy the area of the J plot around the 2 = − 1 diagonal where 1 > ∼ 0 and 2 < ∼ 0. The TVir and SVir structures are approximately evenly distributed on both sides of the diagonal, indicating a balance between centrally condensed and rarified structures. Interestingly, the centrally-condensed structures ( 2 > − 1 ) are predominantly close to circular, whereas the centrally-rarefied structures ( 2 < − 1 ) are predominantly elongated. By contrast, the structures in the Rathborne et al. (2015) observations are almost entirely below the diagonal, indicating that the real Brick is lacking centrally-condensed structures. None of the three data sets exhibit a preferred elongation of the structures, because the orthogonal projections of the data points on the 2 = − 1 diagonal are approximately uniformly distributed across the spanned range of values. Finally, there is no apparent correlation between the size or the colour of a data point and its position on the J plot of the simulated emission maps. In the J plot of the real Brick we see that the largest data points correspond to structures that are elongated and slightly centrally-rarefied. This is due to the fact that on a cloud scale, the Brick is distinctly bean-shaped.
The above described difference in the centrally-condensed regions of the simulated and the real Brick also persists for the other molecular tracers. This difference may be understood in terms of the ongoing process of star formation. The simulations are actively undergoing fragmentation and gravitational collapse, and by the time of the selected snapshots, close to half of the original gas mass has been transformed into sink particles, which are fed from condensations of locally gravitationally-bound gas. These gravity-dominated structures appear as the centrally-condensed data points in the J plot. The reason why we do not see centrally-condensed structures in the J plot of the real Brick is likely because its gravity-dominated regions are less numerous (since the Brick appears primarily quiescent) and smaller than the simulated ones. Indeed, Walker et al. (2021) have detected centrally-condensed star-forming cores in the Brick, with sizes of ∼ 1000 AU, but only few compared to the simulations. In order for these cores to be visible on the J plot, we require a much higher resolution than that of the observations of Rathborne et al. (2015).
The J plot results of Figure 16 reflect differences in the star formation rate and star formation size scale between the simulations and the observations. We have already discussed the differences in length scales relevant for star formation in Section 4.2, so we restrict the brief discussion here to the star formation rate. Federrath et al. (2016) predicted the star formation rate per free-fall time of the Brick to be ff = 0.042 ± 0.030, using the star formation model of Krumholz & McKee (2005) with the parameter fits from . A later measurement by Barnes et al. (2017) found ff = 0.02, which is in agreement with the theoretical prediction. We use the same theoretical models as Barnes et al. (2017) to predict the the star formation rate per free-fall time of the TVir simulation for three different star formation theories. By using our values for M, and (see Section 4.1), and the virial parameter of the simulated cloud at the analysed snapshot ( vir = 1.24, Dale et al. 2019), we obtain ff,TVir = {0.0021, 0.034, 0.298} for the models of Krumholz & McKee (2005), Padoan & Nordlund (2011) and Hennebelle & Chabrier (2013), respectively. This gives us a large spread of possible values due to differences in the star formation theories. In reality, the TVir simulation is actively forming stars, with ff = 0.39 (Dale et al. 2019), and we need to slow this process down in order to match the observations. One way to achieve this is by including magnetic fields. Indeed, if we assume the same turbulent magnetic field as is present in the real Brick (with = 0.34), the predicted values for ff,TVir for the three star formation theories become considerably smaller, with ff,TVir = {0.0001, 0.023, 0.004}. This leads us to the prediction that a magneto-hydrodynamic simulation of our Brick model should better reproduce both the observed star formation rate and the multi-scale structure of the real Brick quantified in Figure 16.

CONCLUSIONS
We have presented a detailed study of the complex, multi-scale structure of numerical simulations and ALMA observations of the CMZ cloud G0.253+0.016 ('the Brick'). To do so, we have created synthetic emission maps of the Brick, using snapshots of two hydrodynamics simulations (Dale et al. 2019;Kruĳssen et al. 2019), which have been post-processed with the line radiative transfer code polaris and subsequently with casa. The chosen snapshots match the current position of the Brick along its orbit. We have adopted the observational setup of Rathborne et al. (2015) for our postprocessing to enable a direct comparison between their observations of the real Brick and our numerically modelled ones.
We have produced emission maps of eight molecular species per simulation snapshot (TVir and SVir), tracing gas of different densities. We find that half of our molecules (CS, HCN, HCO + and HNC) produced emission maps with compact structures, while the other half (HNCO, H 2 CS, HC 3 N and N 2 H + ) traced more diffuse gas. Furthermore, all of our molecular emission maps show at least moderate correlation with the gas surface density. This is in contrast to the real Brick, where Rathborne et al. (2015) found poor correlation between their molecular emission maps and the dust continuum, which is assumed to trace the column density. This discrepancy between the simulations and the observations is likely due to the simplified assumptions of constant temperature and abundances in our modelling, as well as due to potential structural differences between the simulations and the real Brick.
We have studied the structure of our simulated emission maps by constructing brightness and density PDFs and power spectra, and by calculating the fractal dimensions of substructures and their moments of inertia through the J-plots method. The density PDFs of the simulations have a larger width than those of the real Brick, which is likely explained by the lack of magnetic fields in the simulations. Indeed, assuming the observed value of the turbulent plasma beta ( = 0.34) in our simulations compensates for the difference in density dispersion.
We have constructed the spatial power spectra of the column density and the synthetic emission maps to study the turbulent structure of the cloud. We find that the column density power spectrum has a break at ∼ 0.32 pc, which is consistent with the plausible range for the Jeans length in the simulations (0.1−0.5 pc). The large length scale slope of the column density power spectrum is similar to, but slightly lower than the value obtained by Rathborne et al. (2015). The slightly shallower slope may be caused by the lack of magnetic fields in the hydrodynamics simulation. The power spectra of the synthetic emission maps follow a single power law within their reliable fitting range of spatial scales. We find slopes in the range between −4.2 and −2.3, which is approximately consistent with the findings in Rathborne et al. (2015), who reported slopes of −3.
We have studied the fractal structure of the synthetic emission maps, using the perimeter-area method and the box-counting method. We find good agreement between the synthetic ( = 1.2 − 1.4) and the observed emission maps ( = 1.2 − 1.3), when using the perimeter-area method. We obtain ≈ 1.3 − 1.6 with the box-counting method, applied to the synthetic emission maps. These results indicate that the Brick and its simulated counterpart have typical fractal structure for an interstellar gas cloud.
Additionally, we report that we have identified two interpretations of the perimeter-area method in the literature. In one of them individual connected regions within a contour are treated as separate entities, while in the other their areas and perimeters are added up. We advise against the interpretation involving addition, as discussed in Appendix C.
Finally, we have studied the moments of inertia of individual cloud substructures, by applying the J-plots method (Jaffa et al. 2018). Through this method we can distinguish between centrallycondensed (cores), centrally-rarefied (rings) and elongated structures (filaments). We find that the synthetic emission maps contain a mixture of centrally-condensed and centrally-rarefied structures with varying degrees of elongation. In contrast, the real Brick consists of almost exclusively centrally-rarefied structures. We attribute these differences to the ongoing process of star formation and in differences in the star formation rate between the simulations and the real Brick.
Throughout our analysis, we have seen consistent similarities between the TVir and SVir snapshots, even though Dale et al. (2019) find that the tidally-virialised simulations are better models of the CMZ clouds than the self-virialised ones, based on the more rapid gravitational collapse of the latter. This indicates that the small scale properties of the gas are not strongly influenced by the global amount of initial kinematic support of the cloud (at least at the orbital position of the Brick), and it makes our results more robust against variations in the initial conditions of the clouds.
Our results demonstrate that the underlying Brick simulations reproduce some important structural properties of the cloud (in addition to the integrated properties studied by Kruĳssen et al. 2019), such as the hierarchical substructure and the fractal dimension. This implies that many aspects of the cloud's evolution are affected by its orbit through the gravitational potential in the Galactic Centre region. In order to reproduce other observables, such as the power spectrum and the weakly correlated emission maps, we require more complex models. These may include a more thorough treatment of the chemistry, and the addition of magnetic fields in the simulations. Models that include the external gravitational potential of the Galactic Centre take the first necessary step towards our understanding of the complex interplay between turbulence and self-gravity within this region.

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

APPENDIX A: THE INTEGRAL OF A CUBIC SPLINE KERNEL
In order to accommodate for the irregular shape of a Voronoi cell, the integration volume is divided into pyramids, as shown in Figure A1. First the vertices of each cell wall are connected to the particle position, creating a wall pyramid, and then each wall pyramid is divided into vertex pyramids. Each vertex pyramid is characterised by the distance from the particle to the wall, 0 , the distance from the orthogonal projection of the particle position on the plane of the wall to an edge of the wall, 0 , and an angle , as shown in Figure A1.
The integral of over the volume of a vertex pyramid can be written as: In the above = 0 0 , = √︁ 1 − (1 + 2 ) 2 , and = cos . Even though the use of the analytic integral greatly speeds up the density and velocity mapping, there are still multiple expensive calculations that need to be performed. For this reason, we have precomputed and tabulated the values of for various configurations of 0 , 0 and . Any desired value of is then obtained by linearly interpolating between the pre-computed values. The tabulated approach introduces a small error in the mapped SPH parameters, which has been has been shown to be within 0.6% (Petkova 2018).

APPENDIX B: MOMENT MAPS
In addition to the other results, we include the first and second velocity moment maps of HNCO for the two simulation snapshots (see Figure B1 and B2 respectively). In order to construct these maps, we combine the velocity channel data that is directly produced by polaris (without further processing with CASA). Figure B1 shows that both snapshots present clear signs of rotation, with a very similar velocity range. These velocities are an excellent match to the real Brick, which has line-of-sight velocities between ∼ 0 km s −1 and ∼ 45 km s −1 (Rathborne et al. 2015;Federrath et al. 2016;Henshaw et al. 2019). The velocity dispersions in Figure B2 are also similar between the two snapshots, with values up to ∼ 20 km/s for a pixel size of 0.04 pc.

APPENDIX C: INTERPRETATIONS OF THE PERIMETER-AREA METHOD
Let us consider a single, finite, spherically-symmetric blob of gas. Any chosen contour of the blob is perfectly circular, and hence ∝ 1/2 . If we are to measure the fractal dimension using the perimeter-area method, we will obtain the same data set for both interpretations of the method (see Section 4.3.1), since we have a single blob. This will result in a measured fractal dimension of = 1.
If we now introduce a second blob, identical to the first one and non-overlapping with it, the scatter interpretation will yield = 1. The addition interpretation will have a different set of data points {(2 1 , 2 1 ), (2 2 , 2 2 ), ..., (2 , 2 )}, where and are the areas and perimeters of different contours of one of the two blobs. This data can also be fitted with a line of slope 0.5, and a fractal dimension of 1, however the constant term of the liner fit in loglog space will be different than the one in the scatter interpretation (see Figure C1). This agreement between the two methods can be generalised for a cloud consisting of N identical, non-overlapping blobs.
Finally, if we consider two spherically-symmetric but nonidentical blobs, the scatter interpretation once again finds = 1. However, the data points of the addition interpretation are no longer simple multiples of the data points of the scatter interpretation, as they were in the previous case. This leads to a disagreement between the results of the two methods (see Figure C2). Note that in the addition plot we can identify two linear segments with a slope of 0.5, but the liner fit through all of the data points has a steeper slope.
The above two examples both serve to ask a simple question. What is the fractal dimension of a cloud which consists of two (many) independent regions, both (all) of which have the same individual fractal dimension? The scatter interpretation deduces that the fractal dimension of the cloud is always the same as the (identical) fractal dimensions of its substructures. The addition interpretation is inconclusive. In our two examples we get a fractal dimension of either 1 or 1.3, depending on the relative sizes of the cloud substructures. Due to this behaviour of the addition interpretation, we recommend the use of the scatter interpretation in future studies. This paper has been typeset from a T E X/L A T E X file prepared by the author. Each datapoint corresponds to an identified structure within the corresponding dendrogram, and its size is proportional to the area of the structure. The colour indicates the intensity level at which each structure is identified in Jy/beam km/s.   Figure C2. Two non-identical, non-overlapping, spherically-symmetric blobs, which yield = 1 in the scatter interpretation and = 1.3 in the addition interpretation. The blobs were generated using the cubic spline kernel function shown in equation 2, and their sizes are indicated by the dashed circles.