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

With two central galaxies engaged in a major merger and a remarkable chain of 19 young stellar superclusters wound around them in projection, the galaxy cluster SDSS J1531 + 3414 ( z = 0.335 ) offers an excellent laboratory to study the interplay between mergers, active galactic nucleus ( AGN ) feedback, and star formation. New Chandra X-ray imaging reveals rapidly cooling hot ( T ∼ 10 6 K ) intracluster gas, with two “ wings ” forming a concave density discontinuity near the edge of the cool core. LOFAR 144 MHz observations uncover diffuse radio emission strikingly aligned with the “ wings, ” suggesting that the “ wings ” are actually the opening to a giant X-ray supercavity. The steep radio emission is likely an ancient relic of one of the most energetic AGN outbursts observed, with 4 pV > 10 61 erg. To the north of the supercavity, GMOS detects warm ( T ∼ 10 4 K ) ionized gas that enshrouds the stellar superclusters but is redshifted up to + 800 km s − 1 with respect to the southern central galaxy. The Atacama Large Millimeter / submillimeter Array detects a similarly redshifted ∼ 10 10 M e reservoir of cold ( T ∼ 10 2 K ) molecular gas, but it is offset from the young stars by ∼ 1 – 3 kpc. We propose that the multiphase gas originated from low-entropy gas entrained by the X-ray supercavity, attribute the offset between the young stars and the molecular gas to turbulent intracluster gas motions, and suggest that tidal interactions stimulated the “ beads-on-a-string ” star formation morphology


Introduction
Near the centers of galaxy clusters lie the most luminous and massive elliptical galaxies in the Universe-brightest cluster galaxies (BCGs).These galaxies are broadly characterized as "red and dead" due to their smooth ellipsoidal shapes and significant fractions of old stars (Edwards et al. 2016;Kormendy 2016).However, under certain conditions, the hot (>10 6 K), diffuse plasma found between galaxies in clusters, known as the intracluster medium (ICM), cools rapidly, fueling new star formation and black hole activity near the BCG.
Early X-ray observations revealed that in approximately half of galaxy clusters, the intracluster gas harbored dense central regions with temperatures cooler than the cluster virial temperature and cooling times shorter than the age of the Universe (e.g., Fukazawa et al. 1994;Kaastra et al. 2004;Sanderson et al. 2006).These observations led to the development of the simple "cooling flow" model (Cowie & Binney 1977;Fabian & Nulsen 1977), which postulated that in the absence of a heat source to compensate for the rapidly cooling gas, several 100 M e yr −1 of plasma should cool to form large cold gas reservoirs of 5-50 × 10 11 M e near the cluster center (e.g., Fabian 1994).
Follow-up ultraviolet, optical, and infrared observations, however, found less than 1% of the predicted amount of cooled gas and highly suppressed star formation rates (SFRs) of ∼1-100 M e yr −1 near BCGs (e.g., McNamara et al. 2014;Johnstone et al. 1987;Romanishin 1987;McNamara & O'Connell 1989;Crawford & Fabian 1993;Allen 1995;Crawford et al. 1999;Mittaz et al. 2001;Rafferty et al. 2006;O'Dea et al. 2008;Donahue et al. 2015;Mittal et al. 2015;McDonald et al. 2018).Moreover, high-spectral-resolution imaging and spectroscopy from the 1999 launches of the Chandra and XMM-Newton X-ray space observatories found that the supposed "cooling flow" clusters exhibited minimal to no evidence for cooling down to temperatures of ∼0.1 keV and below (e.g., Peterson et al. 2003;Pinto et al. 2014).These discrepancies suggested that either the remaining cooling is hidden from view and/or that some steady form of heating compensates for the bulk of the radiative cooling (for reviews, see McNamara & Nulsen 2007;Fabian 2012;McNamara & Nulsen 2012;Fabian et al. 2022).
Various heating sources have been suggested for quenching cooling flows, such as thermal conduction, gas motions generated by mergers, and cosmic rays.However, feedback from active galactic nuclei (AGN) has since emerged as the most important source.The effects of "radio-mode" AGN feedback are nearly ubiquitous in cool cores, with 19 of the 20 brightest "cooling flow" clusters harboring distinct X-ray cavities (Fabian 2012).These cavities, also known as "bubbles," are often filled with radio emission from the lobes of the BCG's central AGN, suggesting that the lobes displaced the surrounding ICM to create the observed cavities (Bîrzan et al. 2004).As the bubbles buoyantly rise, the work done by the expanding radio lobes heats the surrounding X-ray plasma, offsetting radiative losses from the ICM (Churazov et al. 2002).
Although AGN mechanical feedback is routinely invoked as a mechanism for quenching star formation in simulations of galaxy formation (e.g., Somerville & Davé 2015), it does not completely offset ICM cooling, resulting in residual cooling at either low (<1%) rates (e.g., Tremblay et al. 2012) or in elevated episodes as the AGN varies in power (e.g., O'Dea et al. 2010;Tremblay 2011;McDonald et al. 2015).Observations have revealed massive flows of atomic and molecular gas that appear entrained around the rims of jetblown cavities (e.g., Russell et al. 2017;Tremblay et al. 2018), or closely trailing behind them (e.g., Russell et al. 2016;Vantyghem et al. 2016), suggesting incredibly efficient coupling between the radio jets, the ICM, and the cooled multiphase gas.Given that the cooling hot plasma, warm ionized gas, cold molecular gas, and radio emission from the AGN should each retain imprints of their shared journey within the ICM, multiwavelength observational studies of cool-core clusters have become standard practice for studying the interplay between nebular emission, star formation, and AGN activity, commonly referred to as the "AGN feedback cycle" (e.g., McNamara et al. 2014;O'Dea et al. 2008;McDonald et al. 2015;Russell et al. 2017;Tremblay et al. 2018;Pasini et al. 2019;Ciocan et al. 2021;Calzadilla et al. 2022;Masterson et al. 2023;Tamhane et al. 2023).
The AGN feedback cycle, however, is not without potential disruption.Under the hierarchical model of structure formation, both BCGs and galaxy clusters assemble via a sequence of major and minor mergers driven by gravity (e.g., Peebles & Yu 1970;Rosati et al. 2002;Voit & Donahue 2005;Kravtsov & Borgani 2012).Galaxy-galaxy interactions have long been predicted by theoretical (e.g., Martinet 1995) and numerical analyses (e.g., Toomre & Toomre 1972; Barnes & Hernquist 1992;Mihos & Hernquist 1996) to generate tidal torques that drive gas inflows.These gas inflows may then be accreted by the central supermassive black hole, fueling nuclear activity (Sanders et al. 1988).Multiple observational studies show significant AGN enhancement in merging galaxies (e.g., Alonso et al. 2007;Woods & Geller 2007;Weston et al. 2017).Contrasting findings exist, however, with some studies not observing a higher merger rate in AGNs (Grogin et al. 2005;Kocevski et al. 2011).On larger scales, cluster-cluster mergers can disperse and reheat cooling gas in the cluster core.These clusters are often home to diffuse radio emission, indicating the presence of relativistic particles (i.e., cosmic rays) and cluster-wide magnetic fields amplified by the shocks and turbulence injected into the ICM by the merger.Simulations predict that such turbulence has the potential to stimulate condensation, which then fuels AGN activity (e.g., Gaspari et al. 2017Gaspari et al. , 2018)).
Observational studies of cool-core clusters have traditionally biased "relaxed" systems due to their primary identification via X-ray surveys.Although the detection of mass-selected samples of clusters via the Sunyaev-Zel'dovich effect (Andrade-Santos et al. 2017) has facilitated the study of a larger sample of dynamically disturbed systems (e.g., Olivares et al. 2023), these systems remain comparatively understudied.To contribute to the evolving understanding of AGN feedback in disturbed systems, we present a multiwavelength study of a recently discovered cool-core cluster that features remarkable signatures of cooling-powered star formation amid a major merger between the two central galaxies: the SDSS J1531 +3414 (hereafter SDSS 1531) galaxy cluster.
SDSS 1531 is a strong-lensing cluster of galaxies at z = 0.335 (Hennawi et al. 2008;Oguri et al. 2009;Bayliss et al. 2011;Gralla et al. 2011;Postman et al. 2012).The cluster core was imaged as part of the Hubble Space Telescope (HST) strong-lensing imaging program (Sharon et al. 2020), revealing a remarkable ∼28 kpc scale network of young stellar superclusters (YSCs; or tidal dwarf galaxies) wound beneath and in between two merging giant elliptical galaxies of roughly equal stellar mass in projection (see Figure 1).The superclusters are reminiscent of the "beads-on-a-string" star formation frequently observed in the arms of spiral galaxies, resonance rings, and the tidal arms that bridge interacting galaxies (e.g., Elmegreen & Efremov 1996).Our initial work on this system (1) resolved 19 YSCs in the HST nearultraviolet (NUV) filter (Figure 1; rest-frame NUV inset plot, bottom right), (2) determined that the stellar superclusters and central elliptical galaxies comprising the BCGs are pinned to nearly the same redshift, indicating that the observed features are coplanar and not the result of a projection effect, and (3) estimated an extinction-corrected 5-10 M e yr −1 SFR (Tremblay et al. 2014).In a complementary work that analyzed the cluster's strong-lensing properties, Sharon et al. (2014) concluded that the overwhelming majority of the observed NUV emission is too bright to stem from counterimages of the lensed galaxies or a faint central image of a background source.
In this work, we present new data from the Atacama Large Millimeter and submillimeter Array (ALMA), the Chandra X-ray Observatory, the Gemini-North telescope's Gemini Multi-Object Spectrograph (GMOS), the Low Frequency Array (LOFAR), and the Very Large Array (VLA) to facilitate a multiwavelength view of SDSS 1531ʼs dynamic environment and uncover the origin and evolution of the "beads-on-a-string" star formation complex.This paper is organized as follows: Section 2 describes our procedure for reducing and analyzing the new observations from Chandra, LOFAR, VLA, GMOS, and ALMA.Section 3 presents spatial and spectral results for each gas phase.We use (1) Chandra X-ray observations to create thermodynamic profiles and spectral maps that reveal the general structure and properties of the hot ICM; (2) LOFAR and VLA radio surveys to identify radio emission associated with AGN activity or larger-scale cluster activity; (3) GMOS integral field unit (IFU) observations to create emission-line maps that display the spatial orientation of the warm ionized gas, unveil its kinematics and allow us to explore potential sources of ionization; and (4) ALMA observations to create maps revealing the morphology, kinematics, and mass distribution of the cold molecular gas.Section 4 synthesizes the results and proposes scenarios for the origin and fate of the star formation complex.
Throughout this study we assume H 0 = 70 km s -1 Mpc -1 , Ω M = 0.27, and Ω Λ = 0.73.In this cosmology, 1″ corresponds to 4.773 kpc at the redshift of the southern BCG (z = 0.335), where the associated angular size and luminosity distances are 984.4 and 1756.1 Mpc, respectively, and the age of the Universe is 9.728 Gyr.The spectral index, α, is defined such that the flux density at frequency ν is S ν ∝ ν α .

