The Dynamical State of the Frontier Fields Galaxy Cluster Abell 370

Abell 370 (A370) is a well-studied galaxy cluster at a redshift of $z=0.375$ (e.g., Tyson et al., 1990; Kochanek et al., 1990; Miralda-Escude, 1991; Umetsu et al., 1999; Medezinski et al., 2010; Schmidt et al., 2014; Treu et al., 2015; Lagattuta et al., 2017; Lotz et al., 2017; Lagattuta et al., 2019, and references therein). It was one of the first galaxy clusters in which gravitational lensing was observed (Soucail, 1987). The first spectroscopically confirmed giant arc was also discovered in this cluster (Soucail et al., 1988). A370 was found to be one of the most massive clusters based on weak gravitational lensing, with a total virial mass of $M_{\mathrm{vir}}\approx 3.3\times 10^{15}M_{\odot}$ (Umetsu et al., 2011a, b). Because of its large projected mass and high lensing magnification capability of background galaxies, A370 has been known as a superlens characterized by large Einstein radii (e.g., $\theta_{\mathrm{Ein}}>30\arcsec$ for $z_{s}=2$; Broadhurst et al., 2008), and selected as one of the six Hubble Frontier Fields¹¹1https://frontierfields.org clusters (Lotz et al., 2017). More recently, A370 has been targeted by the BUFFALO survey (Steinhardt et al., 2020) with Hubble Space Telescope (HST), which will expand the existing area coverage of the Hubble Frontier Fields in optical and near-infrared pass bands.

The mass surface density, the galaxy number density, and X-ray surface brightness show a symmetric bimodal distribution along the north-south direction in the core of A370, suggesting that it is a massive major merger with a mass ratio close to unity (Richard et al., 2010; Lagattuta et al., 2017; Strait et al., 2018; Diego et al., 2018). Two brightest cluster galaxies (BCGs) were found in the core of A370. The BCG located in the south of the core of the cluster is very close to the mass peak of the southern cluster component (within a few kpc), while the northern BCG shows a significant offset from the northern cluster mass peak ($\sim 50$ kpc; e.g., Richard et al. 2010; Strait et al. 2018). Botteon et al. (2018) used X-ray observations of A370 to search for surface brightness discontinuities to identify shocks and contact discontinuities. They found evidence for surface brightness edges on the western and the eastern side of the cluster. Although the origin of the surface brightness edges is not clear, the edges might be associated with shocks induced by a merger.

Galaxy cluster scaling relations applied to A370 also point to its dynamical activity. Based on the total mass $M_{500}$ (all relevant symbols are defined in detail at the end of this section) estimated from weak lensing (Umetsu et al., 2011a), the X-ray luminosity in the 0.5 – 2 keV band is a factor of $\sim 2$–$3$ times smaller than that inferred from the $L_{\mathrm{X}}$–$M_{500}$ scaling relation of Pratt et al. (2009). Moreover, $M_{500}$ estimated by Planck Sunyaev–Zel’dovich (SZ) effect observations is a factor of $\sim 2.5$ smaller than the weak-lensing mass (Umetsu et al., 2011a; Coe et al., 2019). All of these pieces of evidence, namely the discrepancies in mass scaling relations, positional offsets between the galaxy and mass surface densities, and multiple X-ray brightness peaks, imply dynamical activity in the cluster.

Motivated by these observational results, here we carry out the first attempt to model A370 as a binary major merger by performing dedicated numerical simulations on CPU clusters using our three-dimensional (3D) $N$-body/hydrodynamical code based on flash (Molnar et al., 2012, 2013a, 2013b; Molnar & Broadhurst, 2015, 2017, 2018) developed at the University of Chicago (Fryxell et al., 2000).

The structure of this paper is as follows. In Section 2, we summarize the results from multi-wavelength observations of A370 including our results based on re-analyzing existing Chandra X-ray observations. We describe our simulation setup to obtain a model of A370 as a binary merger in Section 3. Section 4 presents our results from hydrodynamical modeling of A370. Section 5 contains our summary.

We adopt a spatially flat, cosmological-constant and dark-matter dominated ($\Lambda$CDM) cosmology with $h=0.7$, $\Omega_{\mathrm{m}}=0.3$, and $\Omega_{\Lambda}=0.7$. We use the standard notation $M_{\Delta}$ to denote the mass enclosed within a sphere of radius $R_{\Delta}$, within which the mean overdensity equals $\Delta$ times the cosmic critical density $\rho_{\mathrm{c}}(z)$ at the cluster redshift $z$. We adopt the virial overdensity $\Delta\simeq 128$ (i.e., $M_{\mathrm{vir}}\simeq M_{128}$) based on the spherical collapse model (see Appendix A of Kitayama & Suto 1996). The quoted errors represent the 68$\%$confidence level, unless stated otherwise.

