Simulating and Characterizing Non-Gaussian Small-Scale Polarized Emission

Introduction

Detection of the inflationary $B$-mode signal remains a primary goal for confirming the inflationary model and constraining the tensor properties of the Universe ((M. Kamionkowski et al. 1997; W. Hu & M. J. White 1997; M. Kamionkowski & E. D. Kovetz 2016)). However, this cosmological signal is heavily obscured by Galactic foregrounds. Synchrotron emission dominates the $B$-mode signal at low (radio) frequencies, while thermal dust emission dominates at higher (millimeter/submillimeter) frequencies ((Planck Collaboration I 2020)). Accurately characterizing and removing these foregrounds is therefore critical for unbiased cosmological parameter estimation (see, e.g., (N. Krachmalnicoff et al. 2016; K. Wolz et al. 2024; F. Bianchini et al. 2025)).

While the large-scale foreground emission is relatively well-characterized, small-scale emission remains largely unconstrained by full-sky data. The observational constraints that do exist reveal that the underlying interstellar medium (ISM) and plane-of-sky magnetic field topology are highly complex and structured. For instance, small-scale synchrotron intensity in the Galactic center traces intricate, filamentary magnetic features that are distinctly non-Gaussian ((I. Heywood et al. 2022; S. Vandervoude et al. 2024)). Similarly, magnetic field lines inferred from dust emission exhibit a similarly non-Gaussian, filamentary morphology ((Planck Collaboration Int. XVIII 2015)).

Because these emissions are intrinsically linked to Galactic magnetic field (e.g., G. B. Rybicki & A. P. Lightman 1979; A. Lazarian 2007), physical 3D simulations are often used to generate mock full-sky observations (e.g., A. Waelkens et al. 2009; X. H. Sun & W. Reich 2009; T. Steininger et al. 2018; J. Wang et al. 2020; T. A. Porter et al. 2022). Even these models remain limited in their ability to resolve small-scale fluctuations, thus, synthetic fields can be employed as a surrogate ((J. Delabrouille et al. 2013; B. Thorne et al. 2017; F. Vansyngel et al. 2017; C. Hervías-Caimapo & K. Huf- fenberger 2022; The Pan-Experiment Galactic Science Group et al. 2025)). Standard Gaussian approximations overlook the intrinsic non-Gaussianity of the magnetized ISM (e.g., D. Biskamp 2003; A. G. Kritsuk et al. 2007; G. Kowal et al. 2007) and current stochastic foreground simulations cannot reproduce these filamentary features. There is therefore no existing method that generates @page 1/25 starts_mid=0 ends_mid=1

Introduction

Detection of the inflationary $B$-mode signal remains a primary goal for confirming the inflationary model and constraining the tensor properties of the Universe ((M. Kamionkowski et al. 1997; W. Hu & M. J. White 1997; M. Kamionkowski & E. D. Kovetz 2016)). However, this cosmological signal is heavily obscured by Galactic foregrounds. Synchrotron emission dominates the $B$-mode signal at low (radio) frequencies, while thermal dust emission dominates at higher (millimeter/submillimeter) frequencies ((Planck Collaboration I 2020)). Accurately characterizing and removing these foregrounds is therefore critical for unbiased cosmological parameter estimation (see, e.g., (N. Krachmalnicoff et al. 2016; K. Wolz et al. 2024; F. Bianchini et al. 2025)).

While the large-scale foreground emission is relatively well-characterized, small-scale emission remains largely unconstrained by full-sky data. The observational constraints that do exist reveal that the underlying interstellar medium (ISM) and plane-of-sky magnetic field topology are highly complex and structured. For instance, small-scale synchrotron intensity in the Galactic center traces intricate, filamentary magnetic features that are distinctly non-Gaussian ((I. Heywood et al. 2022; S. Vanderwoude et al. 2024)). Similarly, magnetic field lines inferred from dust emission exhibit a similarly non-Gaussian, filamentary morphology ((Planck Collaboration Int. XVIII 2015)).

Because these emissions are intrinsically linked to Galactic magnetic field (e.g., G. B. Rybicki & A. P. Lightman 1979; A. Lazarian 2007), physical 3D simulations are often used to generate mock full-sky observations (e.g., A. Waelkens et al. 2009; X. H. Sun & W. Reich 2009; T. Steininger et al. 2018; J. Wang et al. 2020; T. A. Porter et al. 2022). Even these models remain limited in their ability to resolve small-scale fluctuations, thus, synthetic fields can be employed as a surrogate ((J. Delabrouille et al. 2013; B. Thorne et al. 2017; F. Vansyngel et al. 2017; C. Hervías-Caimapo & K. Hufenberger 2022; The Pan-Experiment Galactic Science Group et al. 2025)). Standard Gaussian approximations overlook the intrinsic non-Gaussianity of the magnetized ISM (e.g., D. Biskamp 2003; A. G. Kritsuk et al. 2007; G. Kowal et al. 2007) and current stochastic foreground simulations cannot reproduce these filamentary features. There is therefore no existing method that generates many independent realizations while retaining realistic, non-Gaussian ISM structure.

In this work, we present a framework that fills this gap by generating stochastic, non-Gaussian foreground realizations that reflect the structural complexity expected of Galactic emission [4] (C. Hervías-Caimapo & K. Huffenberger 2022; C. Hervías-Caimapo et al. 2025; F. Rahman et al. 2026). The framework introduces anisotropic, multifractal, and advected non-Gaussian fluctuations, and has three properties of particular relevance to CMB foreground modelling. First, it is stochastic: it produces arbitrarily many independent realizations, enabling robust characterization of statistical uncertainties. Second, it inherits information from large angular scales, so that the small-scale structure is anchored to the large-scale emission. Third, it is tunable through a small number of physically interpretable parameters, which allows the model to be matched to new high-resolution observations as they become available and to stress-test and identify break points in CMB analysis pipelines for future missions. We use polarized synchrotron emission as an example use case, but the methods apply equally to other foregrounds such as thermal dust emission. This work lays the foundation for new foreground simulation suites with stochastic realizations of realistic, non-Gaussian structure.

Another central challenge is how to quantify the non-Gaussianity (see, e.g., (B. Burkhart et al. 2009; F. Rahman et al. 2026)). Because the power spectrum (a two-point statistic) is degenerate with respect to various spatial features, it cannot distinguish between a Gaussian random field and a non-Gaussian random field with identical power-law slopes but different phase correlations (e.g., (L.-Y. Chiang & P. Coles 2000; T. Matsubara 2003; D. Falceta-Gonçalves et al. 2014)). In this work, our simulated small-scale fields are extrapolated from large-scale WMAP constraints (C. L. Bennett et al. 2003) and assume power-law correlations (B. Thorne et al. 2017; The Pan-Experiment Galactic Science Group et al. 2025). However, to evaluate the complex phase coherence and morphology introduced by our framework, we must utilize higher-order statistics. Specifically, we evaluate the synthetic fields using Minkowski functionals—which provide an algebraic description of the geometrical and topological properties of a level set (F. Rahman et al. 2024, 2026)—and structure functions, which probe spatial intermittency and scaling by analyzing the moments of scale-dependent probability density functions (A. S. Monin & A. M. Yaglom 1975; D. Russo & W. A. Jury 1987; J. Jewell 2001; P. Padoan et al. 2003; B. Burkhart 2021). By analyzing our synthetic fields with these tools, we demonstrate their sensitivity to varying degrees of non-Gaussian structural complexity.

