Exploration of 3D wavelet scattering transform coefficients for line-intensity mapping measurements

ABSTRACT

The wavelet scattering transform (WST) has recently gained attention in the context of large-scale structure studies, being a possible generator of summary statistics encapsulating non-Gaussianities beyond the reach of the conventional power spectrum. This work examines the three-dimensional solid harmonic WST in the context of a three-dimensional line-intensity mapping measurement to be undertaken by current and proposed phases of the CO Mapping Array Project (COMAP). The WST coefficients demonstrate interpretable behaviour in the context of noiseless CO line-intensity simulations. The contribution of the cosmological z ∼ 3 signal to these coefficients is also detectable in principle even in the Pathfinder phase of COMAP. Using the peak-patch method to generate large numbers of simulations and incorporating observational noise, we numerically estimate covariance matrices and show that careful choices of WST hyperparameters and rescaled or reduced coefficient sets are both necessary to keep covariances well-conditioned. Fisher forecasts show that even a reduced ‘shapeless’ set of ℓ-averaged WST coefficients show constraining power that can exceed that of the power spectrum alone even with similar detection significance. The full WST could improve parameter constraints even over the combination of the power spectrum and the voxel intensity distribution, showing that it uniquely encapsulates shape information about the line-intensity field. However, practical applications urgently require further understanding of the WST in key contexts like covariances and cross-correlations.

methods: statistical, diffuse radiation, large-scale structure of Universe

1 INTRODUCTION

The field of line-intensity mapping (LIM) promises to chart the large-scale structure (LSS) of the Universe as traced by a variety of atomic and molecular spectral lines, and thereby statistically and efficiently constrain high-redshift astrophysics and cosmology (cf. Kovetz et al. 2017, 2019 and Bernal & Kovetz 2022 for comprehensive reviews of the field in recent years). LIM experiments will observe large cosmological volumes for thousands of hours and obtain spatial-spectral data cubes that contain – underneath intimidating layers of noise and other contaminants – line emission from every galaxy in the survey volume faint or bright. The success of LIM science then depends on leveraging fluctuations in this integrated line emission with appropriate statistical methods.

The conventional approach in both simulations and analyses has been to use the spherically averaged three-dimensional power spectrum (or P(k), a function of comoving wavenumber k) as the principal LIM summary statistic. But the processes of structure formation and galaxy formation, and thus the line-intensity map produced by these processes, have significant non-Gaussianities to which the power spectrum is insensitive. Previous works, starting no later than Breysse et al. (2017), have proposed and developed the voxel intensity distribution (VID) as a supplementary LIM summary statistic that is sensitive to these non-Gaussianities.

Outside of LIM contexts, the need for statistical methods that can help extract non-Gaussian information has led to the use of approaches originally developed in computer vision contexts. Works of ambitious scale (e.g. Villaescusa-Navarro et al. 2022) have explored the use of convolutional neural networks, which in essence apply repeated filtering, downsampling, and clipping operations to the input. The resulting multiscale, non-linear nature of the operation strongly motivates use of these networks in LSS studies. However, use of such a network first and foremost requires determination of its weights, typically by optimising performance on a training set as measured by a loss function.

The need for a sizeable and representative training set is problematic in the context of high-redshift astrophysics and cosmology. Prior information at high redshift may well be both incredibly ill-determined and constantly evolving, demanding repeated extensive possibly expensive re-training of the network weights. Pfeffer, Breysse & Stein (2019) considered convolutional neural networks in the context of LIM with carbon monoxide (CO) line emission, and demonstrated strong degradation of the accuracy of inferred astrophysical constraints from unanticipated signal or noise effects. Such degradation reflects a tenuous interpretability (and credibility) of any outputs from these kinds of networks, strongly dependent on the training set and its limitations.

In view of such concerns, the wavelet scattering transform (WST) has recently gained attention in the LSS literature. In comparison to convolutional neural networks, the WST (based on principles developed by Mallat 2012) has a similar multiscale, multilayer nature suiting it to cosmology applications, but potentially does not sacrifice interpretability to the same extent. Cheng et al. (2020) applied the WST to a simulated cosmological observable for the first time, specifically in the context of weak-lensing convergence maps. That work suggested strong constraining power and interpretability compared to both traditional estimators – namely the power spectrum and peak counts – and convolutional neural networks. Further exploration by Cheng & Ménard (2021b) explicitly compared the WST to using the bispectrum, a traditional higher-order statistic that does capture non-Gaussian information, and found the WST to still be more robust and constraining while also being more Gaussian in variance.

This promise has already motivated exploration of the two-dimensional WST in the context of 21-cm cosmology (Greig, Ting & Kaurov 2022), as well as of the three-dimensional WST in the context of simulated matter density fields (Valogiannis & Dvorkin 2022b) and real-world galaxy surveys (Valogiannis & Dvorkin 2022a). However, work has yet to examine the application of three-dimensional WST coefficients (in the manner of Valogiannis & Dvorkin 2022a, b) to noisy observations (as Greig et al. 2022 did but for a two-dimensional WST implementation).

This paper aims to examine the three-dimensional spherical harmonic WST in the context of the CO Mapping Array Project (COMAP; Cleary et al. 2022), and in particular its measurement of the clustering of CO (1–0) line emission at z ∼ 3 (although COMAP ultimately aims to also map CO emission at z ∼ 7). In particular, we will explore the following questions:

• What is the detectability of cosmological CO emission in the WST coefficients derived from a simulated COMAP observation?

• Do parameter constraints from these coefficients improve significantly over the previously studied combination of P(k) with the VID?

The paper’s organization is as follows. We will review the summary statistics to be calculated, including the WST in Section 2 then consider in Section 3, the methods for forecasting inferences using those summary statistics derived from simulated observations. We will examine the resulting forecasts in Section 4 before discussing implications of these results in Section 5 and concluding in Section 6.

Unless otherwise stated, we assume base-10 logarithms, and a ΛCDM cosmology with parameters Ω_m = 0.286, |$\Omega _\Lambda = 0.714$|⁠, Ω_b = 0.047, H₀ = 100 h km s⁻¹ Mpc⁻¹ with h = 0.7, σ₈ = 0.82, and n_s = 0.96 to maintain consistency with previous simulations used by Ihle et al. (2019). Distances carry an implicit h⁻¹ dependence throughout, which propagates through masses (all based on virial halo masses, proportional to h⁻¹) and volume densities (∝ h³).

2 STATISTICS

This section will provide an overview of the three summary statistics applied to our simulated observations: the power spectrum (Section 2.1), the VID (Section 2.2) and the WST (Section 2.3).

2.1 Power spectrum

The full three-dimensional power spectrum of a field |$f(\boldsymbol {x})$|⁠, expressed as a function of a three-dimensional positional vector |$\boldsymbol {x}$| is given through its Fourier transform |$\tilde{f}(\boldsymbol {k})$|⁠, a function of the wavevector |$\boldsymbol {k}$|⁠:

$$\begin{eqnarray*} P(\boldsymbol {k}) \propto |\tilde{f}(\boldsymbol {k})|^2 = \tilde{f}^{*}(\boldsymbol {k})\tilde{f}(\boldsymbol {k}) \end{eqnarray*}$$

(1)

with the proportionality constant being the inverse of the survey volume. The spherically-averaged three-dimensional power spectrum P(k) is then a binned reduction of this |$P(\boldsymbol {k})$| in spherical shells of constant wavenumber k:

$$\begin{eqnarray*} P(k) = \left\langle P(\boldsymbol {k})\right\rangle _{|\boldsymbol {k}|=k}. \end{eqnarray*}$$

(2)

Because the Universe is mostly isotropic and homogeneous on large cosmological scales, P(k), the monopole of the power spectrum contains most of the information that would be contained in |$\tilde{f}(\boldsymbol {k})$| or |$P(\boldsymbol {k})$|⁠, and has the advantage that we can average independently observed Fourier modes across k-shells to obtain a better measurement of P(k). The result is a measurement of the power of field contrast fluctuations associated with different comoving scales L ∼ 2π/k.

Ultimately, the fact that we observe LSS in redshift space and not in real comoving space does introduce anisotropies, which are captured in higher-order multipoles of |$P(\boldsymbol {k})$|⁠. But the plain monopole P(k) remains a veritable workhorse of basic cosmology theory. We commonly describe the distribution of matter across the Universe with the matter power spectrum P_m(k). We then describe the density contrast of a biased tracer of matter with a power spectrum b²P_m(k), given some linear tracer bias b by which the matter density contrast δ_m scales to the tracer density contrast δ = bδ_m. If using a dark matter halo population as a tracer, the halo bias b(M_h) as a function of halo (virial) mass M_h is well studied (cf. Manera, Sheth & Scoccimarro 2010; Pillepich, Porciani & Hahn 2010; Tinker et al. 2010). The bias of galaxies or other associated tracers, including the line intensity contrast in a LIM observation, then follows from the galaxy–halo connection (Wechsler & Tinker 2018). Redshift-space distortions theories of non-linear bias and other complications (as considered by, e.g. Bernal et al. 2019a; Bernal, Breysse & Kovetz 2019b or Moradinezhad Dizgah et al. 2022) operate on top of these fundamental ideas.

