Using X-Ray Morphological Parameters to Strengthen Galaxy Cluster Mass Estimates via Machine Learning

TL;DR

This paper introduces a machine learning method that uses X-ray morphological parameters to improve galaxy cluster mass estimates, reducing scatter by 20% compared to traditional methods, aiding cosmological research.

Contribution

The study demonstrates that incorporating morphological parameters into a random forest model significantly enhances galaxy cluster mass predictions from X-ray data.

Findings

20% reduction in mass estimation scatter with morphological parameters

Morphological features like smoothness and asymmetry are most important

Method effective across idealized and realistic mock observations

Abstract

We present a machine learning approach for estimating galaxy cluster masses, trained using both Chandra and eROSITA mock X-ray observations of 2,041 clusters from the Magneticum simulations. We train a random forest regressor, an ensemble learning method based on decision tree regression, to predict cluster masses using an input feature set. The feature set uses core-excised X-ray luminosity and a variety of morphological parameters, including surface brightness concentration, smoothness, asymmetry, power ratios, and ellipticity. The regressor is cross-validated and calibrated on a training sample of 1,615 clusters (80% of sample), and then results are reported as applied to a test sample of 426 clusters (20% of sample). This procedure is performed for two different mock observation series in an effort to bracket the potential enhancement in mass predictions that can be made possible by including dynamical state information. The first series is computed from idealized Chandra-like mock cluster observations, with high spatial resolution, long exposure time (1 Ms), and the absence of background. The second series is computed from realistic-condition eROSITA mocks with lower spatial resolution, short exposures (2 ks), instrument effects, and background photons modeled. We report a 20% reduction in the mass estimation scatter when either series is used in our random forest model compared to a standard regression model that only employs core-excised luminosity. The morphological parameters that hold the highest feature importance are smoothness, asymmetry, and surface brightness concentration. Hence, these parameters, which encode the dynamical state of the cluster, can be used to make more accurate predictions of cluster masses in upcoming surveys, offering a crucial step forward for cosmological analyses.

Tables

Table 1. Table 1: Pearson correlation between each observable in the model and (i) the true mass, log ⁡ ( M 500 ​ c ) subscript 𝑀 500 𝑐 \log(M_{500c}) , or (ii) the logarithmic mass residual from a mass-luminosity regression, ℛ ℛ \mathcal{R} . These calculations were performed using both the “idealized Chandra ” and “realistic eROSITA ” series of morphological parameters. In both series, S 𝑆 S and A 𝐴 A correlate most strongly with log ⁡ ( M 500 ​ c ) subscript 𝑀 500 c \log(M_{\mathrm{500c}}) . The correlations with mass are generally stronger in the “realistic eROSITA ” series. While log ⁡ ( w ) 𝑤 \log(w) and log ⁡ ( P 10 ) subscript 𝑃 10 \log(P_{10}) correlate most strongly with ℛ ℛ \mathcal{R} in both series, we find that they are not the most important morphological parameters (rather, S 𝑆 S , A 𝐴 A , and c 𝑐 c are).

	Chandra				eROSITA
	Correlation $r$ with				Correlation $r$ with
Feature	$\log \left(M_{500c}\right)$		$ℛ$		$\log \left(M_{500c}\right)$	$ℛ$
$\log \left(L_{\mathrm{ex},z}\right)$	0\@alignment@align.927	0.000	0\@alignment@align.929	0.000
$c$	-0\@alignment@align.101	-0.041	0\@alignment@align.208	-0.146
$e$	0\@alignment@align.060	-0.105	0\@alignment@align.071	-0.111
$\log \left(w\right)$	-0\@alignment@align.096	0.195	-0\@alignment@align.126	0.228
$\log \left(P_{10}\right)$	-0\@alignment@align.170	0.212	-0\@alignment@align.197	0.224
$\log \left(P_{20}\right)$	-0\@alignment@align.140	0.155	-0\@alignment@align.308	0.194
$\log \left(P_{30}\right)$	-0\@alignment@align.139	0.146	-0\@alignment@align.427	0.198
$\log \left(P_{40}\right)$	-0\@alignment@align.145	0.149	-0\@alignment@align.489	0.207
$A$	-0\@alignment@align.294	0.034	-0\@alignment@align.654	0.131
$S$	-0\@alignment@align.493	0.048	-0\@alignment@align.795	0.119
$M_{20}$	0\@alignment@align.032	0.116	-0\@alignment@align.148	0.148

Table 2. Table 2: The 1 ​ σ 1 𝜎 1\sigma percentile intrinsic scatter, defined as half of the 16th–84th percentile range of the mass residuals ℛ ℛ \mathcal{R} , for each model as trained in the text when applied to the test set, computed using the “idealized Chandra ” and “realistic eROSITA ” mock observation series, denoted δ iC subscript 𝛿 iC \delta_{\mathrm{iC}} and δ re subscript 𝛿 re \delta_{\mathrm{re}} , respectively. The multivariable linear models (OLR, RR, LR) improve only marginally relative to the mass-luminosity relation, with regularization yielding no improvement. The RF models, which capture nonlinear relationships between the input features and the mass, are able to further reduce the scatter beyond any linear approach.

Method	${\delta }_{\mathrm{iC}}$		${\delta }_{\mathrm{re}}$
$M-L_{\mathrm{ex},z}$	0\@alignment@align.081	0.081
OLR	0\@alignment@align.078	0.080
RR	0\@alignment@align.078	0.080
LR	0\@alignment@align.079	0.081
RF, defaults	0\@alignment@align.070	0.067
RF with CV	0\@alignment@align.066	0.066
RF, only $L_{\mathrm{ex},z}$, $c$, $A$, $S$	0\@alignment@align.070	0.071

Table 3. Table 3: The feature importances for the cross-validated RF models, computed based on the mean decrease in impurity. The remaining morphological parameters are omitted from the table, as their importances are all ≲ 1 % less-than-or-similar-to absent percent 1 \lesssim 1\% . The smoothness S 𝑆 S , asymmetry A 𝐴 A , and surface brightness concentration c 𝑐 c encode the most additional mass information after the core-excised luminosity L ex , z subscript 𝐿 ex z L_{\mathrm{ex,z}} .

	Importance (%)
Feature	Chandra		eROSITA
$L_{\mathrm{ex},z}$	74\@alignment@align.9	76.2
$S$	10\@alignment@align.0	14.6
$A$	4\@alignment@align.8	2.5
$c$	3\@alignment@align.9	1.5