This paper is organized as follows. In § 3 we describe the generation of small-scale fluctuation fields with anisotropy (§ 3.1) and multifractality (§ 3.2). In § 3.3 we introduce a method to generate coherent filamentary structures. Next, in § 4 we characterize the resulting maps using metrics such as 1D PDFs, Minkowski functionals, structure function moments, and scale filtering. We then extend this framework in § 5 to construct full-sky synchrotron emission maps by incorporating large-scale emission and modulation fields for Galactic inhomogeneities. Finally, in § 6, we conclude with a summary.

The Small-Scale Models

The s6 synchrotron model is available in the PySM{^4} python package. Its synthesis is described in detail in The Pan-Experiment Galactic Science Group et al. (2025), but we will recap some important aspects in this section.

The polarization fraction tensor quantities $i$, $q$, $u$ are used instead of the Stokes $I$, $Q$, $U$ to synthesize non-Gaussian fields. The transformation from Stokes $I$, $Q$, $U$ to $i$, $q$, $u$ is given by

where $P = \sqrt{Q^2+U^2}$ is the usual polarized intensity. We adopt this polarization tensor approach for its compatibility with PySM, which we apply in § 5. Nevertheless, the methodology presented here can be readily adapted to raw small-scale $I$, $Q$, $U$ fields with minimal modification.

The components are generated in the polarization fraction tensor parameterization using

where $i_0$, $q_0$, $u_0$ is the large-scale field on the sky (given by e.g., the WMAP constraint C. L. Bennett et al. (2003)). The fields $i_s$, $q_s$, $u_s$ are the small-scale realizations where typical methodology generates isotropic, homogeneous Gaussian fluctuations. These small-scale fluctuations are the focus of this paper. The fields $i_{\mathrm{m}}$, $q_{\mathrm{m}}$, $u_{\mathrm{m}}$ are modulation fields representing inhomogeneous features on the sky. The result of $i\delta i_{\mathrm{m}}$, $q\delta q_{\mathrm{m}}$, $u\delta u_{\mathrm{m}}$ rebalances the power such that small sub-regions match observational constraints and are fields that have larger amplitude fluctuations where the modulation field is strongest, such as in the Galactic plane.

It is upon transformation back to the Stokes’ $I$, $Q$, $U$ parameters via

where $\xi=\sqrt{q^{2}+u^{2}}$, that small-scale non-Gaussian behavior is obtained.

Towards Filamentary Structures – Anisotropies

Our goal is to generate a field capable of stochastically producing structures with a strongly preferred orientation, thereby mimicking the extreme anisotropies revealed by observations of small-scale synchrotron emission (I. Heywood et al. 2022).

We are interested in generating fields $i_{\delta}$, $q_{\delta}$, $u_{\delta}$ that depart from isotropy/homogeneity. To achieve this, we apply the concepts of generalized scale invariance (see, e.g., S. Pecknold et al. 1997; G. Lewis et al. 1999; D. J. Schertzer et al. 2013). Recall that a scale invariant system needs a way of identifying (say, Fourier wave) vectors of the same scale (size) and it needs a way of measuring the scale of these vectors such that no characteristic scale is introduced to the system. Usually, the length of a vector is opted for, and therefore vectors of equal length are of the same scale. As analogues, in spherical harmonics, isotropic scaling only scales $\ell$ and treats all orders $m$ equally.

Generalized scale invariance moves away from these restrictive cases and states general conditions in which a system can be scale invariant. This is achieved by introducing a spectral-space scale changing operator that warps the harmonic coordinate. Rather than relying on the standard $L_{2}$ norm, generalized scale invariance defines an anisotropic, direction-dependent “effective-scale” metric, which stretches and rotates structures in the harmonic domain.

We define a continuous 2D vector $\ell$ in the flat-sky harmonic plane, where the vertical axis represents the traditional degree $\ell$ (multipole moment) and the horizontal axis represents the (azimuthal) order $m$:

In an isotropic system, the scale or magnitude of this vector is simply the $L_{2}$ norm:

To introduce directional dependence (anisotropy), we construct a transformation matrix $\mathbf{G}$. This matrix is the product of a scaling tensor $\mathbf{S}$, and a spatial rotation tensor $\mathbf{R}$:

where

We map the original isotropic harmonic vector $\ell$ to a warped coordinate vector $\ell'$ via matrix multiplication:

Finally, the “effective scale” is defined as the magnitude of this transformed vector. This represents the distance from the origin in the warped, anisotropic space:

which depends on $\ell,m$. By considering $\ell_{\mathrm{eff}}=(\mathrm{const.})$, a tilted and stretched ellipse in the $(\ell,m)$ plane is defined with the ratio of $s_{\ell}/s_{m}$ determining the overall level of contraction/stretch in the latitudinal and longitudinal directions, and $\theta_{0}$ determines the orientation of the anisotropy. Thus, $\mathbf{G}=\mathbf{S}(s_{\ell},s_{m})\mathbf{R}(\theta_{0})$ defines an anisotropic filter that alters the spherical harmonic coefficients depending on their orientation.

Small-scale anisotropic and self-similar Gaussian fields are then synthesized from harmonic-space following typical procedures with e.g.,

where $\delta_{\ell,m}^{G}$ is white variates whose real and imaginary components are independently Gaussian distributed, and $C_{\ell_{\mathrm{eff}}}^{\mathrm{target}}$ is the power spectrum of the fluctuations interpolated from a target power spectrum $C_{\ell}^{\mathrm{target}}$ onto the effective multipole $\ell_{\mathrm{eff}}(\ell,m)$ for each $\ell,m$. For exploration purposes, we will parameterize the target power spectrum as

with power-law index $\alpha$, and $L$ a large-scale cutoff/correlation scale.

Note that practically, Equation (16) can shift the distribution of power across the $\ell$ modes. Therefore, the estimated isotropic spectrum

typically differs from the target power spectrum $C_{\ell}^{\mathrm{target}}$.

To instead generate fluctuations which are anisotropic in real-space, but whose power is distributed equally across $m$ — so that the isotropic power spectrum equals the target power spectrum — synthesize the coefficients as follows:

where the total weight $W_{\ell,m}$ describes the power scaling for a specific $(\ell,m)$ mode. This weight is constructed from the product of an anisotropic and an isotropic component:

