"Beads-on-a-string" Star Formation Tied to One of the Most Powerful Active Galactic Nucleus Outbursts Observed in a Cool-core Galaxy Cluster

In Figure 6, we present a Pan-STARRS i-band wide-field view of SDSS 1531 (left) and an HST view of the cluster's central region (right) with the background-subtracted, exposure-corrected, wavelet-fit Chandra X-ray surface brightness map of the ICM overlaid. The X-ray emission is centered on the BCGs, and has an angular inclination similar to that of the BCGs in the optical. The emission has a predominantly smooth and circularly symmetric surface brightness distribution, except for two thin, 20–30 kpc "wings" to the southeast and southwest, which create a concave feature.

Zoom In Zoom Out Reset image size

We searched for potential cavities in the unsharp-masked image shown in Figure 7 (see Section 2.1 for full details on the image creation). While helpful for identifying hidden structures, unsharp-masked images can also give the illusion of false cavities. Therefore, we only consider unsharp-masked image fluctuations consistent with those in the surface brightness map. The most plausible cavity candidate lies in the ∼16 kpc radius, concave surface brightness edge created by the two X-ray "wings." However, the radio data presented in Section 3.2 suggest that this feature may instead represent an opening to a large (∼50 kpc radius) avocado-shaped cavity that the current X-ray data are too shallow to resolve.

Zoom In Zoom Out Reset image size

To investigate the region surrounding the "wings" further, we extracted a surface brightness profile from the 0.5 to 7.0 keV Chandra image in the pie regions shown in Figure 8 (left). Figure 8 (center) shows the resulting surface brightness profile. The clear density jump (green dashed line) at ∼12.6 kpc confirms the presence of a surface brightness edge.

Zoom In Zoom Out Reset image size

To test whether the density jump represents a shock or a cold front, we extracted spectra of the ICM inside and outside the discontinuity in the regions shown in Figure 8, right. To obtain projected thermodynamic properties, we fit the 0.5–7 keV band spectra with a phabs * apec model. We measured a clear density jump of n_e,in/n_e,out = 1.40 ± 0.05. We tentatively find continuous pressure at the surface brightness edge (p_out/p_in = 0.94 ± 1.22), and that the ICM temperature close to the cluster center is lower than it is outside the surface brightness edge (kT_in/kT_out = 0.67 ± 0.87). These characteristics are consistent with the properties expected of a cold front. However, it is important to note that these results are not conclusive since there are only ∼200 net counts across the selected regions, resulting in significant errors. Therefore, it is still possible that the observed discontinuity could be better explained by a shock rather than a cold front.

We created high-resolution thermodynamic spectral maps of the ICM using clusterpyxt (see Section 2.1 for details). The temperature map (Figure 9, left) reveals the presence of a central, extended (∼200 kpc), cool (∼2–3 keV) butterfly-shaped structure, with its major axis perpendicular to that of the central surface brightness map. The coolest gas is found in a curved trail extending from the northeast to the southwest. Outside the cool region, temperatures gradually increase to ∼4 keV, with ∼8% uncertainties. The pseudopressure (see Figure 9, right) and entropy maps (not shown) exhibit a uniform circular distribution within the central ∼100 kpc, with the lowest entropy and highest pressure near the center of the cluster.

Zoom In Zoom Out Reset image size

We also created thermodynamic profiles of the larger cluster atmosphere, shown in Figure 10. In Table 2, we list the derived X-ray properties for SDSS 1531, namely the total cluster mass M_Δ, radius R_Δ, central temperature kT₀, central entropy K₀, central pressure P₀, and central cooling time t_cool,0. Each property is consistent with that expected of a strong cool core.

Zoom In Zoom Out Reset image size

In Figure 10 (top left), we plot the projected emission measure profile for ObsID 17218 (green) and ObsID 18689 (black). The emission measure profile is fit with the modified β-model described by Equation (4). Both profiles show excellent agreement between the independent observations and exhibit a central overdensity, suggesting the presence of a cool core. At low redshift, cool-core clusters have characteristic central cusps in their X-ray surface brightness distribution. Vikhlinin et al. (2007) characterizes the central cusp as the power-law index of the gas density profile $\alpha =-2d\mathrm{log}{\rho }_{g}/d\mathrm{log}r$ at r = 0.04 R₅₀₀, with clusters known to host strong cool cores having α > 0.7. From the fit to the emission measure, we obtain α = 1.4, satisfying the cuspiness criteria for the presence of a strong cool core.

In Figure 10 (bottom left), we plot the projected temperature profile, fit with the analytical model given by Equation (5). The density peak in the emission measure profile coincides with a significant temperature decline inside the central 100 kpc, from T = 4.6 ± 0.22 keV at r ∼ 110 kpc to T = 2.4 ± 0.03 keV at r ∼ 10 kpc. We also find a three-dimensional temperature drop ${T}_{\min }/{T}_{0}\sim 0.25$. For comparison, strong cooling flow clusters like A478 and A133 have ${T}_{\min }/{T}_{0}$ ranging from 0.1 to 0.4 (Vikhlinin et al. 2007).

We compared the derived pressure and entropy profiles (see Figure 10, top center and top right), with those from the Archive of Chandra Cluster Entropy Profile Tables (ACCEPT; Cavagnolo et al. 2009; hereafter C09) for 239 galaxy clusters. The profiles are consistent with those of similar cool-core clusters.

According to the cooling time profile (Figure 10, bottom center), the ICM cools within a Hubble time (t_H ) out to a radius of ∼100 kpc, with the innermost ∼30 kpc cooling within <1 Gyr. Another relevant timescale for cool-core clusters is the freefall time t_ff, which describes the characteristic time for the system to collapse. When t_cool > t_ff, the system has sufficient time to reestablish hydrostatic equilibrium. However, t_cool < t_ff signals catastrophic cooling, for the gas is unable to rebound fast enough to account for the rapid loss in pressure (Nulsen 1986; McCourt et al. 2012). Numerical simulations of thermal instabilities in cluster atmospheres have found the t_cool/t_ff ratio is less stringent, with instability occurring when t_cool/t_ff ≲ 10 (Gaspari et al. 2012; Voit & Donahue 2015). The freefall time is estimated as:

$$\begin{eqnarray}&&{t}_{{\rm{ff}}}(r)=\sqrt{\displaystyle \frac{2{r}^{3}}{\left.{GM}(r\right)}},\end{eqnarray} \tag{ 7 }$$

where r is the distance from the center of the cluster and M(r) the radial cluster mass profile. Since the freefall time is most informative near the center of the cluster, we consider two mass profiles: one derived from a strong-lensing analysis (Sharon et al. 2014) and the other from hydrostatic X-ray measurements (see Section 2.1). The strong-lensing analysis yields a cylindrical mass of M( < R_E) = 2.32 ± 0.01 × 10¹³b M_⊙ within the Einstein radius R_E ∼ 52.1 ± 0.5 kpc. On the other hand, the hydrostatic X-ray mass measurement, obtained at R = 52 ±7.1 kpc, gives M_HE( ⪅ R_E) = 9.3 × 10¹² M_⊙. To compare the X-ray mass with the strong-lensing cylindrical mass, we recalculate the hydrostatic X-ray spherical mass within a cylindrical volume of height 10 Mpc as done in Sharon et al. (2015). This yields a cylindrical mass of ∼1.3 × 10¹³ M_⊙. Though small in magnitude on the scale of the total mass, the factor of two difference between the X-ray and strong-lensing mass estimates, coupled with the disturbed morphology of the surface brightness map, indicates that the hydrostatic equilibrium assumption does not accurately describe the cluster's core, as expected given the ongoing major merger.