Equations

$$L_{\mathrm{ex},\mathrm{z}}\equiv L_{X,\mathrm{ex}}E(z{)}^{-7/3}=\frac{N_{\mathrm{tot}}-N_{\mathrm{ce}}}{N_{\mathrm{tot}}}{L}_{X}E(z{)}^{-7/3}.$$

$$\mathcal{R}=\log (M_{500c,\mathrm{pred}})-\log (M_{500c,\mathrm{true}}),$$

Full text

galaxies: clusters: general – X-rays: galaxies: clusters – methods: statistical

††software: SCIKIT-LEARN (Pedregosa et al., 2011), XSPEC (Arnaud, 1996)

1 Introduction

Galaxy clusters are the largest gravitationally bound objects in the universe. They are rare, with masses $\gtrsim 10^{14}\,M_{\odot}$, and their abundance is sensitive to the underlying cosmological model. Cluster counts can be used to constrain cosmological parameters, provided that there is an accurate way to connect the cluster observables (such as X-ray luminosity or temperature) to the underlying dark matter halo mass (for a recent review see Pratt et al., 2019).

Recent cluster-based constraints are in tension with Planck cosmic microwave background (CMB) cosmological constraints. For example, Sunyaev-Zeldovich (SZ; Sunyaev & Zeldovich, 1972) surveys find fewer massive clusters than would be expected from the Planck fiducial cosmology (e.g., Planck Collaboration et al., 2016). This tension could be explained by a mass bias — a systematic under-estimation of X-ray based cluster mass estimates based on the hydrostatic assumption at the level of $30-45\%$ (Planck Collaboration et al., 2016; Bolliet et al., 2018; Zubeldia & Challinor, 2019; Makiya et al., 2019). However, a significant mass bias remains controversial. First, hydrodynamical cosmological simulations predict a hydrostatic mass bias in the range of $15-40$% (e.g. Rasia et al., 2006; Nagai et al., 2007; Lau et al., 2013; Nelson et al., 2014b; Shi et al., 2016; Biffi et al., 2016; Henson et al., 2017) due to non-thermal pressure support provided by bulk and turbulent gas motions (e.g. Lau et al., 2009; Nelson et al., 2014a; Shi et al., 2015) and temperature inhomogeneities in the intracluster medium (ICM) (Rasia et al., 2014). Recent observational results agree that the hydrostatic bias must be small, at least for relaxed systems (e.g. Applegate et al., 2016; Eckert et al., 2019; Ettori et al., 2019; Ghirardini et al., 2019). Second, the hydrostatic mass bias may also arise from the instrument-dependent systematic uncertainties in X-ray temperature measurements (Schellenberger et al., 2015; Israel et al., 2015). Finally, some cluster- and large-scale structure-based efforts put constraints on cosmological parameters that are consistent with those from the CMB (e.g., Mantz et al., 2015b; de Haan et al., 2016; DES Collaboration et al., 2017) while others are in tension with them (e.g., Hildebrandt et al., 2018; Joudaki et al., 2019; Ntampaka et al., 2019a). Given the importance of this problem, concerted efforts are underway to calibrate the cluster mass scales using optical weak lensing measurements of background galaxies (e.g. von der Linden et al., 2014; Hoekstra et al., 2015; Applegate et al., 2016; Dietrich et al., 2019) and CMB lensing (e.g., Raghunathan et al., 2019).

With next-generation observational surveys, such as the eROSITA X-ray survey (Merloni et al., 2012), soon to come online, massive data releases that will offer immense cosmological model constraining power are just around the corner. The eROSITA survey is predicted to identify ${\sim}$93,000 galaxy clusters at or above the 50 photon limit with $M\gtrsim 10^{13.7}\,h^{-1}\,\mathrm{M_{\odot}}$ (Pillepich et al., 2012, 2018). The product of spectral temperature and gas mass, $Y_{X}$, is one of the lowest scatter mass proxies (Kravtsov et al., 2006). However, many of the eROSITA observations will be in the regime of low-photon counts, making $T_{X}$- and $Y_{X}$-based cluster mass estimates inaccessible (Borm et al., 2014). The core-excised luminosity ($L_{X,\mathrm{ex}}$) is another lower-scatter mass proxy that does not require $T_{X}$ measurements; excluding the still poorly understood cluster cores ($r\lesssim 0.15\,R_{\mathrm{500c}}$) reduces the scatter in the $Y_{X}$ mass-$L_{X,\mathrm{ex}}$ (Maughan, 2007; Pratt et al., 2009) and weak lensing mass-$L_{X,\mathrm{ex}}$ (Mantz et al., 2018) relationships, but does so at the expense of drastically reducing the photon statistics.

Methods that provide improvements to $L_{X}$-based mass estimates for these low-photon eROSITA clusters could have a steep payoff. Even in the low-signal regime, there are subtle observable signals that can offer key insights for improving cluster mass estimates. Measures of cluster morphology, including surface brightness concentration (e.g., Santos et al., 2008), centroid shift (e.g., Mohr et al., 1993), and morphological composite parameters (e.g., Rasia et al., 2013), provide additional information about a cluster’s dynamical state (Mantz et al., 2015a), which has been shown to influence the scatter in the mass-$T_{X}$ relationship of simulated clusters (Ventimiglia et al., 2008), the correlated scatter in the relationship between weak lensing mass and integrated SZ Compton parameter $Y_{\mathrm{sph}}$ (e.g., Angulo et al., 2012; Marrone et al., 2012; Shirasaki et al., 2016), and the probability that a cluster is observed (Eckert et al., 2011; Planck Collaboration et al., 2011; Lovisari et al., 2017).

Modern machine learning (ML) techniques have been shown to reduce error in mass estimates of galaxy clusters. The techniques that have been developed use cluster dynamics (Ntampaka et al., 2015, 2016; Ho et al., 2019), X-ray images (Ntampaka et al., 2019b), and multiple wavelength summary statistics (Armitage et al., 2019; Cohn & Battaglia, 2019) as input; similar ML techniques have also been applied to less-massive galaxy groups (Calderon & Berlind, 2019; Man et al., 2019). These methods hinge on using ML to extract additional information from complex correlations in the mass-observable relationships. Here, we use ML to take advantage of the complex correlations among morphological parameters, dynamical state, and cluster mass to improve mass estimates.