where the individual scaling factors are defined as

and $C_{\ell}^{\mathrm{target}}$ represents the target power spectrum. The term $C_{\ell}^{\mathrm{expected}}$ is the expected power introduced by the anisotropic factor, which is obtained by averaging over all $m$ modes:

For exploration purposes, we describe the following transformation matrices

where, as Figure 1 shows, $G_{\mathrm{isotropic}}$ is an identity matrix and returns $\ell_{\mathrm{eff}}=\ell_{\mathrm{isotropic}}$, i.e., the typical $\ell$’s forming circular magnitudes, $G_{\mathrm{meridional}}$ forms $\ell_{\mathrm{eff}}$ with elliptical magnitudes that are diagonal in the $\ell,m$ plane, and $G_{\mathrm{zonal}}$ is elliptical and stretched along the axis of $m$. When considering the true $\ell,m$ space, that for each $\ell> 0$, $m \le\ell$, and for symmetry $m \ge0$, $G_{\mathrm{meridional}}$ favours low-$m$ modes, and $G_{\mathrm{zonal}}$ favours high-$m$ modes.

Vector (PDF)

The resultant fluctuations are shown in Figure 2 with HEALPix² $N_{\mathrm{side}}=512$ (K. M. Górski et al. 2005), a cutoff scale $L=10$, and $\alpha=8/3$ (Equation (17)).³ As the names suggest, the field from $G_{\mathrm{isotropic}}$ consists of visually similar fluctuations, $G_{\mathrm{meridional}}$ generates visually similar fluctuations elongated North-South, and $G_{\mathrm{zonal}}$ generates visually similar fluctuations that are East-West elongated.

Vector (PDF)

An important factor to note is that these modifications generate inhomogeneities in the fluctuations. See, that the “meridional” field fluctuations are dominant in the Galactic plane with fluctuations at the (North-South) poles having less amplitude and different scale whereas the “zonal” field fluctuations are stronger at the North-South poles.

Towards non-Gaussianity – Multifractality

Features generated in Figure 2 are “visually” similar. They contain no abrupt fluctuations, singularities, or statistically anomalous regions of high intensity. In other words, we have been generating fields where the variates (i.e., $\delta^{\mathrm{G}}_{\ell,m}$ in Equation (16)) are statistically independent across all $\ell,m$. However, observations indicate the presence of significant levels of non-Gaussianity characterized by structural features whose intensity values lie deep within the non-Gaussian tails of the probability distribution (S. von Hausegger et al. 2019; Planck Collaboration XI 2020; Planck Collaboration IV 2020).

We now consider a field that is a cascade driven by non-linearities that couple different modes and effectively generate a $\delta_{\ell,m}$ that depends non-linearly on both scale and position (J.-F. Muzy 2019; J.-B. Durrive et al. 2022). In this scenario, rather than a uniform distribution of fluctuations, the field becomes highly intermittent, characterized by localized bursts, sharp spikes and sudden singularities (see, e.g., (Z.-S. She & E. Leveque 1994a; U. Frisch 1995; D. Biskamp 2003)). Given the complex processes that drive and regulate the density and magnetic fields responsible for photon emission, it is entirely justified to stochastically model synchrotron emission as such a cascade (B. G. Elmegreen & J. Scalo 2004; A. Lazarian & D. Pogosyan 2012).

We follow E. Bacry et al. (2001) and S. Lakhal et al. (2023) and generate a multifractal field. We describe the process below:

1. Generate a correlated noise field, $\omega(\hat{n})$, in harmonic space using: $a^{\omega}_{\ell,m} = \sqrt{C^{\mathrm{noise}}_{\ell}}\delta^{\mathrm{G}}_{\ell,m}$. Here, $\delta^{\mathrm{G}}_{\ell,m}$ is complex Gaussian noise and $C^{\mathrm{noise}}_{\ell}$ is the power spectrum of $\omega(\hat{n})$. Whilst $C^{\mathrm{noise}}_{\ell}$ does not determine the power spectrum of the resulting multifractal field, it does determine the strength and correlations of the intense regions.

At scales where $\omega$ is uncorrelated, then field will be monofractal. For example, if $C^{\mathrm{noise}}_{\ell}$ is the same as Equation (17), then the field is monofractal for scales $\ell< L$.

2. Generate log-normal fluctuations by exponentiating the noise field and sampling with noise in physical-space: $h(\hat{n}) = e^{\lambda\omega(\hat{n})}\delta^{\mathrm{v}}(\hat{n})$. By exponentiating the field $\omega(\hat{n})$, the fluctuations become positive and log-normally distributed. The parameter, $\lambda\in(0,1)$ controls the strength and scaling of the multifractal cascade with large $\lambda$ representing stronger non-Gaussianity and for $\lambda\to0$ the field will approach Gaussianity. The second noise, $\delta^{\mathrm{v}}(\hat{n})$, is sampled from a Gaussian

in real-space for every pixel with variance equal to the variance of $e^{\lambda\omega(\hat{n})}$.

3. Finally, synthesize the multifractal fluctuation field in harmonic-space: $a_{\ell,m} = \sqrt{C_{\ell}^{\mathrm{target}}}a_{\ell,m}^{h}$. Here, $a_{\ell,m}^{h}$ is the spherical harmonic coefficients of the field $h(\hat{n})$ that was generated in the previous step. The spectrum $C_{\ell}^{\mathrm{target}}$ describes the scaling of the field.

We show examples of multifractal fields (with $\lambda= 1/4$) in Figure 3. We generate these fields with the same $C_{\ell}$ parameterization from Equation (17) with HEALPix $N_{\mathrm{side}} = 512$ (K. M. Górski et al. 2005), $L = 10$, and $\alpha= 8/3$ for both the correlated noise field, $\omega$, and the final fluctuation field. We also include anisotropies following the procedure we discussed in § 3.1. The visual signature of the multifractal parameter ($\lambda= 1/4$) manifests as strong spatial intermittency. Unlike the random fields of Figure 2 which display a uniform texture, the multifractal fields exhibit heavily localized variance. There are distinct ‘hotspots’ of intense fluctuations interspersed with broad, relatively quiescent ‘voids’. The distribution and spatial structure of the ‘hotspots’/‘voids’ are determined by the correlations of $C_{\ell}^{\mathrm{noise}}$. Despite the change in distribution, the anisotropies from § 3.1 are maintained.

Vector (PDF)

As we will explore in § 4.2, the fields are non-Gaussian. Additionally, the probability density functions depart further from Gaussianity as $\lambda$ increases. Note that the multifractal formalism describes fields whose PDF depends on scale.

Building Coherent Structures – Lagrangian Advection