Observations and Data Reduction
In this section, we describe the new and archival observations from Chandra, VLA, LOFAR, GMOS, and ALMA.All new and archival observations of SDSS 1531 are summarized in Table 1.All Python codes/Jupyter Notebooks used and created to analyze the data are publicly available in an online repository. 18
To reduce and analyze the data from both ObsIDs, we used the Chandra Interactive Analysis of Observations (CIAO) v4.13 red, and yellow, respectively.The cluster features remarkable strong-lensing arcs, numerous elliptical and spiral galaxies, and the focus of this paper: merging elliptical BCGs.From left to right, the three inset panels show a closer view of the merging elliptical nuclei and "beads-on-a-string" star formation in the V band, the BCGs in all bands, and the 19 resolved YSCs in the rest-frame NUV.(Fruscione et al. 2006) and CALDB v4.9.3.We reprocessed the data with chandra_repro, cleaning the ACIS background in VFAINT mode.Flares were identified and filtered using the ChIPS routine LC_CLEAN.To remove the 3σ flares detected during the last 15 ks of ObsID 18689 completely (see Figure 2), we excluded this period from analysis.Point sources were identified through a wavelet decomposition technique (Vikhlinin et al. 1998) and visually inspected before masking.After cleaning, we retained a total of 65 ks (44.1 and 20.9 ks for ObsIDs 17218 and 18689, respectively), yielding a total ∼8200 net counts in a 60″ radius centered on the cluster's X-ray emission peak.The counts present are sufficient to constrain the surface brightness profile within the central region, classify SDSS 1531 as a cool-core cluster, search for tentative X-ray cavities within the cluster center, and perform cursory measurements of key spectral properties, such as the central projected temperature and cooling time (McDonald et al. 2019).
To analyze the morphology of the intracluster gas, we merged the ObsIDs with the MERGE_OBS CIAO script.To highlight surface brightness edges, we used the wavelet decomposition technique (Vikhlinin et al. 1998) to create a reconstructed 0.5-7 keV image.We also created an unsharpmasked image by smoothing the data with a 0 98 Gaussian and subtracting it from the same image smoothed with a 9 8 Gaussian (e.g., Hlavacek-Larrondo et al. 2012).
To derive the spectral properties of the intracluster gas, we first extracted spectra from a series of 20 concentric annuli centered on the X-ray emission peak from 1″ to 60″ in the energy range 0.5-7 keV (see Figure 3) for each ObsID.Next, we defined a 46 5 circle within the same chip but further from the cluster to encompass the appropriate local background, then subtracted it from the annuli.Each annulus contained ∼200-800 net counts across both ObsIDs.Next, we fit and modeled the total 0.5-7.0keV spectrum extracted from each annulus to a PHABS * APEC model in XSPEC.The hydrogen column density and abundance were fixed at N H = 1.79 × 10 20 cm −2 (estimated using NASA's HEASARC N H Column Density tool; Bekhti et al. 2016) and 0.3 Z e (e.g., Panagoulia et al. 2014), respectively, to constrain the fit better.The temperature kT and normalization parameter were allowed to vary.
The normalization factor N(r) of the APEC model provides an estimate for the 3D density profile through the relation:   where D A is the angular distance of the source, n e is the electron density, n p is the proton density, and V is the volume integral performed on the projected annulus along the line of sight (LOS).Assuming spherical symmetry, we can estimate the projected n e as: From the electron density and temperature of the ICM, we obtained the pressure P ≡ 1.83n e kT, the entropy º -K kTn e 2 3 , and the cooling time t cool of the ICM, defined as: e p e H cool where Λ is the cooling function (Sutherland & Dopita 1993), , n p = 0.92n e , n H = 0.83n e and we assume Z = Z e /3.
As a consistency check on our values from the above analysis and to obtain a careful estimate of the cluster mass, we follow a procedure similar to that outlined in Vikhlinin et al. (2006).First, we generated surface brightness profiles for each ObsID by extracting spectra from concentric annuli with finely spaced radii.The annuli were centered on the peak of the X-ray emission and extended from 1″ to 1000″.
The resulting surface brightness profile for each ObsID was corrected for spatial variations in temperature, metallicity, and effective area, and expressed as a projected emission measure integral.We then modeled each calibrated surface brightness profile with the following modified β-model (Cavaliere & Fusco-Femiano 1978;Vikhlinin et al. 2006 where α and β are fit indices, n 0 is the core density, r c is the scaling radius of the core, and r s is the scaling radius of the extended components.We fixed γ = 3 and ò < 5 to exclude unphysically sharp density breaks and set all other parameters free.Projecting this three-dimensional model along the LOS yielded an emission measure profile to fit to the data.An analytic expression for the three-dimensional gas density profile ρ g (r) was obtained by fitting the emission measure profile.The gas density was estimated assuming ρ g = m p n e A/Z, where A = 1.397 and Z = 1.199 represent the mean atomic mass and mean charge for a plasma with 0.3 Z e , respectively (McDonald et al. 2013).From the gas density profile, we computed the total gas mass M g within a spherical volume V(r).From this, we obtain an estimate for the classical cooling rate  M cool , where The temperature profile was fit using nine annuli of equal logarithmic width, r out /r in = 1.4,from the peak of the X-ray emission out to 60″.After subtracting the background, each annulus contained ∼600-1300 net counts across both ObsIDs.The temperature was modeled using a three-dimensional analytic model described by Vikhlinin et al. (2006): r r cool cool .The temperature profile was then obtained by projecting the model along the LOS to fit the observed temperature values to ∼170 kpc, where the temperature is well constrained.
To calculate the total cluster mass within radius r, we utilize a three-dimensional model for the gas density and temperature profiles and employ the hydrostatic equilibrium equation (e.g., Sarazin 1988): where T is in keV and r is in megaparsecs.Given M(r), we then compute the total matter density profile . Since M(r) is most reliably determined in the central region of the cluster where the temperature profile is well constrained, we obtain the total cluster mass M 500 from the Y X -M scaling relation, where Y X approximates the total thermal energy within R 500 , the cluster radius corresponding to a density contrast of δ = 500 (Vikhlinin et al. 2009).We list masses derived from the Y X -M scaling relation as M 500 −Y X and masses derived from the hydrostatic equilibrium assumption as M HE .
To generate high-resolution temperature, pressure, and entropy maps, we used the automated Python pipeline CLUSTERPYXT19 (Alden et al. 2019).We provided the pipeline with the previously described cleaned, merged and exposure- corrected observations and fixed the Galactic hydrogen density and cluster metal abundance to the same values described above.Point sources, identified by visual inspection, were carefully selected and fed to the pipeline to be excluded from the analysis.From here, CLUSTERPYXT utilized the CIAO adaptive circular binning (ACB) algorithm to conduct highresolution spectral fitting and subsequently generate temperature, pressure, and entropy maps.For a more detailed description of CLUSTERPYXTʼs spectral fitting process, we refer the reader to Alden et al. (2019).
We present and analyze the resulting X-ray deep image, spectral fits, profiles, and maps in Section 3.1.

Radio Surveys
To identify any diffuse radio emission associated with AGN or larger-scale cluster activity, as well as provide rough estimates of their spectral indices, we have thoroughly examined all publicly available images of the SDSS 1531 cluster field from the VLA and LOFAR radio facilities.
The VLA has observed SDSS 1531 over multiple epochs.The Expanded VLA conducted the most recent observation of the system on 2014 March 22 for 37 minutes in the 1.5 GHz band in the C configuration as part of project 14A-527 (PI: C. O'Dea).The source J1331+3030 (3C286) was used as a flux and phase calibrator for the observation.The data were successfully reduced using Common Astronomy Software Applications (CASA) version 5.1.0(McMullin et al. 2007), and the observation yielded a nondetection with a 3σ point source upper limit of 135 μJy.The VLA also observed the SDSS 1531 field in the 1.4 GHz band on 1994 June 18 in the B-array configuration (Project ID: AB628) and 1995 May 4 in the D-array configuration (Project ID: AC308).Both observations were calibrated and imaged after several additional rounds of phase-only and phase + amplitude self-calibration using standard procedures in the Astronomical Image Processing System version 31DEC23.The final B-array image achieved an rms sensitivity of ∼0.17 mJy beam −1 and has a restoring beam of size 2 88 × 2 47 at P.A. = −91°.1.The final D-array image achieved an rms sensitivity of ∼0.12 mJy beam −1 and has a restoring beam of size 37 7 × 26 2 at P.A. = 78°.3.
On 2018 September 14, the LOFAR high-band array observed the cluster field at 144 MHz for a total of 8 hr as part of the LOFAR Two-meter Sky Survey (LoTSS; Shimwell et al. 2022) under project code LC10_010.The second LoTSS data release (LoTSS-DR2) provided a public high-resolution 6″ mosaic of the pointing P233+35, which contains the full cluster field near the lower center.The LoTSS-DR2 data reduction and image mosaicing pipelines are described in detail in Shimwell et al. (2022) and references within.

GMOS-N Optical IFU Spectroscopy
In 2014, the GMOS-N instrument on the ground-based Gemini-North telescope observed the SDSS 1531 BCGs in IFU-R mode (Gemini program GN-2014A-DD-3, PI: Tremblay).An 1800 s exposure was taken with spatial dithering at a P.A. of 330°.The R400 grating was used and centered at 7000 Å, providing spectral coverage from 5800 Å < λ < 9000 Å.The resolving power, R, of the R400 instrument is R = λ/ Δλ = 1918 at a blaze wavelength of 764 nm, yielding a velocity resolution of ∼160 km s −1 (FWHM) at 764 nm.The original field of view (FOV) of the instrument is 3 5 × 5″; after dithering, the final image area covered 3 8 × 6 0, corresponding to 18.1 × 28.6 kpc at the target redshift (see Figure 5).
The data was reduced using the Py3D data reduction package for fiber-fed IFU spectrographs (Husemann et al. 2016).To obtain the most accurate emission-line measurements, we decoupled and modeled the stellar and gas components of the galaxy using the package PYPARADISE. 20 PYPARADISE is a Python version of the stellar population synthesis fitting code PARADISE (Walcher et al. 2015).The code models the stellar continuum by iteratively performing nonnegative linear least-squares fitting of the stellar spectrum of each spectral pixel ("spaxel") to a large library of stellar population templates and finds the best-fit LOS velocity distribution with a Markov Chain Monte Carlo method.The best-fit stellar continuum spectrum is subtracted from each spaxel to obtain a gas-only data cube.
Using the publicly available penalized pixel-fitting (pPXF) code (Cappellari 2017), we fit the emission-line model of the gas-only data cube produced by PYPARADISE.We fit the emission lines with single Gaussians to determine their flux, velocity, and velocity dispersion, and also obtained formal uncertainties from the covariance matrix of the fitted parameters.Although the observed gas emission lines were bright enough to fit at the native spatial resolution, the stellar continuum necessitated spatial binning using Voronoi tessellation (Cappellari & Copin 2003) due to the low signal-to-noise ratio (S/N).The GMOS cube was tessellated to achieve a minimum S/N of 20 (per bin) in the line-free stellar continuum.Figure 4 shows the modeled stellar and gas-phase components.
The products from PYPARADISE and PPXF enabled the creation of spatially resolved flux, velocity, and velocity dispersion maps for the following five emission lines: Hα, Hβ, [O III] λ5007, [N II] λ6548, and [N II] λ6583.To analyze the kinematics of the lines, we adopt a threshold of S/N 3. Section 3.3 discusses the maps created.

ALMA
ALMA observed the center of SDSS 1531 in Band 6 as part of ALMA Cycle 3 (ALMA Program ID: 2015.1.01426.S; PI: G. Tremblay) across two scheduling blocks between 2016 April 22-27 and September 8-13.One spectral window was centered at 259.02322 GHz (rest-frame 345.8 GHz at z = 0.335) to target the 12 CO (J = 3 − 2) molecular line transition, an excellent tracer of cold molecular hydrogen (H 2 ).Three 1875 MHz spectral windows, centered at rest frequencies of 347.1, 327.075, and 324.405GHz, were also used to detect line-free continuum.The total integration time on source was 105 minutes.
The BCGs were completely mapped within ALMA's ∼28″ (∼130 kpc) primary beam, but this array configuration is only sensitive to emission on scales up to ∼17″ (∼80 kpc; see Figure 5).Notably, no line or continuum emission was detected in the extended configuration observations when imaged alone.Although line emission was detected in the compact configuration, no corresponding continuum emission was found.
The raw ALMA visibilities were processed into calibrated measurement sets using the ALMA automated pipeline reduction script in CASA version 6.2.1.7.The calibrated visibilities were then imaged and deconvolved with the CLEAN algorithm.After testing multiple configurations of weightings (natural and Briggs), binning (10, 20, 40, and 80 km s −1 channels), and uv tapering (0 6, 0 8, 1 0, and 1 2), we determined that the natural-weighted 20 km s −1 cube without uv tapering provided the optimal setup, recovering the most emission and maximizing sensitivity.The final data cube presented in this paper achieves an rms sensitivity and angular resolution of 0.22 mJy beam −1 per 20 km s −1 channel, and a 0 53 × 0 33 (2.51 × 1.59 kpc) synthesized beam at P.A. = −1°.15.All CO(3-2) fluxes and line widths reported are corrected for the response of the primary beam (PBCOR = TRUE).
We present CO(3-2) integrated intensity, LOS velocity, and velocity dispersion maps in Section 3.4.The moment maps were created from the continuum-subtracted, calibrated spectral cube and constructed using the Python package BETTERMO- MENTS, 21 which applies a quadratic fit to the spectral data (Teague & Foreman-Mackey 2018).To recover as much flux as possible while suppressing noise, we spectrally smoothed the data with a top-hat kernel two channels wide, applied a Savitzky-Golay filter using a polynomial of order zero, spatially smoothed the data by 2.6 pixels, and applied a sigma clip to all pixels with S/N < 3.5σ.
Our observations, presented in Section 3.4, affirm the presence of a cold star-forming reservoir of gas located just beneath the central elliptical galaxies.

Adopting a Systemic Velocity
The SDSS 1531 BCGs host various moving components: the two central galaxies, as well as their stellar and nebular components.Tremblay et al. (2014) used archival SDSS spectroscopy centered on the southern BCG's nucleus to pin the stellar redshift to z = 0.3350 ± 0.0002, in agreement with prior redshift measurements for the BCG (Hennawi et al. 2008) and the cluster (Bayliss et al. 2011).Follow-up spectroscopy from the Nordic Optical Telescope indicated that the velocity of the northern elliptical is blueshifted by 280-300 km s −1 with respect to the southern BCG.Though this offset is negligible within the context of the Hubble flow (i.e., the galaxies are certainly interacting), it becomes significant when placed in context with the motions of the cold molecular and warm ionized gas.To interpret the motions of the stellar components of the BCGs with respect to the motions of the gas, we must select a systemic velocity to be used as a "zero-point," marking the transition from blue-to redshift.In this paper, we adopt a systemic velocity of cz = 100,430 km s −1 , where c is the speed of light and z = 0.3350 ± 0.0002 is the redshift of the southern nucleus.We select the redshift of the southern nucleus instead of that of the northern nucleus because its value is aligned with prior redshift measurements, detailed above.All velocity maps presented in this paper are projected at this value.
Instead of using the stellar redshift, other studies of BCGs in cool cores have set the zero-point relative to where the observed molecular CO emission peaks (e.g., Tremblay et al. 2016;Vantyghem et al. 2019).However, we do not adopt this approach because the bulk of the molecular gas detected lies to the east of the BCGs, rather than concentrated within their central few kiloparsecs (e.g., Vantyghem et al. 2021).Regardless of the reference point chosen, the data reveals a significant velocity offset between the YSCs and the gas.

Results
In this section, we present Chandra X-ray observations of the cooling hot (∼10 6 K) ICM, followed by LOFAR and VLA radio observations tracing emission from AGN and dynamic cluster-scale activity.Next, we present GMOS-IFU optical and ALMA millimeter observations of the warm (∼10 4 K) ionized and cold (∼10 2 K) molecular gas near the BCGs.For a detailed discussion of the results, please refer to Section 4. 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.

Unveiling the Cool Cluster Core
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.

Discontinuity near the Putative Cavity Edge
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.
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.

Spectral Maps and Profiles of the ICM
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) butterflyshaped 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.
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 0 , central entropy K 0 , central pressure P 0 , and central cooling time t cool,0 .Each property is consistent with that expected of a strong cool core.
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 a r =d d r 2 log log g at r = 0.04 R 500 , 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 T 0.25 min 0 . For comparison, strong cooling flow clusters like A478 and A133 have T T min 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.

Cooling and Freefall Timescales
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: 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 13 b M e 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 12 M e .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 13 M e .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).

Unclassified Extended Radio Emission
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.
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 avocadoesque 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 200 ).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.