Despite the loss of non-Gaussian information and anisotropic information, the greatest strength of P(k) as a summary statistic is in the ease of understanding its variance and covariance. Power spectra of statistically independent signal and noise components simply add; the variance of the total P(k) for a given k-shell containing N_m(k) Fourier modes is simply P²(k)/N_m(k); and the covariance matrix of P(k) (especially in a noise-dominated observation) is largely diagonal,¹ easing statistical considerations in measurement and inferences.

Finally, it is straightforward to consider the cross-correlation power spectrum between two fields |$f(\boldsymbol {x})$| and |$g(\boldsymbol {x})$| through their Fourier transforms |$\tilde{f}(\boldsymbol {k})$| and |$\tilde{g}(\boldsymbol {k})$|⁠:

$$\begin{eqnarray*} P_{f\times g}(\boldsymbol {k})\propto \tilde{f}^{*}(\boldsymbol {k})\tilde{g}(\boldsymbol {k}) \end{eqnarray*}$$

(3)

with the proportionality constant and the spherically-averaged P_{f × g}(k) obtained in the same way as for the auto-correlation power spectrum. In fact, it is fair to say that the auto-correlation P(k) is merely a special case of the cross-power spectrum where f(x) = g(x). Given the complexity of galaxy formation and the need for multitracer astrophysics on cosmological scales through different LIM and LSS observations, not to mention the need for mitigation of different foregrounds and systematics present in these observations, cross-correlations between different experiments surveying common comoving volumes are highly desirable. At present, the cross-power spectrum remains the best studied and thus best understood statistic to address these needs.

2.2 VID

The VID was first put forward in LIM contexts by Breysse et al. (2017) who characterize it as an extension of the |$\mathcal {P}(D)$| technique in radio astronomy put forward by Scheuer (1957). The central idea is to characterize the distribution of intensities (or, in the time and language of Scheuer (1957), deflexions on a power recording chart – hence |$\mathcal {P}(D)$| for probability of deflexion) in a confused map of observed emission from a population of sources, and relate this to the underlying luminosity function of sources. LIM experiments tend to operate in a confused regime to focus on resolving the large-scale clustering of line emission, so the VID technique holds extreme relevance. The study of Breysse et al. (2017), as well as follow-up by Ihle et al. (2019) among others, showed that the VID truly encapsulates non-Gaussian small-scale information beyond the power spectrum related to the shape of the source luminosity function and to the stochasticity of line emission.

Being a histogram of voxel intensities (or temperatures in the convention of a cm-wavelength experiment like COMAP) with some bin width ΔT and range (T_min, T_max), the VID is straightforward to calculate from any data cube. Its interpretation, however, is significantly less straightforward due the nature of the VID of a confused map. The VID is given by a series of recursive convoluted probability distributions of voxel intensities conditioned on voxels containing specific numbers of emitters, complicated by the clustering of sources. Observational effects like instrumental resolution and line widths further affect the VID. None the less, the behaviour of the VID – including its covariance – can still be understood through simulations as in Ihle et al. (2019), and even through concrete mathematical expressions (Sato-Polito & Bernal 2022).

Furthermore, while not as straightforwardly as for the power spectrum, one can design one-point cross-correlation statistics between multiple tracers to reject disjoint contaminants. Breysse, Anderson & Berger (2019) propose a one-point analysis for a 21-cm intensity map with VIDs conditioned on the presence or absence of a galaxy in an external catalogue surveying the same volume, showing that the ratio of Fourier transforms of these conditional VIDs provides a one-point cross-correlation statistic unbiased by the presence of foregrounds. Similar joint VID analyses can be applied in the context of LIM observations and internal or external cross-correlations between LIM data sets (Breysse et al., in prep.).

2.3 Solid harmonic WST

As with P(k) and the VID, the WST has a concrete mathematical description and strong interpretability, reviewed pedagogically (albeit in a chiefly two-dimensional context) by Cheng & Ménard (2021a). The essence is to convolve the input field with a family of wavelets, and the complex modulus of this convolved field is the ‘scattered’ output. The coefficients of the WST come from spatial averages of the output after different levels of scattering. The resemblance of the WST to convolutional neural networks arises from the possibility of recursive application of this scattering operation (including the non-linear ‘clipping’ introduced by the modulus operation and the downsampling coming from averaging the scattered outputs) with wavelets of different scales and anisotropies.

In practice, works exploring the WST only deal with the first- and second-order scattering coefficients, i.e. the field average after one or two scattering operations. This is already sufficient to probe scale interactions in ways that depend on the wavelets used and their frequency response.

Following the implementation of the three-dimensional WST by Andreux et al. (2020) in the Kymatio software package, we specifically will consider a family of solid harmonic wavelets. This solid harmonic WST was originally considered by Eickenberg et al. (2018)² in the context of molecular energy regression. Molecular properties like the ground state energy depend on the electronic density, popularly calculated via density-functional theory. Although this is not necessarily a prohibitive calculation, the atomistic state – i.e. the description of a molecule in terms of the constituent atoms and their positions – is still a simpler description than the true electronic density, and can be geometrically represented through different pseudo-electronic densities. The aim of Eickenberg et al. (2018) was to obtain multilinear regressions to map the WST of this molecular geometry directly to various molecular properties in a predictive fashion, which the work demonstrated with excellent performance metrics.

In both the molecular energy regression problem and our case study, the goal is to constrain parameters describing the structure underlying the input field, be it an organic molecule or a sample of high-redshift galaxies. However, in our case the starting point is not an approximate, noise-free description of particle density, but rather an empirical, noisy measurement of true CO volume emissivity. None the less, the central idea behind the WST as a summary statistic probing multiscale interactions remains in all contexts, and even at the time of the work of Eickenberg et al. (2018), solid harmonic wavelets had already been used in three-dimensional image processing contexts outside of quantum chemistry.³

The solid harmonic wavelets are indexed by the logarithmic scale parameter j, the angular frequency index ℓ, and the orientation index m ∈ [− ℓ, ℓ]. Furthermore, we will express these wavelets as a function of the position vector |$\boldsymbol {u}$| in image coordinates (so in units of pixels or grid width) rather than physical or comoving coordinates to mirror the description given in other works.

Before introducing the scale j, the ‘mother’ wavelets are given by the Laplacian spherical harmonic |$Y_\ell ^m$| rescaled by |u|^ℓ to make a (regular) solid harmonic, then multiplied by a Gaussian envelope of width σ:

$$\begin{eqnarray*} \psi _{\ell }^m(\boldsymbol {u}) = \frac{\exp {\left[-|\boldsymbol {u}|^2/\left(2\sigma ^2 \right)\right]}}{(2\pi \sigma ^2)^{3/2}}|\boldsymbol {u}|^\ell Y_\ell ^m(\boldsymbol {u}/|\boldsymbol {u}|). \end{eqnarray*}$$

(4)

The wavelets dilate to different scales as described by j:

$$\begin{eqnarray*} \psi _{j,\ell }^m(\boldsymbol {u}) = 2^{-3j}\psi _{\ell }^m \left(2^{-j}\boldsymbol {u} \right). \end{eqnarray*}$$

(5)

We will fix σ = 1 for the remainder of this work, which is also what Eickenberg et al. (2018) do in effect for their work.⁴

With the family of mother and dilated wavelets defined, we can explicitly write down the first- and second-order wavelet scattering coefficients S₁ and S₂ given an input field |$I(\boldsymbol {u})$| (again, as a function of image coordinates). We also introduce an extra parameter q, by which we exponentiate the complex modulus operation:

$$\begin{eqnarray*} S_1(j,\ell ;q) = \left\langle \left(\sum _{m=-\ell }^\ell |I(\boldsymbol {u})*\psi _{j,\ell }^m|^2\right)^{q/2}\right\rangle _\boldsymbol {u}. \end{eqnarray*}$$

(6)

We only evaluate second-order coefficients across different scales, and not across different ℓ. Then for scales j and |$j^{\prime }$| > ,

$$\begin{eqnarray*} S_2(j,j^{\prime },\ell ;q) = \left\langle \left(\sum _{m^{\prime }}\left|\left(\sum _m|I(\boldsymbol {u})*\psi _{j,\ell }^m|^2\right)^{1/2}*\psi _{j^{\prime },\ell }^{m^{\prime }}\right|^2\right)^{q/2}\right\rangle _\boldsymbol {u}. \nonumber \\ \end{eqnarray*}$$

(7)

Since all coefficients involve sums in m, the calculation discards rotational information, but the variation with ℓ still encodes some level of angular or morphological information. Given that our wavelet basis is given by the solid harmonics, we should expect the variation with ℓ to specifically encode the different amounts of isotropic power contained in field monopoles, dipoles, quadrupoles, octopoles, hexadecapoles, and so on.