We show in Figures 3 and 4 the X-ray surface brightness image of A370 derived from Chandra data and the images of mock Chandra observations based on our flash simulations. Here we applied the same smoothing (5$\arcsec$) on the observed and simulated images for comparison. The images are 1.14 Mpc $\times$ 1.56 Mpc. We used the same color code for all panels and excluded the area containing the foreground elliptical galaxy in the north for clarity. For each mock observation, the epoch and the viewing angle were chosen such that the relative LOS velocity matches $V_{\mathrm{rad}}\sim 750$ km s^-1 and the locations of the two peaks of the mass surface density are aligned with those observed.²²2This can be achieved with a proper choice of the polar angle, $\theta$, if the 3D separation is $\geq 206$ kpc at a given epoch.

In general, a viewing angle that produces a projected mass peak separation of $D_{\mathrm{proj}}\sim 206$ kpc and a relative LOS velocity of $V_{\mathrm{rad}}\sim 750$ km s^-1 may be chosen at a few epochs before and after the core passages, as long as the 3D velocity due to the accelerating or decelerating infalling cluster is greater than $\sim 750$ km s^-1. We found that projections with a LOS relative velocity of $V_{\mathrm{rad}}\sim 750$ km s^-1 before/after the first core passage and before the second core passage appear to result in a single bright peak in the X-ray surface brightness, independently of the infall velocity and impact parameter of the collision and the viewing angle. This is illustrated in Figure 4 (see the first panels in the first and second rows). The relative velocities in the simulations at epochs after the third core passage are usually too low to produce the required relative radial velocity. Therefore, we shall focus on simulations at epochs between the second and third core passages.

The first panel of Figure 3 shows the smoothed Chandra X-ray image based on our analysis (see Section 2.1). Images shown in the other panels are based on our best simulation run (RP300V3.5) at the best epoch, 0.57 Gyr (6.8 Gyr) after the second (first) core passage with a polar angle of $\theta=72.4^{\circ}$ and three different roll angles, namely $\varphi=0^{\circ},-120^{\circ}$, and $-158^{\circ}$. The X-ray morphology with a roll angle of $\varphi=-158^{\circ}$ provides the best match with the observations. Contrary to the observations, views with roll angles of $\varphi=0^{\circ}$ and $-120^{\circ}$ exhibit two X-ray peaks near their associated mass peaks, but no X-ray peak around half way between them.

The best model is provided by our RP100V3.5 run at an epoch of 0.57 Gyr after the second core passage and just before the third core passage, with a polar angle of $\theta=72.4^{\circ}$ and a roll angle of $\varphi=-158^{\circ}$ (see the third panel in Figure 3). The main X-ray peak is located between the two mass centers, about half way with a small offset to the west. It shows an extension toward the southern mass peak and the secondary peak is close to the northern mass center, with a small offset toward south-west. With this viewing angle, the two cluster components of RP100V3.5 has a relative radial velocity of $V_{\mathrm{rad}}=960$ km s^-1 and an integrated SZ amplitude of $Y_{2500}=1.2\times 10^{-10}$ ster.

We show in the right panel of Figure 2 (dashed histogram) the distribution of relative radial velocities constructed from our best model. The corresponding Gaussian fit to our simulated histogram is shown with the dashed curve. We obtain a velocity dispersion of $\sigma_{\mathrm{sim}}\approx 1654$ km s^-1 based on our best model, which agrees well within $10\%$ with the observed value ($\sigma_{\mathrm{obs}}\approx 1795$ km s^-1; see Section 2.3).

We thus conclude that, our simulation (RP100V3.5) with total virial masses of $M_{1}=1.7\times 10^{15}\,M_{\odot}$ and $M_{2}=1.6\times 10^{15}\,M_{\odot}$, an initial impact parameter of $P\sim 100$ kpc, a 3D infall velocity of $V_{\mathrm{in}}\sim 3500$ km s^-1, and concentration parameters of $c_{1}=6$ and $c_{2}=6$ is the best dynamical model for A370. Our best model is in agreement with an earlier interpretation of multiwavelength observations of A370, that is, it is a massive post major merger with $M_{1}/M_{2}\sim 1$. Moreover, the distribution of radial velocities derived from our best model is very similar to that observed in A370 (compare the dashed and solid histograms in the left panel of Figure 2). This agreement provides an independent confirmation of the large mass derived from weak lensing by Umetsu et al. (2011a, b).