Our new X-ray cluster mass measurement technique utilizes cluster dynamical state information, encoded in X-ray morphological parameters, to provide improved, lower-scatter mass estimates relative to a mass-luminosity linear regression. In addition, we demonstrate that this improvement is obtained even in low-photon count eROSITA observations, which makes the inclusion of dynamical state information a promising avenue for future cosmological analyses that depend on robust cluster mass estimates. In Section 2, we introduce the Magneticum simulations and mock Chandra and eROSITA X-ray observations of simulated galaxy clusters used in this work. In Section 3, we provide an overview of the X-ray morphological parameters employed as features in our models. In Section 4, we describe the preprocessing of the mock catalog data and several regression methods used to build our models. We summarize the results of our models in Section 5, followed by our conclusions and proposed follow-up work in Section 6.

Throughout this paper, the WMAP7 $\Lambda$CDM cosmology (Komatsu et al., 2011) is used: $\Omega_{\mathrm{m}}=0.272$, $\Omega_{\Lambda}=0.728$, $\Omega_{\mathrm{b}}=0.046$, $h=0.704$, $\sigma_{8}=0.809$, and $n_{\mathrm{s}}=0.963$. The base-10 logarithm is denoted by $\log$. All errors are quoted at the 68% level. The majority of this work is performed using the scikit-learn (Pedregosa et al., 2011) Python package.

2 Hydrodynamical Simulations

2.1 The Magneticum Simulations

Our cluster catalog is built from the Magneticum111www.magneticum.org/ (Dolag et al., 2015, 2016; Ragagnin et al., 2017) suite of cosmological hydrodynamical simulations. Magneticum uses a WMAP7 cosmology (Komatsu et al., 2011) with a range of baryonic physics included. For additional details about the simulations and the included baryonic physics, see, e.g., Biffi et al. (2013); Steinborn et al. (2015); Teklu et al. (2015); Steinborn et al. (2016); Bocquet et al. (2016); Remus et al. (2017).

We select clusters from the Magneticum Box2 and Box2b high-resolution simulations, selected for having sufficient resolution and volume to produce a suitable cluster catalog. Box2 has cubic side length of $352h^{-1}\,\mathrm{Mpc}$ with a dark matter particle resolution of $M_{\mathrm{dm}}=6.9\times 10^{8}\,h^{-1}\,\mathrm{M_{\odot}}$ and halo catalogs at $z=0.10$, $0.14$, $0.17$, $0.21$, $0.25$, and $0.29$ (as well as higher $z$, but these are not included in our analysis). Box2b is larger in volume ($640h^{-1}\,\mathrm{Mpc}$ on a side), has identical mass resolution, and has cluster catalogs at $z=0.25$ and $0.29$.

We initially select all clusters according to their spherical overdensity masses, $M_{\mathrm{500c}}$,222We define $M_{\mathrm{500c}}$ as the mass enclosed within a sphere of (comoving) radius $R_{\mathrm{500c}}$ whose mean density is 500 times the critical density of the universe at $z=0$. determined using the SUBFIND algorithm (Springel et al., 2001; Dolag et al., 2009). All clusters above $10^{13.5}\,h^{-1}\,\mathrm{M_{\odot}}$ are initially included and then subsampled in order to limit the sample to $\leq 230$ clusters per $0.1$ dex mass bin. The resulting training catalog has a flat mass function at lower masses, which helps to eliminate mass dependence in the scatter. Above ${\sim}10^{14.2}\,h^{-1}\,\mathrm{M_{\odot}}$, the mass function of this sample falls off, following the mass function of the simulation (Bocquet et al., 2016). Hence, the sample has a flat mass function in the range $10^{13.5}\leq M_{\mathrm{500c}}/(h^{-1}\,\mathrm{M_{\odot}})\leq 10^{14.2}$ and a falling mass function in the range $10^{14.2}\leq M_{\mathrm{500c}}/(h^{-1}\,\mathrm{M_{\odot}})\leq 10^{14.8}$ (see Fig. 1). Redshifts in the range $0.1\leq z\leq 0.21$ are roughly equally represented, with ${\sim}$300 clusters per redshift. However, our sample contains ${\sim}$450 clusters at $z=0.25$ and $z=0.29$ due to the addition of Box2b clusters.

The final cluster sample includes a total of 2,041 clusters in the redshift range $0.1\leq z\leq 0.29$, consisting of clusters from both Box2 and Box2b. Within this sample, there are 984 unique clusters, many of which are observed at multiple redshifts. Based on the assumption of self-similarity (Kravtsov & Borgani, 2012), we verify that the distributions of all relevant features included in the model exhibit minimal redshift evolution, justifying our inclusion of multiple snapshots for a particular cluster.

2.2 Mock Observations

From the cluster catalog, we create mock Chandra and eROSITA observations, employing the PHOX algorithm (Biffi et al., 2012, 2013). PHOX models the ICM thermal emission from gas particles by computing the expected number of photons given a fiducial (and large) exposure time and collecting area. The photon energies are then projected onto the sky plane along a chosen line of sight and cosmologically redshifted. A foreground galactic absorption model is applied, and models for Chandra ACIS-I and eROSITA are used to simulate the actual detections. Further details of the Magneticum implementation of PHOX can be found in Biffi et al. (2012), Biffi et al. (2013), and the publicly available Magneticum Cosmological Web Portal (Ragagnin et al., 2017).333https://c2papcosmosim.uc.lrz.de/

This implementation of PHOX allows the user to select from a number of parameters. For all observations, we select the ICM-only setting (i.e., AGN are not included as point sources in this work) and employ a 10 Mpc image line-of-sight size to include all relevant correlated structure. We seek to quantify the level of improvement in cluster mass estimates that can be made possible by incorporating dynamical state information, first in an idealized scenario and then in a realistic case that will be consistent with the observations made in upcoming large, high-throughput surveys such as eROSITA. To this end, our analysis features two different mock observation series: (i) mocks with Chandra-like angular resolution (“idealized Chandra” for short) with a Chandra ACIS-I instrument area and field of view ($2071\times 2071$ pixels, 16.9’ FoV, 0.49” pixel) and a 1 Ms observing time, in the idealized regime of a flat effective area with respect to photon energy ($600\,\mathrm{cm}^{2}$) and no point spread function (PSF) smearing, as well as (ii) “realistic eROSITA” observations with an eROSITA instrument area and field of view ($384\times 384$ pixels, $1.03^{\circ}$ FoV, 9.7” pixel) with a 2 ks observing time (Merloni et al., 2012) and instrument response and PSF modeled (see Ragagnin et al. (2017) for further details regarding eROSITA instrument modeling). To more closely imitate the conditions of the upcoming eROSITA observations, the “realistic eROSITA” mock images also include background noise. The process by which this noise is added is described below. The eROSITA mock observations have a median photon count of ${\sim}$2000 for clusters observed at $z=0.1$ and ${\sim}$100 for clusters observed at $z=0.29$. In contrast, the Chandra mocks have a median photon count of ${\sim}6\times 10^{5}$ for clusters observed at $z=0.1$ and ${\sim}3\times 10^{4}$ for clusters observed at $z=0.29$; clearly, derivative quantities computed from the “idealized Chandra” observations will be affected much less by Poisson noise.

