2 Fractal Cubes

2.1 General Properties

We refer to “fractal cubes” as randomly generated scalar fields in three-dimensional cartesian space, \(\rho(x_1, x_2, x_3)\), from some probability distribution function (PDF), \(P(\rho)\), and fractal correlation in space. The two-point fractal distribution is established by a power-law spectrum \(D(k)\) in Fourier space, and characterized by the slope of the power-law, \(\beta\), the Nyquist limit \(k_\mathrm{max}\), and a lower cutoff wavenumber, \(k_\mathrm{min}\).

Let \(F(\mathbf{k})\) be the Fourier transform of the spatial density distribution, \(\rho(\mathbf{r})\), with \(\mathbf{k}\) and \(\mathbf{r}\) as wave vector and position vector, respectively.

\[D(k) = \int k^2 F(\mathbf{k}) F^{*}(\mathbf{k}) d \Omega \propto k^\beta\;,\]

where the integral of the spectral density, \(F(\mathbf{k}) F^{*}(\mathbf{k})\), is over all solid angle, \(\Omega\).

The target spectrum is generally chosen to have a slope representative of the density or velocity field of turbulent flow, e.g. a Kolmogorov spectrum \(\beta=-5/3\), which is the default in the classes.

It is sometimes useful to restrict the largest or smallest spatially correlated scales contained in the fractal cube. The lower cutoff wave number restricts sampling, and therefore fractal structure, to scales smaller than \(L/k_\mathrm{min}\), where \(L\) is the width of the cube. This is a useful parameter in controlling the size scale (or column depth) of structure in the fractal cube. The Nyquist limit, \(k_\mathrm{max,i} = (n_i + 1)/2\), sets the smallest correlated scales, and is used as the upper sampling wavenumber limit.

The main functionality of the pyFC module is to generate a lognormal fractal cube for a scalar field \(\rho\), but there is also a class that represents gaussian fractal cubes. Note that, although we refer to fractal cubes as cubes, they may in fact be cuboidal, that is \(n_i \ne{} n_j \ne{} n_k\), where these are the number of cells in the \(i\), \(j\), and \(k\) directions of the fractal cube.

Since fractal cubes are generated purely by manipulations in the fourier space of the scalar field, the fractal cube is by construction periodic and can be trivially tesselated.

2.2 Lognormal Fractal Cubes

Lognormal fractal cubes obey a single-point lognormal PDF characterized by the mean, \(\mu\), and variance, \(\sigma^2\) of the distribution.

\[P(\rho) = \frac{1}{s\sqrt{2\pi}\rho}\exp\left(\frac{-(\ln\rho - m)^2}{2s^2}\right)\;.\]

2.2.1 Construction

A lognormal fractal cube is generated by an iterative procedure that attempts to simultaneously ensure that the scalar field obeys a single-point lognormal PDF and a two-point fractal distribution, that is, a power-law in Fourer space. Technically, the aim of the iterative procedure is to find the apodization spectrum for the Gaussian random field, which gives a lognormal fractal cube. The steps of the iteration procedure is as follows:

1. Construct in real space a cube of random variables from a Gaussian PDF, with mean, \(m\), and variance, \(s^2\).

2. Use a Fast Fourier Transform (FFT) to obtain the power spectrum of the Gaussian field.

3. Apodize (multiply) the power spectrum with the target power-law power spectrum, which is first normalized so that the mean and variance of the field is preserved.

4. Use an inverse Fast Fourier Transform (iFFT) to obtain a spatially correlated random Gaussian field.

5. Exponentiate the cube to obtain a spatially correlated random lognormal field.

6. FFT of the lognormal field shows, however, that its power spectrum is not a power-law. The deviation from a power-law is calculated and a correction is applied to the apodization spectrum.

7. The procedure from 3. is repeated and successive corrections are applied until the lognormal field attains a power-law spectrum whose deviation from the target spectrum is within an acceptable tolerance.

The mean, \(m\), and variance, \(s^2\) of the underlying Gaussian distribution are related to the mean, \(\mu\), and variance \(\sigma^2\) of the lognormal distribution by:

\[m = \ln\frac{\mu^2}{\sqrt{\sigma^2 + \mu^2}} \;, \quad s = \sqrt{\ln\left(\frac{\sigma^2}{\mu^2} + 1\right)}\,.\]