With our choice of solid harmonic wavelets, only the isotropic power is encoded as orientation information is lost in the complex modulus operation, but a different choice of wavelet family could encode further information in the WST, which may be desirable given the expected redshift-space distortion of the CO signal. However for this work, we continue with the solid harmonic WST as it has a readily available implementation in Kymatio, and we leave exploration of oriented wavelets in the context of CO LIM to other future work.

The choice of q determines the relative weighting of overdensities against underdensities in determining the WST coefficients. Given the noisy nature of LIM data cubes, a value of q > 1 to excite a moderate but not excessive emphasis on overdensities is likely desirable. We can calculate the WST with q = {0.5, 0.75, 1.0, 1.5, 2.0} to examine detection significance, but beyond that we will consider only the optimal value of q out of those five choices, rather than consider multiple values of q at once. This is to avoid further worsening the conditioning of the covariance matrix, as we expect many of the WST coefficients to be strongly correlated across different values of q.

Note that Andreux et al. (2020) describe a low-pass filter convolution at each step after the modulus operation. For the solid harmonic wavelet transform as implemented by Kymatio at time of writing, this is done with a simple summing operation (equivalent to integration over the volume as described in Eickenberg et al. (2018), given |$\boldsymbol {u}$| is in image coordinates). We simply divide the output from Kymatio by the number of voxels to obtain the spatial average described above. The difference between averaging and summing is immaterial to the key results of this work like detection significance and parameter constraints.

As we have discussed previously, the Universe is mostly homogeneous and therefore one might suppose that there should be greater information contained in variation of the coefficients with scale j than in variation with multipole index ℓ. Given this expectation, we propose the following reduced three-dimensional WST coefficients that discard angular information:

$$\begin{eqnarray*} s_1(j;q)&\equiv 2^j\langle S_1(j,\ell ;q)\rangle _{\ell \ge 0}; \end{eqnarray*}$$

(8)

$$\begin{eqnarray*} s_2(j,j^{\prime };q)&\equiv 2^{j^{\prime }-j}\langle S_2(j,j^{\prime },\ell ;q)/S_1(j,\ell ;q)\rangle _{\ell > 0}; \end{eqnarray*}$$

(9)

$$\begin{eqnarray*} s_{20}(j;q)&\equiv \langle S_2(j,j^{\prime },\ell =0;q)/S_1(j,\ell =0;q)\rangle _{j^{\prime } > j}. \end{eqnarray*}$$

(10)

The calculation of s₁ averages entirely over angular information, while the reduced second-order coefficients separate out ℓ = 0 as a degenerate case (since |$Y_0^0$| is a constant) to obtain s₂₀, but then still average over ℓ > 0 to obtain s₂. This should provide information about peaks across different scales and interactions between scales, but discard all shape information. Therefore we will dub this the ‘shapeless’ WST coefficient set.

We also consider a less aggressively reduced set of WST coefficients that does preserve variation with ℓ and therefore shape information as well as scale-dependence:

$$\begin{eqnarray*} \tilde{S}_1(j,\ell ;q)&\equiv 2^jN_1(j,\ell)S_1(j,\ell ;q); \end{eqnarray*}$$

(11)

$$\begin{eqnarray*} \tilde{S}_2(j,j^{\prime },\ell ;q)&\equiv 2^{j^{\prime }-j}N_2(j,j^{\prime },\ell)S_2(j,j^{\prime },\ell ;q)/S_1(j,\ell ;q). \end{eqnarray*}$$

(12)

We excise |$\tilde{S}_2(j,j^{\prime },\ell =0)$| for all > j + 1 from this coefficient set due to the fact that these correlate almost perfectly with |$\tilde{S}_2(j,j^{\prime }=j+1,\ell =0)$| in practice. With that excision, we have fully defined what we will refer to as the ‘rescaled’ set (summarized with the raw and shapeless WST coefficients in Table 1). Here N₁ and N₂ are arbitrary normalization factors to be defined later that will not depend on S₁, S₂, or the input field in general. Rather N₁ and N₂ will be defined to scale the data in a way that optimizes the condition number of the estimated covariance matrix, which we will discuss further in Section 4.2. Like the shapeless coefficients, the second-order coefficients are rescaled by appropriate first-order counterparts to significantly reduce correlations between coefficients, which is also key for covariance matrix conditioning.

Table 1.

Definitions of the different sets of coefficients derived from the WST, as used in this work.

Label	Definition
Raw	• S1(j, ℓ; q) as defined in equation (6), for allj ∈ {0, …, J} and ℓ ∈ {0, …, L}
	• S2(j, |$j^{\prime }$|, ℓ; q) as defined in equation (7), for all j ∈ {0, …, J − 1}, |$j^{\prime }$| > j, and ℓ ∈ {0, …, L}
Rescaled	• |$\tilde{S}_1(j,\ell ;q)$| as defined in equation (11), for allj ∈ {0, …, J} and ℓ ∈ {0, …, L}
	• |$\tilde{S}_2(j,j^{\prime },\ell ;q)\propto S_2(j,j^{\prime },\ell)/S_1(j,\ell)$| as defined in equation (12), for all j ∈ {0, …, J − 1}, |$j^{\prime }$| > j, and ℓ ∈ {0, …, L} – but omit if |$j^{\prime }$| > j + 1 and ℓ = 0
Shapeless	• s1(j; q) ∝ 〈S1(j, ℓ)〉ℓ as defined in equation (8), for all j ∈ {0, …, J}
	• s2(j, |$j^{\prime }$|; q) ∝ 〈S2(j, |$j^{\prime }$|, fℓ)/S1(j, ℓ)〉ℓ > 0 as defined in equation (9), for all j ∈ {0, …, J − 1} and |$j^{\prime }$| > j
	• |$s_{20}(j;q)\propto \left\langle S_2(j,j^{\prime },\ell =0)/S_1(j,\ell =0)\right\rangle _{j^{\prime }\gt j}$| as defined in equation (10), for all j ∈ {0, …, J − 1}

The choice of these forms is broadly informed by the orientation reductions and normalizations discussed by Cheng & Ménard (2021a), but also by the correlations of the raw WST coefficients in simulations that we will show later. The prefactors of 2^j and |$2^{j^{\prime }-j}$|⁠, present in both the rescaled and shapeless coefficients, are immaterial to correlations but serve to normalize the coefficients against general trends of the coefficient values in presentation, which should improve conditioning of the covariance matrix. In particular, while we will need to define optimal N₁ and N₂ formulae for the full rescaled coefficients based on the behaviour of the signal and noise, the results of Section 4.2 will show that the shapeless coefficients tend to result in well-conditioned covariance matrix estimates without further fine-tuning.

The shapeless coefficients also significantly reduce the dimensionality of any linear algebra involved. If j ranges from 0 to J and ℓ ranges from 0 to L – where for this work we will use J = 4 and L = 3 – the number of raw WST coefficients for fixed q is (J + 1) (J + 2) (L + 1)/2 = 60, of which we only omit J (J − 1)/2 = 6 when considering the rescaled coefficient set. Meanwhile, we have |$J+1\, s_1$| coefficients, |$J\, s_{20}$| coefficients (since s₂₀ (J; q) is not a valid calculation), and |$J(J+1)/2\, s_2$| coefficients for a total of J (J + 5)/2 + 1 = 19 shapeless WST coefficients.

3 FORECAST METHODOLOGY

With our summary statistics defined, we move to outlining methods for modelling CO line emission (Section 3.1), simulating COMAP observations (Section 3.2), choosing bins for summary statistics (Section 3.3), and projecting parameter constraints in a Fisher matrix analysis (Section 3.4).

3.1 CO line model

CO line emission in truth involves a variety of small-scale environmental and dynamical factors that are well beyond the scope of this paper,⁵ whose focus is on recovery of large-scale non-Gaussian information. For our purposes, it is sufficient to use a simple halo model that associates the virial mass M_h and cosmological redshift z of a dark matter halo with a CO luminosity L. Specifically, we use the empirical CO (1–0) model of Chung et al. (2022), which assumes no redshift dependence and relates M_h to L via a double power law:

$$\begin{eqnarray*} \frac{L(M_h)}{L_\odot } = \frac{10^C}{(M_h/M)^A+(M_h/M)^B}. \end{eqnarray*}$$

(13)

Here, C represents an overall normalization, A and B define the slopes of the double power law at extreme masses, and M is a characteristic mass scale that defines the double power-law turnabout point along with C. Compared to the presentation in Chung et al. (2022), note the slight tweak to how C is defined in relation to the normalization of the double power law, which is now 10^C rather than C, and against units of L_⊙ rather than K km s⁻¹ pc². The model also prescribes a mean-preserving lognormal scatter around the above relation with the scatter width σ expressed in units of dex.

Note that in the context of this model, we define A and B such that A < B, so that the L(M_h) relation follows a power law with exponent −A at low mass M_h ≪ M and a different power law with exponent −B at high mass M_h ≫ M. Overall, the double power-law parametrization of L(M_h) reflects the expectation (considered in more detail in the original presentation of the model by Chung et al. 2022) that gas content grows and thus CO emission brightens with halo mass for low-mass objects, but feedback and quenching stunt this growth at sufficiently high mass.