The cluster bolometric luminosities $L_{X}$ are calculated by PHOX using the publicly available X-ray package XSPEC (Arnaud, 1996). Core-excised luminosities $L_{X,\mathrm{ex}}$ are computed as follows: (i) compute the total observed photon count $N_{\mathrm{tot}}$ within $R_{\mathrm{500c}}$, (ii) compute the observed photon count within $0.15R_{\mathrm{500c}}$, denoted $N_{\mathrm{ce}}$, and (iii) scale the bolometric luminosity by the ratio of the photon count observed outside of the core to the total photon count, i.e., $L_{X,\mathrm{ex}}=\frac{N_{\mathrm{tot}}-N_{\mathrm{ce}}}{N_{\mathrm{tot}}}L_{X}$. In this work, the core-excised luminosity is used since it has been shown to have lower intrinsic scatter with the cluster mass (Maughan, 2007; Mantz et al., 2018) and is less sensitive to the details of the complicated core physics models used in the simulations. We note that, for the “realistic eROSITA” observations, the core-excised photon count ratios are computed prior to the addition of background noise. This likely makes our mock luminosities more accurate than in the case of real eROSITA observations. However, this choice puts the core-excised luminosities from our two observation series on equal footing, such that the model performance differences between the idealized and realistic cases will be dominated by the quality of the morphological parameters.

Redshift is not explicitly included as a feature to train the regression models, however the redshift is used in scaling the luminosity. Thus, the core-excised luminosity used in this work is always appropriately scaled by the redshift evolution factor, assuming self-similarity, such that we use

[TABLE]

In the $0.5-2.0\,\mathrm{keV}$ energy band, eROSITA anticipates an average photon plus particle background of $2.19\times 10^{-3}\,\mathrm{counts}\,s^{-1}\,\mathrm{arcmin}^{-2}$ (Clerc et al., 2018). Thus, for eROSITA, the background is given by a Poisson distribution with rate $\lambda=0.113/(2\,\mathrm{ks})$ (Merloni et al., 2012). A unique Poisson background is generated for and added to each eROSITA mock observation.

3 Morphological Parameters

In order to encode information about the dynamical state of the cluster into the model, we incorporate various morphological parameters as features, all of which can be directly calculated from the mock X-ray images. In the following, we define each of these parameters. We refer the reader to Lotz et al. (2004), Rasia et al. (2013), and Lovisari et al. (2017) for more in-depth discussion on each of the parameters. Unless otherwise specified, the aperture used to compute the morphological parameters has a radius of $R_{\mathrm{ap}}=R_{\mathrm{500c}}$ and is centered on the cluster X-ray peak; we discuss the implications for this choice at the end of this section.

First, the concentration parameter $c$ quantifies how centrally concentrated the X-ray emission is within the cluster, and has been shown by Santos et al. (2008) to be an indicator for the presence of cooling-core systems at high $z$. Concentration is defined to be the ratio of the flux within two circular apertures: $0.1R_{\mathrm{ap}}$ and $R_{\mathrm{ap}}$ (Lovisari et al., 2017).

The centroid shift parameter $w$ is defined as the variance of the projected separation between the X-ray peak of the image and the emission centroid obtained within 10 circular apertures of increasing radius up to $R_{\mathrm{ap}}$ (Lovisari et al., 2017).

The power ratios, introduced by Buote & Tsai (1995), use the ansatz that the X-ray surface brightness profiles are a good tracer of the cluster’s projected mass distribution. The “power” is encoded in the coefficients of a 2D multiple decomposition of the cluster X-ray image, where higher-order components probe increasingly smaller scales. The $n$th-order power ratio $P_{n0}=P_{n}/P_{0}$ is, in essence, the ratio between the $n$th multipole moments and the [math]th multipole moment. In this work, we consider $P_{10}$, $P_{20}$, $P_{30}$, and $P_{40}$. The latter two probe large- and small-scale substructures present within the cluster, and thus further convey dynamical information.

The second power ratio $P_{20}$ provides a measurement of the cluster ellipticity. Another ellipticity parameter, denoted $e$, is also calculated, defined as the ratio between the semiminor and semimajor axis (Lovisari et al., 2017).

The asymmetry parameter $A$ quantifies the rotational symmetry of the cluster X-ray emission (Lotz et al., 2004). $A$ is calculated by rotating by $180^{\circ}$ and self-subtracting the background-subtracted cluster image from itself, summing the values of the pixels in this image difference and normalizing by the summed pixels in the original image (Abraham et al., 1996).

The smoothness $S$ quantifies the degree of small-scale substructure within the cluster (Lotz et al., 2004). $S$ is calculated by boxcar-smoothing and self-subtracting the background-subtracted cluster image from itself, again summing the values of the pixels in this image difference and normalizing by the summed pixels in the original image (Conselice, 2003).

Lastly, the $M_{20}$ parameter is an analog of concentration (Lotz et al., 2004). The total second-order moment of the light is a distance-to-center-weighted sum of the flux $f_{i}$ within all pixels $i$ in the cluster, $M=\sum_{i}f_{i}[(x_{i}-x_{cc})^{2}+(y_{i}-y_{cc})^{2}]$, where $cc$ denotes the cluster center. Then, $M_{20}$ is computed as the ratio of the partial second moment $M_{p}$, which sums over only the brightest pixels that contain 20% of the cluster light, divided by the total second moment, written as $M_{20}=\log(M_{p}/M)$.

These morphological parameters encode dynamical state information. For example, disturbed clusters tend to be asymmetric (high $A$), clumpy (high $S$), and not concentrated (low $c$). All of the parameters introduced above are calculated for each mock cluster observation, and are used as features in our regression model.