Similarities between the X-Ray and Radio Morphologies in the
Cluster Core 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.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.

Preliminary Source Classification from Integrated Radio Spectra
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: where S 1 is the flux density of the source at ν 1 = 1.4 GHz, and S 2 is the flux density of the source at ν 2 = 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 −1 ), setting an upper limit of 0.052 mJy beam −1 .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.
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.

Emission-line Maps
We present the Hα flux, velocity, and velocity dispersion maps of the warm (<10 4 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.
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 −1 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.

Rate of Star Formation
From the GMOS data, we measure a total Hα flux of 1.9 ± 0.03 × 10 −16 erg s −1 cm −2 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 4 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 9 cm 3 (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 e yr −1 , which we interpret as a lower limit.Tremblay et al. (2014) estimated an extinction-corrected SFR in the range of 5-10 M e yr −1 .They measured an SFR of ∼5 ± 2 M e yr −1 from correcting the ALFOSC Hα luminosity for extinction using the Balmer decrement from SDSS DR10ʼs reported Hα and Hβ fluxes 22 and an SFR of 9.55 ± 3.4 M e yr −1 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 −15 erg cm −2 s −1 , which corresponds to an SFR of 3 M e yr −1 .We thus estimate the  2.78 ± 0.04 K 0 (keV cm 2 ) 18.1 ± 0.26 P 0 (10 −10 dyn cm −2 ) 5.36 ± 0.076 t cool,0 (Gyr) 0.5 ± 0.015 Note.Summary of the key X-ray properties of SDSS 1531.M 500 and R 500 are calculated assuming hydrostatic equilibrium and the Y X -M relation from Vikhlinin et al. (2009).Central quantities (kT 0 , P 0 , K 0 , and t cool, 0 ) are measured at a radius of 11.9 kpc and extracted from the linearly spaced annuli.
22 SDSS DR10 reported fluxes derived from the Princeton specBS pipeline, accessed at http://das.sdss.org/SpecBS(Adelman- McCarthy et al. 2008).An unphysical (Hα/Hβ < 2.3) Balmer decrement is obtained when using SDSS DR18 fluxes from the Portsmouth catalog, which is based on stellar kinematics as evaluated by PPXF (Thomas et al. 2013).This result agrees with the unphysical value we obtained from the GMOS data cube.Discrepancies between emission-line fluxes reported by SpecBS and the Portsmouth catalog were similarly noted in Lyu & Liu (2016), although less drastic.Without higher-resolution optical data to confirm which result is correct, we elect to adopt the physical value from SDSS DR10.
total SFR within the beaded star formation complex to be between 1 and 10 M e yr −1 .

Ionization Sources
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.
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 2 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.

The Cold Molecular Gas
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 −1 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 ∼0 6 (∼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 −1 channel matches the morphology of the southern string of beads.The beads and the molecular gas are separated by ∼0 3 (∼1.4kpc), 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.
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.

Velocity Structure: Redshifted Gas Flows
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 −1 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 −1 .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 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.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 0 5-1 5).
The velocity of the gas in the northern filament, central clump, and southern filament ranges mostly from +60-260 km s −1 , +280-460 km s −1 , and +480-680 km s −1 , 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.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 2 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 20 cm ( ) --K km s 1 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 13 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.), consistent with the velocity calibration used for the ALMA maps in Section 3.4.2.The gas is redshifted up to +800 km s −1 with respect to the southern BCG and is most disturbed near the southern nucleus.We directly compare these kinematic maps with the ALMA data in Section 3.4.5.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 −1 channels (the spectrum is shown in the left panel of Figure 18).The fit, from −690 to 900 km s −1 , yielded an emission integral of ∼1.7 ± 0.3 Jy km s −1 .Noting the aforementioned caveats, these results give an H 2 gas mass of (3.7 ± 0.7) × 10 10 M e .With natural weighting, we obtain integrals ranging from 0.88 to 2.5 Jy km s −1 depending on the different binning selected, yielding mass estimates that range from 1.9 to 6.5 × 10 10 M e . 23he 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).