We adopt fiducial values based on the UM + COLDz model of Chung et al. (2022), which combines priors from the Behroozi et al. (2019) empirical model of the galaxy–halo connection and observational CO luminosity function constraints from Riechers et al. (2019). The parameter values thus chosen are A = −2.75, B = 0.05, C = 6.3 (after accounting for the change in normalization compared to Chung et al. 2022), log (M/M_⊙) = 12.3, and σ = 0.42. This is also the fiducial model of Chung et al. (2021).

To reflect the finite size of CO line profiles due to peculiar velocities within each dark matter halo, we also use the model of Chung et al. (2021) to associate a CO line width with each halo based on its virial velocity, and apply appropriate line broadening when generating CO signal cubes. We follow the implementation of this model by Chung et al. (2022), and refer the reader to that work for further details. We also assert that small-scale effects outside of this line broadening (e.g. emission from satellite galaxies around each central dark matter halo) will be negligible in the context of a LIM survey, which probes large-scale line-intensity fluctuations rather than individual sources.

3.2 Simulations

Table 2 defines survey parameters for two phases of Ka-band (26–34 GHz) observations with COMAP. The first ‘Pathfinder’ phase is nominally a five-year campaign with the COMAP Pathfinder instrument, which has already been in mostly continuous operation for over three years at the time of writing. The second phase considered is the COMAP Expanded Reionisation Array (COMAP-ERA) concept outlined by Breysse et al. (2022), which involves deployment of additional Ka-band feeds and extended survey operations.

Table 2.

COMAP Ka-band survey parameters for the Pathfinder and COMAP-ERA phases, as used for the simulations in this work.

Observational parameter	COMAP phase:
	Pathfinder	COMAP-ERA
Frequency coverage (GHz)	26–34	26–34
Channel bandwidth (MHz)	31.25	31.25
Solid angle per field (deg2)	4	4
Beam FWHM (arcmin)	4.5	4.5
Map pixel size (arcmin) for:	–	–
• P(k) and WST analyses	2	2
• VID analysis	4	4
Number of fields	3	3
Noise per WST input voxel (µK) at:	–	–
• 26–28 GHz	21.2	4.5
• 28–30 GHz	18.7	4.0
• 30–32 GHz	18.2	3.9
• 32–34 GHz	19.8	4.2

In the context of this work, the only difference between the Pathfinder and COMAP-ERA phases is the expected noise per voxel, which have been extrapolated from noise levels in COMAP Field 1 observations rather than directly calculated from survey parameters.

• For the Pathfinder experiment, we make the same assumption as Chung et al. (2022) that the noise per voxel will decrease relative to COMAP Early Science data by a factor of |$\sqrt{40}$|⁠. The resulting Pathfinder projections broadly agree with the 17.8 µK value quoted for simulations in Chung et al. (2022).

• For COMAP-ERA, we use the assumption of Breysse et al. (2022) that COMAP accrues 1000 hr of integration time per field per year (implying that the COMAP Pathfinder survey should be equivalent to 5000 hr of integration time per field), and their projection that the COMAP-ERA survey should accrue the equivalent of 110 000 Pathfinder hours. This implies a further improvement in noise per voxel by a factor of |$\sqrt{22}$| when going from COMAP Pathfinder data to COMAP-ERA data. The assumed hours per field per year (and therefore the factor obtained here) is credible to within a factor of order unity, certainly a much smaller factor than the level of noise improvement being forecast between COMAP Early Science data and COMAP-ERA data.

Breysse et al. (2022) also define an intermediate COMAP Epoch of Reionisation (COMAP-EoR) experiment, but for our purposes it is sufficient to consider the two extremes of near-term and far-future CO LIM data sets, with their corresponding noise levels.

We will want to generate large numbers of simulations of such data sets, as numerical estimation of WST covariances is necessary in the absence of their analytical understanding. Fortunately, the peak-patch method (Stein, Alvarez & Bond 2019) grants us dark matter halo catalogues of sufficient accuracy for thousands of simulated COMAP observations. In particular, we begin with the same set of peak-patch simulations used by Ihle et al. (2019), which comprises 161 independent light-cone simulations across a comoving box of volume |$L_\text{box}^3=(1140\,$| Mpc)³ with a grid resolution of N_cells = 4096³. As the light-cone extent is equivalent to 9.6° × 9.6° in transverse dimensions and 26–34 GHz in CO (1–0) observing frequency, we split the halo catalogue associated with each lightcone into 16 mock COMAP fields to simulate 2576 semi-independent COMAP observations at z ∼ 3.

We translate the halo catalogue in each field into a simulated CO signal using limlam_mocker⁶ – including line broadening as mentioned in Section 3.1 – and apply angular smoothing by a Gaussian beam with full width at half maximum of 4.5 arcmin, matching the COMAP instrumental resolution. We show an example realization of the CO signal at this stage in Fig. 1.

Figure 1.

Example realization of the CO signal, including line broadening and angular smoothing but excluding further high-pass filtering and additive Gaussian noise. We show the signal at two faces of the simulated volume.

On top of these small-scale effects, we can appropriately account for the removal of large-scale modes from the signal expected from the COMAP pipeline. Appendix D of Chung et al. (2022) suggests an approximation for the COMAP power spectrum transfer function in early science data with the following component in particular likely associated with high-pass filtering and subtraction of large-scale spectral modes:

$$\begin{eqnarray*} \mathcal {T}_\text{hp}\left(k_\perp ,k_\parallel \right) = \frac{1}{\left(1+e^{5-100\text{ {}Mpc}\cdot k_\perp } \right)\left(1+e^{5-200\text{ {}Mpc}\cdot k_\parallel } \right)}, \end{eqnarray*}$$

(14)

with k_∥ denoting the magnitude of the line-of-sight component of the wavevector |$\boldsymbol {k}$| for which we evaluate the transfer function, and k_⊥ denoting the magnitude of the transverse part of the wavevector (such that |$k_\perp ^2+k_\parallel ^2=k^2$|⁠).

We apply an equivalent Fourier filter to our simulated data cubes by taking the Fourier transform, multiplying each point of the Fourier-space data cube by |$\mathcal {T}_\text{hp}^{1/2}$| (since |$\mathcal {T}_\text{hp}$| is effectively the transfer function for the squared magnitude of the Fourier transform), and inverting the Fourier transform to obtain the high-pass filtered CO temperature cube. The final simulated COMAP observation includes uniform Gaussian noise added to this filtered signal with noise per voxel as specified in Table 2 for each COMAP phase.

We will consider each 2 GHz subband of the total frequency coverage separately, and in particular we will excise four frequency channels at each edge of each sub-band. This leaves a data cube spanning 60 voxels across each angular dimension – coarsened to 30 voxels on each side for the VID analysis – and 56 voxels along the line of sight. This results in a roughly isotropic cube in comoving space, corresponding to a ≈220 × 220 × 240 Mpc³ cube at the central redshift of z ≈ 2.8.

3.3 Bins for statistics

We have already significantly detailed the power spectrum, the VID, and the WST in Section 2. Here we briefly clarify the scale or temperature bins used for each of these statistics.

As stated above, for P(k) and WST analyses, we use a temperature cube of 60 × 60 × 56 voxels for each 2 GHz sub-band with each voxel spanning roughly 3.6 × 3.6 × 4.2 Mpc³. The extent and grid size of the cube suggests k_min ∼ 0.03 Mpc⁻¹ and k_max ∼ 0.8 Mpc⁻¹. This is reflected in the k-bins used for the power spectrum calculation, which are 29 linearly spaced bins with Δk ≈ 0.03 Mpc⁻¹ (where the exact value depends on the sub-band) and the lower edge of the first bin (which contains as little as one Fourier mode) placed at k = 0.

We use WST parameters previously outlined in Section 2.3, including a maximum scale index of J = 4 and a maximum spherical harmonic multipole index of L = 3. Each wavelet spans a band of Fourier wavenumbers, but if we associate a scale index j roughly with a span of 2^j voxels (meaning it is essentially senseless to set J above the base-2 logarithm of the extent of the data cube), the scales of j = {0, 1, 2, 3, 4} correspond roughly to {4, 8, 16, 32, 64} Mpc scales. If this is roughly the half-wavelength of a Fourier mode, the corresponding wavenumbers would be k ∼ {0.8, 0.4, 0.2, 0.1, 0.05} Mpc⁻¹. Therefore, the range of scales being probed is quite similar between P(k) and the WST.

To better match the angular size of the CO map voxels to the COMAP beamwidth, we use a coarser pixel size for VID calculations, meaning we work with a cube of 30 × 30 × 56 voxels. We generate voxel temperature histograms using 12 log-spaced bins whose edges span 10¹–10^1.6 µK (or approximately 10–40 µK, limiting ourselves to probing >4σ excesses above COMAP-ERA noise). This is likely a conservative choice, given that Ihle et al. (2019) used bins going up to 100 µK, but justified given our somewhat dimmer CO emission model and the lower noise level of COMAP-ERA simulations.