In the subsequent analysis, we utilize two distinct series of morphological parameters, which are computed from our two mock observations series, described above. The “idealized Chandra” series is computed from the background-free Chandra observations, using $R_{\mathrm{ap}}=R_{\mathrm{500c}}$. The “realistic eROSITA” series is computed from the eROSITA observations with added background, also using $R_{\mathrm{ap}}=R_{\mathrm{500c}}$; in this case, the mean background is subtracted prior to computing the parameters. The former series is intended to give an upper bound on the expected improvement in cluster mass estimates made possible by including dynamical state information present in idealized, high spatial resolution (0.5”), and high-photon count observations. The latter series is intended to give a more realistic estimate of the expected improvement that will be possible in upcoming cosmological analyses that will be performed with low-photon count cluster observations. We acknowledge that by using the exact $R_{\mathrm{500c}}$ for our aperture when computing the morphological parameters, we are neglecting additional scatter that will be present due to this effect. Also, in the case of the “idealized Chandra” series, properly including PSF effects would introduce additional smoothing to these observations. Hence, our subsequent results will remain as optimistic estimates. As we find below, the most important morphological parameters are smoothness, asymmetry, and concentration. In Fig. 2, we plot the distributions of these three parameters, comparing the “idealized Chandra” and “realistic eROSITA” series. While we find generally good agreement between the two series, it is clear that the eROSITA cluster observations result in systematically lower concentrations and higher smoothness parameters.

The lack of high-concentration objects in the eROSITA mocks is due to the broader PSF of eROSITA with respect to Chandra. Photons originating from the central regions of the observed systems are redistributed over a wider area, which reduces the concentration with respect to the true value. Since we do not attempt to correct for PSF smearing by applying PSF deconvolution, our procedure for reconstructing $c$ values from eROSITA mocks underestimates the concentration of highly-peaked objects. While also impacted by the broader PSF, the shift to larger $S$ (i.e., less smooth) in the eROSITA mocks is additionally due to both (i) the lower exposure time, which results in a less “filled in” photon distribution due to Poisson noise, and (ii) the presence and subtraction of background, which introduces additional Poisson noise.

Additionally, in Fig. 3, we show several example Chandra cluster images to demonstrate the morphological parameters, in particular the concentration $c$, asymmetry $A$, and smoothness $S$. All clusters shown have roughly the same mass, lying in the range $10^{14.3}\leq M_{\mathrm{500c}}/(h^{-1}\,\mathrm{M_{\odot}})\leq 10^{14.6}$, and are all at the lowest redshift of $z=0.1$. The images are all scaled by the cluster $R_{\mathrm{500c}}$. Clearly, a cluster with a larger concentration has a substantially larger fraction of its flux coming from its core. Furthermore, the asymmetry parameter is successfully able to capture disturbances or substructure in the cluster that result in reduced symmetry. The smoothness parameter is capable of quantifying small-scale structures; note that a cluster with a larger value of $S$ is overall less smooth, and more likely to contain substructures.

4 Analysis Methods

4.1 Data Preprocessing

As stated previously, our sample consists of 2,041 mock cluster observations across six redshifts in the range $0.1\leq z\leq 0.29$. Each observation consists of many features, including core-excised $L_{\mathrm{ex,z}}$ and all of the morphological parameters described in the previous section. The logarithm of the power ratios $P_{i0}$, centroid shift $w$, and luminosity are used due to their large dynamic ranges, whereas the remaining features are not transformed. The regression target for each observation is $\log(M_{\mathrm{500c}})$ of the cluster.

The sample is then split into a training set that comprises 80% of the observations (1,617 clusters) and a test set that comprises the remaining 20% of the observations (425 clusters); this train-test split is a common rule of thumb based on the Pareto principle. The split is performed such that all redshift observations of each unique cluster are assigned either to the training or test set, but not split between the two. For optimization of hyperparameters, $k$-fold cross-validation is employed on the training set, with $k=10$. The folds are generated such that all observations of each unique cluster are confined to only one fold.

Many regression algorithms require the distribution of each observable to be scaled to have roughly zero mean and unit variance. In order to scale in such a way that is robust to outliers, we subtract the median and divide by the $1\sigma$ (16th/84th) percentile range computed over the training set in order to standardize each feature. The medians and $\sigma$ are stored from the training set such that an identical transformation is applied to the test set.

4.2 Regression Methods

The work of Armitage et al. (2019) found that ordinary linear regression (OLR) and ridge regression (RR; Hoerl & Kennard 1970) models were able to produce the least scatter in cluster mass estimates using a variety of X-ray, spectroscopic, and photometric datasets. The authors also tested an ordinary decision tree model (Quinlan, 1986), as well as AdaBoost (Freund & Schapire, 1996) and gradient boosted regression (Friedman, 2001), but did not test the popular random forests (RF) regression algorithm (Breiman, 2001). Motivated by their work, and considering that our feature set of morphological parameters contains different information, we will focus our analysis on various linear regression methods and expand by applying random forest regressors. We train different regression models on our mock catalogs, including a standard mass-luminosity power law ($M-L_{\mathrm{ex,z}}$), an OLR model, several regularized linear regression models (including RR and Lasso regression [LR; Tibshirani 1996]), and RF regression models.

Ordinary linear regression is performed as follows. For $n$ clusters, each of which are described by $p$ features (i.e., observables), one has a data matrix $\mathbf{X}=\{\mathbf{x}_{0}^{T},\mathbf{x}_{1}^{T},\mathbf{x}_{2}^{T},...,\mathbf{x}_{n}^{T}\}$, where each $\mathbf{x}$ is a vector of length $p$. Each vector of observables $\mathbf{x}_{i}$ is associated with a true logarithmic mass $y_{i}$. The mass is predicted as a linear function of the observables, $y_{i}=\mathbf{x}_{i}^{T}\bm{\beta}+\epsilon_{i}$, where $\bm{\beta}$ is the model parameter vector of length $p$ and $\epsilon_{i}$ is the random error in the model for cluster $i$. The best-fit model parameters are chosen by minimizing the cost function, which is selected to be the sum of squared residuals, $E_{\text{OLR}}(\bm{\beta})=\sum_{i=1}^{n}(y_{i}-\mathbf{x}_{i}^{T}\bm{\beta})^{2}$.

