Decoding the Early Universe: Exploring a Merger Scenario for the High-Redshift Cluster JKCS041 using Numerical Models

The cluster-scale dark matter haloes in the simulation are described by the super NFW (sNFW) profile (Lilley et al., 2018a, b), a model very similar to NFW (Navarro et al., 1996, 1997) but with finite mass ensuring a smooth cutoff at larger radii. The equations for the sNFW profile given below are taken from Lilley et al. (2018a). NFW haloes are usually parameterized by a mass $M_{200}$ or radius $r_{200}$ and mass concentration parameter $c$. $M_{200}$ is the mass contained inside radius $r_{200}$, at which the mean enclosed density is 200$\rho_{\rm crit}$. The critical density, $\rho_{\rm crit}=3H(z)^{2}/8\pi G$, where $H(z)$ is the Hubble parameter at the redshift $z$. The NFW scale radius is defined as $r_{s}=r_{200}/c$. The sNFW density profile is given by:

where $a$ is the sNFW scale radius which is related to the half-mass radius $R_{e}$ by $a=R_{e}/5.478$. The $M_{s}$ parameter is determined by $a$ and an sNFW concentration parameter $c_{sNFW}$. Fitting sNFW haloes with NFW density profiles, Lilley et al. (2018a) found that the concentrations $c_{sNFW}$ and $c$ are related by:

The scale radius of an NFW profile $r_{s}$ is related to the sNFW $a$ parameter by Lilley et al. (2018a):

Essentially, the parameters we study are NFW parameters $M_{200}$ and $c$ that map to the sNFW profile (equation 2) parameters $M_{s}$ and $a$ as described in Lilley et al. (2018a).

We set up the initial clusters with gas density following the modified $\beta$-profile as described in (Vikhlinin et al., 2006), along with requiring that the gas and dark matter are in hydrostatic equilibrium. This assumption is valid only for relaxed clusters and will be discussed in Section 6. This is achieved by giving the dark matter particles velocities after the mass profiles are set, as per Kazantzidis et al. (2004) in which the energy distribution is calculated using the Eddington formula (Eddington, 1916).

Most clusters in the local Universe are characterized by a peak in the X-ray surface brightness and a lower temperature in the centre denoting a ‘cool core’. The distribution of the gas in cool core clusters is described by a modified-$\beta$-profile (Vikhlinin et al., 2006),

where $n_{e}$ and $n_{p}$ are the electron and proton number densities. $\alpha$ is the inner density slope which determines the strength of the cool core and $r_{c}$ is the core radius. Cool cores are well-represented by parameters similar to $r_{c}=0.05a$, $r_{s}=0.6a$ and $\alpha=2$ (Chadayammuri et al., 2022). We use $\beta=2/3$ which was found to fit X-ray surface brightness of clusters in Jones & Forman (1984) and Vikhlinin et al. (2006). The parameter $\epsilon$ describes a steepening of slope near the radius $r_{s}$ and the parameter $\gamma$ controls the width of the transition. The default values that will be adopted for the remaining plasma cloud parameters are $\gamma=3$ and $\epsilon=2$, consistent with studies of multiple clusters in Vikhlinin et al. (2006).