In Figure 3, the fluctuations still look like clumpy, disconnected clouds. Even though the power spectrum may be “correct”, and the higher-order moments successfully generate localized bright spots, is this purely anisotropic multifractal field truly what we would expect the small-scale fluctuations of the interstellar medium to look like? The answer is likely no. The model fails to reproduce thread-like filaments seen in, e.g., I. Heywood et al. (2022). The multifractal and anisotropic model lacks this complex spatial coherence. In a purely anisotropic Gaussian model, the phases of the spherical harmonics are random and independent. In the multifractal and anisotropic non-Gaussian model, the phases are indeed coupled — which is what generates the heavy-tailed intermittency and localized ‘clumps’ — but they lack the directional coherence needed to form elongated structures. To synthesize filaments, sweeping arcs, and

thread-like structures we must generate deeply correlated phases.

The ISM is a highly turbulent, magnetized plasma. In a highly conductive plasma, the magnetic field lines are “frozen-in” to the fluid. As turbulent eddies spin, shear, and cascade energy from large to small scales, they physically stretch and fold the magnetic field lines. It is the filaments that are the physical signatures of magnetic field lines being stretched and sheared by turbulent gas velocities. To emulate this procedure, we adapt the Lagrangian mapping approach (C. Rosales & C. Meneveau 2006; P. Subedi et al. 2014; J. Lübke et al. 2024). We generate two fields: an anisotropic multifractal field, and a turbulent velocity field which then advects the multifractal field along the velocity streamlines scale-by-scale. The advection process takes the clumpy, disconnected multifractal “clouds” and physically stretches and smears them into filaments. It naturally correlates phases across scales and conserves the heavy-tailed statistics.

We follow the Lagrangian mapping (LM) approach of J. Lübke et al. (2024) and generate a multifractal and geometrically complex random field. We describe the process below:

1. Generate a random scalar field: $A(\hat{n})$. Following the anisotropic and multifractal procedures detailed above, we synthesize a random scalar field $A(\hat{n})$.

2. Generate a random second field: $B(\hat{n})$. This second field can also be synthesized using the anisotropic and multifractal procedures.

3. Derive the turbulent velocity field, $\mathbf{u}(\hat{n})$. The spatial derivatives of the second field, $B(\hat{n})$, take the place of a turbulent velocity field on the sphere:

where $\mathbf{u} = (u_{\theta},u_{\phi})$ with global maximum $u_{\max} = \max\left\{\sqrt{u_{\theta}^{2}+u_{\phi}^{2}}\right\}$.

4. Generate coordinate displacements: $\Delta\theta(\hat{n})$, $\Delta\phi(\hat{n})$. We define the scale-dependent deformation timescale as

where $c$ is a free parameter and $\ell_i$ is a representative multipole moment. Using this timescale, the

velocity fields perturb the coordinates of $A(\hat{n})$ by

It is important to note that this step is scale-dependent, governed by the multipole moment $\ell_i$ in Equation (28). Under the minimal multiscale Lagrangian mapping approach (P. Subedi et al. 2014; J. Lü bke et al. 2024), large scales successively advect smaller scales. To implement this practically, the initial random scalar fields are decomposed using a sequence of band-pass filters from large to small scales. For each scale iteration $i$, the filter isolates a specific band of spherical harmonics characterized by a representative multipole moment $\ell_i$. This characteristic $\ell_i$ dictates the appropriate deformation timescale for that specific advection step.

In this step, $c$ is a free parameter that must be chosen carefully. Conceptually, $c$ represents the number of eddy turnover times that the “turbulence” is allowed to stir/distort the synthetic field at a given scale. In other words, it determines how long the artificial advection is allowed to run before the next scale. Increasing $c$ allows the field to become more twisted and distorted. Initially, this generates prominent coherent structures, though excessive distortion will eventually decorrelate the field entirely.

5. Advect the random scalar field.

Generate the LM field by advecting the random scalar field, $A(\hat{n})$, by sampling/interpolating from itself by the displacements of $\Delta\hat{n}=(\Delta\theta(\hat{n}),\Delta\phi(\hat{n}))$:

6. Correct the power of the displaced field.

Finally, correct the power of the displacement field in harmonic space:

where $W^{\mathrm{isotropic}}_{\ell,m}$ is as described in Equation (22).

Examples of advected fields are shown in Figure 4. In these examples, the base field, $A(\hat{n})$ is synthesized with: isotropic, meridional, and zonal fluctuations, and the advection field $B(\hat{n})$, is synthesized with only isotropic fluctuations. The advection parameters, for all fields are set to $c=5$ and $\lambda=1/4$. Visually, the influence of

Vector (PDF)

the Lagrangian advection ($c=5$) is immediately apparent in the form of sheared, filamentary structures that stretch across the fields. Rather than an isotropic distribution of peaks and troughs, the fields exhibit distinct directional streaking and whorl-like patterns indicative of a turbulent flow. The applied anisotropies further elongates these features along a preferred axis.

Polarization

We generate the linear polarization $q_\delta$, $u_\delta$ using the $e_\delta$, $b_\delta$ pseudo scalar fields (where $e$ and $b$ are the polarization fraction quantities analogous to $E$ and $B$ for $Q$ and $U$). For a covariance matrix described by:

the Cholesky decomposition ($\mathbf{A}=\mathbf{L}\mathbf{L}^{\mathsf{T}}$) is used to generate correlated noise (when non-zero cross-spectra are required). Generally, three independent fluctuation fields will be required, which are then scaled by combinations of the power-spectra described in $\mathbf{A}$ to produce the correct cross-correlations. These three noise fields can be synthesized using the above procedures thus generating polarization fluctuation fields with anisotropies, multifractality, and advected structures.

For a simple example, we provide power spectra as:

with target $C_\ell^{\mathrm{target}}$ given by Equation (17) with $\alpha=8/3$ and $L=10$. We synthesize polarization fields with multifractal parameter $\lambda=1/4$ and advection strength $c=5$. The resulting polarization fields: $e_\delta$, $b_\delta$, and polarized intensity $\xi_\delta=\sqrt{q_\delta^2+u_\delta^2}$ are shown in Figure 5. Comparing the fields, the higher relative power in the $E$-mode ($C_\ell^{ee}=3C_\ell^{\mathrm{target}}/4$) compared to the $B$-mode ($C_\ell^{bb}=2C_\ell^{\mathrm{target}}/5$) translates into visibly stronger, higher-contrast structures in $e_\delta$. The polarized intensity map ($\xi_\delta$) reveals a complex structure of highly polarized ridges and depolarization canals where $q_\delta$ and $u_\delta$ approach zero. The polarized intensity becomes small at the poles due to the mathematics of the meridional anisotropic fluctuations. In § 5 we will explore spatial combinations of maps to correct this artifact. Overall, the visual pattern of anisotropies, non-Gaussianity, and advected structures are displayed in the polarization fields.

Vector (PDF)

Summary