We note that the shape of the outer X-ray brightness distribution from our best model is rounder than the observed one, which is a consequence of a shallower fall off of the gas density distribution than that of A370. Here we did not aim to improve the fit of the outer regions, because it would require an even more extensive search in an extended parameter space including the outer density slope of the initial clusters, which would increase the computer time significantly. However, on the basis of our extensive simulations, we do not expect a significant difference in the cluster core regions from changing the slope in the outer regions with much lower gas densities.

We have also performed a few simulations using different values of $P$, $V_{\mathrm{in}}$, and $c_{\mathrm{vir}}$ to derive crude constraints on our best initial parameters (see Table 1 for the list of parameters of simulations discussed in our paper). Specifically, we ran simulations with the infall velocity $V_{\mathrm{in}}$ in the range of $2000$ to $4000$ km s^-1, the impact parameter $P$ from 100 to 400 kpc, and the concentration parameter $c_{\mathrm{vir}}$ from 5 to 8 to cover a representative range of parameter space. We did not perform simulations with zero impact parameter because that would result in an X-ray morphology with axial symmetry contrary to the observations. We chose $V_{\mathrm{in}}=2000$ km s^-1 for the lower limit of the infall velocity because mergers with smaller infall velocities result in an integrated SZ effect ($Y_{2500}$) that is too large to match the Bolocam measurements. An infall velocity of $V_{\mathrm{in}}=4000$ km s^-1 already results in an unacceptably long time interval between the first and second core passages (longer than the age of the universe) for highly massive clusters. Since such massive clusters have lower concentration parameters on average, we chose the lower limit of 5 and ran simulations within the range $c_{\mathrm{vir}}\in[5,8]$. A full parameter search is not feasible with CPU clusters due to the large demand on computer resources.

In Figure 4, we show mock Chandra images based on our simulations between the second and third core passages that do not provide the best match with the observations. Most of the initial parameters were the same as in our best model (RP100V3.5), but we changed one or two parameters for each run to show how the changes impact the X-ray morphology (for the initial parameters, see Table 1). The first panel in the first row shows our simulation with an impact parameter changed to $P=200$ kpc (RP200V3.5) at an epoch 1.1 Gyr (8.7 Gyr) after the second (first) core passage. The X-ray morphology shows only one bright peak. Thus, it does not match the observations although it has an integrated SZ amplitude of $Y_{2500}=1.1\times 10^{-10}$ ster, which agrees with the Bolocam constraint.

The second panel displays our simulation with the infall velocity changed to $V_{\mathrm{in}}=3000$ km s^-1 (RP100V3.0) at an epoch 0.68 Gyr (3.4 Gyr) after the second (first) core passage. The X-ray surface brightness has two peaks with similar amplitudes (but no dominant peak) and no extension to the southern mass center. Moreover, it also has a somewhat large integrated SZ amplitude, $Y_{2500}=1.4\times 10^{-10}$ ster.

In the third and fourth panels, we show X-ray images based on our simulation with an infall velocity of $V_{\mathrm{in}}=3700$ km s^-1 (RP100V3.7) at two different epochs, 0.32 Gyr (11 Gyr) and 1.4 Gyr (12 Gyr) after the second (first) core passage. Again, the X-ray morphologies do not match with the observations: the image in the third panel shows a single peak, while the fourth panel shows two peaks with similar amplitudes and no extension toward the southern mass peak. The integrated SZ amplitudes are $Y_{2500}=2.6\times 10^{-10}$ ster and $9.4\times 10^{-11}$ ster at these epochs, and thus the first epoch can be excluded based on its large integrated SZ amplitude as well.

The first and second panels in the second row of Figure 4 display our simulations with concentration parameters $c_{1}=c_{2}=6$ (RP100V3.5b) and $c_{1}=c_{2}=8$ (RP100V3.5c), while the other initial parameters are the same as those of our best simulation. These images are shown at epochs of 0.16 Gyr (7.8 Gyr) and 0.92 Gyr (6.4 Gyr) after the second (first) core passage. The X-ray morphologies of these images do not provide a good match with the observations, even though they satisfy all other requirements including their integrated SZ amplitude, $Y_{2500}=1.0$ and $1.1\times 10^{-10}$ ster, which agree with the Bolocam constraint.