Undetected Continuum
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 −4 mJy beam −1 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: 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 2 kg 1 , 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 8 M e (and thus a molecular gas-to-dust ratio of 100).

Comparison with the Warm Ionized Gas
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. 2016Tremblay et al. , 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.
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 , intensity-weighted velocity (moment 1), and intensity-weighted velocity dispersion (moment 2) maps for SDSS J1531.The molecular gas extends ∼24 kpc from north to south.There is a notable offset between the YSCs and the molecular gas distribution, which is discussed further in Section 3.4.1.Like the Hα gas, the molecular gas is similarly redshifted up to +600 km s −1 with respect to the southern BCG.The filled dark red circle represents the size of the ALMA beam.
ionized and molecular gas exhibit similar overall velocity structure within ±160 km s −1 , which is within the GMOS data's ∼160 km s −1 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.

Gas Depletion Timescale
To estimate the amount of time needed to deplete the entire ∼10 10 M e reservoir of molecular gas, we adopt an SFR of 1-10 M e yr −1 .We calculate the depletion time as ).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.

Discussion
The new X-ray, optical, and radio data presented in this paper reveal that the merging central ellipticals and the associated beaded strings of star formation in SDSS 1531 are situated within a rapidly cooling, highly magnetic ICM.The two most striking results of this study are (1) the remarkable spatial alignment between the diffuse radio emission from Source C and the concave X-ray surface brightness discontinuity, and (2) the ∼10 × 10 10 M e molecular nebula's mostly offset spatial location from the YSCs.We propose that the alignment between the X-ray discontinuity and the radio source indicates the presence of a large X-ray cavity, likely blown during an older epoch of AGN feedback.We posit that the molecular gas originated from low-entropy gas entrained by the massive cavity and is slowly encircling and settling into the gravitational potential well of the BCGs.A combination of ram  pressure, tidal interactions, and ionization by star formation likely contribute to the projected separation between the young stars and the molecular gas.In the following section, we provide a detailed description of this scenario based on the results presented in the previous section, which are summarized in Figure 20.We also evaluate alternative mechanisms that could have contributed to the observed molecular gas supply, such as its potential capture from previous encounters with gasrich galaxies and/or an origin from the central ellipticals.A summary of our expectations for each scenario and whether or not the current observations support them is summarized in Section 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 clusterscale 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.
If the merging ellipticals in SDSS 1531 result from a clusterscale 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.(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 2 ), 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).