The process of convergence in the generation of a lognormal fractal cube is shown in the figure below for successive iterations. About 3 to 7 iterations are required. The left panel shows a midplane slice of the lognormal field, the central panel shows the PDF of \(\log_{10}(\rho)\), which is Gaussian for a lognormal field, and the right panel shows the isotropic power spectrum of the gaussian field, which is apodized with the succesively corrected power spectrum, the isotropic spectrum of the lognormal field (which is simply the exponentiated Gaussian field), and the target power-law spectrum (black solid line). The convergence is measured with respect to the power-spectrum. In the process of modifying the apodization spectrum, although it remains normalized, the single-point PDF may deviate slightly from the target spectrum.

Left: Midplane slice, Middle: single-point lognormal PDF, and Right: power-law isotropic power spectrum.

2.3 Gaussian Fractal Cubes

A class and functions also exist to generate fractal cubes with single-point Gaussian statistics, described by the mean and variance. These are useful to generate velocity dispersions.

Future work includes combining lognormal and Gaussian fractal cubes in a consistent way to be used in hydrodynamic simulations.

2.4 Implementation in hydrodynamic simulations and generating porosity

The fractal cube can be included in hydrodynamical simulations to represent density inhomogeneities. In astrophysical simulations, these are typically dense atomic or molecular clouds in the interstellar medium, which are embedded in a more tenous hotter medium. Fractal cubes are thus useful to create two or multi-phase gaseous initial conditions for simulations.

To construct a dense cloud distribution with the fractal cube data for hydrodynamical simulations, it is often useful to set the mean of the fractal cube \(\mu=1\). The fractal cube data can then be apodized (multiplied) by the desired mean density distribution of the clouds, \(\bar{\rho}_\mathrm{clouds}\). This could, for example, be a quasi spherically symmetric and isotropic gas distribution, or a galactic disc. Once apodized, the inhomogeneities will usually have a physical lengthscale, or scaling associated with them, that is defined by the relative sizes of the simulation box and of \(1/k_\mathrm{min}\).

The porosity of the warm phase arises by imposing a lower density threshold for the clouds. Individual clouds with the associated scaling emerge. In a two-phase medium, the lowest realistic density threshold is the density of the tenuous hot background gas, but one may choose an arbritrary lower density threshold, below which the simulation cell is set to contain hot background values. Across phase transitions, if we consider a quasi-static system, the total pressures of neighboring phases are usually in equilibrium, and a lower density threshold may be prescribed according to a pressure matching condition.

2.4.1 Simple example: Clouds in the interstellar medium

Let us consider the construction of a two-phase interstellar medium (ISM) using the fractal cube data for the dense, colder phase (correctly, known as the “warm” phase), embedded in a uniform, tenous, hot phase (e.g. as in Sutherland & Bicknell 2007). For simplicity, let us assume both phases are in pressure equilibrium, both are described by an ideal equation of state, and have the composition with the same mean mass per particle, \(\mu\).

We generate porosity by imposing an upper temperature cutoff for the existence of clouds at \(T_\mathrm{crit}=3\times10^4\,\mathrm{K}\), beyond which clouds are deemed thermally unstable. The upper temperature cutoff corresponds directly to a lower density cutoff, \(\rho_\mathrm{crit}=\mu p/(k T_\mathrm{crit})\), if the pressure, \(p\), is defined. Since the clouds are in pressure equilibrium with the surrounding hot phase, \(\rho_\mathrm{crit}=\rho_\mathrm{hot} T_\mathrm{hot}/T_\mathrm{crit}\), where \(\rho_\mathrm{hot}\) and \(T_\mathrm{hot}\) are the hot phase density and temperature, respectively. The volume filling factor of the warm phase is:

\[\begin{split}f_\mathrm{vol}&=&\int_{\rho_\mathrm{crit}}^\infty P(\rho) d \rho \\ &=&\frac{1}{2}\left[1 + \mathrm{erf}\,\left(\frac{\ln\left\{(\rho_\mathrm{crit}/\mu)\sqrt{\sigma^2/\mu^2 + 1}\right\}}{\sqrt{2\ln\left(\sigma^2/\mu^2 + 1\right)}} \right)\right]\end{split}\]

Note that after imposing a lower density threshold for the dense phase, the mean of the PDF becomes higher and \(\bar{\rho}_\mathrm{clouds}\) is effectively changed.

Future work includes formulating a prescription for the velocity dispersion of the clouds, which should have the same spatial correlation as a given density distribution.