1 Introduction↩︎

The present-day Universe exhibits a complex large-scale structure known as the cosmic web, consisting of voids, filaments, and nodes [1], [2]. It is well known that the properties of galaxies depend on the environments in which they reside (e.g., [3]–[14]). For example, the well-established morphology-density relation shows that early-type galaxies characterized by elliptical morphology, quenched star formation, and old stellar populations tend to reside in denser environments while late-type galaxies with disk-like morphology, active star formation, and young stellar populations are likely to be in less dense environments (e.g., [15]–[18]). This galaxy segregation suggests that galaxy evolution is significantly influenced by its surrounding environment.

Galaxy clusters are the most massive (\(M_\mathrm{h}\gtrsim 10^{14}\,M_{\odot}\)) virialized systems in the Universe (e.g., [19], [20]). They host hundreds to thousands of galaxies, as well as hot X-ray-emitting gas, known as the intracluster medium (ICM). Galaxy clusters are also the densest structures in the Universe, located at the nodes of cosmic filaments (e.g., [21]–[23]), and are dominated by quiescent galaxies (QGs), providing a unique laboratory for studying the role of dense environments in quenching galaxy star formation.

In galaxy clusters, various quenching mechanisms unique to massive host halos are proposed. These include direct stripping of cold gas from galactic disks through ram-pressure (ram-pressure stripping, RPS: e.g., [24]), halting gas replenishment due to removal of the gas reservoir (strangulation: e.g., [25]), and multiple close encounters with other satellite galaxies (harassment: e.g., [26], [27]). In addition, a fraction of galaxies accreted into clusters may have already been quenched in larger-scale filamentary structures before infall through environmental processes operating there (preprocessing: [28]–[33]). Understanding the origin of the quenched population in clusters, therefore, requires disentangling the relative contributions of these different quenching mechanisms.

It is known that the distribution of satellite galaxies in clusters is anisotropic: satellites tend to be aligned to the major axis of the central galaxy (i.e., the brightest cluster galaxy: e.g., [34]–[37]). Interestingly, this alignment is more significant for quenched satellites than star-forming ones [34], [38]–[40]. Recently, this representative alignment of QGs, dubbed anisotropic quenching or angular conformity, has attracted attention in the field of galaxy quenching pathways, and several possible interpretations of its physical origin have been proposed.

[41] have detected anisotropic quenching in a low-redshift (\(z < 0.08\)) galaxy group/cluster and compared their results with a cosmological simulation. They have argued that the anisotropy could originate from directional active galactic nucleus (AGN) feedback from the central galaxy, which weakens RPS along the feedback direction (i.e., minor axis of the central galaxy). Such directional AGN feedback is supported by observations [42] and simulations [43]. However, [44] have presented a counterargument to this interpretation by reproducing similar anisotropy in a simulation that does not include AGN feedback. [45] has also detected anisotropic quenching for massive clusters at \(0.2<z<0.5\) and proposed another scenario in terms of cluster morphology: a cluster halo is elongated along the major axis of the central galaxy, and some environmental effects are more efficient along that direction at fixed cluster-centric radii. More recent work has proposed an alternative explanation: [46] and [47] suggest that anisotropic quenching may instead arise from the preferential infall of preprocessed galaxies from larger-scale filaments. In this picture, the observed anisotropy reflects the anisotropic accretion of quenched populations rather than in-situ quenching within cluster halos. While the physical origin of anisotropic quenching remains debated, distinguishing between the AGN-feedback and preprocessing scenarios would provide important insights into the role of AGN feedback in dense environments and the impact of filaments on galaxy evolution prior to cluster infall.

Another important question about anisotropic quenching is its redshift dependence. While many works in the literature have mainly focused on the relatively low-redshift Universe (e.g., \(0<z<0.5\)), [48] have explored higher redshift using \(>5000\) clusters derived from the Hyper-Suprime Cam Subaru Strategic Program (HSC-SSP: [49], [50]) and detected anisotropic quenching at \(0.25<z<1\). This finding implies that anisotropic quenching occurs universally throughout a wide range of cosmic history and that the underlying mechanism operates even at relatively early times. If the physical origin of anisotropy at \(z\gtrsim1\) can be determined, it would place important constraints on quenching pathways in clusters at the relatively early stage of evolution.

In this study, we extend the exploration of anisotropic quenching beyond \(z=1\), using the spectroscopically confirmed high-redshift cluster catalog from the Gemini Observations of Galaxies in Rich Early Environments (GOGREEN) survey and the Gemini CLuster Astrophysics Spectroscopic Survey (GCLASS). We report the tentative detection of anisotropic quenching even at \(z>1\). Then we use a cosmological simulation to discuss the physical origins of the observed anisotropy and their implications for the role of preprocessing.

The structure of this paper is as follows. In section 2, we describe the observational data and the selection of the sample. In section 3, we show the detected anisotropic quenching signal beyond \(z=1\). In section 4, we conduct a semi-analytic model analysis based on a cosmological simulation and discuss the implications of the observed anisotropy. Overall discussion is in section 5, and section 6 is devoted to a summary and conclusion. Throughout this paper, we assume a flat \(\mathrm{\Lambda}\)CDM cosmology with \((\Omega_\mathrm{m},\, \Omega_\mathrm{\Lambda},\, h)=(0.3,\, 0.7,\, 0.7)\) and a [51] initial mass function. We use AB magnitudes [52] and physical distances.