Thermally Unstable Cooling
The dynamic environment within SDSS 1531 offers several viable routes to condensation from the hot atmosphere.Defining the "classical" cooling rate as: cool M e yr −1 .Assuming a consistent cooling rate over 3 Gyr, the hot atmosphere could easily supply up to ∼10 12 M e of cool material, with only 100 Myr needed to fill the entire 10 10 M e 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, Figure 18.3σ ALMA CO(3-2) spectra extracted from different components of the molecular gas.From left to right, the panels display multi-Gaussian fits to the CO(3-2) spectra extracted from the region encompassing the entire extent of the molecular gas, the northern filament, the central clump, and the southern filament.The black curves represent the raw spectra extracted from each region, while the purple curves depict the fits to each component of the molecular gas within the velocity extent defined for each filament/clump.The residuals are shown in light pink.The multi-Gaussian fit of the total structure's spectrum estimates a cold gas mass of approximately 10 10 M e , with the majority concentrated in the northern filament and a roughly even distribution between the central clump and southern filament.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: 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 −1 .Accounting for the conversion from the LOS to a three-dimensional velocity dispersion with a factor of 3 , we adopt σ v,L = 80-300 km s −1 .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.
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 6 M e 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 turbulenceinduced 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.

An Older Epoch of AGN Activity?
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: 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 39.6 min 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σ 2 /R.SDSS reports two values for the stellar velocity dispersion of the BCG: (1) σ = 444 ± 58 km s −1 from the SDSS pipeline's SPECOBJ table and (2) σ = 356 ± 34 km s −1 from the Portsmouth catalog.We adopt the value reported from the SDSS pipeline as the stellar velocity dispersion and use the ∼100 km s −1 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 −1 ).
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): 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 2 , k is the ratio of energy density in nonradiating particles to that in synchrotron-emitting particles, V is the source volume, C 12 is a constant depending on the spectral index and frequency cutoffs, f 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.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 5 kpc 3 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): 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 −6 (1 + z) 2 G.With = B min m 2.9 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 −13 dyne cm −2 , which is three orders of magnitude lower than the external gas pressure derived from the X-ray observations (∼2 × 10 −10 dyne cm −2 ).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 f 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 −1 to reach the projected ∼120 kpc distance away from the southern BCG.This value is well below the measured velocity dispersion of -+ -