In this work, we consider several OLR models. First, we train a simple OLR model with only one feature, the core-excised luminosity $L_{\mathrm{ex,z}}$. This allows us to set a baseline for performance and compare the results of our mock observations in terms of the scatter to an observed mass-luminosity relationship. Then, we train an OLR model on the full feature set, including core-excised luminosity and all morphological parameters.

In an effort to reduce the feature dimensionality and highlight the most important features in the model, one can use regularized linear regression, where an additional term is added to the cost function that introduces a penalty for models with large $\|\bm{\beta}\|$. In ridge regression, the new cost function is of the form $E_{\text{RR}}(\bm{\beta};\alpha)=E_{\text{OLR}}(\bm{\beta})+\frac{\alpha}{2}\|\bm{\beta}\|_{2}^{2}$, where $\|\cdot\|_{2}$ denotes the Euclidean norm. Similarly, in Lasso regression, the cost function is instead $E_{\text{LR}}(\bm{\beta};\alpha)=E_{\text{OLR}}(\bm{\beta})+\alpha\|\bm{\beta}\|_{1}$, where $\|\cdot\|_{1}$ denotes the Manhattan norm.

Regularization acts to reduce the weights of unimportant features in the model, reducing the capability of the model to overfit the training data. Lasso regression is a more strictly regularized model than ridge regression. The hyperparameter $\alpha$ is selected via a grid-search cross-validation (CV) of logarithmically-spaced $\alpha$ values, where the model performance is evaluated via $k$-fold CV for each $\alpha$. In $k$-fold CV, the training set is split into $k$ random subsets (split according to unique cluster ID; see Sec. 4.1). Then, $k-1$ of the subsets are used to train the model, and model predictions are made on the remaining subset. This process is iterated $k$ times such that predictions are made for all clusters in the training set. The model performance is quantified by the mean squared error (MSE) of all of the predictions. The CV process is repeated for all $\alpha$ in the grid, and the $\alpha$ that minimizes the MSE is selected for the training of the final model, which is trained on the full training set and subsequently applied to the test set.

For our last set of models, we use the non-parametric random forest regression model, which is an ensemble technique based on decision trees. RF models reduce the issues of overfitting that are endemic to decision trees by randomly growing an ensemble of trees, each trained on a different subset of the total training data, and taking the average of their predictions. Furthermore, RFs increase the tree diversity relative to a standard decision tree ensemble by splitting each node according to the best feature in a random subset of the features, instead of the full feature set. This increased tree diversity results in a more generalizable model that is less prone to overfitting the training set. RFs have several important hyperparameters: (i) the number of trees in the forest, (ii) the maximum number of features that can be included in one node splitting condition, (iii) the maximum depth allowed for a tree (i.e., number of decisions that must be made to reach an output), (iv) the minimum number of samples in the training set at a particular node that are required in order for the node to split, (v) the minimum number of samples in the training set required to form a leaf, and (vi) whether or not to use bootstrap resampling (i.e., using “bagging” vs. “pasting”). Reducing the “maximum” hyperparameters (i.e., [ii] and [iii]) or increasing the “minimum” hyperparameters (i.e., [iv] and [v]) is an effective way to regularize the model and reduce the tendency for overfitting. The interested reader should refer to Géron (2017) for additional details of various machine learning regressors, including ensemble and tree-based regression.

In this work, we consider several RF models with different sets of hyperparameters and different input feature sets in order to demonstrate the level of sensitivity that RF models have to the hyperparameters and to tune an optimal model for future mass predictions. The first is a RF model with the default hyperparameters from the scikit-learn implementation. The hyperparameters of the second model are optimized using grid search CV over the six-dimensional parameter space of hyperparameters described above. The third model includes a reduced set of features ($L_{\mathrm{ex,z}}$, $S$, $A$, and $c$), but the hyperparameters are also tuned via grid search CV.

After selecting hyperparameters for the various models using CV on the training set, the final models are each trained on the entire training set. The models are then applied to predict the masses of the test set, which we emphasize was never used for either hyperparameter selection or model training, and thus should represent a true example of the generalization capability of the models. The entire preprocess-split-cross-validate-train-test procedure is performed separately for each of the two series of morphological parameters, i.e., those from the “idealized Chandra” and “realistic eROSITA” observations. In the next section, we report the results for these final models as applied to the test sets.

5 Results

For both series of observations, we compute the Pearson correlation coefficient between each observable and the cluster mass, shown in Table 1. Additionally, the best fit linear regression model between $\log(M_{\mathrm{500c}})$ and $\log(L_{\mathrm{ex,z}})$ is used to make mass predictions, and the corresponding mass residuals are then correlated against the observables, also shown in Table 1. The mass residuals $\mathcal{R}$ are defined as

[TABLE]

where we again emphasize that the base-10 logarithm is used throughout.

After luminosity, the observables that correlate or anti-correlate most strongly with mass are smoothness $S$ and asymmetry $A$. The centroid shift $w$ and first power ratio $P_{10}$ correlate most strongly with the mass residuals, although these correlations are still quite weak ($|r|\approx 0.2$). Thus, the naive expectation is that $w$ and $P_{10}$ should be the most important additional features (i.e., after luminosity) in a multivariable model of the mass. However, as we will show below, this ends up not being the case. We note that while the ranking of the morphological parameters in terms of their correlation strengths remains close to the same between the two mock observation series, the strengths are systematically stronger in the “realistic eROSITA” observations. In particular, the high-order power ratios, smoothness, and asymmetry (i.e., the parameters that quantify substructure) correlate much more strongly with mass in the eROSITA observations, which is likely a result of the deviation from a smooth profile driven by Poisson noise in the low-photon count regime. On the other hand, the correlations between the Chandra morphological parameters and the cluster masses are in good qualitative agreement with Lovisari et al. (2017), which, using XMM-Newton cluster observations, found no significant correlation between the total mass and any of $c$, $w$, or the power ratios.

The primary model of interest is our cross-validated random forest regressor, which, as we will show below, performs the best among all of the regression methods tested for both series of morphological parameters. The mass predictions generated by the random forest model for the 426 clusters in the test set are shown in Figure 4, with the two separate panels corresponding to the models trained and tested on the two different series of mock observations. In both cases, it is clear that the model begins to systematically underpredict the masses of the high-mass clusters with $M_{\mathrm{500c}}\gtrsim 10^{14.4}\,h^{-1}\,\mathrm{M_{\odot}}$, which roughly corresponds to the regime where our sample transitions from a flat to falling mass function. In order to employ this method to predict the masses of observed clusters, it is crucial that the training sample consists of a flat mass function that covers the entire range of masses of interest. The performance of machine learning models, such as RFs, will greatly improve as larger training samples that are uniform in the prediction (in this case, the mass) become available, for example from state-of-the-art cosmo-hydrodynamical simulations.