2 Data and sample↩︎

2.1 Cluster sample↩︎

We use a high-redshift cluster catalog from the first public data release of the GOGREEN survey and the GCLASS survey [53]–[55], which includes 17 spectroscopically confirmed clusters at \(0.8<z<1.5\). Of these, 14 clusters are selected from the Spitzer Adaptation of the Red Cluster Sequence (SpARCS) Survey ([56]–[58], which are identified as overdensities of red-sequence galaxies [59]. Others are discovered by the South Pole Telescope (SPT) survey through the detection of the Sunyaev–Zeldovich (SZ) signature. [60]–[62]. These clusters are observed in multi-wavelength imaging from the near ultraviolet (NUV) to the near infrared (NIR), as well as in optical spectroscopy. In addition, high-resolution NIR images F160W taken by the Hubble Space Telescope (HST) are available for 12 clusters. These HST images are vital to avoid blending issues when measuring the major axes of central galaxies in dense cluster cores. Therefore, we only use these 12 clusters with HST images in our analysis, which are listed in table ¿tbl:tab:cluster?.

The cluster mass (\(M_{200}\)) and radius (\(R_{200}\)) have been reported in [63], which are obtained from the scaling relation between these quantities and the line-of-sight velocity dispersion of the spectroscopically confirmed cluster members. We summarize these values in table ¿tbl:tab:cluster?. For SpARCS1613, whose information is not available in [63], we refer instead to the values presented in [64].

2.2 Galaxy sample↩︎

We use a photometric redshift (photo-z) catalog provided by the GOGREEN/GCLASS data release. We mainly use photo-z, stellar mass, and rest-frame color information. Photo-z is estimated using the eazy code [67] based on multi-wavelength (1113 bands) photometry from NUV to NIR [55], [66]. The other galaxy properties are measured using the fast code [68] with the stellar population synthesis model of [69]. See [66] for the details.

[65], [66] have presented the limiting \(K_{s}\)-band magnitudes for each cluster at which 80 percent of the injected mock sources to the observed images can be detected. They have also calculated the corresponding limiting stellar mass by assuming relatively old stellar populations (i.e., QGs). The derived mass completeness in \(\log(M_{*}/M_{\odot})\) for the 12 clusters ranges from \(9.35\) to \(10.07\), with a median of \(9.67\). With a limited number of cluster samples, we decide not to apply a single stellar mass threshold to our galaxy sample (e.g., \(\log(M_{*}/M_{\odot})>10\)). Instead, we adopt the individual mass limits for each cluster, enabling us to maximize the sample galaxy size. The limiting magnitudes and stellar masses are summarized in table ¿tbl:tab:cluster?.

To classify galaxies as QGs and star-forming galaxies (SFGs), we adopt a rest-frame U-V and V-J color criterion [70], [71]: galaxies are QGs if

\[\begin{align} (U-V) > 1.3\;\cap \;(V-J)<1.5) \;\cap \\ (U-V) > 0.88(V-J)+0.59, \end{align}\] and others are SFGs.

2.3 Cluster membership↩︎

We define cluster membership by considering line-of-sight and projected distances as: \[\begin{align} |z_\mathrm{phot}-z_\mathrm{cl}|&<&w_{z} (1+z_\mathrm{cl}), \\ r_\mathrm{proj} &<& 1.5 R_{200}, \end{align}\] where \(z_\mathrm{cl}\) and \(w_{z}\) are the cluster redshift and the redshift window size for membership selection, respectively, and \(r_\mathrm{proj}\) is the projected distance from the cluster center at \(z_\mathrm{cl}\). We use photo-z rather than spec-z to avoid any biases due to the incompleteness of spectroscopic observation.

In general, the redshift window size should be determined by a compromise between contamination rate and incompleteness in cluster member selection: narrower (wider) windows result in lower (higher) contamination and higher (lower) incompleteness. Using the GOGREEN galaxy subsamples with spec-z’s, [66] have compared the number counts of ‘false negatives’ and ‘false positives’ (i.e., galaxies selected as cluster members by spec-z but not selected by photo-z, and vice versa), finding these two are almost the same at \(w_{z}\sim0.08\). In other words, contamination and incompleteness do not significantly affect the net galaxy number counts under this window size. Given these considerations, we adopt \(w_{z}\sim 0.05\), a slightly stricter selection criterion than [66], to keep our cluster memberships purer. We summarize the number of selected total and quiescent members in table ¿tbl:tab:cluster?.

2.4 Measuring position angle↩︎

We first select the galaxies nearest to the cluster center as central galaxies, and others are satellites [55]. With the HST F160W images, we measure the position angle (PA) of each central galaxy by fitting a single Sersic profile [72] using galight, a Python-based 2D model image fitting package [73]. galight can not only fit the primary object (i.e., the central galaxy) but also simultaneously fit all detected sources around it (i.e., satellite galaxies) using multiple Sersic models, thereby reducing blending effects in the dense cluster core. The Markov Chain Monte Carlo (MCMC) algorithm determines the best-fit model parameters, and PAs are very accurately determined: the uncertainties are less than \(2°\).

For each satellite, we then calculate the orientation angle, \(\theta\), the projected azimuthal angle relative to the PA of its central galaxy. We note that, in principle, all quantities related to this orientation angle should be symmetric in every \(90°\).