km s
194 120 1 for 11 galaxies in the cluster, and also below the sound speed c s ∼ 1100 km s −1 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 −1 (e.g., Roettiger et al. 1993Roettiger et al. , 1996Roettiger et al. , 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 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.
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. 2013Giacintucci et al. , 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).

Cooling Stimulated and Mitigated by AGN Feedback
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 11 M e ) of hot gas.The molecular gas has a lower mass of ∼10 10 M e , 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 cosmicray 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 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 43 erg s −1 .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: where E cav is the enthalpy of the rising bubble and t cav is the age of the cavity (∼150 Myr).Note.A summary of the expectations for each scenario explored in Section 4 to explain the origin of the cold molecular gas and whether the current observations support these expectations (Y-yes), find them unlikely (Nno), or are inconclusive (Inc.).
For a surrounding pressure of ∼2.7 × 10 −10 dyne cm −2 , we obtain a cavity enthalpy of ∼1 × 10 61 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 45 erg s −1 , 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. 2019Rose et al. , 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, lowellipticity 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 −1 .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 −1 .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: 2 0 2 0 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 0 ) is the initial velocity of the cloud, which is the same as that of the ICM: v(r 0 ) = 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 0 , 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 ), and M, the total gravitating mass of the BCGs.The BCGs' total mass interior to 30 kpc is ∼10 13 M e (Sharon et al. 2014).The K-band luminosity from the SDSS observations provides a rough estimate of 3.8 × 10 11 M e 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 −1 .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. 2001Jaffe et al. , 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.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 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 −1 from the cluster mean velocity, and by +374 km s −1 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 −1 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. 2016Tremblay et al. , 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. 2019Zabel et al. , 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 −1 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: where M total is the total mass of the BCGs, M gas is the ICM mass, R 0 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 12 M e and M gas = 5.9 × 10 10 M e within ∼30 kpc, using the total mass and gas density profiles outlined in Section 2.1.We also set v gal = 300 km s −1 , the relative velocity between the merging galaxies.Within a radius of 30 kpc, ρ ICM ∼ M gas /V 30 kpc ∼ 5.4 × 10 −26 g cm −3 .Using the simplifying equations and assumptions above, we find that ram pressure can strip up to 5.5 × 10 10 M e of gas within ∼15 kpc.This is comparable to the derived ∼10 10 M e 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 = 0 km s −1 .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: 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 2 and an acceleration of 1.8 × 10 −14 km s −2 .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: 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 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 3 , 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.