The PDFs of the mass residuals for these cross-validated random forest regression models are shown in Figure 5. Additionally, the $1\sigma$ intrinsic scatter (i.e., half of the 16th–84th percentile range of $\mathcal{R}$) in the test set for each of the trained models and both of the observation series are shown in Table 2. For our test sample, the mass residuals of the standard mass-luminosity relationship have a bias of $\mu=-0.017$ dex and $1\sigma$ scatter of $\delta=0.081$ dex. Interestingly, for both the “idealized Chandra” and “realistic eROSITA” observations, the mass residuals have virtually negligible biases and $1\sigma$ intrinsic scatter of $\delta=0.066$ dex, which amounts to a 20% reduction in scatter relative to the mass-luminosity relationship. Table 2 demonstrates that ordinary linear regression with a combined input feature set that includes the luminosity and all morphological parameters improves only marginally over the single variable mass-luminosity regression. The incorporation of regularization (i.e., the RR and LR models) does not result in an improved model. The lack of improvement in these linear models after the inclusion of morphological parameters illustrates that the relationship between cluster morphology and mass is nonlinear and justifies the use of nonlinear approaches, such as a RF regressor.

In the “idealized Chandra” observations, 140 of the most massive clusters have $R_{\mathrm{500c}}$ that extend beyond the instrument field of view. This results in changes to the morphological parameters calculated for these clusters; for example, the concentrations will systematically increase (although only slightly since the cluster outskirts have the lowest surface brightness) and parameters that quantify substructure ($S$, $P_{30}$, $P_{40}$) may deviate from the correct value if substructures lie outside of the field of view. We verified that the reported scatters are insensitive to the presence or removal of these clusters from the dataset. However, this effect, in addition to the dearth of high-mass clusters in the training sample, is likely responsible for the less accurate mass predictions for high-mass clusters (and lower predicted-to-true slope $a$) when using the Chandra observation series.

The intrinsic scatter for core-excised luminosity-based mass estimates in the observational literature ranges from $\lesssim 15\%$ for the weak lensing mass-$L_{X,\mathrm{ex}}$ relationship (Mantz et al., 2018) to 16–21% for the $Y_{X}$ mass-$L_{X}$ relationship (Maughan, 2007). While the work of Mantz et al. (2018) finds a lower scatter than we do for our $M-L_{\mathrm{ex,z}}$ relation, we note that they employ a mass cutoff of $M\geq 3\times 10^{14}\,M_{\odot}$, which results in a much smaller mass range than covered by our dataset. The scatter in the $M-L_{\mathrm{ex,z}}$ mass residuals for the subset of our clusters with $M_{\mathrm{500c}}\geq 3\times 10^{14}\,M_{\odot}$ is 17%, which is rather close to that of Mantz et al. (2018). The sample used by Maughan (2007) includes clusters down to $8\times 10^{13}\,M_{\odot}$, which is more consistent with our cluster sample. Thus, the scatter in the $M-L_{\mathrm{ex,z}}$ mass residuals of our full sample, $\delta=0.081$ dex (20.5% scatter), is consistent with similar such calculations performed using observations of either the $Y_{X}$ mass-$L_{X,\mathrm{ex}}$ or weak lensing mass-$L_{X,\mathrm{ex}}$ relationships. The current state-of-the-art mass estimation methods require high-resolution, long-exposure cluster observations with good spatial and spectral resolution, with the $Y_{X}$ approach resulting in 5–7% scatter (Kravtsov et al., 2006). The observational conditions necessary for utilizing the $M-Y_{X}$ method will simply not be present for the vast majority of clusters observed in upcoming surveys such as eROSITA. However, we have demonstrated that our method achieves a 20% improvement in cluster mass estimates over $M-L_{\mathrm{ex,z}}$ even in the low-spatial resolution, short-exposure (2 ks) conditions of eROSITA observations ($\delta=0.066$ dex, 16%). Since we also find the same level of improvement for our Chandra observations, which likely places an upper bound on the method performance, this suggests that the mass-encoding dynamical state information, as quantified by our set of morphological parameters, remains present even in the short-exposure eROSITA observations.

One metric available for interpreting the results of the random forest model is the feature importance ranking. For example, in an OLR model, the feature importances are roughly quantified by the magnitudes of the regression coefficients. The standard metric for RF feature importance, and the one that is implemented in scikit-learn, is the mean decrease in impurity, which measures how effective a feature is at reducing the variance of the predictions. Based on this importance measurement, the most important features and their Chandra importances are (after $L_{\mathrm{ex,z}}$, 75%), in decreasing order: smoothness $S$ (10%), asymmetry $A$ (5%), and concentration $c$ (4%). The ranking of these features are the same for both observation series, with importance magnitudes being similar, and all other morphological parameters have negligible importance ($\lesssim 1\%$) in both cases (see Table 3).

Another measurement of feature importance, known as permutation importance, quantifies the drop in the $R^{2}$ score when the values of a feature are permuted over the samples. Thus, the larger the drop in $R^{2}$ when a feature is permuted, the more important the feature. This metric is considered less biased, as it is not sensitive to the dynamic range of the input variables; this detail is likely irrelevant since our features are scaled. When this importance metric is employed, we find the same feature importance ranking as before for the “realistic eROSITA” parameter series. However, for the “idealized Chandra” series, we find that, following luminosity, $c$ and $S$ are nearly tied for being most important, followed by $A$. These three features are not the highest correlators with mass residual, which goes against the naive expectation that the most important features after luminosity should correlate the strongest with mass residual. However, $S$ and $A$ do correlate the most strongly with mass after $L_{\mathrm{ex,z}}$, as seen in Table 1. Thus, there must be some nonlinear relationship between $L_{\mathrm{ex,z}}$, $S$, $A$, and $c$ (slightly supplemented by the combination of all the other morphological parameters) that the RF model identifies in order to make the improved mass estimates.