The third panel in the second row shows our simulations with an impact parameter of $P=300$ kpc and an infall velocity of $V_{\mathrm{in}}=3000$ km s^-1 (RP300V3.0), at an epoch 0.20 Gyr (4.3 Gyr) after the second (first) core passage. The X-ray morphology of this run shows two X-ray peaks, whereas the morphology does not provide a good match with the observations. The fourth panel displays our simulations with an impact parameter of $P=400$ kpc, an infall velocity of $V_{\mathrm{in}}=3000$ km s^-1, and concentration parameters of $c_{1}=6$ and $c_{2}=7$ (RP400V3.0) at an epoch 0.32 Gyr (4.2 Gyr) after the second (first) core passage. Again, we find two X-ray peaks, but the morphology does not match the observations. These two runs can be excluded based on their large integrated SZ amplitudes, $Y_{2500}=2.1$ and $2.4\times 10^{-10}$ ster, which are more than twice as large as the Bolocam constraint, $Y_{2500}=0.9\times 10^{-10}$ ster.

On the basis of our simulations with fixed initial masses ($M_{1}=1.7\times 10^{15}M_{\odot}$ and $M_{2}=1.6\times 10^{15}M_{\odot}$) and different impact parameters, infall velocities, and concentration parameters, we can provide a crude estimate on the uncertainties of our best model parameters as $P=100_{-100}^{+100}\,$kpc, $V_{\mathrm{in}}=3500_{-500}^{+200}$ km s^-1, and $c_{1}=6\pm 1$ and $c_{2}=6\pm 1$. However, we should consider these initial conditions only as a guide to set up the conditions for the second core passage, since our simulations have high fidelity only from around the second core passage until the best epoch. That is why we displayed the time elapsed from the first core passage for our simulations only in parentheses. The reason for this is that, as usually done in controlled binary merger simulations, the expansion of the universe and mass accretion from the surrounding large-scale structure are ignored (e.g., Ricker & Sarazin 2001; Ritchie & Thomas 2002; Poole et al. 2007; McCarthy et al. 2007; ZuHone et al. 2011; for a discussion, see Molnar 2016). This is a good approximation for a few Gyrs, but not for a significant fraction of the age of the universe, such as 5–10 Gyrs, which is the time scale of our simulations from the first to the second core passage, due to the large masses and infall velocities we had to assume. On a time scale of a significant fraction of the age of the universe, the expansion of the background and the mass accretion of clusters cannot be ignored. A study of the dynamics of A370 through its evolution would need a cosmological $N$-body/hydrodynamical simulation, which would incorporate these effects naturally. This could be done, for example, using constrained realizations of Gaussian random fields (Hoffman & Ribak, 1991), or one of its variants (e.g., Roth et al. 2016; Rey & Pontzen 2018; Sawala et al. 2020). However, this is out of the scope of our paper.

In general, during each core passage, outgoing shockwaves are generated in the central region of merging galaxy clusters. We show in Figure 5 the outgoing shocks generated after the second core passage based on our best simulation (RP100V3.5) at the best epoch. We identified the shocks using the magnitude of the pressure gradient. The locations of large gradients in the pressure delineate the shock surfaces. In these images, the red curves show the location of the shocks. The red crosses mark the dark-matter centers and the horizontal red bars indicate the physical scale of 1 Mpc.

In the left panel of Figure 5, we show the magnitude of the pressure gradient in a two-dimensional (2D) cut through the centers of the two components, before the rotation of the system out of the plane of the sky. This image shows the physical (not projected) locations of the shocks. The two shocks generated after the second core passage can be clearly seen. The shocks are moving outward from the center. Using this 2D pressure cut, we can derive the Mach numbers of these shocks to obtain ${\cal M}\sim 3.2$. At this epoch, the merging cluster is in an infalling phase, where the two dark-matter components already turned over and are moving toward each other, approaching the third core passage.

In the right panel of Figure 5, we display the magnitude of the projected pressure gradient using our best model, which is rotated out of the plane of the sky with a polar angle of $72.4^{\circ}$. This panel shows where the shocks could be observed. The two well-separated shock waves form entangled circles in projection. We find that the closest shock surface to the cluster center is on the east side, about 890 kpc away from the center. We note that, owing to projection effects, the change in the projected pressure is not equal to the physical pressure change due to the shock. The Mach number derived from the pressure drop in this projection is ${\cal M}\sim 1.2$, much less than the physical Mach number. We expect that observations of this outgoing shock would measure a Mach number close to this value (for a discussion on how much bias is caused by projection effects in the Mach numbers derived from X-ray observations, see, e.g., Molnar & Broadhurst 2017). The position of the closest shock in the east is in rough agreement with the position of the X-ray surface brightness edge we found by reanalyzing the Chandra observations of A370 (at $\sim 690$ kpc; Section 2.1; see Figure 5b). Our simulations suggest that the X-ray surface brightness edge found in A370 on the east from the cluster center is an outgoing shock generated after the second core passage. According to our best model, the outgoing shock on the west of the cluster center is located much farther out than the X-ray surface brightness edge found by Botteon et al. (2018). Reanalyzing the Chandra data, we found no significant X-ray edge on the east from the cluster center, which is in agreement with our best model.