In summary, this section has introduced methods to generate complex morphological and statistical fields on the sphere for use in modelling cosmological foreground contamination. Synthesis takes $\sim2$ min at HEALPix $N_{\mathrm{side}}=1024$ (K. M. Górski et al. 2005) on a desktop computer.⁴ The resulting fields visibly demonstrate the effects of their underlying parameterization. Specifically, the transformation of $(\ell,m)$ by $\mathbb{G}$ (parameterized by $s_{\ell}$, $s_{m}$, $\theta_{0}$) introduces anisotropies that establish a global preference of spatial scale, allowing for the modelling of preferred directional alignments. Additionally, the multifractal nature ($\lambda$) breaks the spatial uniformity of standard Gaussian noise; this produces strong intermittency characterized by clustered high- and low-intensity regions governed by $C_{\ell}^{\mathrm{noise}}$. Furthermore, advection from the Lagrangian mapping (LM) generates coherent structures such as filaments, where the advection strength $c$ and the parameters for the flow field determine the spatial displacements. Although this framework introduces several new variables per field, such as $C_{\ell}^{\mathrm{noise}}$, $c$, and $\lambda$, alongside the advecting flow field parameters, in practical applications it is reasonable to assume these parameters are identical across the fields to simplify the model.

Figure 6 shows additional examples of small-scale temperature fluctuations, $i_{\delta}$, for varying multifractal and advection strength parameters: $\lambda$ and $c$, respectively. This grid serves as a visual map of this parameter space. Moving across the grid highlights how the spatial clustering driven by $\lambda$ and the coherent structural displacement driven by $c$ couple together to break Gaussian uniformity, providing a robust toolkit for simulating complex foregrounds.

Vector (PDF)

Caveats: We note the following caveats for those wishing to employ these methods:

• The anisotropies (§ 3.1) generate stronger amplitude fluctuations in the Galactic plane (for meridional models) and at the Galactic pole (for zonal models). If analyzing sufficiently small sky patches in these regions, the localized power can deviate from that of isotropic fluctuations.

• For multifractal fields (§ 3.2), the noise power-spectrum ($C_{\ell}^{\mathrm{noise}}$) controls the distribution and scaling of the intense regions. Therefore, an excessively large effective correlation scale ($L$) in the noise power spectrum ($C_{\ell}^{\mathrm{noise}}$) creates large

regions of intense fluctuation. Flat-sky analyses confined to these specific regions yield local power spectra that differ from localized patches in other regions of the sky.

• Large values of the multifractal parameter ($\lambda\gtrsim0.75$) inherently introduce singularities that can cause numerical instabilities when transforming between $E/B$ and $Q/U$ bases. Furthermore, the resulting strong non-Gaussianity can alter the specified variance ($C_{\ell}^{\mathrm{target}}$), thereby distorting the intended full-sky power-spectrum.

• Setting the advection strength parameter $c$ too high causes the field to become entirely decorrelated.

Analysis – Methods of Measuring Non-Gaussianity

All of the fields in Figure 2, Figure 3, and Figure 4 have exactly the same power spectrum. This highlights that the power spectrum is degenerate with respect to various spatial features. Given the above developments of anistropies, non-Gaussianity, and coherent structures, we are faced with a challenge: how do we actually measure and characterize their states and significance. In this section, we review and compare several methodologies.

1D Distributions

In terms of quantifying non-Gaussianity, it would seem obvious to first assess whether the probability distribution is indeed Gaussian or non-Gaussian. Figure 7 shows the full-sky distributions for isotropic, meridional, and zonal fluctuations for scalar and polarized intensity fluctuations (top and bottom, respectively). It is clear that the distributions become non-Gaussian, and increasingly so for increasing multifractal parameter $\lambda$. For the polarized intensity ($\xi_\delta= \sqrt{u_\delta^2 + q_\delta^2}$) the expected $\chi$ distribution with 2 degrees of freedom becomes increasingly heavy-tailed, departing significantly from this standard Rayleigh profile as the multifractal parameter $\lambda$ increases.

Vector (PDF)

Scale-dependent Distributions

The multifractal formalism of § 3.2 implies a scale-dependence to the distribution. To investigate this feature, we use the structure function.

The structure function (SF) is another popular tool for analyzing the statistical properties of fields (see, e.g., D. Russo & W. A. Jury 1987; D. Biskamp 2003; N. Cressie 1993; U. Frisch 1995). The “angle-averaged $p$-th order structure function” is defined as the moments of the increments, $\Delta f(\boldsymbol{x},\boldsymbol{r}) = f(\boldsymbol{x}) - f(\boldsymbol{x} + \boldsymbol{r})$:

where the angle-brackets denote ensemble averaging as well as performing an average over angles (spherical/circular shells) where $r = |\boldsymbol{r}|$. Similar tools have been used in previous analyses for the detection of remaining non-Gaussian signatures in component separated CMB maps (e.g., J. Jewell 2001; J.-L. Starck et al. 2004; M.-A. Miville-Deschênes et al. 2007).

In the Gaussian case, we have been using fields where the noise (viz. $\delta_{\ell,m}^{G}$ in Equation (16)) has the same distribution for all $x$, and the moments of the increments, the structure function (SF), scale as:

where $H$ is the Hurst parameter (a power-law index). Note that the exponent ($pH$) is linear in $p$ and thus $f(x)$ is monofractal. We now consider a field, $g(x)$, where its increments are instead described by

where $\psi(r)$ is an arbitrary, random, and statistically independent flux-scaling function. This flux-scaling function controls the scaling of $g(x)$, which is no-longer purely self-similar. The associated SF of $g(x)$ is

If $\psi(r)$ has no scale dependence, then $\langle|\psi(r)|^p \rangle= \text{(const.)}$ and we retrieve the monofractal behavior described above. However, consider a $\psi(r)$ scaling such that its statistics follow

where $K(p)$ is some moment-dependent function, and $L$ is the largest “outer-scale” of the scaling regime. Scaling of $K(p) \ne0$ implies that the fluctuations of various amplitudes change scale with different exponents (multiscaling) and each realization of such a process is a multifractal. Now, for $\overline{S}^{(p)}(r)$ of $g(x)$:

has scaling $\zeta_p = pH - K(p)$ and is non-linear for convex and non-zero $K(p)$.

The multifractal method of § 3.2 provides a model for synthesis of a non-Gaussian and multifractal field with log-normal flux model:

with moments scaling as

and

Note that theoretical multifractal scaling is strictly defined over ensemble averages, however, practically the ensemble averages are replaced with spatial averages. Because higher-order moments are dominated by rare, extreme events (intermittency), the spatial average of an individual, finite realization may lack these extremes, resulting in an empirical scaling profile that deviates from ensemble theory ([7] (T. Dudok de Wit 2004; W.-X. Zhou 2012)). Additionally, the LM method (§ 3.3) generates coherent structures which alters the probability distribution of extreme increments and fundamentally modifies the multifractal spectrum. This can cause higher-order moments to saturate or grow much slower than the purely multiplicative log-normal model (C. Rosales & C. Meneveau 2008; P. Subedi et al. 2014; J. Lübke et al. 2024).