3 Result↩︎

We followed the method used in the literature to quantify the anisotropic quenching signal (e.g., [41], [45], [46], [48]). Satellite galaxies are divided into arbitrarily chosen angular bins, \(\theta_{i}\), and the quiescent fraction is calculated for each bin as:

\[f_{\mathrm{q},i}=\frac{N_{\mathrm{q},i}}{N_{\mathrm{all},i}}, \label{eqn:f95q}\tag{1}\] where \(N_{\mathrm{all,}i}\) and \(N_{\mathrm{q,}i}\) are the numbers of all galaxies and QGs in the i-th angular bin, respectively. The angular dependence is then fitted with a sinusoidal function, and the sinusoid amplitude is used to quantify the degree of anisotropy. We use the cosine function as a fiducial model:

\[f_\mathrm{q ,model}(\theta)=A_\mathrm{q}\cos(2\theta)+f_\mathrm{q,0}, \label{mtgxfyrz}\tag{2}\] where \(A_\mathrm{q}\) and \(f_\mathrm{q,0}\) are the modulation amplitude and the baseline, respectively. Given the small cluster sample size (\(N=12\)), we adopt six angular bins spanning \(15°\). The black points in figure 1 are the derived quiescent fraction with error bars calculated assuming the binomial distribution. There is a decreasing trend of \(f_\mathrm{q}\) from the major axis (\(\theta=0°\)) to the minor axis (\(\theta=90°\)).

The parameter fitting for equation 1 is performed by the MCMC procedure using emcee code [74]. We adopt uniform prior distributions. The log-likelihood function is: \[\mathcal{L}=-\frac{1}{2} \sum_{i}\left(\frac{f_{{\mathrm{q}},i}-f_\mathrm{q,model}(\theta_{i})}{\sigma_{\mathrm{q},i}} \right)^{2}, \label{eqn:likelihood951}\tag{3}\] where \(\sigma_{\mathrm{q},i}\) is the error of the quiescent fraction in the i-th angular bin.

In figure 1, we show the median \(f_\mathrm{q}\) values (black dashed line) in the MCMC sampling with the 68% (dark yellow) and 95% (light yellow) intervals. The cosine function fits the data well, and we derive the median value and 68% credible interval of each parameter from the posterior distribution (figure 2) as \(A_\mathrm{q}=0.045^{+0.021}_{+0.021}\) and \(f_\mathrm{q,0}=0.639^{+0.015}_{-0.015}\). We summarize these values in table ¿tbl:tab:fitting?. The solution of \(A_\mathrm{q}=0\) (i.e., no angular dependence) is outside the 98.4% credible interval, suggesting anisotropic quenching is detected at \(2.14\sigma\).

Although this calculation follows a standard approach, it is still possible that using an arbitrary angular bin may affect the anisotropy amplitude. As a cross-check, we also fit a model without binning. by adopting the log-likelihood calculated assuming a binomial distribution as:

\[\mathcal{L}=\sum_{j}\left[ q_{j}\cdot \log(f_\mathrm{q,model}(\theta_{j})) + (1-q_{j})\cdot \log(1-f_\mathrm{q,model}(\theta_{j}))\right], \label{eqn:likelihood952}\tag{4}\] where \(q_{j}\) takes \(1\) if the j-th satellite galaxy with orientation angle \(\theta_{j}\) is QG, and \(0\) otherwise. The median value and 68% credible interval for each parameter are \(A_\mathrm{q}=0.041^{+0.021}_{-0.021}\) and \(f_\mathrm{q,0}=0.639^{+0.015}_{-0.015}\), respectively. The detection significance is \(1.95\sigma\).

Regardless of the calculation method, we derive the non-zero amplitudes of anisotropic quenching. As a complementary check, we fit the data points with a constant value (i.e., an angular-independent) model, obtaining \(f_\mathrm{q,0}=0.643^{+0.018}_{-0.019}\). To examine whether the cosine function model is better to fit the data, we calculate the Akaike Information Criterion (AIC: [75]) for each model: \[\mathrm{AIC_{model}} = -2\ln \mathcal{L}_\mathrm{model} +2k,\] where \(\mathcal{L}_\mathrm{model}\) is a likelihood under a given model, and \(k\) is the number of parameter of the model. We then derive their difference of \(\Delta \mathrm{AIC}=\mathrm{AIC_{fiducial}}-\mathrm{AIC_{const}}=-2.58\), suggesting the fiducial cosine model is mildly favored. Based on these pieces of evidence, we conclude that the signature of anisotropic quenching in the clusters at \(z>1\) is confirmed for the first time, although still tentative (\(2\sigma\) level).

4 Origin of anisotropic quenching↩︎

There are two major possible origins of anisotropic quenching. One is an intrahalo quenching scenario presented by [41]. In this scenario, satellites are less likely to be quenched along the minor axis: AGN feedback from a central galaxy preferentially blows out the ICM along this direction, making quenching driven by RPS less effective. The other is the preprocessing effect [46], [47]. Clusters are frequently connected to large-scale filamentary structures where galaxies are more frequently quenched than in the field, and the major axis of the central galaxy and such filaments are likely to be aligned [40], causing the anisotropic distribution of QGs in a cluster. Indeed, the anisotropic quenching signature has been reported beyond the cluster virial radius [46]. Since anisotropic quenching beyond \(z=1\) is confirmed in this study, investigating its origin may have significant implications for understanding the role of AGN feedback, RPS by the ICM, and preprocessing in the emergence of the quenched populations in high-redshift clusters.