3.4 Fisher formalism

We undertake an analysis in the Fisher matrix formalism⁷ to consider potential constraints on the CO emission model outlined in Section 3.1. Given the necessary but tenuous assumption of Gaussian covariances, these should not be taken as quantitative forecasts, but qualitative comparisons of constraining power for different summary statistics should still be of interest.

For this work, we fix B and M at their fiducial values. In the basis of the remaining parameters {λ_r} = {A, C, σ}, the Fisher matrix for an observable vector |$\boldsymbol {O}_\alpha$| with covariance matrix |$\bf{C}_{\alpha \beta }$| is given by

$$\begin{eqnarray*} \bf{F}_{rs} = \sum _{\alpha ,\beta }\frac{\partial {\boldsymbol {O}}_{\alpha }}{\partial \lambda _r}\bf{C}^{-1}_{\alpha \beta }\frac{\partial {\boldsymbol {O}}_{\beta }}{\partial \lambda _s}, \end{eqnarray*}$$

(15)

approximating the derivatives with central difference quotients based on median values for the observables across 2576 simulated COMAP survey volumes (as described in Section 3.2) obtained at parameter values 0.05 above and below the fiducial value for each parameter. We impose no additional priors.

In fact, we have observable vectors for each one of four COMAP sub-bands, which we will index with i. To minimize effects of spurious artefacts in numerically estimated covariances, we evaluate |$\bf{F}_{i;rs}$| separately for each sideband and obtain a total Fisher matrix |$\sum _i\bf{F}_{i;rs}$| by simple summation. We thus treat the four COMAP bands as independent observations, which seems justifiable in practice as correlation coefficient matrices for all observables show negligible correlation (⁠|${\lesssim}1\,\hbox{per cent}$|⁠) across frequency sub-bands. Given the large comoving length spanned by the COMAP frequency coverage, this is perhaps not surprising, certainly for statistics at smaller scales.

Recall also that COMAP observes three independent survey fields, which should decrease covariances by a factor of three. We equivalently multiply the Fisher matrix by 3 to reflect this fact, then invert this Fisher matrix to obtain the parameter covariance matrix for the full COMAP Pathfinder or COMAP-ERA survey.

4 RESULTS

This section will eventually answer the questions, we posed at the outset of this paper about the WST in the context of COMAP. However, it will be instructive to first qualitatively examine how the WST changes with our CO model parameters, which we do in Section 4.1. We must also examine the covariance of the WST coefficients in Section 4.2 before considering detection significances in Section 4.3, and both will play a role in our choice of q for the WST for finally forecasting parameter constraints in Section 4.4.

4.1 Qualitative examination of WST coefficients in CO line-intensity simulations

A key motivation for using the WST is the fact that the coefficients summarize scale interactions in an interpretable fashion, characterizing not only total power at different scales but also information about morphology and sparsity. To demonstrate this, we ran an additional set of simulations for each one of our 2576 mock COMAP volumes with parameter values shifted ±0.1 away from fiducial values. These are separate from the simulations used for our Fisher analyses, which shift parameter values only by ±0.05. The larger shift is to better illustrate how the WST coefficients change in a qualitative comparison, although we can (and successfully do) use these simulations to also verify that our estimates of |$\partial \boldsymbol {O}_\alpha /\partial \lambda _r$| are relatively stable.

Fig. 2 shows WST coefficients for simulations with beam smoothing and line broadening but without the high-pass transfer function |$\mathcal {T}_\mathrm{hp}$| or random Gaussian noise per voxel as described in Section 3.2. We specifically show the first-order WST coefficients S₁(j, ℓ; q = 1.5) and the second-order WST coefficients by their first-order counterparts S₂(j, |$j^{\prime }$|, ℓ; q = 1.5)/S₁(j, ℓ; q = 1.5), such that the latter are partially rescaled (noting that what we have termed the ‘rescaled’ |$\tilde{S}_2$| differs by additional factors of |$2^{j^{\prime }-j}$| and the as-yet-undetermined N₂).

Figure 2.

Variation of the median raw first-order WST coefficients S₁(j, ℓ; q = 1.5) (left-hand panels) and the partially rescaled second-order coefficients S₂(j, |$j^{\prime }$|, ℓ; q = 1.5)/S₁(j, ℓ; q = 1.5) (right-and panels), given noiseless CO line-intensity simulations with changes in model parameters A (upper panel sets), C (middle panel sets), and σ (lower panel sets). We show both the original values (upper panel in each panel set) and the change relative to the results for fiducial parameter values (lower panel in each panel set), as well as 68 per cent sample intervals (red shaded areas) around the median given the fiducial parameter values.

We observe a few interesting general trends for the WST coefficients even considering their values for the fiducial parameter values alone. First, note that both S₁(j, ℓ) and S₂(j, |$j^{\prime }$|⁠, ℓ) tend to decrease with higher multipole index ℓ. Thinking about the variation of these coefficients with ℓ as encoding how much brightness is in different kinds of shapes that fill out the cosmic web, this suggests more brightness in simpler shapes – i.e. the majority of observed CO luminosity comes from CO brightness ‘monopoles’ rather than ‘dipoles’, with more brightness contained in dipoles than ‘quadrupoles’, and so on.

Note also the amount of relative variation in the WST coefficients, which appears to trend upwards with higher values of j, i.e. larger scales. This is not too surprising, and is at least qualitatively analogous to how the relative variation in the power spectrum σ[P(k)]/P(k) scales simply as the inverse of the square root of the number of modes ∼k⁻¹ due to there simply being fewer large-scale modes to sample in a finite volume.

Most interesting, however, is how the first- and second-order coefficients respond differently to changes in different parameters. Changing C, the overall normalization of L(M_h), affects S₁ uniformly across all j and ℓ but essentially has zero effect on S₂/S₁. Meanwhile, shifts in the other two parameters affect both S₁ and S₂/S₁. The first-order coefficients S₁ are somewhat uniformly sensitive to changes in A – which we recall is the negative of the faint-end L(M_h) power-law exponent – and the reduced second-order coefficients S₂/S₁ are less sensitive. Meanwhile, changes in the lognormal scatter σ around the mean L(M_h) relation result in a strongly scale- and ℓ-dependent change in S₁(j, ℓ), and significant changes in S₂/S₁.

All of this highlights, the simultaneously rich and interpretable information that the WST coefficients contain. Since C is an overall normalization, shifting its value does not affect different scales or shapes of CO emission in any non-uniform way. Therefore, only S₁ is sensitive to it, while sparsity as measured by S₂/S₁ is unaffected. Meanwhile, shifting A down will decrease (and shifting it up will increase) the contribution of faint CO emitters to the total cosmological signal, which makes the signal both dimmer overall and sparser (or brighter and less sparse when shifting A up). This leads to both S₁ and S₂/S₁ changing, although perhaps with S₂/S₁ less sensitive.

Finally, shifting σ up increases the scatter in the halo mass–luminosity relation, leading potentially to significantly brighter but also more sparsely distributed emitters. In the power spectrum, this effect manifests as an increase in the shot noise contribution to the total P(k), which chiefly affects higher k or smaller scales where clustering becomes subdominant. We see the same thing here with the WST coefficients, which are more sensitive to changes in σ for lower j and smaller scales. But beyond shifts in S₁, we also see shifts in S₂/S₁, demonstrating the effect that the change in σ has not merely on the typical scale of CO fluctuations but on their morphology and sparsity. We also see increased sensitivity to σ with higher ℓ, apparently reflecting how increased scatter results in more luminosity being contained in more complex shapes.

All of this, however, is in the absence of observational effects. We must now consider the fact that real-world observations of CO line intensity will always be noisy, as well as being subject to various observational transfer functions like the |$\mathcal {T}_\mathrm{hp}$| that we described in Section 3.2. The WST of a noisy data cube is very different from that of a noiseless data cube, because spaces that would have been line-intensity voids will be filled with spurious noise fluctuations, and features that would have been line-intensity peaks are effectively smeared out in a WST by the background noise. This is in fact what we see with noisy simulations in Fig. 3. Where before noise only the first-order coefficients were sensitive to changes in C, now the reduced second-order coefficients S₂/S₁ are also sensitive. This makes sense because having brighter (or dimmer) CO contrast versus the noise level will in fact increase (or decrease) the sparsity of the observed temperature field.

Figure 3.

Same as Fig. 2, but only for the model parameter C, and given noisy CO line-intensity simulations with |$\mathcal {T}_\mathrm{hp}$| applied to the CO signal as described in Section 3.2. We show results for both COMAP Pathfinder simulations (upper panel sets) and COMAP-ERA simulations (lower panel sets). We show a noise-only baseline (dotted) in addition to WST coefficients in the presence of signal with different values of C.

Furthermore, in the case of COMAP Pathfinder simulations, the difference in S₂/S₁ from adding signal on top of noise is negligible compared to the variance due to the overall damping of sparsity in the observed temperature field. Therefore, in a first-generation CO line-intensity survey, only the first-order coefficients are likely to allow inferences about the CO field.