Dissipation structures may not be entirely scale-invariant, particularly, in magnetohydrodynamic turbulence. For example, current sheets can dominate dissipative structures and challenge the scale invariance assumptions (P. D. Mininni & A. Pouquet 2009; P. Lesaffre et al. 2024; A. Pouquet et al. 2025). Nevertheless, the log-normal and advected models are suitable, since, in the first place, there is little observational constraint for $\zeta_p$ on 2D projections of multifractal three-dimensional fluctuations in the ISM. Additionally, the aim is to introduce a model with parameters to control the level of non-Gaussianity for the purposes of testing foreground removal techniques for future CMB missions.

Structure Functions

Since there is an angular power spectrum, and an angular correlation function, it would be natural to extend the relationship to derive some form of angular SF for both scalar and tensor quantities, e.g.,

where the average is taken over all pairs $(\hat{\mathbf{n}}_1,\hat{\mathbf{n}}_2)$ separated by Geodesic distance $d$ and oriented at angle $\phi$.

In a flat plane, the increment $|\mathbf{r}|$ is a straight line. On a sphere, there are two choices for measuring the distance between two locations $\hat{\mathbf{n}}_1$ and $\hat{\mathbf{n}}_2$:

• Geodesic distance: the actual distance traveling along the curved surface, $d = \arccos\left(\hat{\mathbf{n}}_1 \cdot\hat{\mathbf{n}}_2\right)$.

• Chord distance (Euclidean distance): The straight-line distance piercing through the sphere, $r = |\hat{\mathbf{n}}_1 - \hat{\mathbf{n}}_2| = 2\sin(d/2)$.

For “small-scale structure” (fluctuations of scale $d \ll1$ radian), the sphere locally looks flat, and these two definitions converge ($\sin(d/2) \approx d/2$, so $r \approx d$).

We estimate⁵ the angular SF of order $p$, parameterized by the Geodesic distance $d$ and the local orientation angle $\phi$ relative to Galactic North. The resulting second-order ($p = 2$) angular SFs for the isotropic and anisotropic fields (introduced in § 3.1) are presented in Figure 8. Each field exhibits a distinct directional signature. As expected, the isotropic field displays no dependence on $\phi$, maintaining a constant amplitude for a given scale $d$. Conversely, the anisotropic fields reveal strong angular dependence. The meridional field, characterized by North-South elongated fluctuations, maintains higher spatial correlations along the longitudinal axis. Because the SF measures the variance of increments, this higher correlation manifests as a suppression of the SF amplitude near $\phi\sim0^\circ$ (and by symmetry $\phi\sim180^\circ$). Similarly, the zonal field is highly

correlated East-West, resulting in minimized SF values near $\phi\sim90^\circ$.

Computing the complete anisotropic SFs on the sphere with HEALPix $N_{\mathrm{side}} \ge512$ can be computationally expensive ((K. M. Górski et al. 2005)). Thus, an analogous approach is to examine the statistics of the field at a given scale by band-pass filtering it in harmonic space, as is standard in wavelet-based analyses of multifractals (see, e.g., [8] (J. F. Muzy et al. 1991; J. F. Muzy et al. 1993; A. Turiel et al. 2006)). We multiply the spherical harmonic coefficients by an isotropic window $h_{\ell}^{(j)}$, which selects a narrow, roughly octave-wide range of multipoles around a scale $\ell_j$:

where $Y_{\ell,m}$ are the spherical harmonics and $a_{\ell m}$ are the spherical harmonic coefficients of the field. The windows are smooth and overlap only with their neighbours, so each $\Delta f_j$ isolates the fluctuations of $f$ at scale $\ell_j$. The filtering is done in harmonic space with healpy.almxfl() ((A. Zonca et al. 2019)), at the cost of one healpy.map2alm() and one healpy.alm2map() per band. Note that this approach ignores anisotropies (similar to averaging Equation (45) over fixed angles $\phi$), because $h_{\ell}^{(j)}$ is independent of $m$.

Figure 9 highlights the dependence of the PDF on scale. We band-pass filter the field(s) using Equation (46). For a monofractal (Gaussian, $\lambda= 0$) field, the fundamental shape of the distribution is scale-invariant; it remains strictly Gaussian at all scales, entirely characterized by its changing variance. In contrast, the multifractal fields ($\lambda> 0$) exhibit strong scale-dependent intermittency. As the physical scale decreases (larger $\ell$), the distributions progressively depart from a Gaussian profile, developing a sharper central peak and increasingly heavy tails.

Vector (PDF)

Figure 10 shows the scaling parameter $\zeta_p$ (Equation (41)) for the isotropic fields across a range of multifractal parameter $\lambda$. For the monofractal Gaussian field ($\lambda= 0$), $\zeta_p$ depends linearly on the moment order $p$. In contrast, the multifractal fields ($\lambda> 0$) exhibit a distinct non-linear dependence on $p$. The degree of this non-linearity grows with $\lambda$, quantitatively capturing the scale-dependent intermittency and increasingly heavy tails — the transition to non-Gaussianity — at smaller scales as shown in Figure 9.

Vector (PDF)

The statistics discussed thus far provide a robust characterization of the field’s scale dependence and intermittency, but they inherently discard spatial phase information. Consequently, while the isotropic high-pass filtering approach provides a computationally efficient means to evaluate scale-dependent intermittency, it fundamentally integrates out this angular phase information. The 2D angular SF therefore serves as a complementary metric, essential for isolating the global anisotropic signatures that are otherwise lost in 1D scale-filtered statistics. Additionally, the localized anisotropies of coherent structures generated via Lagrangian mapping (§ 3.3) tend to average out globally.

Minkowski Functionals

Another method to measure the non-Gaussianity in the maps is the Minkowski functionals (MF; J. Schmalzing & K. M. Górski 1998), which are a common tool to quantify the map-space higher-order statistical properties (e.g., P. Chingangbam & C. Park 2009; Planck Collaboration XVI 2016). Hadwiger’s theorem implies that, for any $D$ dimensional excursion set defined with

a threshold value $\nu$:

there exist $D + 1$ MFs that geometrically and topologically describe the morphology of the excursion set (see H. Hadwiger 1957). In the case of maps on the sphere $S^2$, such as temperature or polarized intensity, there are three kinds of MFs: $V_0$, $V_1$, $V_2$, which correspond to the area, the perimeter, and the connectivity of the excursion set, respectively. To compute the MFs we use pynkowski⁶ (A. Carones et al. 2024; J. Carrón Duque et al. 2024).

We first compute the MFs for the isotropic scalar and polarized intensity fluctuations: $i_\delta$, $\xi_\delta= \sqrt{q_\delta^2 + u_\delta^2}$, respectively. The resulting $V_0$, $V_1$, $V_2$ are shown in Figure 11 for Gaussian ($\lambda= 0$) and non-Gaussian ($\lambda> 0$) fields. The multifractal parameter, $\lambda$, exhibits a clear effect on the MFs for both types of fields (scalar, and polarization). The departure from Gaussianity increases with increasing $\lambda$.