Figure 10 (b'ottom right) illustrates the enclosed t_cool/t_ff profiles estimated from the hydrostatic mass (thick pink line and black points) and the strong-lensing mass reprojected to a spherical mass (brown points). Both estimates lead to a minimum t_cool/t_ff ∼ 20 for r ≤ 10 kpc. Although the theoretically expected t_cool/t_ff ratio often falls near or below unity (Gaspari et al. 2012; McCourt et al. 2012; Sharma et al. 2012), the observed behavior in cool-core clusters typically deviates from this threshold, with most cool-core clusters harboring cold filaments when t_cool/t_ff ≤ 20 (Olivares et al. 2019). Other works have suggested that thermally unstable cooling can occur when this ratio lies above unity, in the range of 10−40 (Gaspari et al. 2012; McCourt et al. 2012; Sharma et al. 2012).

The LOFAR 144 MHz image, shown in Figure 11, reveals several sources of diffuse radio emission throughout the cluster, few of which are marginally visible in the VLA images. From this image, we have identified and labeled six distinct radio sources affiliated with the cluster.

Zoom In Zoom Out Reset image size

Of the six radio sources, Source C is most closely aligned with the BCGs. Its total angular size is 28'', corresponding to a physical size of 133 kpc. The emission takes on an avocado-esque shape and increases in intensity as it extends southeastward from the BCGs. The emission peaks ∼60 kpc away from the southern BCG, forming a distinct "hotspot." Notably, the source is undetected in all available VLA images.

Source D is a hybrid blob/extended source with an angular size of 22'', corresponding to a physical length of ∼105 kpc. The ≥6σ emission resembles the avocado-like shape of Source C, with a similar hotspot near the edge of the source. The >3σ emission also allows for an interpretation where Source C resembles a miniature head–tail source whose tail curves around the "head," where the emission peaks. The source has no confirmed counterpart but is close to a galaxy at z = 0.3357. Directly north of Source D is Source E, the faintest source detected. Source E also does not have a direct optical counterpart, but its proximity to Source C and the affiliated galaxy suggests it could belong to the cluster. The source has a similar "avocado-shaped" morphology and a very small hotspot detected at 6σ.

Source B, the largest source within the cluster, has a distinct head–tail morphology. The "head" of the source comprises a concentrated emission region, with two objects in the center, one of which has z = 0.341. The tail extends outward up to 45'', corresponding to a physical length of ∼215 kpc. Only the head of the source is also detected in the VLA B-array image.

Source A is the most luminous radio source detected and located very close to the cluster's virial radius (r₂₀₀). The emission has an oval morphology and peaks in the center. The source is brightly detected in all VLA images. Several objects appear associated with the emission, one of which has a redshift of z = 0.3376. Conversely, Source F is the faintest source associated with the cluster and is located opposite Source A near the virial radius. The round, compact source is juxtaposed with two foreground objects, one of which is located at z = 0.3298.

Figure 12 displays ≥3σ LOFAR radio contours of Source C superimposed on the Chandra X-ray surface brightness image of the cluster's central field. The most concentrated radio emission, detected at ≥12σ and highlighted in blue, aligns remarkably well with the region delineated by the X-ray "wings" and is connected by a bridge of radio emission to the southern BCG. The spatial coincidence is reminiscent of radio AGN lobes that push aside the X-ray gas to create cavities, as observed in several other cool-core clusters with active AGN, such as Perseus (Fabian et al. 2000) and A2597 (McNamara et al. 2001). If the concave region is indeed a cavity, the Chandra observations are likely too shallow to resolve the full extent uncovered by the high-resolution LOFAR image. However, the fact that the radio emission is undetected by the VLA suggests that we may be observing an aged AGN lobe, which we discuss further in Section 4.1.3.

Zoom In Zoom Out Reset image size

None of the remaining sources identified in the cluster have readily identifiable X-ray counterparts, as they are located in the less dense portions of the ICM.

Various models have been used to describe the origin of diffuse radio emission in clusters, a list of which is discussed in Kempner et al. (2004). To make preliminary classifications of the radio sources in the cluster, we estimate the spectral index, α, from 144 MHz to 1.4 GHz for each source using the equation:

$$\begin{eqnarray}&&{\alpha }_{144{\rm{M}}{\rm{H}}{\rm{z}}}^{1.4{\rm{G}}{\rm{H}}{\rm{z}}}=\displaystyle \frac{\left.{\rm{l}}{\rm{o}}{\rm{g}}({S}_{1})-{\rm{l}}{\rm{o}}{\rm{g}}({S}_{2}\right)}{\left.{\rm{l}}{\rm{o}}{\rm{g}}({\nu }_{1})-{\rm{l}}{\rm{o}}{\rm{g}}({\nu }_{2}\right)},\end{eqnarray} \tag{ 8 }$$

where S₁ is the flux density of the source at ν₁ = 1.4 GHz, and S₂ is the flux density of the source at ν₂ = 144 MHz.

To obtain integrated flux densities in the LOFAR and VLA B-array images, rectangular regions were defined around each source. For sources undetected in the VLA image, we assume that the flux densities of the undetected sources were no more than 3σ_rms at 1.4 GHz (σ_rms = 0.135 mJy beam⁻¹), setting an upper limit of 0.052 mJy beam⁻¹. The rms noise at 144 MHz and 1.4 GHz were used to estimate the spectral index error. Table 3 lists all integrated flux densities and derived spectral indices for each source.

Table 3. Integrated Spectral Indices

Source	ν	Flux Density	${\alpha }_{144\mathrm{MHz}}^{1.4\mathrm{GHz}}$
	(MHz)	(mJy)
A: bright round source	144	73.1	−0.15 ± 0.005
	1400	52.5
B: head–tail source	144	13.1	−0.19 ± 0.03
	1400	8.6
C: central avocado-shaped source	144	26.0	−1.7 ± 0.4
	1400	<0.5
D: off-center amorphous source	144	9.9	−1.3 ± 0.4
	1400	<0.5
E: faint amorphous source	144	2.4	−0.7 ± 0.4
	1400	<0.5
F: faint round source	144	1.9	−0.6 ± 0.4
	1400	<0.5

Note. Integrated spectra for each radio source identified in Figure 11. The integrated spectral index is calculated from 144 MHz to 1.4 GHz. For sources undetected with the VLA, the 3σ = 0.5 mJy upper limit is used to estimate the spectral index.

Download table as:  ASCIITypeset image

Sources A and B have relatively flat spectra with α < −0.20. Given Source B's tailed morphology, optical counterpart at z = 0.336, and flat spectral index, the source is likely a radio jellyfish (e.g., Roberts et al. 2021). Source A may be an early stage radio jellyfish.

The integrated spectral index estimated for Source C is −1.7 ± 0.4. Given the diffuse emission's projected location in the cluster's central region, steep spectral index, striking overlap with the concave X-ray surface brightness edge, and emission bridge connecting the southern BCG to the emission peak, the bulk of the emission in Source C likely stems from an AGN relic. Source D also has a steep spectrum, with α ∼ −1.3 ± 0.4. If Source D is not a separate radio source, it could be Source C's missing relic lobe counterpart, a scenario we explore further in Section 4.1.3.

Sources E and F have moderately steep spectra. Their small angular sizes make them difficult to classify, so we refrain from doing so.

Higher-resolution observations at 1.4 GHz and other frequencies are required for precise spectral analysis, accurate classification, and detailed study of each newly identified radio source.

We present the Hα flux, velocity, and velocity dispersion maps of the warm (<10⁴ K) ionized gas derived from the continuum-subtracted GMOS gas cube in Figure 13. The ionized gas surrounds the YSCs, with the brightest knots of emission coincident with the young stars. The gas is marginally detected throughout the nuclei of the BCGs. The emission likely extends beyond the GMOS FOV (depicted as a rectangular box), but the emission near the edges lies within the noise of the data. It is unclear how much further out the emission extends and whether it envelops the full extent of the molecular gas, which we discuss later on in Section 3.4.5.

Zoom In Zoom Out Reset image size

The velocity map of the gas reveals that the ionized gas is significantly redshifted with respect to the stellar systemic velocity of the southern nucleus. The gas enshrouding the southern string of beads is redshifted up to +800 km s⁻¹ with respect to the systemic velocity, while the gas tracing the northern strings closely approaches the systemic velocity of the southern nucleus. The FWHM map shows that the gas is most disturbed near the southern nucleus.

From the GMOS data, we measure a total Hα flux of 1.9 ± 0.03 × 10⁻¹⁶ erg s⁻¹ cm⁻² across all spaxels. Although each "bead" of star formation is likely embedded within dust supplied from the aging stars of the BCGs (Donahue & Voit 2022), we are unable to estimate the internal extinction for each spaxel reliably via the Balmer decrement (Hα/Hβ). Although both Hα and Hβ emission lines are clearly detected, the Balmer decrement yields a median value of ∼2.2, an unphysical value given the intrinsic line ratio of 2.86 appropriate for low-density gas with a temperature of 10⁴ K under the standard Case B recombination scenario (Osterbrock & Ferland 2006). Such a flat Hα/Hβ ratio is expected for environments where there is no dust attenuation or n_e > 10⁹ cm³ (Adams & Petrosian 1974), neither of which are commonly expected in cool-core clusters. The unphysical value may instead be due to difficulty in accurately subtracting the relatively weak stellar continuum, as described in Section 2.3. Without correcting for extinction, we use Equation (2) in Kennicutt (1998) to obtain a total SFR of 0.60 ± 0.009 M_⊙ yr⁻¹, which we interpret as a lower limit.

Tremblay et al. (2014) estimated an extinction-corrected SFR in the range of 5–10 M_⊙ yr⁻¹. They measured an SFR of ∼5 ± 2 M_⊙ yr⁻¹ from correcting the ALFOSC Hα luminosity for extinction using the Balmer decrement from SDSS DR10's reported Hα and Hβ fluxes ²² and an SFR of 9.55 ± 3.4 M_⊙yr ⁻¹ when matching spectral energy distribution templates to extinction-corrected SDSS ugriz magnitudes. We take the SFRs uncorrected for extinction as a lower limit on the SFR in the system. Without correcting for extinction, ALFOSC obtained a Hα flux of 1.1 × 10⁻¹⁵ erg cm ⁻² s⁻¹, which corresponds to an SFR of 3 M_⊙ yr⁻¹. We thus estimate the total SFR within the beaded star formation complex to be between 1 and 10 M_⊙ yr⁻¹.

To identify potential sources of ionizing radiation in the star formation complex, we utilize the broad wavelength coverage of the optical spectra to compare key diagnostic emission-line ratios. We created a spatially resolved Baldwin–Philips–Terlevich (BPT) diagram (Baldwin et al. 1981), namely the [N ii] BPT diagram (O iii λ5007/Hβ versus [N ii] λ6583/Hα), shown in Figure 14 left). The classic [N ii] BPT diagram is used to distinguish sources of ionizing radiation, namely between H ii regions and star formation, Seyferts, or a combination of a Seyfert and star formation. The solid line represents the demarcation line of pure star formation (Kauffmann et al. 2003), the dashed line represents the demarcation line of extreme star formation (Kewley et al. 2001), and the dashed–dotted line represents the empirical division between low-ionization nuclear emission-line region (LINER) and Seyfert-like sources (Schawinski et al. 2006). We have color coded the data points based on the regions in which they sit.

Zoom In Zoom Out Reset image size