Decreasing the noise level from that expected for the COMAP Pathfinder survey to that expected for COMAP-ERA, we observe a corresponding decrease in the noise-only baseline (apparently by a factor of 22^q/2, which we might expect from the factor-of-22^1/2 decrease in the noise level), and potentially detectable deviations in S₂/S₁ from the noise-only baseline. Note, however, that the interaction of noise and signal is not simply additive, even for first-order coefficients (ostensibly most analogous to P(k)). Comparing the first-order coefficients in Fig. 2 to those in the COMAP-ERA case in Fig. 3, the noise seems to add to S₁ at smaller scales but then drag down the CO signal to the noise-only baseline with the effect being more marked on larger scales. This is also what Greig et al. (2022) appear to see for their two-dimensional S₁ coefficients at scales larger than the instrumental resolution. Given the non-linear, multiscale nature of the WST, we should not expect the interaction of signal and noise to be as simple as in a P(k) analysis, and ultimately differentiating signal from noise in a WST analysis may require an even stronger understanding of both.

4.2 Covariance of full and shapeless WST coefficient sets

Reducing strong correlations between the components of our observable vector is not merely an aesthetic matter. Undertaking likelihood analyses and Fisher forecasts requires credible computation of the inverse of the covariance matrix |$\bf{C}$|⁠, as one may observe from equation (15). High multicollinearity in the covariance matrix or other pathological behaviour will lead to the matrix inversion being an ill-conditioned problem, as measured by the condition number |$\kappa \equiv ||\bf{C}||\, ||\bf{C}^{-1}||$| of the matrix (the norm of the matrix⁸ times the norm of the inverse). We consider in Table 3 representative condition numbers for the full and shapeless WST coefficients across all values of q initially considered.

Table 3.

Base-2 logarithms of condition numbers of WST covariance matrices (which should be as low as possible for credible invertibility) estimated from the COMAP Pathfinder simulations in the 30–32 GHz frequency band including both signal and noise for different reductions of the WST coefficients (or lack thereof). We do not show the condition numbers for rescaled WST coefficients as we only define them properly for q = 1.5, but in that case we find the corresponding value is |$\log _2{(||\bf{C}||\, ||\bf{C}^{-1}||)}=18.0$|⁠.

q	|$\log _2{(||\bf{C}||\, ||\bf{C}^{-1}||)}$| for covariance of:
	raw WST	shapeless WST	s1 and s2 only
0.5	28.4	19.5	16.0
0.75	28.2	19.2	15.7
1.0	74.4	45.6	15.8
1.5	27.9	17.1	17.0
2.0	38.8	22.8	22.8

Ideally the condition number is as close to one as possible as effectively it represents the factor by which inversion (or solutions of other linear algebra problems involving this matrix in general) can result in loss of precision. 32-bit and 64-bit floats, respectively, have 24 and 53 bits of total significand precision, so the base-2 logarithm of the condition number should not exceed these numbers. Most modern CPUs and programming languages use 64-bit float arithmetic by default, but in many cases GPU programming conventions (as might be used in GPU-accelerated WST calculations) still favour 32-bit floats for throughput and performance. With Kymatio at time of writing, even the CPU-based WST implementation uses a 64-bit complex float data type (i.e. 32-bit floats for each part of the complex number) in calculating the spatial ‘integral’ (or sum), so the WST coefficients are encoded as 32-bit floats. This work will leave this choice unmodified, partly to make conservative use of memory and disk space, but partly because large numbers of WST evaluations will likely be necessary in future work and these may use GPU-accelerated calculations based on 32-bit float arithmetic. Therefore, WST coefficient sets that result in covariance matrices with a condition number above 2²⁴ (or even within a factor of 2–3 of this threshold) should be treated with extreme suspicion.

With q = 1, the raw WST covariance is absolutely unacceptable for inversion. It is unlikely that the fact of the condition number exceeding our threshold is a fluke or numerical artefact, given the large number of realizations used to generate the covariance matrix, but we must be cautious. As an additional check, we show in Fig. 4 covariance matrix condition numbers for additional values of q close to 1. The results demonstrate that as q approaches 1 from either direction, the condition number consistently ends up exceeding the acceptable range of values. That said, the condition number improves to an acceptable value when we discard information for ℓ = 0 (e.g. by discarding s₂₀(j) from the shapeless WST coefficients). This suggests that as q → 1, the ℓ = 0 coefficients (which we recall correspond to a degenerate wavelet with constant amplitude) lead to some combination of extreme ranges of covariance values and pathological correlations.

Figure 4.

Plot of the same condition number logarithms given in Table 3 as a function of q with additional calculations for extra values of q that approach 1 from either direction. The condition numbers show a clear trend of rising to unacceptably high values as q approaches 1 from either direction, unless we discard ℓ = 0 WST coefficients. We also show the approximate borderline of κ ∼ 2²⁴ (thick yellow band), above which inversion of the covariance matrix with 32-bit float arithmetic becomes unreliable, either because the matrix is ill-conditioned or because the matrix is outright singular.

With other values of q inversion of the raw WST covariance would be dangerous with 32-bit float arithmetic. By using the ℓ-averaged WST coefficients (although again, in the case of q = 1 we must discard s₂₀(j) altogether), the condition number of the covariance matrix decreases by around a factor of ∼10³, sufficient to bring it within the precision of 32-bit floats (albeit only by the equivalent of a couple of decimal digits). In the case of q = 2, however, the condition number remains at an unacceptably high level.

The high-condition numbers of the covariance matrices |$\bf{C}_{\alpha \beta }$| are not entirely due to correlations alone with the range of covariance values playing a large part also⁹. The latter effect can be discarded by examining the condition numbers of the correlation coefficient matrices |$\bf{R}_{\alpha \beta } = \bf{C}_{\alpha \beta }/\sqrt{\bf{C}_{\alpha \alpha }\bf{C}_{\beta \beta }}$|⁠, and are what the covariance matrices would be if all variables were normalized by their standard deviation, We show correlation coefficient matrices in Fig. 5 for selected q for both raw and shapeless WST coefficients. The condition numbers for these matrices are significantly lower than the covariance matrix condition numbers, but in the case of q = 2 the condition number of the correlation coefficient matrix for the full raw WST coefficient set is still in excess of 2³⁰. Overall, q = 2 appears to result in highly pathological correlations between the WST coefficients (not as easily understood as in the case of q → 1, and potentially requiring careful rescaling of coefficients beyond the scope of this work), and we must discard it as a viable choice of q for the purposes of this work. This is likely a consequence of the q = 2 WST tending to over-emphasize spurious peaks in noisy observations.

Figure 5.

Correlation coefficient matrices R_αβ for the raw WST coefficients (left-hand panels) and the shapeless (ℓ-averaged) WST coefficients (right-hand panels) derived from the COMAP Pathfinder simulations in the 30–32 GHz band including both signal and noise for three different values of the density-weighting exponent q (top to bottom, indicated in annotations above each panel).

4.3 Detection significance

We consider a kind of median expectation for the detection significance as follows. For each 2 GHz sideband i of the COMAP Pathfinder or COMAP-ERA observation, we have a numerical estimate of the covariance matrix |$\bf{C}_{i;\alpha \beta }$| of the power spectrum or the WST coefficients (indexing k-bins or WST coefficients with α and β here to avoid confusion with the indices j and |$j^{\prime }$| used for the WST scales). We also have a median expectation value – write it as 〈ΔP(k)〉 or {〈Δs₁(j)〉, 〈Δs₂₀(j)〉, 〈Δs₂(j, |$j^{\prime }$|)〉}, for example – for the difference between the measured power spectrum or WST coefficient set and the expected statistic values from noise alone (i.e. in the absence of any CO signal). Then for each sideband i, we can effectively consider the Mahalanobis distance as a measure of the significance of the deviation of the measurement from noise, if we consider the expectation from noise alone to be the presumed mean of the distribution.

For the power spectrum,

$$\begin{eqnarray*} d_{M,i} = \sum _{\alpha ,\beta }\langle \Delta P(k_\alpha)\rangle \bf{C}_{i;\alpha \beta }^{-1} \langle \Delta P(k_\beta)\rangle . \end{eqnarray*}$$

(16)

The total significance across all four sidebands is then

$$\begin{eqnarray*} \mathrm{\frac{S}{N}}\equiv \sqrt{\sum _i d_{M,i}}. \end{eqnarray*}$$

(17)

If the power spectrum covariance matrix is diagonal for all i such that |$\bf{C}_{i;\alpha \beta }=\sigma ^2_{P(k);\alpha }\delta _{\alpha \beta }$|⁠, then the detection significance is simply |$\sqrt{\sum _\alpha \left\langle \Delta P(k_\alpha)\right\rangle ^2/\sigma _{P(k);\alpha }^2}$|⁠. Furthermore, since the total observed power spectrum is the sum of the signal and noise power spectra, 〈ΔP(k_α)〉 should be the same as the expected value of the signal power spectrum by itself. All of this is consistent with previous detection significance forecasts for COMAP.