The initial conditions for the merger simulations consist of specification of sNFW parameters (using NFW parameters, $M_{200}$ and $c$ through equations (4) and (3) and gas parameters for the two haloes, initial distance between the haloes, impact parameter and the initial halo velocities (or relative velocity of the haloes).

Fig. 1 shows the geometry of the merging process, and location of the observer. The simulation box is a cube of 14 Mpc side length where the objects are placed 3 Mpc apart initially. The individual clusters are then given an initial velocity towards each other along the X-axis and hence the merger happens primarily along the X-axis. The impact parameter specifies the displacement in the Y direction of the second cluster from the origin at the start of the simulation. Unless otherwise specified, we use the Z-axis projection to study the evolution of the merger with time (equivalent to the observer being on the Z-axis). We run the simulation for 2.5 Gyr, which allows the system to evolve for sufficient time after first pericentre passage in all of our runs. In the highest mass configuration we ran, this time was enough for a second pericentre passage. Considering the age of the Universe at $z=1.803$ ($\approx$3.6 Gyr), a longer time evolution is not favored.

For ease of analysis, for the merger scenario we analyze each of the configurations with data outputs between 1.2 Gyr and 2.5 Gyr at a frequency of every 0.1 Gyr. This time window is chosen as it takes approximately 1.2 Gyr from the start of simulation for the clusters to get sufficiently close so that the behavior starts to deviate from that of isolated clusters. This time range also gives us enough pre-pericentre passage data.

For the main set of merger simulations that we run, and considering the masses we use, pericentre passage occurs between 1.5 Gyr and 1.8 Gyr after the start of the simulation. Therefore, we choose the initial redshift as $z=3.4$ which corresponds approximately to $1.8$ Gyr before the observed system redshift $z=1.8$.

We use the output of the GAMER-2 (Schive et al., 2018; Zhang et al., 2018) simulations in combination with the yt code (Turk et al., 2011; Smith et al., 2022) to obtain synthetic SZ, X-ray and gravitational lensing maps.

First, we explore the evolution of the SZ peak (as defined in Section 2) during a typical merger configuration. Fig. 2 shows the variation of the SZ peak with time for a typical configuration (run $v$ in Table 1). When we refer to SZ peak values, we implicitly mean $|\Delta T|$ peak values, which are in practice directly obtained from SZ observations. At the start of the simulations, the SZ peak remains almost constant until the clusters are sufficiently close together, when it increases to reach a maximum around the pericentre and then decreases after that. After pericentre passage the SZ peak attains a lower value than at the start because of the lower gravitational binding of the gas making the distribution more diffuse compared to the distribution before pericentre passage. Therefore, a merger can produce lower SZ peak than the SZ peak of the cluster components in isolation. These trends for mergers are consistent with Wik et al. (2008) and Krause et al. (2012).

In order to narrow down parameter space of mergers that are consistent with the observables, in this section we start by considering the impact of total mass (first, fixing mass ratio $q=1$) on SZ maps. We then constrain the concentration $c$, mass ratio $q$, Impact parameter $I$ and gas fraction ${f}_{\rm gas}$. This process allows us to examine the sensitivity of the SZ signal to these parameters and rule out regions of parameter space.

Note that for $q=1$, in order to have a single SZ component observed, the clusters would need to be seen in projection directly along the merger axis. However, in that case, no offset would be observed between the SZ and X-ray peaks. Later we consider cases where $q\neq 1$, and where a single SZ component can be observed if the lower-mass component would go undetected. Although the $q=1$ case is ruled out by observations, it is useful to first get a general sense of the impact of total mass $M$ on the evolution of the SZ signal, where a single SZ component would be observed. For each of the runs in Var_M and Var_c in Table 1 ($q=1$ mergers), the simulation box is projected along the merger axis to give a single SZ peak. For the snapshot with an SZ peak closest to the observations, Fig. 3 shows the peak SZ temperature decrement as a function of $M$. The navy blue line and blue shaded region show the same for the observed system and error range, respectively. As expected, increasing the total mass of the merger system, increases the temperature decrement when the other parameters are fixed. In this case where $q=1$, the total mass of $2\times 10^{14}M_{\odot}$ produces a peak SZ temperature decrement closest to the observations, whereas higher mass systems produce too large an SZ decrement. This provides us with a starting point for the simulations and depending on the other parameters that we will study in the subsequent sections, the second lowest mass ($2.4\times 10^{14}M_{\odot}$) may become consistent with the observations. Additionally our preliminary runs confirm that masses above $3\times 10^{14}M_{\odot}$ for $q\neq 1$ (with masses different enough that we can have a single prominent SZ peak independent of projection axis) cannot produce a peak SZ temperature decrement as low as the observed value.

In Fig. 4 we vary the initial cluster concentrations in the range $c=3-5$ and show that different concentrations can result in a temperature decrement consistent with the observations. We identify the earliest post-merger simulation snapshot with a $\Delta T$ peak consistent with the observations (solid circles). To assess the impact of snapshot sampling on this measurement, we also identify the snapshots immediately before (green downward arrow) and after (red upward arrow). Concentration is degenerate with SZ temperature decrement and time from pericentre passage in combination. Therefore, we fix $c=5$ in our subsequent analysis, based approximately on the $M-c$ relation (Child et al., 2018).

Mass ratios, $q$, close to 1:1 are likely to produce the largest offset between SZ and X-ray centres (Zhang et al., 2014, 2015), in comparison with the smaller mass ratios considered here. In Fig. 5, we show the dependence of the SZ temperature decrement and of the offset on mass ratio for a total mass of $2\times 10^{14}M_{\odot}$ (Var_q in Table 1). As in Fig. 4, the coloured circles denote the $\Delta T$ and offset measured on the earliest snapshot that is consistent with the observed range, with the downward arrow and upward arrow indicating the previous and subsequent snapshots. For the mass ratios (denoted by different colours) considered, configurations consistent with the observed temperature decrement can be identified, whereas the SZ - X-ray offset is best matched by a mass ratio $q=2:3$ and hence we use it in the subsequent studies. Mass ratios $q\sim 1$ are inconsistent with the observations; in that case the requirement of a single prominent SZ peak necessitates a projection along the merger axis which would result in zero offset between the X-ray and SZ peaks in the plane of the sky.

Fig. 6, shows the peak SZ temperature decrement vs observed SZ-X-ray offset as a function of impact parameter. We can identify snapshots consistent with the observed SZ temperature decrement using initial impact parameters in the range $0$ to $500$ kpc studied. Although satisfying the observed offset suggests a smaller impact parameter, we note that snapshots with larger impact parameters can yield offsets roughly consistent with the data.

The SZ and X-ray signal depend on the initial ${f}_{\rm gas}$ of the clusters. Fig.  7, shows the SZ-X-ray offset vs peak SZ temperature decrement ( Andreon et al. (2023), Andreon et al. (2009); Newman et al. (2014); Andreon et al. (2014)) for different $f_{\rm gas}$. As expected, increasing ${f}_{\rm gas}$ of the clusters results in increased SZ signal due to increased inverse Compton scattering of the CMB photons. Fig. 7 also strongly suggests that configurations with lower initial ${f}_{\rm gas}$ i.e. $0.05$ - $0.1$ are more consistent with the observations.

The majority of studies measuring the gas mass fractions of galaxy clusters have been at low redshift although recent studies using cosmological simulations find lower ${f_{\rm gas}}$ at higher redshifts, e.g., Altamura et al. (2023). Note that some of these studies use different definitions of ${f}_{\rm gas}$, calculated at $r_{500}$ or $r_{2500}$, while we performed the calculation at $r_{200}$. We have measured the gas mass fractions of the most massive clusters extracted at $z\approx 2$ from the Illustris TNG simulations (Weinberger et al., 2017; Pillepich et al., 2018; Springel et al., 2018), finding that our conclusion of a lower ${f_{\rm gas}}$ for JKCS041 (in comparison with clusters at $z=0$) is consistent with that sample, albeit that the number of high-redshift massive clusters is limited (Shavelle et al., 2023).

Fig. 8 shows the dependence of (X-ray based) gas mass in different apertures on mass ratio (left panel), impact parameter (central panel) and ${f_{\rm gas}}$ (right panel), as well as the observed values (black line with shading). The gas mass is measured from the first snapshot after core passage to match the observed peak SZ temperature decrement, or the snapshot closest to the observed value (solid blue circles in Figs 5-7). At some point after first core passage, all of the impact parameters considered yield gas masses consistent with the observations. Mergers of very comparable mass clusters and gas-rich clusters ($f_{\rm gas}=0.17$) produce gas masses as a function of aperture size that are inconsistent with the observations.

To summarise our results so far: the amplitude of the SZ temperature decrement is sensitive to the total mass, and to ${f}_{\rm gas}$. The total mass is found to be about $2\times 10^{14}M_{\odot}$, given the remaining parameters. The integrated gas mass inside apertures depends on ${f}_{\rm gas}$, and values below, or of the order of 0.1, are consistent with the observations. The SZ-X-ray offset is sensitive to the mass ratio and to (a lesser extent) the impact parameter. Given also the non-detection of a strong second SZ peak (which rules out roughly equal mass mergers, $q\sim 1$), the mass ratio is found to be about 2:3, given the remaining parameters. None of the observables considered here have a great sensitivity to the initial cluster concentrations.

The 220 kpc offset between the SZ and X-ray peaks is measured in the plane of sky, and the velocity gradient in Prichard et al. (2017) and in Andreon et al. (2009) indicates that the offset also has a component along the line of sight. Thus, 220 kpc is a lower limit on the 3D separation of the SZ and X-ray peaks. In the results above, we focused on configurations where the mergers occur in the plane of the sky, and hence the projected and 3D SZ-X-ray offset are equivalent.

From our simulations and subsequent analysis, the initial configuration with a total mass of $M=2\times 10^{14}\ M\odot$, mass concentrations of each cluster $c=5$, mass ratio of $q=2:3$, ${f}_{\rm gas}=0.05$ (up to 0.10) and impact parameter of $I=0$ (head on collision) satisfies all of the observational constraints. The evolution of the SZ map of this particular configuration (run $xviii$ in Table 1) with time during the merger is shown in Fig. 9. The snapshot at +0.3 Gyr after first core passage is our best fit result to the observations, and is shown in Fig. 10. The compact appearance of the SZ emission associated with the main peak in Andreon et al. (2023) is reproduced by our simulations.

It is worth noting that this configuration is not exclusive, and that similar configurations with small variations in parameters around this configuration are also consistent with the observations. Generically though, we emphasize that a merger scenario is required to explain all of the observations of JKCS041. Merger scenarios with different initial relative velocities apart from 1100 km/s were also explored (relative velocities of 750 km/s and 1600 km/s), and configurations that match the observations can be identified. All of these velocities were found to be consistent with the observations, with the other parameters fixed as in the best fit configuration. One of the notable differences is that higher initial relative velocities can produce larger SZ-X-ray offsets (1600 km/s produces the largest offset among the velocities we tried). These larger offsets can be reoriented to produce 2D offsets that are consistent with the observations. When exploring other merger configurations (departing from the best fit configuration for the 1100 km/s scenario), in mergers with lower initial velocities (such as 750 km/s), often the maximum 3D SZ-X-ray offset is less than 220 kpc and hence would not be consistent with the observations. For general viewing $\theta$ and $\phi$ angles as indicated on Fig. 1 (not aligned with the Z-axis where $\theta=0$ and $\phi=0$), some merger scenarios can easily be excluded, as with the Z-axis viewing angle. For example, projected SZ-X-ray offset may be too small, or the peak SZ temperature decrement may be too large. However, in general, merger scenarios consistent with the observational constraints can be identified for multiple viewing angles.

The method we adopt is as follows: for the best fit merger simulation configuration, we obtain the projected density ($\Sigma$) map on a grid, and corresponding lensing convergence ($\kappa$) map, where $\kappa=\Sigma/\Sigma_{\rm crit}$ and $\Sigma_{\rm crit}=c^{2}D_{s}/4\pi GD_{ds}D_{d}$. The angular diameter distances to the source and lens and the lens-source distance are $D_{s}$, $D_{d}$ and $D_{ds}$ respectively. The redshifts of the lens and sources are taken to be $z=1.8$ and $z_{s}=2.1$ respectively, determining the angular diameter distances in the adopted cosmology. Note that we take $z_{s}=2.1$ for simplicity, consistent with the peak in the redshift distribution in Kim et al. (2024). We discuss the impact of this assumption in Section 5.2. A synthetic lensing shear ($\gamma$) map is obtained from the $\kappa$ map using Fourier transform techniques, also yielding a reduced shear map ($g=\gamma/(1-\kappa)$). On the sky, a synthetic distant unlensed galaxy is projected at a location on the $g$ map, sampling a particular value. As outlined below in Eq.6, $g$ determines the distorted lensed images of distant galaxy sources.

In the weak lensing regime, unlensed ($\epsilon^{s}$) and lensed galaxies ($\epsilon$) are described by ellipses, conveniently represented by complex numbers; the modulus of the complex ellipticity is related to the galaxy axis ratio $r=b/a$ via $(1-r)/(1+r)$ and the phase of the complex ellipticity is twice the galaxy position angle. Similarly, $g$ at any location is represented by a complex number, with a modulus giving the strength, and a phase which is twice the position angle. For each galaxy the lensed and unlensed ellipticities are related via:

where * denotes complex conjugation. Equations relating the intrinsic (unlensed) and lensed shapes of galaxies can be found in Schneider et al. (2000), for example.

In order to generate a synthetic catalogue of lensed galaxies, we adopt a galaxy number density of $\approx$ 87 arcmin^-2, randomly distributed and randomly oriented in a square field of view of side 1.5 Mpc for consistency with the data used for the weak lensing analysis of Kim et al. (2024). The sources are taken to be at a single redshift of $z_{s}=2.1$ (the effective redshift of the Kim et al. (2024) study); the consequences of this assumption are discussed in Section 5.2. The dispersion of the galaxy ellipticities, quantifying the departure of galaxy shapes from circular, is taken to be $\sigma_{\epsilon^{s}}=0.25$, also in keeping with Kim et al. (2024) and previous studies, e.g., (Schneider et al., 2000). For all the galaxies in the field, synthetic lensed ellipticities were determined using Equation 6. Although synthetic catalogues can also be determined using different sets of random galaxy locations on the sky, here we focus on the impact of galaxy shape noise.

After creating synthetic lensed galaxy catalogues, we obtained weak lensing mass estimates using the method from Schneider et al. (2000) and King & Schneider (2001), fitting a single NFW component to the lensed galaxies. In keeping with Kim et al. (2024) we use the relation between mass and concentration determined from cosmological simulations, see Diemer & Joyce (2019), during the fitting. Note that our comparison is with their main mass component.

Repeating the synthetic catalogue generation and mass recovery 300 times is sufficient to determine mass estimate uncertainties. Fig. 11 shows the distribution of weak lensing $M_{200}$ fits for synthetic data generated using the best fit configuration. The recovered $M_{200}$ distribution has a mean and standard deviation of $2.3(\pm\,0.9)\times 10^{14}M_{\odot}$, consistent with the observational lensing mass estimate $4.7(\pm 1.5)\times 10^{14}M_{\odot}$ (Kim et al, 2023). In our synthetic lensing catalogue generation and analysis we did not account for random distant galaxy positions, or for a source redshift distribution, which would further broaden the distribution of possible $M_{200}$ values. Further, contrary to the case of low redshift clusters, the lensing efficiency for the sources will vary substantially making the distribution of redshifts important.