Vector (PDF)

In Figure 12, we explore the influence of the full-sky anisotropies on the MFs. There are two clear classes of MF curves: the Gaussian, and the non-Gaussian. The non-Gaussian fields, including those with advected and anisotropic structures, exhibit similar profiles and thus cannot be distinguished. The MFs are not a strong detector for the full-sky anisotropies (§ 3.1), nor can they be used to characterize the small-scale coherent structures of advected fields (§ 3.3).

Vector (PDF)

Although neither structure functions nor Minkowski functionals are explicitly sensitive to translational phase shifts of LM advection itself, MFs provide a rigorous mathematical framework for extracting the morphology, spatial connectivity, and resulting shape of the intermittent structures.

Summary

Neither structure function or Minkowski function metric can reliably detect properties introduced by the Lagrangian mapping (§ 3.3). Further work is required to classify and detect properties of the fields when advection is included (§ 3.3). However, examination of multifractal scaling (J. Lübke et al. 2024), or techniques like wavelet scattering transforms (E. Allys et al. [5]; A. K. Saydjari et al. 2021; C. Auclair et al. 2024), Minkowski tensors (G. E. Schröder-Turk et al. 2013; C. Collischon et al. 2021), and fourth-order field curvature statistics (A. A. Schekochihin et al. 2002) could hold promise.

It is worth highlighting that the log-normal multifractal model (§ 3.2) produces a field that becomes more non-Gaussian with scale (see, § 4.2). This potentially contrasts with the results of, e.g., F. Rahman et al. (2021); F. Martire et al. (2023); F. Rahman et al. (2024, 2026), in which foreground emission exhibits a kurtosis-dominated non-Gaussian signature that is approximately constant with scale. If this observed trend is a robust feature of the interstellar medium (ISM), it represents a behavior that explicitly cannot be captured by the log-normal or advected multifractal models as currently formulated, in which case, setting $\lambda\approx0$ would be appropriate. However, because the framework in this paper specifically targets the currently unconstrained small-scale fluctuations, future high-resolution observations will be necessary to definitively test whether this scaling behavior persists at these scales. Additionally, both Minkowski functionals and structure functions struggle to describe the complete statistical and structural behaviour. This emphasizes the need for further expansion of the available tool box to comprehensively characterize the ISM.

The framework proposed in this paper provides tools to test methods for quantifying non-Gaussianity and coherent structures. Additionally, when combined with PySM to include data-driven large-scales for full-sky synthesis, it enables robust testing of non-Gaussian foreground pipelines under potential observational conditions.

Synthesizing New Full-sky Simulations

In this section, we return to the full-sky (and all-scale) synthesis from § 3. To synthesize a mock synchrotron observation, we apply the modulation field and large-scale constraints. Following the synthesis procedures of the PySM-s6 model, we use the small-scale power spectra ($C_\ell^{\mathrm{target}}$), modulation fields ($i_m$, $q_m$, $u_m$), large-scale fields ($i_0$, $q_0$, $u_0$), and a polarization fraction tensor transform to generate reference $I_{\mathrm{ref}}$, $Q_{\mathrm{ref}}$, $U_{\mathrm{ref}}$ maps ((The Pan-Experiment Galactic Science Group et al. 2025)). For the s6 model, these reference fields correspond to an observation at 23 GHz that can then be scaled with appropriate power-law index maps to realize fields at other specific observational frequencies.

We generate the small-scale temperature and polarization fluctuations using the provided PySM-s6 $C_{\ell}^{\mathrm{target}}$ spectra, which include non-zero $tt$, $ee$, $bb$, and $te$ correlations. We introduce spatial variation to these fluctuations by applying a spatial filter: the Galactic plane—defined by the inverse of a $2^\circ$ apodized mask with $f_{\mathrm{sky}} = 70\%$ from the Planck Legacy Archive⁷—receives meridional anisotropic fluctuations, while the rest of the sky receives isotropic fluctuations. This masking procedure serves a dual purpose: it generates intense, coherent filamentary structures within the Galactic plane, and it rebalances the loss of power at the Galactic poles caused by meridional fluctuations (as seen in e.g., Figure 2 and Figure 3).

We additionally require the small scales to be appropriately coherent with the large scales. By default, there is no phase correlation between the large scales associated with, e.g., $i_0$, $q_0$, $u_0$ and the small-scale realizations $i_\delta$, $q_\delta$, $u_\delta$. To induce this coherence, we introduce a large-scale deterministic component to the advection. The Lagrangian mapping approach uses a linearized advection equation and accumulates displacements on a scale-to-scale basis. The cascade (small-scale) advection components are calculated as described in § 3.3, and the accumulated displacements which include contributions from a large-scale (LS) field, $A_{\rm LS}$, are given by

The advection components derived from the large scales are:

where $c_{\rm LS}$ is an additional free parameter representing the large-scale advection strength. For PySM-s6 derived fields, the large-scale fields $A_{\rm LS}$ correspond to the large-scale constraints divided by the modulation maps: $t_0/t_m$, $e_0/e_m$, $b_0/b_m$ for the $t_\delta$, $e_\delta$, $b_\delta$ fields, respectively (see, § 3). In this manner, the small-scale fields are advected by the large scales and become slightly

phase-coherent across the large and small scales, while the large scales remain unmodified. A comparison of the full-sky $Q$ and $U$ maps and their corresponding $C_\ell$ spectra are provided in Appendix § 9.

Figure 13 shows “high-latitude” sub-regions of the polarized intensity reference field, $P_{\rm ref} = \sqrt{Q_{\rm ref}^{2} + U_{\rm ref}^{2}}$, for small-scale models generated with varying multifractal ($\lambda$) and advection ($c$, $c_{\rm LS}$) parameters. The top-left panel represents the default PySM-s6 isotropic Gaussian fluctuations. Moving across the parameter grid reveals how the spatial clustering driven by $\lambda$ and the coherent structures driven by $c$ couple to break Gaussianity. This demonstrates that the observational imprint on the small scales, as established in Figure 6, is preserved in the full synthesis.

Vector (PDF)

Figure 14 displays the corresponding pseudo-$C_\ell$ of $T_{\rm ref}$, $E_{\rm ref}$, $B_{\rm ref}$ for these regions and parameters (estimated using NaMaster; D. Alonso et al. [6]). The spectra of the modified fields (grey lines) closely track the default PySM-s6 case ($\lambda= 0$, $c = 0$, $c_{\rm LS} = 0$, black line). This reinforces our earlier assertion that two-point correlation functions (power spectra) alone fail to capture the higher-order statistics and phase coherence introduced by the multifractal and advected fields. Furthermore, it confirms that the local power is not significantly altered by the methods introduced in this paper.