The computation for the WST coefficients (raw or reduced) is analogous, except the sum is now over the range of coefficients rather than the k-bins of the power spectrum; we omit an explicit expression for the WST S/N. As expected from the condition numbers for the covariance matrices, this calculation is outright impossible for the raw WST coefficients with q = 1 due to the ill-conditioned covariance matrix, and although we show the results for the raw coefficients with q = {0.5, 0.75, 1.5}, major caveats apply as to their credibility. As for the shapeless WST coefficients, their covariance matrices are credibly invertible and therefore, we show the results without the same caveats. We do not show results for the rescaled WST coefficients as they should be largely equivalent to those for the raw WST (modulo numerical instabilities).

As previously discussed, the covariance matrices could be estimated from either noise-only simulations or simulations that include both signal and noise. Given that the null hypothesis would be that our measurement is consistent with the distribution of P(k) or WST measurements expected from noise alone, it arguably makes sense to calculate detection significance against the noise-only covariance. However, to quantify information content, we should calculate detection significance against the signal-plus-noise covariance. We will present both calculations in Table 4.

Table 4.

Total expected detection significance of P(k) and WST coefficients against noise covariance (values in parentheses are computed against signal-plus-noise covariance) for the COMAP Pathfinder and COMAP-ERA simulations. Dashes indicate that the detection significance could not be calculated due to a singular covariance matrix.

Statistic	Computed S/N for:
	Pathfinder	COMAP-ERA
P(k)	11 (8)	240 (21)
Raw WST, q = 0.5	30 (24)	500 (96)
Shapeless WST, q = 0.5	12 (9)	184 (34)
Raw WST, q = 0.75	30 (22)	598 (87)
Shapeless WST, q = 0.75	12 (9)	197 (32)
Raw WST, q = 1.0	–	–
Shapeless WST, q = 1.0	12 (9)	206 (30)
Raw WST, q = 1.5	39 (32)	625 (61)
Shapeless WST, q = 1.5	13 (9)	265 (30)

Regardless of whether we use the raw or shapeless WST coefficients, the WST detection significance is generallly in excess of the P(k) detection significance. Based on the Pathfinder computations, the optimal value of q would appear to be 1.5. This is still true when looking at COMAP-ERA detection significance against noise alone, which increases monotonically with q. However, for COMAP-ERA, the detection significance calculated against the total (signal-plus-noise) covariance matrices appears to decrease with greater q, perhaps a sign that information in CO underdensities becomes more apparent with lower noise levels. Overall, however, we will proceed with q = 1.5, a configuration that still results in the WST apparently holding several times more information than the power spectrum.

4.4 Parameter constraints

Now fixing q = 1.5 based on the above results, we consider Fisher forecasts for parameter constraints with the formalism established in Section 3.4. The observables considered are the power spectrum P(k) only, the shapeless WST coefficients, the rescaled WST coefficients, and the combination of P(k) with the VID.

We now define the additional normalizations N₁ and N₂ required to fully define the rescaled coefficients |$\tilde{S}_1$| and |$\tilde{S}_2$|⁠. As might be gathered from Fig. 3, the rescaling required to optimize covariance conditioning differs between the noise-dominated COMAP Pathfinder scenario and the significantly less noise-dominated COMAP-ERA scenario. For COMAP Pathfinder simulations, we use the following definitions:

$$\begin{eqnarray*} N_1(j,\ell) &= 1; \end{eqnarray*}$$

(18)

$$\begin{eqnarray*} N_2(j,j^{\prime },\ell) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}2^{J-j}&\text{if }\ell > 0,\\ 1&\text{if }\ell =0.\end{array}\right. \end{eqnarray*}$$

(19)

For the COMAP-ERA scenario we instead use the following:

$$\begin{eqnarray*} N_1(j,\ell) &= (\ell +1/2)2^{-j/2}; \end{eqnarray*}$$

(20)

$$\begin{eqnarray*} N_2(j,j^{\prime },\ell) &= (\ell +1/2). \end{eqnarray*}$$

(21)

Recall also that as part of defining the rescaled WST coefficients, we discard |$\tilde{S}_2$| with |$j^{\prime }$| > j + 1 and ℓ = 0 due to very high correlation of these coefficients with |$\tilde{S}_2(j,j^{\prime }=j+1,\ell =0)$|⁠.

The above definitions of N₁ and N₂ lower the base-2 logarithm of the condition numbers of the covariance matrices from ≈28 to ≈18–19 for COMAP Pathfinder covariances, or a slightly higher range of 19–20 for COMAP-ERA. In all cases, our condition numbers are one order of magnitude better than the 10⁷ maximum imposed by Eickenberg et al. (2022), which corresponds to ≈2²³ (just barely within 32-bit float precision).

We show the Fisher forecast results in two different ways: we show 68 per cent and 95 per cent credibility ellipses in Fig. 6, and marginalized parameter errors in Table 5. These results in the COMAP Pathfinder case should not be directly compared to forecasts for the Pathfinder survey in Chung et al. (2022). That work uses a Monte Carlo method and fully incorporates non-Gaussian and skewed parameter covariances in this way, which likely significant affects errors for A. We also fix B and M and thus it is natural to obtain tighter constraints on the remaining parameters, and this is likely a significant effect for C, which is highly degenerate with M in the modelling and forecasts of Chung et al. (2022).

Figure 6.

68 and 95 per cent credibility ellipses from the Fisher forecasts for summary statistics (indicated in legends) from simulations of the COMAP Pathfinder (upper panels) and COMAP-ERA (lower panels).

Table 5.

Marginalized parameter constraining power predicted by Fisher forecasts using various observables as described in the main text. The final row provides constraints from only |$\tilde{S}_1(q,\ell > 0)$| and |$\tilde{S}_2(q,q^{\prime },\ell > 0)$|⁠.

Statistic	COMAP Pathfinder (COMAP-ERA) 1σ error for:
	A	C	σ
P(k) only	0.130 (0.053)	0.042 (0.025)	0.119 (0.030)
Shapeless WST	0.126 (0.068)	0.041 (0.016)	0.095 (0.019)
P(k) & VID	0.100 (0.052)	0.034 (0.015)	0.091 (0.019)
Rescaled WST	0.014 (0.009)	0.010 (0.008)	0.015 (0.008)
–withℓ > 0only	0.021 (0.017)	0.015 (0.010)	0.027 (0.012)

Overall, we find that the WST has the fundamental sensitivity to allow parameter constraints finer by a factor of 2–6 compared to with the combination of P(k) with the VID. However, note that the constraining power of the shapeless WST coefficients is comparable to the P(k)–VID combination. Also of note is the change in parameter constraints between the COMAP Pathfinder and COMAP-ERA scenarios. Whereas the improvement in marginalized parameter errors does not reach a factor of 2 for the rescaled WST coefficients, it reaches a factor of 2–4 for all other summary statistics used. This appears to match the relative changes in detection significance computed for P(k) and the WST (both raw and shapeless coefficient sets) against the signal-plus-noise covariance in Section 4.3.

We also consider constraints with the full rescaled WST coefficients but only using the subset of 45 |$\tilde{S}_1$| and |$\tilde{S}_2$| values with ℓ > 0, i.e. only using the information content in the WST that derives from non-monopole local CO emission gradients. We find that the loss of ℓ = 0 information results in parameter constraints that are somewhat weaker than with the full rescaled WST coefficient set, but none the less stronger than those derived from the P(k)–VID combination.

5 DISCUSSION

5.1 The shapeless WST coefficient set: key takeaways

The reduction of the full WST coefficient to the shapeless set is somewhat sensible in a few different ways. The key motivation is the naïve expectation that we do not gain signficant information from distinguishing between WST coefficients of different spherical harmonic multipoles ℓ, and that the majority of information content in the WST is in variations of s₁(j), s₂(j, |$j^{\prime }$|⁠), and s₂₀(j) with scales j and |$j^{\prime }$|. We also observe an important numerical advantage – likely helped by reduced dimensionality – that correlation coefficient and covariance matrices are much better conditioned than for the raw WST, and without significant fine-tuning as was required for the rescaled WST coefficients.

However, subsequent calculations of detection significance and parameter constraints in the latter half of Section 4 clearly show that the shapeless WST coefficients do not capture the full information in the WST. In fact, the similarity of constraining power between the P(k)–VID combination and the shapeless WST coefficients should be striking. Neither P(k) nor the VID will be sensitive to the ‘shape’ of CO fluctuations – the addition of the VID provides non-Gaussian but shape-free information. Therefore, it is natural that we do no better than the P(k)–VID combination once we discard variations with ℓ and thus morphological information from the WST.

5.2 The full WST coefficient set: practical considerations