Most of the spaxels lie inside the star formation region. Few spaxels lie near the outer edges, where the signal is weakest, in the "composite" and "LINER" regions. This suggests that the warm nebula is primarily ionized by the YSCs, with some potential contribution from AGN activity. We do not overinterpret the BPT classifications since cool-core BCGs likely have several ionization sources (Ferland et al. 2009; McDonald et al. 2012).

The morphology of the cold (<10² K) molecular gas is shown in Figure 15, which displays the masked CO(3–2) integrated intensity map (moment 0), the intensity-weighted mean velocity field (moment 1), and the intensity-weighted velocity dispersion (moment 2). The integrated intensity map pieces clumps of molecular gas together to reveal a cloud-like structure spanning ∼24 kpc across the sky, with three bright emission peaks to the north, southwest, and southeast. Surrounding the bulk of the emission are small pieces belonging to a smoother gas distribution below the >3.5σ threshold mask. The ALMA observations do not resolve the 19 individual clumps detected in the NUV, hindering the opportunity for a clump-by-clump analysis. Instead, the emission appears concentrated in three filaments and clumps that are spatially interconnected yet distinct in velocity space. These clumps are shown in Figure 18 and further discussed in Section 3.4.2.

Zoom In Zoom Out Reset image size

The molecular gas is peculiarly spatially offset with respect to the central BCGs and the YSCs. To quantify the offset, we compare the morphology of the gas in Figure 16 to that of the YSCs in the NUV (HST F390W). The emission in the 220 km s⁻¹ channel resembles that of the northern YSCs. Although the molecular gas overlaps with the majority of the YSCs in the northern filament, the northernmost string of beads is offset from the peak of the molecular gas by ∼06 (∼2.8 kpc). Furthermore, the bulk of the emission remains shifted to the east. Unlike the northern filament, the central clump is almost completely decoupled from the YSCs, and resembles none of the YSCs morphologically. The emission in the 660 km s⁻¹ channel matches the morphology of the southern string of beads. The beads and the molecular gas are separated by ∼03 (∼1.4 kpc), with the molecular emission lying to the southeast. In the mostly decoupled southern corner, the cold gas appears to fit into the southwestern bay-like feature like two pieces of a puzzle.

Zoom In Zoom Out Reset image size

It is important to note that the ALMA observations reveal faint, smooth emission across most of the YSCs, near the southern nuclei of the BCGs, and in between the BCGs. However, this emission lies below the threshold applied to all CO(3–2) maps presented in this paper (see Section 2.4). The detection of some faint molecular emission in the inner extent of the merging BCGs and across most of the YSCs suggests that if we observed the system at greater depths with ALMA, we might detect CO(3–2) completely enveloping the YSCs and certain areas of the BCGs. However, the bulk of the emission would most likely still lie to the east of the merging system.

Regardless, the coincidence between the morphologies of the YSCs and that of the molecular gas suggests that they are associated. However, the question remains: why is the molecular gas spatially offset from the young stars? One possibility is that the molecular gas that originally enveloped the YSCs collapsed to form stars, which then ionized the surrounding gas. However, this does not explain why the collapse was localized to one side of the gas. We discuss potential scenarios to explain the gas' relation to the YSCs in Section 4.1 and describe its relation to the warm ionized gas in Section 3.4.5.

The first and second moment maps in Figure 15 reveal complex velocity structure spanning the cold molecular gas' extent. The projected LOS velocities across the structure are significantly redshifted by ∼+700 km s⁻¹ with respect to the stellar velocity of the southern BCG. The velocity of the gas peaks near the southern string of beads, with a corresponding velocity dispersion peak of ∼150 km s⁻¹. The gas with the largest velocity dispersion traverses a diagonal path from the end of the southern string of beads to just beneath the northern string of beads.

To analyze the gas flows in greater detail, we define a rectangular region that spans the total extent of the molecular gas and average the emission in each velocity channel over the width of the region to produce a position–velocity (PV) diagram, shown in Figure 17. The PV diagram reveals a smooth decrease in velocity from the southern to the northern nuclei, with distinct velocity components present across the structure. There is one primary peak in the northern filament (Position 2''–5''), one in the central clump (Position 1''–2''), and another in the southern filament (Position 05–15).

Zoom In Zoom Out Reset image size

The velocity of the gas in the northern filament, central clump, and southern filament ranges mostly from +60–260 km s⁻¹, +280–460 km s⁻¹, and +480–680 km s⁻¹, respectively. We present the intensity-weighted velocity map of each structure defined by the aforementioned velocity ranges in Figure 18. The smooth velocity gradient across the total structure could indicate positive or negative radial velocities. This suggests that the molecular gas clumps in this system experienced radial velocity changes at different epochs. Specifically, the gas in the northern filament is closely associated with the stellar systemic velocity, while the gas in the southern filament appears more consistent with the average velocity for seven other galaxies found within ∼600 kpc of the cluster core (Bayliss et al. 2011). This may indicate the presence of molecular gas that is mostly influenced by the bulk motions of the intracluster gas (near the southern filament), slowly entering the gravitational potential well of the BCGs (near the northern filament). We provide a more detailed interpretation of the gas motions in Section 4.1.4.

Zoom In Zoom Out Reset image size

Assuming a CO(3–2)/CO(1–0) line ratio of ∼0.8 (Edge 2001), we estimate the mass of the molecular H₂ nebula by following the relation reviewed by Bolatto et al. (2013):

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{mol}} & = & \left(\displaystyle \frac{1.05\times {10}^{4}}{0.8}\right)\,\left(\displaystyle \frac{{X}_{\mathrm{CO}}}{{X}_{\mathrm{CO},\,\mathrm{MW}}}\right)\\ & \times & \left(\displaystyle \frac{1}{1+z}\right)\left(\displaystyle \frac{{S}_{\mathrm{CO}}{\rm{\Delta }}v}{\mathrm{Jy}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{{D}_{{\rm{L}}}}{\mathrm{Mpc}}\right)}^{2}{M}_{\odot },\end{array}\end{eqnarray} \tag{ 9 }$$