Some of the morphological parameters employed, particularly $w$ and the power ratios, will be more difficult to measure accurately for eROSITA-observed clusters. Motivated by our finding that $S$, $A$, and $c$ are the most important morphological parameters, we consider an additional RF model that is cross-validated and trained with a reduced feature set that includes only these three parameters and $L_{\mathrm{ex,z}}$. As displayed in Table 2, we find that this reduced model yields a $1\sigma$ scatter of $\delta=0.070$ for the “idealized Chandra” series and $\delta=0.071$ for the “realistic eROSITA” series. While this is still a 12–13% improvement over the mass-luminosity relation, this finding highlights the benefit of including the additional morphological parameters, even considering that each of them has an importance of $\lesssim 1\%$. The combined effect of the additional morphological parameters, including $e$, $w$, $M_{20}$, and the power ratios, is ultimately responsible for roughly a third of the overall improvement offered by our approach.

6 Conclusion

In this work, we have presented a method for estimating cluster masses from mock X-ray observations of galaxy clusters. The mock observations are generated from 2,041 clusters with masses in the range of $10^{13.5}\leq M_{\mathrm{500c}}/(h^{-1}\,\mathrm{M_{\odot}})\leq 10^{14.8}$ and over a redshift range of $0.1\leq z\leq 0.29$ from the Box2 and Box2b Magneticum simulations. The mass predictor is based on a random forest regression model, trained with a feature set that includes only the core-excised luminosity and a set of morphological parameters, all of which can be computed directly from an X-ray image.

We demonstrate that this method can be used to estimate the masses of galaxy clusters with negligible bias and a scatter of $\delta=0.066$ dex (16%). This model achieves a 20% reduction in the scatter relative to a more standard core-excised luminosity power law. Importantly, the same level of improvement is present both when using idealized, high-resolution, long-exposure (1 Ms) Chandra mock observations with no background and when using realistic, low-resolution, short-exposure (2 ks) eROSITA mock observations with added background noise. The majority of this improvement comes from three parameters: smoothness $S$, asymmetry $A$, and surface brightness concentration $c$. A more conservative model, which includes only the luminosity and these three parameters, estimates the cluster masses with a scatter of $\delta=0.070$, demonstrating that a third of the overall improvement comes from the inclusion of the additional morphological parameters ($e$, $w$, $M_{20}$, and the power ratios). However, it is yet to be seen how additional sources of error present in real observations will affect the performance of this model; for example, the scatter in $R_{\mathrm{500c}}$ measurements will propagate to increased scatter in the morphological parameters.

While excising the cluster core reduces the scatter in mass estimates, this improvement comes at the cost of lowering the photon counts used in the analysis. However, even at the eROSITA detection threshold of ${\sim}$30 core-excised photon counts, the statistical uncertainty on $L_{\mathrm{ex,z}}$ will be $\lesssim 20\%$; hence, the mass estimate errors from a mass-luminosity relationship will be dominated by intrinsic scatter even in the low-photon limit. Since the dynamical state of the cluster encodes important information that affects mass errors, the inclusion of morphological parameters, which utilize the full photon distribution, in the mass model enable more accurate predictions with reduced intrinsic scatter relative to $M-L_{\mathrm{ex,z}}$. However, we expect that the statistical uncertainty in the mass estimates will still closely follow the corresponding statistical uncertainty of $L_{\mathrm{ex,z}}$. The relationship between the morphological parameters, the luminosity, and the mass is complicated and nonlinear; we have demonstrated that the nonlinear RF regression method offers a substantial improvement over linear models.

Our model was trained to predict the spherical overdensity masses, $M_{\mathrm{500c}}$, of galaxy clusters identified in Magneticum by the SUBFIND algorithm. As demonstrated in the Knebe et al. (2011) halo finder comparison project, this algorithm is able to estimate the $M_{\mathrm{200c}}$ masses of NFW host haloes to within $\lesssim 3\%$; more broadly, all modern halo finders compared in Knebe et al. (2011) are able to determine host halo $M_{\mathrm{200c}}$ to within $\lesssim 10\%$. Hence, we expect that uncertainty introduced due to the mass estimates of our simulated clusters is sub-dominant, but not insignificant, compared to the intrinsic scatter of the $M-L_{\mathrm{ex,z}}$ relationship. The recent FABLE simulations project (Henden et al., 2018), a set of hydrodynamical simulations with similar sub-grid physics to Magneticum, reports that the $M-L_{X}$ of simulated clusters is in excellent agreement with the observed relation based on X-ray hydrostatic masses. Thus, if there is indeed an X-ray hydrostatic mass bias, this would indicate that simulated clusters may be too gas-rich, resulting in $L_{X}$ values that are high relative to weak lensing-calibrated $M-L_{X}$. In addition, current models of AGN feedback result in simulated cluster cores that do not match the observed cool-core and non-cool-core cluster populations. Because of this, the morphological parameters of observed clusters that depend most sensitively on the cluster core (i.e., X-ray surface brightness concentration) may be biased relative to observations. We expect that these simulation sources of uncertainty will improve as the hydrostatic mass bias quandary approaches a resolution and as more sophisticated AGN feedback models are developed.

Random forest models are notoriously bad at extrapolation. We expect that our model’s systematic underprediction of the masses of clusters in the high-mass tail of the halo mass function will improve when trained on a cluster sample that is uniform across the full mass range of interest. When large-volume, high-resolution hydrodynamical simulations become available for creating such a training sample, we expect that this can be used to train a model that predicts well across the entire mass range. This model, once trained on a sufficiently large simulated sample, could then be applied to a set of Chandra-observed clusters, such as the HIFLUGCS sample (Zhang et al., 2011), and the predictions could be tested against accurate mass estimates, such as those based on $Y_{X}$.

ML-based methods of estimating galaxy cluster masses from X-ray observations, including the method presented here as well as others in the literature (e.g., Ntampaka et al., 2019b), offer a promising step towards extracting the maximum information content present in imminent datasets such as eROSITA. Modern ML methods will enable the completion of an unprecedentedly accurate cosmic census and position the halo mass function to be used to place ever stronger cosmological constraints. Ultimately, the continued progress in cosmological hydrodynamical simulations, both in terms of physical realism and size, are rapidly facilitating the coming of age of these techniques, which will soon be ready for deployment on state-of-the-art cluster observation samples.

The authors would like to thank Stefano Ettori for providing the original morphological parameter estimation code. SBG is supported by the US National Science Foundation Graduate Research Fellowship under Grant No. DGE-1752134. LL acknowledges support from NASA through contracts 80NSSCK0582 and 80NSSC19K0116. KD acknowledges the support of the DFG Cluster of Excellence ‘Origins’.