SDSS 1531 in Context
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.5BCG (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 11 M e reservoir of cold molecular gas, offset from the center, likening its characteristics to SDSS 1531.SPARCS 1049ʼs ∼860 M e yr −1 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 e yr −1 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 61 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 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-2430similarly 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 −1 (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 8 M e of H 2 (e.g., ESO 137-001; Jáchym et al. 2019) was stripped and flowed toward the BCGs.This results in ∼10 9 M e gas reservoir, an order of magnitude below the estimated gas mass of ∼10 10 M e .

Molecular Gas Native to the BCGs?
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 10 M e of H 2 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.

Conclusions, Summary, and Future Work
Initially observed by the Wisconsin-Indiana-Yale NOAO ground-based telescope (Hennawi et al. 2008), SDSS 1531 originally appeared to contain one large and disturbed central BCG.However, subsequent observations from HST, SDSS, ALFOSC, and more revealed two elliptical galaxies engaged in a major merger, and a remarkable chain of stellar superclusters-a kiloparsec-scale manifestation of the Jeans instability (Tremblay et al. 2014).Less than a decade later, with new observational data from Chandra, LOFAR, GMOS, and ALMA, this paper presents four main results: 1. SDSS 1531 is a cool-core cluster.SDSS 1531 satisfies all the criteria for a cool-core cluster, including a cuspy emission measure profile (α = 1.4) with a clear central overdensity, t cool < t H ) within the inner 100 kpc, and t ff /t cool ∼ 20 within the central 10 kpc.The gas within r cool = 100 kpc cools at a rate of   M M 185 yr −1 , which corresponds to a cumulative gas mass of 10 12 M e within 3 Gyr if the cooling were uninterrupted.2. There is compelling evidence for an old, extremely powerful AGN outburst.The Chandra X-ray observations reveal a concave surface brightness discontinuity near the edge of the cool core.LOFAR low-frequency observations fill this discontinuity, in alignment with the picture of a radio AGN lobe pushing aside the X-ray gas to create an X-ray cavity.The source's steep spectrum suggests that the emission stems from aged plasma from an AGN outburst that occurred during an earlier epoch in the cluster's history.The energy required to excavate the cavity is ∼10 61 erg, making it one of the most powerful AGN outbursts observed.If the missing symmetrical lobe has not faded into a less dense portion of the ICM, we propose that the nearby radio Source D, which bears a similar morphology and spectral index to Source C, could be the missing lobe.3. Comoving cold and warm gas are tangentially connected to the cavity opening.In projection, a bright cloud-like structure of cold molecular and warm ionized gas lies to the north of the X-ray cavity opening, suggesting that the origin of the gas is tied to the older AGN outburst.The multiphase gas likely cooled in the wake of the buoyantly rising cavity and is now infalling back to the BCGs.The warm gas envelops the YSCs and is redshifted up to ∼800 km s −1 .The massive ∼10 10 M e reservoir of cold gas is mostly comoving with the warm gas, with the central regions differing in velocity by up to ∼160 km s −1 .The cold gas lies mostly to the southeast of the YSCs, offset by ∼1-3 kpc. 4. The "beads-on-a-string" star formation complex is likely a product of the dynamic cluster environment.The beaded star formation and the edges of the cold molecular gas suggest that the gas played a critical role in powering the observed star formation.A cooling wake from the central galaxies moving through the ICM and/or strong ram pressure forces may have caused the observed offset between the YSCs and the molecular gas.Tidal interactions resulting from the major merger between the central ellipticals may have stimulated the beaded star formation via gas compression and contributed to the observed morphology of the stars and gas.
Potential origin scenarios for star-forming gas are summarized in Table 5.Further constraining its origin will require follow-up observations across multiple wavelengths in addition to numerical simulations.A forthcoming spectroscopic survey of the cluster with the Multiple Mirrors Telescope will be crucial for identifying whether a subcluster merger is mitigating ICM cooling in SDSS 1531, confirming our estimate of galaxies with high gas fractions within the cluster core that could have been stripped to power the observed cold gas reservoir, and constraining the velocity offset between the central galaxies and the average cluster member.
To confirm the presence of an extremely powerful, old AGN outburst in the cluster and investigate its link to the multiphase gas, we would need deeper, more sensitive radio observations that sample several lower and higher frequencies.Additional molecular line observations could probe the extent of ram pressure's impact on the spatial offset between the molecular and ionized gas, as we would expect to see spatial offsets between different CO line maps in order of critical density.The molecular line data could also allow the creation of line ratio maps to explore potential excitation mechanisms within the molecular filament, which we predict to be most excited near the western border.Deeper Chandra observations will be critical for confirming the presence or absence of X-ray cavities, while future X-ray microcalorimetric observations could establish if large-scale motions of the surrounding ICM dominate the motion of the cold and warm gas.This scenario could also be investigated using numerical simulations of cluster mergers that include the effects of radiative cooling and star-forming gas to determine if merger-induced gas motions can produce the observed separations between the cooled gas and newly formed stars.

Figure 1 .
Figure 1.Hubble's Wide Field Camera 3 (WFC3) view of SDSS 1531, the focus of this paper.NUV, V-band, H-band, and I-band emission are shown in blue, cyan, red, and yellow, respectively.The cluster features remarkable strong-lensing arcs, numerous elliptical and spiral galaxies, and the focus of this paper: merging elliptical BCGs.From left to right, the three inset panels show a closer view of the merging elliptical nuclei and "beads-on-a-string" star formation in the V band, the BCGs in all bands, and the 19 resolved YSCs in the rest-frame NUV.

Note.
Summary of all new and archival observations presented in this paper.The observations are presented in descending order of wavelength, from X-ray to radio.Columns: (1) wave band; (2) facility name; (3) instrument (and aperture/detector/band) used for observation; (4) imaging filter or spectroscopic configuration; (5) wave band, filter/grism central wavelength, or emission lines covered by observation; (6) on-source exposure time; (7) date of observation; and (8) comment specific to observation.

Figure 2 .
Figure 2. Count rate vs. time in Chandra ObsID 18689.3σ flares are persistently detected from ∼45-60 ks, requiring the complete removal of the last 15 ks, one-quarter of the total exposure time.

Figure 3 .
Figure 3. Chandra X-ray image of SDSS 1531 in the 0.7-2 keV band.We primarily analyze the spectra extracted from 20 linearly spaced annuli within 60″ (red solid circle) of the cluster's central region, marked by the peak of the X-ray emission (black cross).The derived value for R 500 is marked by the white dashed circle.

Figure 4 .
Figure 4.An illustration of the best-fit stellar and gas-phase components obtained using PYPARADISE for one of the brightest spaxels (x = 8, y = 13) in the GMOS-N IFU data cube of the BCGs.The top panel shows the original spectrum with both stellar and gas emission (dark blue) and the best-fit stellar continuum by PYPARADISE (red).The error spectrum is shown in green.The spectrum of the stellar continuum-subtracted gas-only component (dark blue) is depicted in the bottom panel, while the fit performed by the PPXF routine is shown in green.

Figure 5 .
Figure 5. Top: FOVs covered by ALMA (blue) and GMOS (orange), overlaid on an HST F606w image of SDSS 1531.The Chandra, VLA, and LOFAR observations (FOVs not shown) cover the entire cluster field.Bottom: GMOS Vband view of the BCGs, with HST contours (black) overlaid.The photo-centroids of the optical nuclei of the galaxies in the GMOS cube match those of the HST data.

Figure 6 .
Figure 6.Left: Pan-STARRs i-band image (blue) of SDSS 1531 and its surrounding ∼1720 × 1720 kpc environment.Overlaid in red is a wavelet-fit broadband Chandra X-ray surface brightness map in the 0.7-3.0keV range, revealing the extended emission from the intracluster gas.The white box outlines the central 50 kpc region shown to the right.Right: a closer 100 × 100 kpc view of the SDSS 1531 cluster in the HST F606W band (blue), with the same X-ray surface brightness map overlaid.The peak of the X-ray gas (black cross) coincides with the location of the southern BCG.Two 25 kpc (right) and 30 kpc (left) "wings" below the peak mark a potential cavity opening.

Figure 7 .
Figure 7. Unsharp-masked 0.5-7.0keV image of the central region of SDSS 1531 with white contours of the surface brightness map from Figure6overlaid in white contours.The opening to a potential massive cavity is located between the X-ray "wings."

Figure 8 .
Figure 8. Unfiltered 0.5-7 keV Chandra X-ray image of SDSS 1531 shown in Figure 7 with the regions used to extract the surface brightness profile (left) and spectra (right) overlaid.The dashed green line marks the location of the tentative surface brightness edge.The center plot shows the resulting surface brightness profile and depicts a clear discontinuity in surface brightness at the site of the potential shock/cold front.Due to a limited number of counts, we can only conduct spectral fitting within the two large pie regions shown.

Figure 9 .
Figure9.Wide-field temperature (left) and pressure (right) map of SDSS 1531 generated using the ACB technique described in Section 2.1.Both maps are overlaid with contours of the X-ray surface brightness map shown in Figure6.Uncertainties range from 7% to 10%.The temperature map reveals an extended cooling trail.The distribution of the pressure and entropy (not shown) are roughly circularly symmetric, with the highest pressure and lowest entropy near the X-ray peak.

Figure 10 .
Figure10.Top left: projected X-ray emissivity per unit area for both Chandra ObsIDs (black and green points), fit with a modified β-model (thick pink line).Bottom left: projected X-ray temperature profile using linearly spaced annuli (black points) vs. log-spaced annuli (orange points), which are used to fit the temperature model described inVikhlinin et al. (2006; V+06 hereafter).Middle column: pressure (top) and entropy (bottom) profiles.Right column: cooling time (top) and t cool /t ff (bottom) profiles.Brown points represent the freefall time estimate using the strong-lensing mass profile fromSharon et al. (2014).The upturn of the V+06 t cool /t ff profile at low radii is a numerical artifact resulting from preventing the derivative from approaching zero.In the middle and right columns, gray points represent the thermodynamic profiles of nearby clusters from the ACCEPT sample(Cavagnolo et al. 2009).All profiles consistently indicate the presence of a strong cool core in SDSS 1531.

Figure 11 .
Figure 11.Top: LOFAR 144 MHz image of SDSS 1531 with 4″ × 6″ beam (orange ellipse) and noise σ rms = 0.15 mJy beam −1 .The image unveils six prominent sources of diffuse, low-frequency radio emission within the cluster, each labeled with a pink letter.There are six prominent sources of radio emission associated with the cluster, each labeled with a pink letter.A white dashed circle outlines the central 200 kpc of the cluster.Bottom: a closer view of each identified source.The redshift of each potentially associated galaxy is labeled.Source C is most closely aligned with the BCGs.In both the top and bottom panels, the contour levels start at 3σ rms and increase sequentially with steps of 6σ, 12σ, 20σ, and 37σ.

Figure 12 .
Figure12.LOFAR contours of Source C overlaid on the Chandra surface brightness map shown in Figure6.The 12σ contour (cyan) aligns remarkably well with the X-ray "wings," suggesting that the wings represent an opening to a giant X-ray cavity, filled by the relic AGN lobe.

Figure 13 .
Figure 13.Left to right: GMOS maps of Hα flux, LOS velocity, and velocity dispersion, extracted after subtracting the stellar continuum as described in Section 2.3.White (left) and gray (center, right) contours depict the HST F606W view of the BCGs.Blue circles indicate the spatial location of the 19 YSCs.The Hα-emitting gas fully engulfs the YSCs and may extend beyond the eastern edge of the GMOS FOV.The velocities shown are projected around a zero-point at z = 0.335 (e.g., cz = 100,430 km s −1), consistent with the velocity calibration used for the ALMA maps in Section 3.4.2.The gas is redshifted up to +800 km s −1 with respect to the southern BCG and is most disturbed near the southern nucleus.We directly compare these kinematic maps with the ALMA data in Section 3.4.5.

Figure 14 .
Figure 14.Emission-line diagnostic diagrams for spaxels with S/N 3 in each emission line.Left: BPT (Baldwin et al. 1981) diagnostic plot.The spaxels are color coded based on their location relative to boundaries between well-known empirical and theoretical classification schemes (seeKewley et al. 2001;Kauffmann et al. 2003;Schawinski et al. 2006) shown in gray dashed and solid lines.The majority of the spaxels lie within the star-forming region.Right: sky distribution of the spaxels, color coded according to their BPT diagram classification.

Figure 15 .
Figure15.Left to right: ALMA integrated line intensity flux density (moment 0), intensity-weighted velocity (moment 1), and intensity-weighted velocity dispersion (moment 2) maps for SDSS J1531.The molecular gas extends ∼24 kpc from north to south.There is a notable offset between the YSCs and the molecular gas distribution, which is discussed further in Section 3.4.1.Like the Hα gas, the molecular gas is similarly redshifted up to +600 km s −1 with respect to the southern BCG.The filled dark red circle represents the size of the ALMA beam.

Figure 16 .
Figure 16.ALMA CO (3-2) channel maps displaying the three distinct segments of molecular gas in velocity space, and how the peak emission is offset by 1-3 kpc from the YSCs.The left and right panels show emission morphologically similar to the YSCs at 220 and 660 km s −1 , respectively.The central panel at 400 km s −1 shows emission comprising the central clump, mostly spatially distinct from the YSCs.The blue contours highlight 3σ emission.

Figure 17 .
Figure 17.PV diagram created from the CO(3-2) data.The top image displays the moment 1 velocity map with the PV extraction apertures overlaid in gray.The lower image shows the PV diagram extracted from the gray rectangular region.The molecular gas exhibits a smooth, continuous velocity distribution.

Figure 19 .
Figure19.Left to right: maps of the Hα and CO(3-2) flux, velocity difference, and velocity dispersion ratio, respectively, corrected for differences in spatial resolution (see Section 2.3).Whether the Hα emission extends beyond the molecular gas is uncertain due to the limited field covered by GMOS (black rectangular box).The edges of the velocity difference and dispersion ratio maps should be disregarded as they are artifacts of the subtraction/division. Though the ionized and molecular gas are not fully cospatial, they are largely comoving.

Figure 20 .
Figure20.Large-scale and small-scale (inset plot) views of the multiphase gas in SDSS 1531.Purple denotes the Chandra X-ray surface brightness map tracing the hot intracluster gas, orange represents the warm ionized gas traced by Gemini/GMOS-N IFU Hα flux, and blue shows the cold molecular gas traced by the ALMA CO (3-2) line.LOFAR contours of Sources C and D are outlined in white with 12σ contours, which best outline the extent of each potential relic lobe, highlighted in cyan.With the caveat that projection effects complicate interpretation, the 12σ radio contours perfectly fill the putative X-ray cavity opening, and the warm and cold gas phases are oriented perpendicular to the north of the cavity opening.The warm ionized gas fully envelops the YSCs, which are mostly offset from the molecular gas by ∼1-3 kpc.

Figure 21 .
Figure21.The redshift distribution of galaxies in the SDSS 1531 field, based on spectroscopic and photometric data from 14 galaxies.The distribution shows hints of two distinct peaks, suggesting a possible bimodal nature of the galaxy population, which may indicate a subcluster merger.

Figure 22 .
Figure 22.Profiles of the cooling time to eddy turnover time ratio (t cool /t eddy ) in the ICM.The dashed and solid lines represent injection length scales L = 30 and 70 kpc, respectively, corresponding to the radial extent of the CO nebula and the radio cavity.The dark blue and black lines correspond to velocity dispersions σ = 90 and 260 km s −1 , respectively.

Figure 23 .
Figure 23.A freefall model depicting the motion of a parcel of molecular gas released from a height of ∼18 kpc along the gravitational potential of the SDSS 1531 BCG.The velocity profile of the parcel is represented by a dashed purple line, which is overlaid on the contours of the PV diagram presented in Figure 17.The model successfully describes the observed PV distribution of the molecular gas.

Table 1
Summary of Observations

Table 3
Integrated Spectral Indices

Table 4
Radio Properties Assuming Equipartition and/or Pressure Equilibrium