where S_COΔv is the integrated CO(3–2) intensity, z is the galaxy redshift (z = 0.335), and D_L is its luminosity distance (1767 Mpc in our adopted cosmology). The Galactic CO-to-H₂ conversion factor X_CO dominates the uncertainty in this relation (e.g., Bolatto et al. 2013). For SDSS 1531, we adopt the average value for the disk of the Milky Way, X_CO =X_CO,MW = 2 × 10²⁰ cm ${\left({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{-1}$, which has a ∼30% uncertainty. The true value of the conversion factor largely depends on the gas metallicity and whether the CO emission is optically thick. It is unclear whether the X_CO factor measured in the Milky Way and nearby spiral galaxies should match that of BCGs, as it often varies wildly in ultraluminous infrared galaxies (Bolatto et al. 2013). However, Vantyghem et al. (2017) recently reported one of the first detections of ¹³CO in a BCG (RX J0821+0752) and found an X_CO factor only a factor of two lower than that of the Milky Way. To account for our limited information about the system's metal abundance and to facilitate direct comparison with the molecular gas masses derived in other studies of BCGs (e.g., Russell et al. 2016; Tremblay et al. 2018; Olivares et al. 2019; North et al. 2021), we adopt X_CO = X_{CO, MW} as our most reasonable choice. A caveat to this selection is that we may overestimate the total molecular mass by a factor of a few. This should be interpreted as the overriding uncertainty on all mass estimates quoted below.

To estimate S_COΔv, we fit a four-component Gaussian line to the CO (3–2) spectrum extracted from an elliptical aperture containing all ≥3σ emission in the calibrated cube binned to 20 km s⁻¹ channels (the spectrum is shown in the left panel of Figure 18). The fit, from −690 to 900 km s⁻¹, yielded an emission integral of ∼1.7 ± 0.3 Jy km s⁻¹. Noting the aforementioned caveats, these results give an H₂ gas mass of (3.7 ± 0.7) × 10¹⁰ M_⊙. With natural weighting, we obtain integrals ranging from 0.88 to 2.5 Jy km s⁻¹ depending on the different binning selected, yielding mass estimates that range from 1.9 to 6.5 × 10¹⁰ M_⊙. ²³

The remaining panels in Figure 18 show the spectra of the gas across the northern filament, central clump, and southern filament regions. The fit to each structure's spectrum indicates that most of the gas mass is concentrated in the northern filament (∼45%), and distributed relatively evenly throughout the central clump and southern filament (∼32.5% each).

ALMA detected no continuum within the vicinity of the cluster's BCGs in any of the three line-free spectral windows placed in Band 6 (see Section 16). From the standard deviation of the low-resolution data, we place 2σ upper limits of 4.87 × 10⁻⁴ mJy beam⁻¹ on the continuum. If detected, the continuum would likely have originated from thermal emission from dust and/or synchrotron and hot dust emission from the central AGN. The nondetection is, therefore, consistent with the picture of an AGN that is currently inactive.

To place an upper limit on the mass of dust present, we modeled the emission S_ν as a modified blackbody using the following equation:

$$\begin{eqnarray}&&{S}_{\nu }=\displaystyle \frac{\left.{\kappa }_{\nu }{M}_{d}B(\nu ,{T}_{d}\right)}{{D}^{2}},\end{eqnarray} \tag{ 10 }$$

where M_d is the dust mass, B(ν, T_d ) is the Planck function, which depends on frequency ν and dust temperature T_d , and D is the distance to the galaxy. κ_ν is the dust absorption coefficient, described by a power law with dust emissivity index β such that κ_ν ∝ ν^β . Here, we utilized an empirical κ_ν , where κ_500μm = 0.051 m² kg¹, and β = 1.8 (Clark et al. 2016). Assuming a dust temperature of 25 K (e.g., Davis et al. 2017), we obtain a maximum dust mass of 3.7 × 10⁸ M_⊙ (and thus a molecular gas-to-dust ratio of ≳100).

Given that young stars form in molecular clouds, one would expect the stellar superclusters to be deeply embedded within the cold gas from which they formed and that the cold gas is cospatial with the warm ionized gas, as previous analyses of cool-core clusters have observed (e.g., Tremblay et al. 2016, 2018; Vantyghem et al. 2016; Olivares et al. 2019). A direct comparison between the integrated ALMA CO(3–2) intensity map and GMOS Hα observations (Figure 19, left) reveals that the molecular and ionized gas have similar physical extents, with the ionized gas extending slightly further than the molecular gas within the GMOS FOV. However, it remains unclear whether the ionized gas fully enshrouds the entire molecular structure with the current observations. It is, nonetheless, evident that the molecular gas lies ∼3 kpc to the southwest of the ionized gas (Figure 19, left), with the offset being too large to be attributed to uncertainty in either ALMA or Gemini-North's respective astrometric reference frame.

Zoom In Zoom Out Reset image size

The aforementioned studies also found the two gas components tend to be comoving. To explore this possibility, we plot the difference in velocity and the velocity dispersion ratio between the overlapping CO(3–2) and Hα gas in Figure 19 (center, right). Using the REPROJECT module from the ASTROPY package, the ALMA maps were resampled onto GMOS pixel grids. The velocity difference map shows that the ionized and molecular gas exhibit similar overall velocity structure within ±160 km s⁻¹, which is within the GMOS data's ∼160 km s⁻¹ velocity resolution. This indicates that the cold and warm ionized gas phases are likely comoving, consistent with prior studies.

The velocity dispersion ratio map shows that the Hα velocity dispersion is mostly consistent with the CO(3–2) or slightly broader, consistent with multiwavelength studies of cool-core clusters (Olivares et al. 2019). Tremblay et al. (2018) suggested that lines of sight are more likely to intersect warm gas than cold clouds as an explanation for the broader velocity distribution of the warm ionized gas compared to the cold gas.

To estimate the amount of time needed to deplete the entire ∼10¹⁰ M_⊙ reservoir of molecular gas, we adopt an SFR of 1–10 M_⊙ yr⁻¹. We calculate the depletion time as ${t}_{\mathrm{dep}}=\tfrac{{M}_{\mathrm{gas}}}{\mathrm{SFR}}$ and recover timescales ranging from ∼1 to 10 Gyr, which are approximately within a Hubble time (t_H = 9.728 Gyr). The majority of cool-core BCGs typically have gas depletion timescales of ∼1 Gyr (O'Dea et al. 2008), except starburst BCGs, which can have depletion timescales as short as ≤30 Myr (McDonald et al. 2014). Both ends of the depletion timescale range for SDSS 1531 lie above ∼1 Gyr, implying that the star formation efficiency in this system is lower than in other BCGs. This is likely because only the BCG-facing border of the molecular gas has collapsed to form stars, which we discuss further in Section 4.1.5.

Although the two central galaxies in SDSS 1531 are definitively engaged in a major merger, determining whether this event is part of a larger-scale subcluster merger demands more nuance. Sharon et al. (2014) find that the mass distribution of the cluster is well described by two cluster-scale halos. Moreover, each central elliptical is of roughly equal stellar mass, and there is an apparent bifurcation in the optical line redshift distribution for 14 of the nearby cluster members for which there is existing optical spectroscopy/photometry (see Figure 21). Although the current redshift sample is too limited to draw any definitive conclusions, the presence of diffuse radio sources throughout the cluster signals broader systemic dynamical activity. Such cluster-scale synchotron emission is almost exclusively found within systems dynamically disturbed by mergers.

Zoom In Zoom Out Reset image size

If the merging ellipticals in SDSS 1531 result from a cluster-scale merger, it is likely in the last stages of it. The large-scale X-ray gas distribution is relatively smooth and relaxed, with no secondary galaxy concentrations. The exact nature of the merger—whether it involved two clusters of equal size or a massive cluster and a smaller subcluster, or whether it was a head-on or off-center collision—is not discernible from the available data. In the context of the impact on the cool core, however, the literature provides varied insights. Chadayammuri et al. (2021), Valdarnini & Sarazin (2021), and Rasia et al. (2015) find that head-on major mergers significantly disrupt the structure of cool cores, at least for a few gigayears. On the other hand, Poole et al. (2008), ZuHone et al. (2010), and Hahn et al. (2017) find that although cool cores are mostly disrupted by head-on major cluster mergers, low-angular-momentum mergers do not always destabilize the core. Among these varied findings, the consensus seems to be that while mergers can heat cluster cores, the parameters of a merger, such as mass ratios and impact parameters, likely determine the outcome. Currently, SDSS 1531 aligns with the portrait painted by the latter studies, for despite the dynamic activity within the cluster, it has maintained its cool core. Furthermore, the central 30 kpc of SDSS 1531's hot atmosphere still exhibits low entropy (S ≤ 30 keV cm²), thermally unstable (t_cool/t_ff ∼ 20) gas with short cooling times (t_cool < 1 Gyr), characteristics often found in relaxed galaxy clusters with large molecular gas reservoirs, star formation, and nebular emission in their cores (Cavagnolo et al. 2008).

The dynamic environment within SDSS 1531 offers several viable routes to condensation from the hot atmosphere. Defining the "classical" cooling rate as:

$$\begin{eqnarray*}&&{\dot{M}}_{\mathrm{cool}}=\displaystyle \frac{{M}_{\mathrm{gas}}(r\lt {r}_{\mathrm{cool}})}{{t}_{\mathrm{cool}}},\end{eqnarray*}$$

where r_cool ∼ 40 kpc is the radius within which the cooling time is less than 3 Gyr (McDonald et al. 2018), we obtain ${\dot{M}}_{\mathrm{cool}}\sim 185$ M_⊙ yr⁻¹. Assuming a consistent cooling rate over 3 Gyr, the hot atmosphere could easily supply up to ∼10¹² M_⊙ of cool material, with only 100 Myr needed to fill the entire 10¹⁰ M_⊙ gas reservoir.

If the multiphase gas observed directly cooled out of the hot atmosphere, the "circumgalactic precipitation" model (Voit et al. 2017) posits that the growth of small, local thermal instabilities can lead to unstable cooling when t_cool/t_ff < 10. For SDSS 1531, however, t_cool/t_ff > 20 within the cluster core.

On the other hand, the chaotic cold accretion (CCA) model (Gaspari et al. 2018) proposes that the ICM condenses through turbulence driven by AGN outflows, mergers, and on smaller scales, supernovae and stellar winds. In this model, Gaspari et al. (2018) suggest that the ratio of t_cool to the eddy turnover timescale better indicates thermal instability. The eddy turnover time t_eddy is defined as:

$$\begin{eqnarray}&&{t}_{{\rm{eddy}}}=2\pi \displaystyle \frac{{r}^{2/3}{L}^{2/3}}{{\sigma }_{v,L}},\end{eqnarray} \tag{ 11 }$$

where L is the injection length scale of the turbulence and σ_v,L is the velocity dispersion of the turbulence at the injection scale in the ICM. Given that neither of these parameters are directly observable for SDSS 1531, we use the velocity dispersion of the Hα and CO nebulae as a proxy for the velocity dispersion of the ICM. The velocity dispersion of the Hα and CO clouds range from ∼50 to 200 km s⁻¹. Accounting for the conversion from the LOS to a three-dimensional velocity dispersion with a factor of $\sqrt{3}$, we adopt σ_v,L = 80–300 km s⁻¹. We infer the injection length scale L from the extent of the CO emission and the length of the X-ray/radio cavity, ranging from L = 30–70 kpc. Figure 22 shows that the t_cool/t_eddy profile approaches unity throughout the region where CO is observed, supporting the idea that the ICM is thermally unstable and rapidly cooling out to at least 30 kpc.

Zoom In Zoom Out Reset image size

If the thermal instability originated from the hot atmosphere, possibly due to turbulence from the merger, we would expect the cooled gas to be relatively dust free (Donahue & Voit 1993). This is because dust grains sputter rapidly (∼1 Myr) in the ICM and can only form when the gas is shielded from UV and X-ray irradiation (Draine & Salpeter 1979). The GMOS observations found minimal extinction across the extent of the warm ionized component of the multiphase gas, and the ALMA observations placed an upper limit of ∼10⁶ M_⊙ on the dust mass present. These observations, however, do not completely rule out the presence of dust within the cooled gas, in which case would favor an origin from within the BCGs rather than directly condensing out of the hot atmosphere.

In Section 4.1.4, we argue that the observed turbulence-induced thermal instability was most likely engendered by uplift from the massive X-ray/radio cavity, a viable route to condensation in all three precipitation, CCA, and stimulated feedback models.

Figure 20 shows clear spatial alignment between a potential X-ray cavity opening and the low-frequency radio emission, hinting at AGN activity within the cluster. A similar large concave surface brightness discontinuity is also observed in the Ophiuchus cluster (Giacintucci et al. 2020), and on smaller scales in the Perseus Cluster, A1795, A2390, and the Centaurus Cluster, where they were all interpreted as the inner walls of cavities resulting from AGN activity (Walker et al. 2014; Sanders et al. 2016). The concave X-ray surface brightness discontinuity likely does not fully wrap around the radio lobe due to its low surface brightness contrast, thus requiring deeper X-ray observations for successful detection.

The steep radio emission, strongly detected at lower frequencies but absent at 1.4 GHz, likely originates from aged, relic plasma from a powerful AGN outburst in an older stage of the cluster's history. The velocity dispersion of the warm ionized gas is most disturbed near the southern BCG's nucleus. Moreover, a radio emission bridge connects the lobe to the BCG, suggesting the AGN outburst originated there.

If Source C began its life as a buoyant bubble injected near the southern BCG and rose buoyantly to its current location in the plane of the sky at the terminal velocity v_t , we can estimate its age as:

$$\begin{eqnarray}&&{t}_{\mathrm{cav}}=R/{v}_{t}=R\cdot \sqrt{{SC}/2{gV}},\end{eqnarray} \tag{ 12 }$$

where R is the projected distance from the southern BCG to the cavity, V is the volume of the bubble, S is the cross-sectional area of the bubble, C = 0.75 is the drag coefficient (Churazov et al. 2001), and g is gravitational acceleration. We model the potential cavity as an ellipsoidal volume with an axis perpendicular to the plane of the sky and equal to the projected major axes (r_maj = 43.6 kpc) and minor axes (${r}_{\min }=39.6$ kpc).

Following Bîrzan et al. (2004), we calculated the gravitational acceleration using the stellar velocity dispersion of the southern BCG under the approximation that the galaxy is an isothermal sphere, as g ≈ 2σ²/R. SDSS reports two values for the stellar velocity dispersion of the BCG: (1) σ = 444 ±58 km s⁻¹ from the SDSS pipeline's SpecObj table and (2) σ = 356 ± 34 km s⁻¹ from the Portsmouth catalog. We adopt the value reported from the SDSS pipeline as the stellar velocity dispersion and use the ∼100 km s⁻¹ difference between the two measurements as a proxy for the systematic uncertainty in the velocity dispersion of the host galaxy. It would take at least ≥150 Myr for the cavity's hotspot to have traveled to its current distance of ∼67 kpc moving at the terminal velocity (v_t ∼ 750 km s⁻¹).

We can also obtain a rough age of the cavity by assuming that synchrotron radiation and inverse-Compton losses were the only mechanisms by which the relativistic electrons lost energy. However, this method is less reliable because we only have two low-resolution spectral data points. Assuming equipartition of energy between the relativistic particles and the magnetic field (Burbidge 1959), we calculate the minimum magnetic field, pressure, and particle energy of the lobe using the following relations in O'Dea & Owen (1987):

$$\begin{eqnarray*}\begin{array}{rcl}{L}_{\mathrm{rad}} & = & 1.2\times {10}^{27}\,{D}_{L}^{2}S{\nu }_{\mathrm{br}}^{-\alpha }{\left(1+z\right)}^{-(1+\alpha )}({\nu }_{u}^{1+\alpha }-{\nu }_{l}^{1+\alpha }){\left(1+\alpha \right)}^{-1},\\ {P}_{\min } & = & {\left(2\pi \right)}^{-3/7}(7/12){[{C}_{12}{L}_{\mathrm{radio}}(1+k)\phi {V}^{-1}]}^{4/7},\\ {B}_{\min } & = & {[2\pi (1+k){C}_{12}{L}_{\mathrm{radio}}{V}^{-1}{\phi }^{-1}]}^{2/7},\\ {E}_{\min } & = & V\phi {\left(2\pi \right)}^{-1}{\left({L}_{\mathrm{radio}}(1+k){C}_{12}\right)}^{4/7},\end{array}\end{eqnarray*}$$

where L_rad is the radio luminosity in ergs per second, D_L is the luminosity distance in megaparsecs, z is the source redshift, S is the total flux density in Jy, ν is the frequency at which S is measured in hertz, ν_u and ν_l are the upper and lower frequency cutoffs in hertz, respectively, ${P}_{\min }$ is the minimum pressure in dynes cm², k is the ratio of energy density in nonradiating particles to that in synchrotron-emitting particles, V is the source volume, C₁₂ is a constant depending on the spectral index and frequency cutoffs, ϕ is the volume filling factor, ${B}_{\min }$ is the magnetic field at minimum pressure in gauss, and ${E}_{\min }$ is the particle energy (electrons and protons) at minimum pressure in ergs. All equipartition-derived values are given in Table 4.

Table 4. Radio Properties Assuming Equipartition and/or Pressure Equilibrium

Source	Lrad	ϕ	k	${P}_{\min }$	${E}_{\min }$	${B}_{\min }$	te
	(erg s−1)			(dyn cm−2)	(erg)	(μG)	(Myr)
C	1.5 × 1041	1	1	7.7 × 10−13	1.2 × 1058	2.9	54
⋯	⋯	0.1	5499	2.7 × 10−10	4.2 × 1059	54	4.7
D	7.1 × 1040	1	1	1.2 × 10−12	2.8 × 1057	3.5	56
⋯	⋯	0.1	2699	2.7 × 10−10	6.5 × 1058	54	5

Note. Summary of the radio properties derived for Sources C and D assuming equipartition only, and assuming equipartition and pressure equilibrium between the external gas pressure and the magnetized gas.

Download table as:  ASCIITypeset image

The total radio luminosity was estimated assuming the radio spectrum extends from ν_l = 10 MHz to ν_u = 100 GHz with a spectral index of α = 1.7 ± 0.5, and S=0.025 Jy at ν = 144 MHz. Over the cavity's 3.2 × 10⁵ kpc³ volume, we assume the particle energy is equally divided between the electrons and protons (k = 1) and that the lobes are filled with radio plasma, which is justified by the fact that cavities must be mostly empty of thermal gas or they would not be evident in X-ray images (Fabian et al. 2000; McNamara et al. 2000; Blanton et al. 2001). We then estimate the lifetime of electrons in the radio component undergoing both synchrotron radiative and inverse-Compton losses on cosmic microwave background (CMB) photons using the following relation from van der Laan & Perola (1969):

$$\begin{eqnarray}&&{t}_{e}=2.6\times {10}^{4}\displaystyle \frac{{B}^{1/2}}{{B}^{2}+{B}_{R}^{2}}\displaystyle \frac{1}{{\left[\left(1+z\right){\nu }_{b}\right]}^{1/2}}\,{\rm{y}}{\rm{r}},\end{eqnarray} \tag{ 13 }$$

where B is the equipartition magnetic field ${B}_{\min }$ and the break frequency, ν_br, is assumed to be 144 MHz, the lowest available observing frequency. We estimated the magnetic field equivalent to the radiation, which is assumed to be primarily CMB photons, as B_r ≈ 4 × 10⁻⁶(1 + z)² G. With ${B}_{\min }\,=2.9\,\mu {\rm{G}}$, we derive an electron lifetime of ∼50 Myr, which is about 3 times below our initial estimate.

The above calculation, however, yields a minimum pressure of p ∼ 7 × 10⁻¹³ dyne cm⁻², which is three orders of magnitude lower than the external gas pressure derived from the X-ray observations (∼2 × 10⁻¹⁰ dyne cm⁻²). This scenario is clearly unphysical, for the bubble would collapse under the pressure of the external medium. Furthermore, the location of cool X-ray and multiphase gas near the rim surrounding the radio lobe argues strongly against supersonic motion near the lobe boundaries (Nulsen et al. 2002). To maintain equipartition and achieve local pressure equilibrium between the magnetized gas in the radio lobe and unmagnetized gas in the outer rim of the X-ray cavity, we must assume 1 + k ≥ 5500 and reduce the filling factor to ϕ ≤ 10%. This, however, generates a much stronger magnetic field of B_eq ∼ 50 μG, and therefore an unrealistically short plasma age of ∼5 Myr. This either means that equipartition does not hold in the relic lobe, or that there is additional pressure support, potentially from the lobe entraining hot thermal gas. This additional pressure support has been proposed to compensate for the significant pressure differences observed in other sources, such as the intermediate Fanaroff–Riley (FR) type I and II sources in MS0735.6+7421 (Biava et al. 2021) and Hydra A (Croston & Hardcastle 2014).

There are additional puzzles with interpreting Source C as a relic AGN lobe. First, AGN jets typically come in symmetric pairs, producing a pair of radio lobes in the ICM on two sides of the AGN. Assuming Source C represents one relic lobe, its counterpart is missing. Facing a similar dilemma in the Ophiuchus Cluster, which also hosts a massive relic AGN lobe, Giacintucci et al. (2020) speculated that the counterpart may have propagated into a less dense ICM on the other side of the cluster and completely faded away.

While this may also be the case for SDSS 1531, we also explore the possibility that Source D is the missing counterpart given its similar morphology to Source C, and similarly steep spectral index (α ∼ −1.3). Adopting 150 Myr as the minimum age for both relic lobes, Source D would need to have been displaced at a minimum speed of ∼780 km s⁻¹ to reach the projected ∼120 kpc distance away from the southern BCG. This value is well below the measured velocity dispersion of ${998}_{-194}^{+120}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for 11 galaxies in the cluster, and also below the sound speed c_s ∼ 1100 km s⁻¹ at 4.7 keV, an upper limit on the motions of gas within the ICM. Merging clusters, which SDSS 1531 likely is, can generate large-scale (∼100 kpc), long-lived (∼1 Gyr) motions of the ICM, called sloshing (for a review, see Markevitch & Vikhlinin 2007). Numerical simulations of cluster mergers predict gas velocities up to ≥1000 km s⁻¹ (e.g., Roettiger et al. 1993, 1996, 1998; Ricker & Sarazin 2001). Therefore, it is plausible that turbulent gas motions in the ICM displaced the lobe.

If Source D is not the missing counterpart, and the lobe has not faded away in a less dense portion of the ICM, we also explore whether Source C could be a wide-angle-tail (WAT) radio source whose tails are unresolved by the LOFAR beam. WATs are powerful, bent radio sources thought to be produced via interaction of outflowing radio jets with an magnetized ICM (for a review, see O'Dea & Baum (2023). As a result, WATs are preferentially found in merging clusters, which we believe to be the case for SDSS 1531. Although WATs are rarely found in cool-core clusters because mergers strong enough to produce WATs are expected to disrupt the cool core (Ritchie & Thomas 2002; ZuHone et al. 2010), the few WATs observed in cool cores may not have had their cores disrupted yet due to a potential time delay between the merger and disruption. The main detractor from the WAT explanation lies in the morphology of Source C, which significantly deviates from typical WAT characteristics. WATs usually have tails of radio emission extending from the BCG, but Source C has an "avocado" shape, with its major axis perpendicular to the expected tail orientation.

If Source C is not a relic from a past AGN outburst, it could be a radio minihalo. The ≥3σ emission fully traces the X-ray surface brightness distribution within the central 60 kpc of the cool core. Such large-scale emission is a hallmark of radio minihalos, which are typically confined to the cores of relaxed clusters (∼50–500 kpc) and characterized by higher emissivity (e.g., Giacintucci et al. 2013, 2017). While ram pressure might have altered the minihalo's symmetrical profile to extend more to the southeast, this scenario is unlikely due to the striking concurrence between the X-ray cavity and the bulk of the radio emission. Furthermore, the image may artificially inflate the size of the source due to smearing from LOFAR's larger synthesized beam. In the event a minihalo is present, its emission is likely intertwined with that from the AGN relic, but the low-resolution observation makes it difficult to distinguish them clearly, as observed in A2029 (Govoni et al. 2009).

As reviewed in the introduction, recent multiwavelength observations have found evidence for multiphase gas draped around the rims of X-ray cavities in projection, supporting the idea that the gas condenses in the wake of the rising bubbles. Introduced by McNamara et al. (2016), the "stimulated feedback" model proposes that X-ray bubbles blown by AGN outbursts lift cold, low-entropy gas away from the location where the heating rate balances the cooling rate. If the overdensity does not return to its original position within its cooling time, the gas becomes thermally unstable, leading to the condensation of small clouds that eventually rain back onto the BCG.

Despite the current inactivity of the southern BCG's AGN, the southern end of the multiphase gas still retains a tangential connection to the cavity opening in projection, linking it to its potential origin from the previous epoch of AGN activity. During the bubble's inflation to its current projected location, it would have displaced (M_disp = μ(m_p + m_e )n_e V ∼ 10¹¹ M_⊙) of hot gas. The molecular gas has a lower mass of ∼10¹⁰ M_⊙, implying a hot-to-cold accumulated mass ratio of 0.1 and indicating that the bulk of the cold reservoir could have been accumulated over the most recent cycle of jet activity.

To determine whether it is energetically possible for the cavity to have displaced the gas and mitigate the rapid ICM cooling, we can calculate the work done by the buoyantly rising cavity. Bîrzan et al. (2020) examined the low-frequency radio emission in 19 nearby (z < 0.3) and six higher-redshift cool-core clusters (z > 0.3). When detected, they found that the low-frequency radio-emitting plasma rarely extended beyond the X-ray cavity edges, suggesting limited evidence for cosmic-ray electrons leaking. Thus, the ≥12σ radio contour of Source C likely provides a better approximation for the actual extent of the cavity compared to the opening uncovered by the shallow X-ray observations.

To estimate the bolometric X-ray cooling luminosity, L_cool, which describes the total luminosity within r_cool = 40 kpc, we use ${L}_{\mathrm{cool}}=4\pi {d}_{L}^{2}\cdot {f}_{{\rm{SB}}}$, where f_SB is the flux from the surface brightness profile and d_L is the luminosity distance. This gives an X-ray luminosity of 2.18 ± 0.14 × 10⁴³ erg s⁻¹. Using the hybrid X-ray-radio method put forward in Timmerman et al. (2022), which uses the volume measurement derived from radio observations and pressure measurement derived from X-ray observations, we can derive the cavity power P_cav as:

$$\begin{eqnarray*}&&{P}_{\mathrm{cav}}=\displaystyle \frac{{E}_{\mathrm{cav}}}{{t}_{\mathrm{cav}}}=\displaystyle \frac{4{pV}}{{t}_{\mathrm{cav}}},\end{eqnarray*}$$

where E_cav is the enthalpy of the rising bubble and t_cav is the age of the cavity (∼150 Myr).

For a surrounding pressure of ∼2.7 × 10⁻¹⁰ dyne cm⁻², we obtain a cavity enthalpy of ∼1 × 10⁶¹ erg, which is the same order of magnitude as the energetic outbursts observed in the few clusters with "supercavities," namely Hydra A (Nulsen et al. 2005b; Wise et al. 2007), MS 0735+7421 (Gitti et al. 2007), Hercules A (Nulsen et al. 2005a), and the more recent Ophiuchus Cluster (Giacintucci et al. 2020). The mean cavity power is 2.2 ± 0.5 × 10⁴⁵ erg s⁻¹, yielding P_cav/L_cool ∼ 86. This indicates that AGN feedback alone would be sufficient to extinguish the cooling flow and raises questions regarding this energetic outburst's impact on the cool core, which we discuss further in Section 4.1.6.

Since the radio/X-ray cavity is likely a relic of past, rather than ongoing AGN activity, the cavity is likely no longer uplifting the large reservoir of molecular gas, as one would expect in its youth. Instead, the uplifted gas is likely now raining back down onto the BCGs. The gas, however, exhibits smooth radial velocity gradients, which could suggest either an inflow toward or outflow away from the gravitational potential of the BCGs.

Establishing the direction of gas motion is challenging without directly observing gas clouds in absorption. Absorption lines have been detected against the submillimeter nuclear continuum emission in a few systems such as NGC 5044, A2597, and Hydra A (David et al. 2014; Tremblay et al. 2016; Rose et al. 2019, 2023). In NGC 5044 and A2597 (David et al. 2014; Tremblay et al. 2016), the apparent motion of the gas relative to the AGN indicates inflowing clouds that could serve as fuel for the central AGN. In contrast, the motion of the gas in Hydra A (Rose et al. 2019) suggests that it is on a stable, low-ellipticity orbit around the central galaxy. For SDSS 1531, however, ALMA detected no millimeter continuum emission, likely because the gas has yet to be accreted by the AGN.

If the projected separation between the young stars and molecular gas is caused by ram pressure stripping (discussed in-depth in Section 4.1.5), it can provide insight into the direction of gas motion. Section 4 of Li et al. (2018) investigated this using high-resolution hydrodynamical simulations (∼244 pc for the smallest cell) to examine the effect of ICM ram pressure on the cold clouds in the centers of cool-core clusters for different AGN feedback models. In their execution of the precipitation model, most cold gas condenses out of the ICM due to local thermal instabilities, causing the gas only to move inward. They find that when there is a detectable separation between the young stars and cold gas, the locations of the young stars are always closer to the cluster center than the cold gas. In SDSS 1531, the young stars lead the molecular gas relative to the merging galaxies, and the gas is entirely redshifted, consistent with the scenario where the gas is falling toward the center from between the cluster center and the observer. This observation, coupled with the location of the molecular gas by the cavity rim, supports the idea that we are observing the gas at a stage where it has decoupled from the AGN-blown bubble.

According to precipitation models, when the uplifted gas decouples from the ICM, it follows a drag-limited ballistic orbit, achieving speeds of up to ∼300–1000 km s⁻¹. However, observations of condensed clouds in cluster cores often show that clouds drift with subvirial velocities, in agreement with the CCA model. In the CCA model, the condensed gas displays bulk velocities up to a few 100 km s⁻¹. To test whether the bulk motions of the molecular gas follow a semiballistic trajectory, we follow Lim & Ao (2008) and assume that the gravitational potential can be modeled by a Hernquist profile (Hernquist 1990). Thus, a freely falling gas cloud should accelerate to a velocity v(r) of:

$$\begin{eqnarray}&&\displaystyle \frac{v{\left(r\right)}^{2}}{2}-\displaystyle \frac{{GM}}{\left(r+a\right)}=\displaystyle \frac{v{\left({r}_{0}\right)}^{2}}{2}-\displaystyle \frac{{GM}}{\left({r}_{0}+a\right)},\end{eqnarray} \tag{ 14 }$$

with respect to the ICM. Following Vantyghem et al. (2016), we modify the velocities in the above equation, such that v(r) is v(r) − v_ICM, where v_ICM is the velocity offset between the BCGs and the cooling gas, and v(r₀) is the initial velocity of the cloud, which is the same as that of the ICM: v(r₀) = v_ICM. The inclination angle of the cloud's trajectory and v_ICM are free parameters in this model.

The remaining parameters in Equation (14) are r₀, the initial radius from which the cloud originally formed, r, the current distance of the cloud from the center of the gravitational potential, a, the scale length related to the half-mass radius r_1/2 ($a={r}_{1/2}/(2+\sqrt{2}$), and M, the total gravitating mass of the BCGs. The BCGs' total mass interior to 30 kpc is ∼10¹³ M_⊙ (Sharon et al. 2014). The K-band luminosity from the SDSS observations provides a rough estimate of 3.8 × 10¹¹ M_⊙ for the stellar mass of the BCG (Maraston et al. 2009). Assuming the half-mass radius r_1/2 is approximately half of the ∼30 kpc Petrosian radius of the BCGs, the scale length is a ∼ 4.4 kpc. The amplitude of the velocity profile is primarily controlled by the total gravitating mass, which is degenerate with the inclination angle. Therefore, our results are not strongly affected by the adopted total mass value, since the inclination angle is uncertain and can be adjusted to compensate.

In Figure 23, we present the velocity trajectory (purple) of a freely falling clump of gas dropped from an initial height of 16 kpc on top of the ALMA PV diagram of the molecular gas. The trajectory assumes an inclination angle of 72° and a v_ICM of +600 km s⁻¹. This simple gravitational freefall model broadly reproduces the observed velocity gradient and velocities spanned by the molecular gas as it flows from south to north. Though the bulk motions of the cloud support gravitational infall, it is important to note that the molecular gas is not akin to a monolithic slab. Given the low volume filling factor of CO, the structure is more like a "mist" of smaller individual clumps and filaments seen in projection (e.g., Jaffe et al. 2001, 2005; Wilman et al. 2006; Emonts et al. 2013; McCourt et al. 2018; Tremblay et al. 2018). The freefall model fails to account for the fact that the individual clumps are likely not smoothly connected in velocity. Furthermore, the surface density of giant molecular clouds is such that the clouds reach terminal velocity and accelerate freely within the ICM (Combes 2018). It is possible that some of the individual clumps are in approximate freefall, while others are animated by the bulk velocities of the surrounding ICM velocity, as evidenced by the redshifted velocities characterizing the molecular filament.

Zoom In Zoom Out Reset image size

The observed "beads-on-a-string" star formation complex in SDSS 1531 is likely a product of the dynamic cluster environment. Although the YSCs are separated from the molecular gas by ∼1–3 kpc, the similarity in the morphologies of the beaded star formation and the edges of the cold molecular gas suggests that the gas played a critical role in fueling the star formation. Similar offsets between star formation and cold gas have been found in the Perseus Cluster (NGC 1275), A1795, and in simulations of AGN feedback (Canning et al. 2014; Li et al. 2015; Tamhane et al. 2023). Below, we discuss the potential contributions of a cooling wake, ram pressure, and tidal interactions to the observed separation between the star formation, the BCGs, and the molecular gas.

A cooling wake. The young stars and the cold molecular gas may be offset from the BCGs due to gas cooling in the wake of the central galaxies as they move through the cluster atmosphere. Fabian et al. (2001) proposed this scenario for A1795, a cool-core cluster where the BCG is offset by +150 km s⁻¹ from the cluster mean velocity, and by +374 km s⁻¹ within the central 270 kpc (Oegerle & Hill 1994). A1795 features a ∼50 kpc trail of multiphase gas that extends to the southeast of the BCG, with the gas motions reflecting that of the cluster instead of the BCG. In SDSS 1531, the central ellipticals are offset by ∼+100–400 km s⁻¹ from the average cluster velocity, and the motions of the multiphase gas are similarly offset from the BCGs. Without a direct measurement of the ICM velocity, we cannot confirm whether the BCGs are in motion with respect to the cluster. If they are, the molecular gas was likely deposited above the BCGs and progressively slowed by dynamical friction and/or ram pressure as it approached the BCGs' gravitational potential well (Combes 2018). This would explain why the gas is redshifted with respect to the southern nucleus and why the gas is not close enough to be accreted and fuel an AGN response (e.g., Tremblay et al. 2016, 2018).

Ram pressure. As discussed in Section 4.1.1, the presence of a major merger between the central galaxies, a bifurcated cluster redshift distribution, and the presence of diverse sources of diffuse radio emission within the cluster collectively suggest a turbulent ICM environment in SDSS 1531. A significant pointer toward the effects of the turbulent ICM motions is the morphology of Source B, reminiscent of a radio "jellyfish" galaxy. When such galaxies move through the dense ICM, they experience ram pressure forces strong enough to strip gas out of the disk directly and leave behind a wake of material trailing the galaxy. Complementing this, the relic AGN lobe extends southeastward from the southern BCG, mirroring the orientation of the molecular gas and the YSCs. Their parallel displacement suggests that ram pressure also plays a significant role in creating the observed projected offset between the YSCs and the cold molecular gas.

Since ram pressure primarily acts upon gas rather than stars, in order of critical density, it should slow down the ionized gas more efficiently than the molecular gas. In SDSS 1531, it is possible that the stars initially formed in the infalling molecular clouds and subsequently decoupled from the gas. As a result, the stars are now moving in the gravitational potential of the central galaxies without significant resistance, while the remaining molecular gas is still acted upon by ram pressure from the ICM.

To assess ram pressure's impact on the multiphase gas, we look to the few spatially resolved, high-resolution observations of molecular gas in jellyfish galaxies (e.g., Jáchym et al. 2019; Zabel et al. 2019, 2020; Cramer et al. 2020; Moretti et al. 2020). Each study shows that ram pressure from the ICM strips interstellar gas from infalling galaxies while the stars in the galaxy remain unaffected. Although the majority of these studies have found that the bulk of Hα-emitting gas is cospatial in projection with CO, Figure 3 in Moretti et al. (2020) shows an offset on the order of ∼1 kpc between Hα and CO, similar to what is observed in SDSS 1531.

Jellyfish galaxies experience great ram pressure forces as they fall through the ICM to the center of the cluster potential well. The central galaxies in SDSS 1531, however, presumably sit at the center of the cluster potential well, absent a more detailed redshift survey to prove otherwise (e.g., A1795; Oegerle & Hill 1994). Although the galaxies are not falling into the center at high >1000 km s⁻¹ speeds, the high gas density ρ_ICM in the center of the cluster increases the strength of the ram pressure forces.

To estimate the total ICM mass stripped M_strip within a stripping radius R_strip, we use the analytically determined relations from Gunn & Gott (1979) and Domainko et al. (2006), modified to account for the fact that the gravitational potential in BCGs is dominated by the total dark matter and gas mass and that elliptical galaxies do not have a stellar disk:

$$\begin{eqnarray}&&{M}_{\mathrm{strip}}=(x+1){e}^{-x}\cdot {M}_{\mathrm{gas}},\end{eqnarray} \tag{ 15 }$$

where M_total is the total mass of the BCGs, M_gas is the ICM mass, R₀ is the Petrosian radius of the BCGs, and v_gal is the relative velocity between the BCGs and the ICM. We adopt a value of 30 kpc for the radius, based on the Petrosian radius (Tremblay et al. 2014). This radius encloses the BCGs, the beaded string of star formation, and the molecular gas. From our X-ray observations, we obtain M_HE,total = 3.5 × 10¹² M_⊙ and M_gas = 5.9 × 10¹⁰ M_⊙ within ∼30 kpc, using the total mass and gas density profiles outlined in Section 2.1. We also set v_gal = 300 km s⁻¹, the relative velocity between the merging galaxies. Within a radius of 30 kpc, ρ_ICM ∼ M_gas/V_{30 kpc} ∼ 5.4 × 10⁻²⁶ g cm⁻³. Using the simplifying equations and assumptions above, we find that ram pressure can strip up to 5.5 × 10¹⁰ M_⊙ of gas within ∼15 kpc. This is comparable to the derived ∼10¹⁰ M_⊙ molecular gas mass.

To estimate the time needed for ram pressure to separate the YSCs from the molecular gas, we consider a simple model. In this model, we assume that the infalling molecular gas is subject to ram pressure and that the YSCs formed at t = 0. At t = 0, the YSCs and gas have an initial relative velocity of v₀ = 0 km s⁻¹. From this time onwards, we assume that only ram pressure causes the separation between the stars and gas, with no additional acceleration or deceleration due to other factors such as turbulence, and that both components experience roughly equivalent gravitational acceleration. We can estimate the gas' acceleration a_mol due to ram pressure as:

$$\begin{eqnarray}&&{a}_{\mathrm{mol}}={F}_{\mathrm{ram}}/{m}_{\mathrm{mol}},\end{eqnarray} \tag{ 16 }$$

where F_ram = P_ram A_mol, where A is the cross-sectional area of the molecular gas, and m_mol is the mass of the molecular gas. Modeling the projected extent of the molecular gas as an ellipse with major and minor axes of ∼11.1 and ∼6.8 kpc, respectively, we obtain a cross-sectional area of 235 kpc² and an acceleration of 1.8 × 10⁻¹⁴ km s⁻². Under the simplifying assumption that the molecular gas is only moving away from the young stars, we can estimate the time t_sep needed for the gas to travel a projected distance d = 1–3 kpc as:

$$\begin{eqnarray}&&{t}_{\mathrm{sep}}=\sqrt{\displaystyle \frac{2d}{{a}_{\mathrm{mol}}}}.\end{eqnarray} \tag{ 17 }$$

We obtain a timescale of ∼60–100 Myr for the YSCs to separate from the molecular gas, which is consistent with the rough <300 Myr age of the YSCs (Tremblay et al. 2014).

Tidal interactions. The beaded star formation may have been stimulated by tidal forces resulting from the major merger between the BCGs. When galaxies merge, tidal forces can pull and distort the stars and gas within them, moving stars from the disk to the spheroid component (Toomre 1977; Kaviraj et al. 2012). Tidal forces can also compress and shock gas into rapidly forming stars, resulting in a tidal-induced "starburst" (Wang et al. 2004). Gas-rich (wet) mergers typically host such starbursts due to the abundant fuel available for star formation (e.g., Lin et al. 2008; Perez et al. 2011; Athanassoula et al. 2016). In contrast, gas-poor (dry) mergers have less fuel, making starbursts less common in these systems (e.g., Bell et al. 2006; Naab et al. 2006; Lin et al. 2008).

Although gas-rich mergers between elliptical galaxies have been observed in some systems (e.g., Kaviraj et al. 2012; George 2017), the massive ellipticals in SDSS 1531 are more likely to be gas poor given the lack of significant molecular or ionized gas detected near their nuclei. However, tidal interactions can still impact the large reservoir of molecular gas that cooled from the ICM, providing the gas needed to stimulate star formation as expected in a wet merger. Since tidal forces scale with distance as r³, they are strongest in the region between the merging galaxies and decrease in strength outwards. These forces may have compressed the gas along the western border of the gas, stimulating a burst of star formation. Since tidal interactions impact both stars and gas, they could also explain why some of the strings of star formation appear tightly wound beneath the nuclei of the two central galaxies in projection. Tidal forces may have been too weak to compress the shielded gas further east, accounting for the absence of star formation in that region. Further investigations with hydrodynamical simulations would help further understand the role of tidal interactions in SDSS 1531.

Another beads-on-a-string system. To our knowledge, SDSS 1531 is one of only two known cool-core clusters that hosts a "beads-on-a-string" star formation complex near its center. The second, SPARCS 104922.6+564032.5 BCG (hereafter SPARCS 1049), resides among several merging cluster members and features diffuse emission in a tidal tail–like structure adorned with ∼66 kpc long beads-on-a-string star formation (Webb et al. 2015). The star formation is similarly offset from the central galaxy, though to a much larger extent of ∼25 kpc.

While initial investigations linked the distinct morphology of SPARCS 1049 to starbursts induced by gas-rich mergers, and a gas-rich BCG (Webb et al. 2015), more recent observations Chandra (Hlavacek-Larrondo et al. 2020) and the IRAM interferometer NOEMA (Castignani et al. 2020) have uncovered rapid ICM cooling and a significant ∼10¹¹ M_⊙ reservoir of cold molecular gas, offset from the center, likening its characteristics to SDSS 1531. SPARCS 1049's ∼860 M_⊙ yr⁻¹ SFR matches the rapid ICM cooling rate, making it one of the few known "true" classical cooling flow systems (Hlavacek-Larrondo et al. 2020). Observations also suggest that the AGN is currently inactive (Trudeau et al. 2019).

SDSS 1531, with a ∼1–10 M_⊙ yr⁻¹ SFR, has likely diverged from SPARCS 1049's current evolutionary path due to an earlier epoch of powerful AGN activity suppressing its cooling rate. Although the redshift distribution of 27 cluster members in SPARCS 1049 provides no evidence for significant subcluster structure (Webb et al. 2015), the large-scale star formation and merger activity within its core suggests a dynamic cluster environment similar to that of SDSS 1531. In both clusters, the ongoing mergers and residual cluster motions likely stimulated the "beads-on-a-string" star formation morphology. As a result, SDSS 1531 could represent a future evolutionary stage of SPARCS 1049, after AGN feedback severely limits the observed runaway cooling and star formation.

Another fossil AGN outburst detected in a cool-core system. Although AGN feedback is a common occurrence in cool-core clusters, the potential relic AGN lobes found in SDSS 1531 and Ophiuchus (Giacintucci et al. 2020) stand out due to the massive amount of energy (∼10⁶¹ erg) required to excavate each lone supercavity. Neither cavity is fully resolved by the X-ray data, and both clusters show limited signs of young AGN activity. Both are also expected to have undergone cluster-scale merger activity, with the merger in Ophiuchus strong enough to shift the peak of the cooling gas off from the BCG. In theory, such large-scale mergers and powerful AGN outbursts should heat and displace a large amount of the ICM in the cooling cluster cores, contributing significantly to their eventual disruption. In MACS J1931.82634, which also hosts a powerful AGN outburst and cluster-scale merger, Ehlert et al. (2011) argued that feedback from both processes has significantly disrupted the cooling core, evidenced by the cluster's flattening entropy profile, and the metallicity profile's consistently flat slope, which suggest that large masses of metal-rich gas were stripped from the center of the cluster and dispersed to the surrounding regions. In Ophiuchus, the metallicity and entropy profiles are still centrally peaked, though Werner et al. (2016) argued that these profiles were truncated to smaller radii by the activity within the core. In SDSS 1531, the entropy profile appears typical for a cool-core cluster. However, deeper X-ray observations are needed to determine the extent to which AGN feedback and the merger have potentially disrupted the cool core. Nevertheless, the description of AGN feedback as "gentle" (McNamara & Nulsen 2012) still holds for all three systems. Though feedback may have disrupted their cool cores, they do not appear to be transitioning to a non–cool core state, and if they are, the journey is far less chaotic than it could be.

Kinematics of the molecular gas. Molecular gas observed in cool-core BCGs typically exhibits extended, filamentary morphologies with complex, chaotic velocity structures (e.g., 13 out of 15 clusters studied in Olivares et al. 2019). Moreover, the majority of the molecular gas mass is typically concentrated in extended filaments; e.g., only two out of 12 BCGs studied by Russell et al. (2019) have ≤10% of their mass located in filaments. In contrast, SDSS 1531's molecular gas is mostly distributed in compact clumps. Additionally, the entire molecular gas structure exhibits a linear velocity gradient. RXC J2014.8-2430 similarly features a clumpy spatial distribution of molecular gas but does not exhibit a similar velocity structure.

Like SDSS 1531, the 2A 0335+096 galaxy cluster harbors two merging central galaxies with velocity offsets of ∼300 km s⁻¹ (Vantyghem et al. 2021). Although the spatial distribution of the molecular gas in 2A 0335+096 is more filamentary than clumpy, it exhibits a similar linear velocity gradient that can be replicated by a simple freefall model. However, the gas is situated in a filament between the central galaxies and spans a significantly smaller spatial extent than the gas in SDSS 1531.

Other galaxy clusters that have been observed to possess molecular gas distributions and velocity structures similar to SDSS 1531 include the MACS 1931 BCG (Fogarty et al. 2019), Hydra A (Olivares et al. 2019), and A262 (North et al. 2021). The molecular gas in Hydra A is thought to have a profile consistent with a rotating molecular disk, given the linear gas gradient, double-peaked spectra, and a central broader component of the velocity dispersion. The molecular gas in A262 is believed to either have an outflow or cooling filament of gas stimulated by CCA. In contrast, the molecular gas in SDSS 1531 does not have a velocity profile and structure consistent with a disk.

A commonly proposed alternative for the origin of multiphase gas in cluster cores is that the gas was accreted from interactions with gas-rich spiral or dwarf galaxies. This theory is consistent with the high luminosity, dusty, and small spatial extent of many molecular filaments observed in cool-core clusters. It could also provide a unified picture of AGN and radio galaxies, wherein the material responsible for activating the supermassive black hole is thought to be driven by a stochastic merger process (Sparks et al. 1989).

We cannot completely rule out this option for SDSS 1531, as BCGs are also predicted to grow through mergers (Ostriker & Hausman 1977), and there are at least two gas-rich spiral galaxies within the central 50 kpc. One of the aforementioned spiral galaxies has a redshift of z = 0.329 and shines brightly in the HST NUV filter, indicating the presence of young stars. We detect no molecular gas on the galaxy with the current ALMA observations, suggesting that it could have been stripped or that our observations are not deep enough to resolve the molecular gas content. The merger hypothesis is also viable because captured gas would be expected to rotate around the central galaxies, which we presume to be the case here.

The mass and velocity distributions of the observed molecular gas are the most significant challenges to the proposed scenario. Although the velocity distribution of the molecular gas shows signs of distinct velocity peaks, as expected in a major merger event (e.g., Gao et al. 2001; Greve et al. 2005; Schulz et al. 2007), the overall structure is still remarkably coherent, and cannot be easily accounted for by the contributions of multiple galaxies on different orbits. Furthermore, assuming the eight galaxies within the central 100 kpc of the cluster core that shine brightly in the HST NUV filter contain significant young stellar populations and thus high gas fractions, we can estimate an average of ∼10⁸ M_⊙ of H₂ (e.g., ESO 137-001; Jáchym et al. 2019) was stripped and flowed toward the BCGs. This results in ∼10⁹ M_⊙ gas reservoir, an order of magnitude below the estimated gas mass of ∼10¹⁰ M_⊙.

Lastly, we consider the possibility that the molecular gas is native to the merging central galaxies. In this scenario, the cold gas was expelled during the ongoing major merger, and the beaded star formation results from the gas dynamically responding to the gravitational torques and shear induced by the strong tidal field created by the merger. Of all the scenarios proposed, we find this to be the least likely. The molecular gas does not appear significantly disturbed morphologically, as one would expect if from the result of a major merger (e.g., SP423; McDonald et al. 2019). Moreover, major mergers involving elliptical galaxies have been found to trigger minimal or no star formation due to the small amounts of gas present (Cattaneo et al. 2008; Lin et al. 2008). Therefore, it is unlikely that the gas from the merger could account for the ∼10¹⁰ M_⊙ of H₂ observed, particularly considering minimal amounts of molecular and ionized gas were detected on the nuclei themselves.

However, the presence of YSCs between the merging BCGs offers some merit to this scenario, albeit on a smaller scale. If the merger enhanced gas compression and stimulated star formation (Wang et al. 2020), we would expect to find young stars between the merger participants, which is clearly observed. The bulk of the molecular gas lies east of the BCGs, so the little gas native to the BCGs could be responsible for the YSCs in between them. It would naturally follow that the YSCs and the ongoing merger entirely ionized the molecular gas in this region.