As a result of the outgoing shocks generated after the second core passage, the gas density decreases in the central region. In a later phase, the gas falls back toward the center. We display in Figure 6 the ratio $\Sigma_{\mathrm{best}}/\Sigma_{\mathrm{2CP}}$ between the gas mass surface density $\Sigma_{\mathrm{best}}$ at the best epoch (i.e., our best model) and that $\Sigma_{\mathrm{2CP}}$ at the epoch of the second core passage based on our best run (RP100V3.5). These projected maps are obtained using the viewing angle of our best model within 2.5 Mpc from the cluster center. The red horizontal bar in the lower-left corner marks a scale of 1 Mpc, which is about the size of the BOLOCAM extraction region around the cluster center ($R\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}R_{2500}$; central blue region). At the best epoch after the second core passage, the cluster passed the second turnover, while the gas is still moving outward as a result of the shock, which is located about $1$–$2$ Mpc from the cluster center (see Figure 5). Accordingly, the gas is depleted in the central area where $\Sigma_{\mathrm{best}}/\Sigma_{\mathrm{2CP}}\sim 0.7$ (dark blue region). This depletion of the gas from the central part of the system explains why the X–ray emission and SZ amplitude of this cluster is well below those suggested by cluster mass scaling relations (Section 1).

We have used our existing library of binary merging cluster simulations performed in our previous work and carried out self-consistent $N$-body/hydrodynamical simulations with flash to study the dynamical state of the massive merging cluster A370. The cluster is a superlens system characterized by its large Einstein radius. Our simulations were constrained by X-ray, SZ-effect, gravitational lensing, and optical spectroscopic observations. Specifically, we utilized the locations of the two mass peaks derived from strong gravitational lensing (Strait et al., 2018), the X-ray morphology (Section 2.1), the relative LOS velocity (Section 2.3), and the amplitude of the integrated SZ effect (Section 2.2; Czakon et al. 2015).

We preformed flash simulations by fixing the masses of the two cluster components based on weak- and strong-lensing measurements (Umetsu et al., 2011a). We used different impact parameters, infall velocities, and concentration parameters to find the best match with the multi-wavelength observations. We found that our best model with masses of $M_{1}=1.7\times 10^{15}\,M_{\odot}$ and $M_{2}=1.6\times 10^{15}\,M_{\odot}$, an impact parameter of $P=100$ kpc, an infall velocity of $V_{\mathrm{in}}=3500$ km s^-1, and concentration parameters of $c_{1}=6$ and $c_{2}=6$ can explain the main features of the X-ray morphology and the mass surface density and simultaneously satisfy the constraints from the observed relative LOS velocity and the integrated SZ amplitude. Moreover, our best model reproduces the observed velocity dispersion of cluster member galaxies, which supports the large total mass of A370 derived from weak lensing. Our results strongly suggest that A370 is a post major merger observed 0.57 Gyr after the second core passage, just before the third core passage.

It is intriguing to note that, similar to the case of A370, Cl0024$+$1654 is another superlens cluster at $z=0.395$ that is very faint in X-ray and SZ signals (Umetsu et al., 2011a, b). A careful interpretation of X-ray, lensing, and dynamial data based on flash simulations suggests that Cl0024$+$1654 is also a post major merger occurring along the line of sight, viewed approximately 2–3 Gyr after impact before the gas has recovered (Umetsu et al., 2010).

Our simulations represent the first attempt to model one of the most massive dynamically active merging galaxy clusters, A370, using dedicated self-consistent $N$-body/hydrodynamical simulations. In order to improve our dynamical model for A370, deeper X-ray observations would be necessary to verify the positions of the mean and secondary peaks, obtain a spatially resolved temperature distribution of the intracluster gas, and verify the shock feature in the east from the cluster center. Moreover, a detailed optical spectroscopic study would be needed to precisely determine the relative LOS velocity of the two main merging components. With these proposed new observations a more extensive parameter search would be worth while to pursue, which is more likely feasible using GPU clusters.