The fundamental sensitivity of the WST in a COMAP Pathfinder or COMAP-ERA observational scenario is extremely promising. In principle, deviations of the WST from the noise-only baseline can be detected at considerably higher significance than the CO P(k), and parameter constraints with the full rescaled WST coefficients are correspondingly tighter by an order of magnitude in the COMAP Pathfinder case. They are only tighter by a factor of 3–5 in COMAP-ERA simulations, but this is likely conservative as making use of additional values of the density-weighting exponent q could further improve WST-based constraints. With lower noise in the COMAP-ERA scenario, CO voids should become more readily apparent, and evaluations of WST coefficients with q = 0.5 or q = 0.75 could provide additional information content. Indeed, Eickenberg et al. (2022) are apparently able to use two or even five different values of q at once, although their work only calculates summary statistics including the WST on noiseless matter density fields.

However, before even considering using multiple values of q at once, LIM experiments should consider whether they can make practical use of the WST with even just one value of q. Our results are strongly contingent on the ability to understand the interaction between signal and noise in the WST. Our use of a beam-smoothed, high-pass filtered CO field summed with a uniform Gaussian field mirrors the approximate input for a pseudo-P(k) analysis as undertaken by COMAP (Ihle et al. 2022). However, non-uniform noise residuals from removal of systematics and foregrounds may eventually manifest above uniform noise with greater integration times. The interaction between the CO signal and uniform noise is already quite complex, and only further complexity can arise from systematics like standing waves in the instrument or other pathological variations in the noise level across frequency, making the WST analysis more tenuous in practical terms even beyond mode mixing considerations.

Making matters even more tenuous is the fact that the WST does not have an established extension to cross-correlations as the power spectrum does. The pseudo-cross P(k) analysis of Ihle et al. (2022) is made more robust by the fact that the result is not generated from a single co-added map, but from cross-correlation power spectra between various pairs of pseudo-temperature maps generated from different data splits. This methodology intrinsically rejects disjoint systematics and uncorrelated noise between the splits, thus improving the estimate of the underlying CO power spectrum. Such cross-correlation power spectra are mathematically well motivated since the auto-correlation is merely a special case of the cross-correlation spectrum. Meanwhile, the key insight behind generalization of the VID to cross-distribution statistics (Breysse et al. 2019) is the straightforward convolution of signal and noise distributions and the ability to deconvolve uncorrelated noise out of distribution shapes in Fourier space. Due to the non-linear nature of the WST, gaining similar insights and defining cross-field extensions may prove difficult.

The ease with which P(k) by contrast can be applied to cross-correlations is key not only for mitigation of systematics within a single experiment, but also for cross-correlations between independent measurements of the same cosmic web for multitracer astrophysics on cosmological scales. The COMAP survey has been designed for cross-correlations with spectroscopic galaxy surveys in particular the Hobby–Eberly Telescope Dark Energy eXperiment (HETDEX; Hill et al. 2008, 2021; Gebhardt et al. 2021) surveying Lyman-alpha emitters (LAE) at redshifts overlapping with COMAP target redshifts. Fisher forecasts already predict approximately a factor-of-two improvement in determination of the mean CO temperature-bias product with the CO–LAE cross-power spectrum versus the CO auto P(k). Certainly the CO ‘auto’ WST promises even further improvement but at the cost of reduced robustness compared to a pseudo-cross P(k) analysis and no readily apparent path to a cross-field extension.

On top of all of this, the WST covariance does not appear to be analytically understood, given that it is numerically estimated in every work that has been cited in this paper (and in this paper as well). The non-linear, multiscale nature of the WST appears to be both an advantage in the resulting rich morphological and non-Gaussian information content, and a severe disadvantage in complicating covariance considerations and thus affecting the ability to extract desired constraints from all of that information content.

6 CONCLUSIONS

With our results, we have answered the questions that we set out in the Introduction (reproduced below):

• What is the detectability of cosmological CO emission in the WST coefficients derived from a simulated COMAP observation? The fundamental sensitivity of the COMAP Pathfinder should allow for the detection of deviations of the full WST coefficient set from noise. The predicted overall signal-to-noise ratio of 30–40 with the Pathfinder survey is several times what is expected for the CO P(k). The detectability increases even further with the lower noise level expected from COMAP-ERA, although information content only increases by the same relative amount as for P(k).

• Do parameter constraints from these coefficients improve significantly over the previously studied combination of P(k) with the VID? Fisher forecasts suggest that the full (rescaled) WST coefficient set allows for parameter constraints at least 2–5 times finer than the P(k)–VID combination. This potential advantage is specifically due to the additional shape information provided across the basis of solid harmonic wavelets with varying ℓ, and could potentially be magnified by use of different values of the density-weighting exponent q.

These exploratory results are highly encouraging and seem to corroborate that wavelets are a more optimal basis than Fourier modes for probing LSS, even as traced by CO. However, as we have discussed above, the covariance of the WST still needs to be better understood. Illustrating this need is our discovery that for the solid harmonic WST applied to simulated noisy data, the covariance matrix of the full coefficient set becomes singular as the WST hyperparameter q goes to 1. For the forecasts in this work, appropriate choices of q and rescaled or reduced coefficient sets kept covariances well-conditioned and avoided such pathological behaviour. None the less, this singularity, which appears to be newly found in this work, has clear important implications even outside the context of CO LIM. A full understanding of this and other aspects of WST covariances is vital to any real-world likelihood analyses that would make use of the WST, given strong inherent correlations between WST coefficients and given the complex interaction of signal and noise.

Perhaps even more important is the design of a cross-correlation extension of the WST, without which it is not possible to design an analogue to a pseudo-cross P(k) analysis that rejects uncorrelated systematics and noise in LIM data. As previously mentioned, Breysse et al. (2019) already formulate a VID extension to cross-correlate galaxy counts with 21-cm intensity measurements, and other work is already extending this to make use of joint PDFs between different line-intensity fields (Breysse et al., in prep.). Hopefully, similar future extensions to the WST will allow robust auto- and cross-correlation measurements of the morphological and non-Gaussian information content that seems uniquely accessible to a WST-type LIM analysis.

ACKNOWLEDGEMENTS

Thanks to Dick Bond, Jonathan Braden, and Patrick Breysse for various wavelet- and/or WST-related conversations leading up to and during the writing of this paper. Thanks also to Veronica Stewart for a conversation that helped clarify the physical chemistry context underlying the original development of the solid harmonic WST by Eickenberg et al. (2018).

It is a pleasure to acknowledge the COMAP collaboration for allowing use of approximate noise levels and transfer functions to inform simulation parameters. Many thanks to George Stein for running and making available the original peak-patch simulations for Ihle et al. (2019) that are reused in this work (and in Chung et al. 2022). Thanks also go to an anonymous referee whose constructive comments improved this manuscript.

Research in Canada is supported by NSERC and CIFAR. Parts of these computations were performed on the GPC supercomputer at the SciNet HPC Consortium, and on workstations hosted at CITA. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Foundation Fund – Research Excellence; and the University of Toronto.

DTC is supported by a CITA/Dunlap Institute postdoctoral fellowship. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. The University of Toronto operates on the traditional land of the Huron-Wendat, the Seneca, and most recently, the Mississaugas of the Credit River; DTC is grateful to have the opportunity to work on this land. This research made use of Astropy,¹⁰ a community-developed core Python package for astronomy (Astropy Collaboration et al. 2013, 2018). This research also made use of NASA’s Astrophysics Data System Bibliographic Services.

DATA AVAILABILITY

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

Footnotes

1

The covariance in practice will not be diagonal as observational effects could certainly introduce mode coupling. Furthermore, off-diagonal terms absolutely must be considered against higher-order multipoles of |$P(\boldsymbol {k})$|⁠.

2

See also Eickenberg et al. (2017) for a similar overview of the solid harmonic WST in its original motivating context.

3

Specifically, Eickenberg et al. (2018) provide the example of Reisert & Burkhardt (2009), whose work considers harmonic filters in the context of three-dimensional feature detection with a demonstrative example using microscopic imagery of pollen grains.

4

This will also happen to sidestep potential issues with normalization in implementation – note for example that the functional form of the mother wavelet as defined by Valogiannis & Dvorkin (2022a), Valogiannis & Dvorkin (2022b) does not include division by σ³ as we do, despite their use of σ = 0.8.

5

For further details, we refer the interested reader to the review of Carilli & Walter (2013); the works of Chung et al. (2021) and Chung et al. (2022), from which we take the form of our model; and references in each work. Again, the state of physical understanding and observational constraints around high-redshift CO line emission is unfortunately well beyond the scope of the present work.

6

https://github.com/georgestein/limlam_mocker

7

For more on the topic of Fisher analyses, Section 11.4 of Dodelson (2003) for a pedagogical overview and Coe (2009) for additional guidance in practice.

8

The numpy implementation and thus this work both use the 2-norm in this context. With the 2-norm, the condition number is simply the ratio of the largest singular value to the smallest singular value. For a symmetric positive definite matrix, this is simply the ratio of the largest eigenvalue to the smallest eigenvalue. For further insights we refer the reader to the text of Golub & van Loan (1996).

9

This is clear when recalling the statement from the previous footnote that with the 2-norm, the condition number is simply the ratio of the covariance matrix’s largest and smallest singular values (or eigenvalues if the covariance is positive definite).

10

http://www.astropy.org