Vector (PDF)

Figure 15 focuses on “low-latitude” cases of the polarized intensity reference fields. Here, we compare the default s6 isotropic Gaussian fluctuations with a model with anisotropic and non-Gaussian fluctuations ($\lambda= 1/4$, $c = 0$, $c_{\rm LS} = 0$), as well as one adding advection (“Anisotropic + NG + LM”: $\lambda= 1/4$, $c = 5$, $c_{\rm LS} = 1$). Near the Galactic plane, meridional anisotropies impose a strong directional preference on the fluctuations. However, the purely “Anisotropic + NG” model retains a somewhat diffuse, cloudy morphology. In contrast, the advected model exhibits filamentary structures, breaking the global homogeneity seen in the purely anisotropic case. At higher latitudes (e.g., Figure 13), the spatial masking smoothly replaces anisotropic fluctuations with isotropic ones. Yet, the advected fields continue to display filamentary features with power consistent with those regions of the sky.

A distinct advantage of this methodology is its stochastic nature. In other words, an ensemble of constrained realizations can be generated as opposed to synthesizing small-scale fields from singular observations or simulations. This provides a significant advantage for testing cosmological foreground subtraction methodologies and constraints for future CMB $B$-mode searches (e.g., [2, 1, 3] (S. Hanany et al. 2019; Y. Liu et al. 2025). Figure 16 shows a sample of five realizations for small-scale $i_\delta$ and PySM-s6 polarized intensity reference fields $P_{\mathrm{ref}}$. As illustrated, the five realizations exhibit visually consistent textures and coherences for both $i_\delta$ and $P_{\mathrm{ref}}$, confirming they are drawn from the same statistical distribution. However, each individual panel maintains a uniquely randomized spatial configuration. Once again, we emphasize the need to expand the analytical tool box beyond the Minkowski functionals and structure functions in order to adequately characterize such complex, spatially coherent fields.

Vector (PDF)

In summary, the synthesis framework presented here produces all-scale mock synchrotron skies. By integrating the large-scale observational constraints of PySM with our multifractal and advected small-scale models, we bridge the gap between simple statistical power spectrum matching and generating morphological properties in simulated skies.

Conclusion

In this paper, we present parameterized mock models of small-scale Galactic intensity and polarization fluctuations, designed specifically to provide rigorous testing grounds for component separation algorithms.

Our framework provides tunable parameters to govern both the morphology and statistical distributions of these fluctuations. Inspired by generalized scale invariance, directional preferences – or anisotropies – are introduced via a transformation matrix $G$ governed by the parameters $s_{\ell}$, $s_{m}$, $\theta_{0}$. Intermittency and non-Gaussianity are controlled using a log-normal multifractal model, which introduces scale-dependent distributions tuned by the multifractal parameter $\lambda$. Finally, phase-coherent structures are generated utilizing a Lagrangian mapping approach, which linearly advects the base field scaled by an advection parameter $c$.

To characterize the impact of these parameters, we evaluated the simulated fields using two complementary statistical metrics: structure functions and Minkowski functionals. While Minkowski functionals successfully capture geometric differences driven by non-Gaussianity in the multifractal fields, they fail to detect directional anisotropies or the small-scale coherent structures introduced by advection. Conversely, structure functions effectively isolate full-sky anisotropies and naturally quantify the multifractal intermittency through the fractal scaling parameter $\zeta_{p}$. However, when integrated over the full-sky, standard structure functions inherently average out spatial phase information, making them largely blind to localized advected structures.

By integrating our small-scale synthesis pipeline with the large-scale constraints and modulation fields provided by PySM, we successfully generate realistic, full-sky mock observations exhibiting both non-Gaussian intermittency and coherent filamentary features. The Python implementations for both the fluctuation synthesis and structure function estimators are publicly available at https://github.com/mab68/small_scale_synch.

Future work will focus on delivering a simulation suite with appropriate realism for CMB foreground subtraction, as well as establishing tighter constraints on the introduced free parameters by comparing our mock maps with high-resolution physical simulations and observational data. Additionally, the introduced methodologies can be similarly applied to, e.g., thermal dust emission.

Acknowledgments

Mark A. Bishop is a recipient of the New Zealand Space Scholarship administered by the New Zealand Space Agency on behalf of the New Zealand Government.

Mark A. Bishop utilized Google Gemini 3.1 Pro for minor language editing and Python syntax debugging.

The authors reviewed and verified all content and assume full responsibility for the manuscript.

This research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

Author Contributions

MAB: contributed conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing – original draft, writing – review & editing.

BSH: contributed conceptualization, formal analysis, investigation, methodology, validation, visualization, writing – review & editing.

KMG: contributed conceptualization, methodology, writing – review & editing.

CL: contributed conceptualization, methodology, writing – review & editing.

GR: contributed conceptualization, methodology, writing – review & editing.

Software: Astropy ([9]), HEALPix (K. M. Górski et al. 2005), healpy (A. Zonca et al. 2019), PySM (B. Thorne et al. 2017; A. Zonca et al. 2021; The Pan-Experiment Galactic Science Group et al. 2025), NumPy (C. R. Harris et al. 2020), Pynkowski (A. Carones et al. 2024; J. Carrón Duque et al. 2024), Matplotlib (J. D. Hunter 2007), SciPy (P. Virtanen et al. 2020), Seaborn (M. L. Waskom 2021), CMasher (E. van der Velden 2020), Kea (M. A. Bishop et al. 2026), NaMaster ([6]).

Full-sky Comparison

Example full-sky fields, synthesized following the procedure discussed in § 3 and § 5, for the $Q$ and $U$ polarization fields are shown in Figure 17 and Figure 18, respectively. We compare the base PySM-s6 model to a model with anisotropies, multifractality, and coherent advected structures. The large-scale structures in both models exhibit the same patterns and power. However, the small-scale fluctuations of the default s6 appear ‘cloudy’ whereas the “Anisotropic + NG + LM” model exhibits remarkably more correlated and filamentary structures. Additionally, the anisotropies in the Galactic plane can be seen, with the fluctuation sizes becoming elongated along lines of fixed longitude.

The auto-power spectra for the fields are shown in Figure 19. Both models have the same large-scale power ($\ell\lesssim50$). On the small scales ($\ell\gtrsim50$), the “Anisotropic + NG + LM” systematically exhibits slightly more power than the default PySM-s6 model by a factor of $\lesssim3$. This is likely due to the concentration of the fluctuations into the Galactic plane under the meridional anisotropic procedure (§ 3.1) and by using the same modulation map as the isotropic model, excess power could be generated at these scales near the Galactic plane. When masking out the Galactic plane using a mask of $f_{\mathrm{sky}} = 60\%$, the “Anisotropic + NG + LM” fluctuations become consistent with the s6 model.

Vector (PDF)