To understand the physical origin of anisotropic quenching, cosmological simulations are useful because they provide the accretion histories of cluster-member galaxies. Here we use the IllustrisTNG (TNG), a state-of-the-art cosmological simulation [76] to examine where and when galaxies are quenched around clusters. We first describe the simulation data in section 4.1. We then show how satellites’ orientation angles are changed after accretion in section 4.2. Finally, in section 4.3, we perform semi-analytic modeling of the observed anisotropic quenching signal to infer the origin of this phenomenon at \(z\sim1\).

4.1 Simulation data↩︎

TNG is a series of cosmological magnetohydrodynamical simulation suites with various baryon physics implemented: metal-line radiative cooling, stochastic star formation, stellar population evolution with associated chemical enrichment and mass loss, stellar feedback, seeding and growth of supermassive black hole and its feedback [77], [78]. There exist simulation runs with different box sizes and resolutions. We focus on the TNG300-1 (TNG300) run, which has the largest box size with a side length of \(300\,\mathrm{comoving\;Mpc}\) to include as many massive halos as possible in the analysis.

In TNG, a friends-of-friends algorithm (FoF: [79] is adopted to identify groups (halos), and the Subfind [79] is run to detect substructures in each FoF group as subhalos (galaxies). Subfind defines central galaxies as those containing the most gravitationally bound particle of the parent FoF group. In addition, Sublink [80] is used to link subhalos in different snapshots and generate merger trees (i.e., histories of halo assembly).

From the group catalog, we first select massive halos of \(\log(M_\mathrm{200}/M_{\odot})>13.5\) at \(z=1.25\) (snapshot #44) and \(\log(M_\mathrm{200}/M_{\odot})>14\) at \(z=0\) (snapshot #99). To ensure a sufficient number of groups for statistical analysis, we adopt the halo mass threshold for the \(z=1.25\) snapshot slightly lower than the typical cluster mass regime (e.g., \(\log(M_\mathrm{200}/M_{\odot})>14\)). Among subhalos belonging to each massive group, we restrict our analysis to those with stellar masses of \(\log(M_\mathrm{*}/M_{\odot})>10\) at given redshifts (\(z=0\) or \(z=1.25\)). We have 307 (280) groups and 8104 (10500) galaxies at \(z=1.25\) (\(z=0\)) snapshot.

To derive the 3-D positions of progenitors of group members at a given snapshot, we trace back merger trees of each galaxy along the main branch. However, merger trees only record the 3-D position at snapshots taken at typical intervals of \(150\,\mathrm{Myr}\). Since velocity dispersion in a massive cluster is \(>1000\,\mathrm{km\,s^{-1}}\), we may miss sudden changes in the orientation angle of satellites, especially near the cluster centers. Therefore, we perform a cubic Hermite spline (cf. [81]) to interpolate the object’s motion from their 3-D positions and velocities. In the end, the motion of each satellite is traced back over \(5\,\mathrm{Gyr}\) with a time step of \(10\, \mathrm{Myr}\).

Since observation uses projected position, the orientation angles of satellite subhalos are measured in the x-y, y-z, and z-x planes, tripling the effective sample size. The PA of a central subhalo on the i-j plane (\(i,j=\{x,y,z\}\)) is computed from the inertia tensor \(I_{ij}\) using stellar particles associated with that galaxy: \[\begin{align} I_{ij}&=&\sum_{n} m_{n} x_{n}^{i}x_{n}^{j}, \\ \mathrm{PA}_{ij}&=&\frac{1}{2}\arctan\left(\frac{2I_{ij}}{I_{ii}-I_{jj}}\right), \end{align}\] where \(m_{n}\) is the mass of the n-th stellar particle, and \(x_{n}^{i(j)}\) is its position on the coordinate i (j) relative to the subhalo center. Not to include stars extended far beyond the subhalo main body (i.e., intracluster light), stellar particles beyond the minimum of \(2r_{*,1/2}\) and \(50\,\mathrm{kpc}\) are excluded, where \(r_{*,1/2}\) is a half stellar mass radius of a given subhalo.

Finally, for each satellite, we calculate the time since it first reached \(R_{200}\) of its host group (hereafter, time since infall, or TSI). This quantity is commonly used as a measure of how long satellites have been subjected to the extreme environmental effects within a cluster halo (e.g., [82]). In Figure 3, we show the normalized distribution of TSI at \(z=0\) and \(z=1.25\) snapshots. At \(z=0\) snapshot, the satellites have a flat TSI distribution, and one third of them exhibit \(\mathrm{TSI}>5\,\mathrm{Gyr}\), which are not shown in the plot. In contrast, those at \(z=1.25\) snapshot are skewed toward shorter TSI (\(\mathrm{TSI}\lesssim2\,\mathrm{Gyr}\)), suggesting that our higher-redshift groups in TNG300 are dominated by the recent infalling satellites.

4.2 Orientation angle↩︎

Given large velocity dispersions in massive clusters, the orientation angles of satellite galaxies may change quickly. Therefore, the observed orientation angles of satellites can differ significantly from those at which satellites are under the strong quenching effects. [48] have assumed this timescale to be \(1-2\,\mathrm{Gyr}\), comparable to a cluster crossing time, and constrain quenching timescale to be \(\lesssim 1\,\mathrm{Gyr}\) if intrahalo quenching is the main driver of the anisotropy. However, the actual timescale over which the orientation angle changes significantly from its initial value has not been investigated.

To quantify this timescale, we trace the satellites’ orientation angles as a function of TSI. Regardless of snapshots, we always measure orientation angle from the major axes of the central galaxies at \(z=1.25\) (\(z=0\)). In figure 4, we show the distributions of orientation angles relative to the initial orientation against TSI, \(\Delta \theta= \theta(\mathrm{TSI})-\theta(\mathrm{TSI}=0)\), at \(z=0\) (top) and \(z=1.25\) (bottom). At \(z=0\), a clear peak is seen around \(0.6\, \mathrm{Gyr}\), probably representing the peri-center passage. After that, the median orientation angle approaches its initial value once (\(\mathrm{TSI}\sim1\, \mathrm{Gyr}\)) and is gradually randomized in \(5\,\mathrm{Gyr}\). This suggests that, if preprocessed galaxies are preferentially accreted onto clusters along the major axis of the central galaxy, their initial orientations are on average preserved at least several \(\mathrm{Gyr}\) given the symmetry of the orientation angles in every \(90°\), making the anisotropic quenching signature. This fact also suggests that the main contributors to the anisotropic quenching signal are recently infalling satellites. In the TNG \(z=0\) snapshot, [47] have calculated the amplitude of anisotropic quenching for the young (\(\mathrm{TSI}<3\,\mathrm{Gyr}\)) and old (\(\mathrm{TSI}>3\,\mathrm{Gyr}\)) population separately. They have shown that the young population exhibits anisotropy with an amplitude of \(0.035\pm0.019\), while the old population shows no anisotropy (\(0.002\pm0.009\)). This is consistent with our expectations.

At \(z=1.25\), the median \(\Delta \theta\) shows a peak around \(0.7\, \mathrm{Gyr}\), and that peak is much less clear than the \(z=0\) case. In addition, the orientation angles are randomized more immediately (\(\sim2.0\, \mathrm{Gyr}\)). This difference may be due to a difference in cluster crossing time. Since the crossing timescale is approximately proportional to \(H(z)^{-1}\), the crossing time at \(z\sim1\) is roughly half that at \(z=0\). A shorter timescale for randomizing the initial orientation angle compared to \(z=0\) may indicate that recently accreted galaxies play a more significant role in arising anisotropic quenching signals. We revisit this point in later sections.

4.3 Semi-analytic quenching modeling↩︎

To further understand the physical origin of the observed anisotropic quenching, it is necessary to consider not only when and where each galaxy is quenched, but also that the orientation angles of satellites can change over time, as discussed in the previous section. One possible approach is to investigate the quenching mechanisms implemented in hydrodynamical simulations such as TNG300. While this approach has the advantage of fully accounting for complex hydrodynamical processes, it can be challenging to isolate the physical mechanism driving the anisotropic quenching. More importantly, it has been reported that even the latest simulations, including TNG300, fail to reproduce the quiescent fraction observed in the GOGREEN clusters [83]. Therefore, instead of examining the quenching mechanisms implemented in TNG300, we adopt a semi-analytic approach that combines a simple anisotropic quenching model with the satellites’ infall histories obtained from TNG300. Although a semi-analytic approach cannot fully capture the complex hydrodynamical processes, it allows for a more flexible description of the data and facilitates physical interpretation.

We incorporate anisotropic quenching effects into the so-called ‘delay-then-rapid’ quenching model. Our model has four parameters (\(f_\mathrm{q,pre}\), \(\Delta f_\mathrm{q,pre}\), \(\tau_\mathrm{major}\) and \(\tau_\mathrm{minor}\)) and determines whether a given satellite is quenched depending on its initial orientation angle \(\theta_{0}\equiv\theta(\mathrm{TSI = 0})\) and TSI at \(z=1.25\). First, we randomly classify galaxies as QGs with probabilities of \(f_\mathrm{q,pre}+\Delta f_\mathrm{q,pre}\) and \(f_\mathrm{q,pre}\) for those accreted along the major (\(0°<\theta_{0}<45°\)) and minor axis (\(45°<\theta_{0}<90°\)), respectively; all remaining galaxies are classified as SFGs. This represents quenching that occurs before the cluster infall (i.e., preprocessing), and the anisotropy is implemented as \(\Delta f_\mathrm{q,pre}>0\). Remaining star-forming satellites also get ‘quenched’ if they spend time longer than delay times after accretion: \(\mathrm{TSI}>\tau_\mathrm{major}\) and \(\mathrm{TSI}>\tau_\mathrm{minor}\) for those accreted along the major and minor axis, respectively. This represents quenching within the cluster halo, including any environmental effects such as RPS, starvation, harassment, etc. In this model, \(\tau_\mathrm{major}\) and \(\tau_\mathrm{minor}\) can take different values if intrahalo quenching works anisotropically. After all satellites are classified into SFG or QG, we calculate the quiescent fraction in arbitrary angular bins on the x-y, y-z, and z-x planes. Note that we do not consider any stellar- or halo-mass dependence for simplicity; this will be implemented in future work.

To constrain four parameters, we fit this model to the observationally derived quiescent fraction (see figure 1) via an MCMC procedure adopting the log-likelihood function of equation 3 and uniform prior distributions. Figure 5 shows the derived median posterior predictions of quiescent fraction in the MCMC sampling. Our semi-analytic model (black dashed line) successfully reproduces the observed anisotropy in the quiescent fraction (black points). In table ¿tbl:tab:toymodel?, we summarize the median and 68% credible intervals of four model parameters derived from the posterior distributions (Figure 6) as fiducial results. In addition, we also report the mode with a 68% highest probability density interval of the posterior distribution and a representative posterior sample whose parameter set is closest to the median posterior predictions in a least-squares sense. Although the delay times \(\tau_\mathrm{major}\) and \(\tau_\mathrm{minor}\) are not constrained well, at least the rapid quenching (\(\lesssim1.5\,\mathrm{Gyr}\)) is not preferred. The baseline quiescent faction \(f_\mathrm{q,pre}\) is successfully determined as \(0.46^{+0.08}_{-0.13}\). In addition, the excess quiescent fraction along the major axis, \(\Delta f_\mathrm{q,pre}\), is also constrained as \(0.20^{+0.12}_{-0.12}\). Taking the median value, this suggests that by \(\sim 20\) percent more galaxies are likely to be quenched along the major axis than those along the minor axis. This is evidence that preprocessing occurs anisotropically and is a primary driver of anisotropic quenching signals, even at \(z>1\).

Conceptually, \(f_\mathrm{q,pre}\) and \(f_\mathrm{q,pre}+\Delta f_\mathrm{q,pre}\) should be consistent with the quiescent faction observed in the outskirts of the clusters. We estimate the angular-averaged quiescent fraction to be \(\langle f_\mathrm{q,pre}\rangle+\langle \Delta f_\mathrm{q,pre}\rangle/2=0.56\), where brackets denote median. Although we do not use the quiescent fraction in the cluster outskirts as a constraint in the fitting, this value is roughly consistent with that in the infalling region (i.e., \(1<R/R_{200}<3\)) of the GOGREEN clusters [32], \(f_\mathrm{q}\sim 0.5\) at \(\log(M_{*}/M_{\odot})=10\text{-}11\). It is possible that we constrain our quenching model more effectively if we directly use the quiescent fraction in infall regions as a constraint (cf. [84]), which is left for future work. One caveat is that the quiescent fraction depends strongly on stellar mass, which we do not account for in our model. To test the validity of our model, including its dependence on stellar mass, a larger cluster sample with deeper imaging data is required. We will address these points in future work.

Finally, we examine which quenching pathways dominate for each ‘observed’ angular bin. We divide the quenched galaxies in the model output into four subcategories according to their quenching path: (1) quenched via preprocessing along the major axis, (2) quenched via preprocessing along the minor axis, (3) intrahalo quenching with a timescale of \(\tau_\mathrm{major}\), and (4) intrahalo quenching with a timescale of \(\tau_\mathrm{minor}\). Figure 7 shows the fraction of quenched satellites in each subcategory (\(N_{\mathrm{q},i}^{k}/N_{\mathrm{all},i}\), where \(k\) represents the subcategories) with their median values and 68% intervals during the MCMC sampling.

The preprocessed populations that quenched along the major (orange triangles) and minor (green squares) axes are clearly aligned with the axes along which they were accreted. The sum of these two populations accounts for the majority of the total quenched population (gray circles). In contrast, intrahalo quenched populations with \(\tau_\mathrm{major}\) (magenta diamonds) and \(\tau_\mathrm{minor}\) (blue downward triangles) each contribute at most \(\sim10\) percent and show no apparent anisotropy. Given that the rapid quenching is not favored by the MCMC analysis, group members at \(z=1.25\), most of which are recent infallers, do not spend sufficient time within the halo to become quenched. Overall, our semi-analytic model indicates that preprocessing is the dominant quenching channel, and that the anisotropy in preprocessing between the major- and minor-axis detections is the primary driver of anisotropic quenching at \(z>1\).

5 Discussion↩︎

5.1 Literature comparison↩︎

5.1.1 Anisotropic quenching at \(z>1\)↩︎

Using the GOGREEN/GCLASS cluster sample, we detect anisotropic quenching at \(z>1\), extending such measurements to the highest redshift yet achieved. This finding implies that anisotropic quenching has been operating in clusters in their early evolutionary stages. We compare our results with [48], who have pioneered the exploration of anisotropic quenching across various redshifts.

[48] have investigated anisotropic quenching using a photometric redshift catalog of HSC-SSP [49], [50], [85], [86] and the CAMIRA cluster catalog at \(0<z<1.25\) [87], [88], detecting the anisotropic signature at \(>5\sigma\) significance in the \(0.25<z<1\) redshift range. They have also examined the anisotropy in the \(1<z<1.25\) redshift bin and found the very weak anisotropy with the amplitude of \(0.0175\pm0.0082\) (\(\sim2.1 \sigma\)). [48] have also argued that they detected the hint of the anisotropy at \(z>1\) in the cluster center (\(<0.75 R_{200}\)) if the quiescent fractions only near the major and minor axes (i.e., \(\Delta \theta \pm 15°\)) are compared. However, due to the weakness of the anisotropy at this redshift bin, they have left the detection of anisotropic quenching at \(z>1\) inconclusive. The detection of anisotropic quenching in this work confirms what [48] have inferred. Compared to them, we detect a stronger anisotropic signal by a factor of \(\sim 2.6\), enabling us to overcome significant statistical errors due to a small cluster sample size. However, the origin of the large difference in the observed anisotropy amplitudes is unclear.

While all the clusters used in this work have been spectroscopically confirmed, only a part of the CAMIRA clusters used in [48] have been confirmed, which may introduce contamination into the cluster sample from less massive systems and/or filament-like structures aligned with the line-of-sight by chance, possibly weakening the anisotropy. Another possibility is a cluster-to-cluster variation. Compared to [48], whose sample cluster size is large (\(N=592\)), we only use 12 clusters. It is possible that our cluster sample is biased toward those with a large anisotropy by chance. This effect is included in the relatively large uncertainties in the anisotropy amplitude. In addition, the definition of QGs is different: [48] have defined QGs by the specific star formation rate (sSFR), while we do so by the rest-frame color. Since [48] have used sSFR inferred from the SED fitting only with optical photometry [86], [89], their QGs can be contaminated by dusty star-forming galaxies, which may not significantly affect the color-selected QGs in this work [70]. In any case, a larger confirmed cluster sample at \(z > 1\) and the analysis of these clusters in a self-consistent manner with IR data are crucial for conclusively detecting anisotropic quenching in this high-redshift Universe.

5.1.2 The physical driver of anisotropic quenching↩︎

In the literature, anisotropic quenching has mainly been investigated at \(z\lesssim 0.5\) (e.g., [41], [45]). Some recent works have reported that anisotropic quenching occurs not only inside the cluster virial radius but also beyond it. [46] have used 11 massive clusters at \(0.2<z<0.5\) from the CLASH survey and detected anisotropic quenching out to \(3 R_{200}\) at \(>2\sigma\) level. [47] have investigated anisotropic quenching at \(z=0\) in TNG. They have detected the anisotropy out to \(5 R_{200}\). These studies have explained how anisotropic quenching arises as follows. Clusters are often connected to large-scale cosmic filaments. Because filaments are relatively dense environments and sometimes even host small groups, galaxies within them are quenched more frequently than those outside filaments (e.g., [30]). The major axis of a central cluster galaxy tends to be aligned with the filaments, along which these preprocessed galaxies are accreted. Consequently, quiescent satellites are preferentially aligned with the major axis of the central galaxy. As shown in section 4.3, our semi-analytic model infers that preprocessing is a plausible origin of anisotropic quenching at \(z>1\), consistent with this explanation. This suggests that preprocessing in filaments along the major axis direction of central galaxies is a main physical driver of anisotropic quenching throughout a wide range of cosmic history.

[41] have originally proposed that the diluted ICM density due to AGN feedback from a central galaxy, and consequently, the weakened RPS causes anisotropic quenching. This scenario is not preferred by our semi-analytic model, which predicts only a limited contribution of intrahalo quenching to anisotropic quenching at \(z>1\). In our cluster sample, 3 out of 12 clusters are detected by SZ effects, providing strong evidence for the presence of an ICM. While the other cluster samples are massive enough to host ICM, it is unknown whether they do so. Therefore, there is no clear evidence that supports [41]’s scenario yet. We note that our quenching model incorporates the anisotropy in intrahalo quenching through the initial accretion angles, but a realistic treatment of the ICM distribution would be required to test this scenario more rigorously.

[48] have supported [41]’s explanation rather than the preprocessing scenario because they have not detected the significant anisotropy at \(1<R/R_{200}<2\). This non-detection is opposite to the detection of the anisotropy out to \(3 R_{200}\) by [46]. One possible explanation of the non-detection by [48] is that their sample consists of less massive clusters with \(\log(M_\mathrm{h}/M_{\odot})\sim 14.2\) compared to those of [46] with \(\log(M_\mathrm{h}/M_{\odot})>14.5\). This may lead to lower contrast between filament galaxies and background/foreground galaxies, resulting in the non-detection of anisotropy beyond \(R_{200}\).

Using the TNG data at the \(z=0\) snapshot, [44] have pointed out that satellites located along the major axis are, on average, accreted onto the halo earlier (i.e., have larger TSI) than those along the minor axis, potentially reproducing the anisotropic quenching signal when a fixed quenching timescale in group halo is assumed. We check whether the same trend persists in our quenching model with the TNG dataset at \(z=1.25\). Figure 8 shows the relation between orientation angle and TSI for the satellites. Taking angular bins with \(\Delta \theta=5°\), we calculate the median values and scatters (68% and 90% intervals) of TSI (blue points). The median TSI is almost constant at \(\theta<60°\) but slightly declines towards the minor axis (\(\theta=90\,\mathrm{deg}\)), and the difference between the leftmost and rightmost data point is \(\Delta \mathrm{TSI}=0.345\,\mathrm{Gyr}\). However, the upper sides of 68% and 90% intervals show no clear trend with orientation angles. Considering only the satellites with \(\mathrm{TSI}\gtrsim1.5\,\mathrm{Gyr}\) can contribute to the quiescent fraction, this slight difference in TSI distributions is not likely to contribute significantly to anisotropic quenching.

[32] have pointed out that the majority of massive galaxies in the GOGREEN clusters are likely to be quenched before accretion, either through self-quenching or preprocessing, rather than after accretion. This is consistent with the preprocessing scenario. To the contrary, it is possible that the intrahalo quenching is more effective for less massive galaxies (e.g., \(\log(M_{*}/M_{\odot})<10.5\)). With deeper imaging data that can detect low-mass galaxies, we can test this hypothesis in the future.

5.2 Redshift dependence of anisotropic quenching↩︎

By fitting the cosine function to the angular-binned quiescent fraction, we find the amplitude of the anisotropy to be \(A_\mathrm{q,0}=0.045\pm0.021\). This value is slightly higher than those reported in the lower-redshift observations ([41], [48]). But to make a fair comparison, we need to take care of halo masses of the samples since it has been suggested that groups hosted by more massive halos show the more significant anisotropy [40], [41], [45], [46], [90]. Therefore, when discussing the redshift dependence of anisotropic quenching, group samples with similar halo masses should be compared. Here, we compare our results with those from [48], whose sample halo masses roughly match those in this work.

[48] have detected anisotropic quenching at \(0.25<z<1.0\). Similar to this work, they have used the cosine function to measure the amplitudes of the anisotropic signals. They have derived the amplitudes of the anisotropy in three redshift bins at \(5>\sigma\) level: \(0.0167\pm0.0048\) (\(0.25<z<0.50\)), \(0.0259\pm0.0030\) (\(0.50<z<0.75\)), and \(0.0303\pm0.0041\) (\(0.75<z<1.0\)). With our result \(0.045\pm0.021\) (\(0.9\lesssim z\lesssim1.4\)), the amplitudes of the anisotropic signature seem to increase with redshift monotonically. One possible explanation is that the efficiency of the mechanism responsible for anisotropic quenching changes with redshift. For example, the amplitudes of the anisotropy in preprocessing may increase toward higher redshifts. Another possible explanation is that the fraction of recently accreted satellites is higher for higher-redshift clusters. As discussed in section 4.2, anisotropic quenching is primarily driven by recently accreted satellites. With a small number of satellites that have spent a long time within the cluster, the fraction of recently accreted galaxies relative to the total population will increase in the high-redshift clusters compared to the lower-redshift ones, resulting in stronger anisotropic signals.

To further discuss the redshift dependence of anisotropic quenching, one approach is to extend the semi-analytic modeling in this work to lower redshifts. If redshift-dependent parameters are implemented in the model, we discuss the redshift evolution of preprocessing in a self-consistent manner.

6 Summary and conclusion↩︎

In this paper, we investigate anisotropic quenching in the \(z\gtrsim1\) clusters observed by the GOGREEN and GCLASS surveys. We measure the quiescent fraction as a function of the orientation angle relative to the major axis of the central galaxy and quantify the anisotropy amplitude using a cosine fit. We derive the anisotropy amplitude to be \(A_\mathrm{q}=0.045\pm0.021\), suggesting that anisotropic quenching already occurs even at \(z>1\).

To further investigate the physical origin of anisotropic quenching, we use the Illustris TNG cosmological simulation. We then construct a semi-analytic quenching model based on a delay-then-rapid quenching framework, which has four parameters that account for anisotropy in the intrahalo quenching (\(\tau_\mathrm{major}\), \(\tau_\mathrm{minor}\)) and preprocessing (\(f_\mathrm{q,pre}\), \(\Delta f_\mathrm{q,pre}\)). What we find in the simulation analysis and the semi-analytic modeling is as follows:

1. At the \(z=0\) snapshot, the distribution of TSI is almost flat back to \(5\,\mathrm{Gyr}\). On the other hand, in \(z=1.25\) clusters, most satellite are recently infalled galaxies (\(\mathrm{TSI}<2\,\mathrm{Gyr}\)), and only a small fraction of satellite have \(\mathrm{TSI}>3\,\mathrm{Gyr}\).
   

2. In the \(z=0\) clusters, it takes \(\sim5\,\mathrm{Gyr}\) for the orientation angles of satellite galaxies to be randomized. On the other hand, in the \(z=1.25\) clusters, orientation angles are more immediately randomized in \(\sim2\,\mathrm{Gyr}\) after accretion.
   

3. Our semi-analytic quenching model well reproduces the observed angular dependence of the quiescent fraction. We find that a short quenching timescale (\(\lesssim1.5\,\mathrm{Gyr}\)) is not preferred, and that intrahalo quenching is not significant, since most satellites are recent infallers. Instead, preprocessing is a dominant source of quenching.
   

4. The quiescent fraction along the major axis due to preprocessing exceeds that along the minor axis by \(\sim20\%\), resulting in the observed anisotropic quenching. This suggests that preprocessing is a main driver of anisotropic quenching in high-redshift clusters.
   

Given the detection of anisotropic quenching beyond \(z=1\), this phenomenon is found to be universal across a wide range of cosmic ages. Moreover, if the anisotropic quenching signal reported in this study arises from preprocessing within filaments, this link in turn opens the possibility of using cluster galaxies as a probe of environmental processes in the cosmic web in the early Universe. The large samples of high-redshift galaxy clusters expected from upcoming surveys with the Subaru Prime Focus Spectrograph (PFS), the Rubin Legacy Survey of Space and Time (LSST), Euclid, and the Nancy Grace Roman Space Telescope will enable detailed investigations of environmental effects working in filaments through anisotropic quenching measurements.

References↩︎