# The GOGREEN Survey: Evidence of an Excess of Quiescent Disks in Clusters at 1.0

Jeffrey C. C. Chan, Gillian Wilson, Michael Balogh, Gregory Rudnick, Remco F. J. van der Burg, Adam Muzzin, Kristi A. Webb, Andrea Biviano, Pierluigi Cerulo, M. C. Cooper, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Jeffrey C. C. Chan, Gillian Wilson, Michael Balogh, Gregory Rudnick, Remco F. J. van der Burg, et al.. The GOGREEN Survey: Evidence of an Excess of Quiescent Disks in Clusters at 1.0. The Astrophysical Journal, 2021, 920, 10.3847/1538-4357/ac1117. insu-03720856

HAL Id: insu-03720856

https://insu.hal.science/insu-03720856v1

Submitted on 11 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The GOGREEN Survey: Evidence of an Excess of Quiescent Disks in Clusters at $1.0 < z < 1.4$

Jeffrey C. C. Chan¹, Gillian Wilson¹, Michael Balogh²,³, Gregory Rudnick⁴, Remco F. J. van der Burg⁵, Adam Muzzin⁶, Kristi A. Webb²,³, Andrea Biviano⁷,⁸, Pierluigi Cerulo⁹, M. C. Cooper¹⁰, Gabriella De Lucia⁷, Ricardo Demarco⁹, Ben Forrest¹, Pascale Jablonka¹¹,¹², Chris Lidman¹³, Sean L. McGee¹⁴, Julie Nantais¹⁵, Lyndsay Old¹⁶, Irene Pintos-Castro¹⁷, Bianca Poggianti¹⁸, Andrew M. M. Reeves²,³, Benedetta Vulcani¹⁸, Howard K. C. Yee¹⁷, and Dennis Zaritsky¹⁹

¹ Department of Physics & Astronomy, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA; jeffreyrdcs@gmail.com

² Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

³Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, Ontario, N2L3G1, Canada

⁴ Department of Physics and Astronomy, The University of Kansas, Malott Room 1082, 1251 Wescoe Hall Drive, Lawrence, KS 66045, USA

⁵European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748, Garching, Germany

⁶ Department of Physics and Astronomy, York University, 4700 Keele Street, Toronto, Ontario, ON M3J 1P3, Canada

⁷INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143, Trieste, Italy

⁸IFPU—Institute for Fundamental Physics of the Universe, via Beirut 2, I-34014 Trieste, Italy

⁹ Departamento de Astronomía, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Concepción, Chile

¹⁰ Department of Physics and Astronomy, University of California, Irvine, 4129 Frederick Reines Hall, Irvine, CA 92697, USA

¹¹Laboratoire d'astrophysique, École Polytechnique Fédérale de Lausanne (EPFL), 1290 Sauverny, Switzerland

¹²GEPI, Observatoire de Paris, Université PSL, CNRS, Place Jules Janssen, F-92190 Meudon, France

¹³The Research School of Astronomy and Astrophysics, Australian National University, ACT 2601, Australia

¹⁴School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK

¹⁵ Departamento de Ciencias Físicas, Universidad Andres Bello, Fernandez Concha 700, Las Condes 7591538, Santiago, Región Metropolitana, Chile

¹⁶European Space Agency (ESA), European Space Astronomy Centre (ESAC), E-28691 Villanueva de la Cañada, Madrid, Spain

¹⁷ Department of Astronomy and Astrophysics, University of Toronto 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada

¹⁸ INAF-Padova Astronomical Observatory, Vicolo dell'servatorio 5, I-35122 Padova, Italy

¹⁹Steward Observatory and Department of Astronomy, University of Arizona, Tucson, AZ, 85719, USA

Received 2021 April 30; revised 2021 June 29; accepted 2021 July 2; published 2021 October 8

## Abstract

We present the results of the measured shapes of 832 galaxies in 11 galaxy clusters at $1.0 < z < 1.4$ from the GOGREEN survey. We measure the axis ratio ($q$), the ratio of the minor to the major axis, of the cluster galaxies from near-infrared Hubble Space Telescope imaging using Sérsic profile fitting and compare them with a field sample. We find that the median $q$ of both star-forming and quiescent galaxies in clusters increases with stellar mass, similar to the field. Comparing the axis ratio distributions between clusters and the field in four mass bins, the distributions for star-forming galaxies in clusters are consistent with those in the field. Conversely, the distributions for quiescent galaxies in the two environments are distinct, most remarkably in $10.1 \le \log(M/M_\odot) < 10.5$ where clusters show a flatter distribution, with an excess at low $q$. Modelling the distribution with oblate and triaxial components, we find that the cluster and field sample difference is consistent with an excess of flattened oblate quiescent galaxies in clusters. The oblate population contribution drops at high masses, resulting in a narrower $q$ distribution in the massive population than at lower masses. Using a simple accretion model, we show that the observed $q$ distributions and quenched fractions are consistent with a scenario where no morphological transformation occurs for the environmentally quenched population in the two intermediate-mass bins. Our results suggest that environmental quenching mechanism(s) likely produce a population that has a different morphological mix than those resulting from the dominant quenching mechanism in the field.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Galaxy clusters (584); High-redshift galaxy clusters (2007)*

## 1. Introduction

It is well-established that the environment of a galaxy plays a crucial role in its evolution. In the local Universe, the galaxy population in high-density environments comprises mainly galaxies that have ceased forming stars. The dominance of quiescent galaxies in groups and clusters, as reflected by the higher quiescent fraction at fixed stellar mass compared to the field (e.g., Balogh et al. 2004; Baldry et al. 2006; Wetzel et al. 2012), suggests that there are physical processes that correlate with the environment to suppress star formation. These quiescent galaxies are composed of mostly early-type objects, as opposed to the late-type morphologies that are seen in star-forming galaxies (e.g., Dressler 1980; Postman et al. 2005;

Holden et al. 2007; Bassett et al. 2013). This implies that a morphological change must have taken place at a certain evolutionary stage. Despite the focused effort in recent decades, the physical processes that drive the quenching of star formation and the morphological transformation of the galaxies in dense environments are not yet fully understood.

Detailed studies of the properties of the galaxy population in clusters and groups in the local Universe and low redshifts have revealed various mechanisms that can contribute to environmental quenching (see, e.g., Boselli & Gavazzi 2006, 2014, for reviews). For example, the cut-off of the cold gas accretion from the cosmic web during infall into a massive halo can gradually quench the star formation of a galaxy as the fuel slowly runs out ("strangulation" or "starvation", Larson et al. 1980; Balogh et al. 1997, 2000).
Quenching can also occur due to rapid removal of the cold gas in the galaxies when it passes through the intracluster medium (ICM) ("ram pressure stripping", Gunn & Gott 1972) or due to interactions between galaxies with other group or cluster members ("galaxy harassment", e.g., Moore et al. 1998). Nevertheless, the relative importance of each of these mechanisms is still not well understood, in part because the efficiency of these mechanisms depends on both the properties of the galaxies (e.g., gas content or star formation rate) and the environment that they are in (e.g., halo mass or ICM density). Environmental quenching at low redshift is shown to be largely separable from quenching driven by mechanisms that act internally in the galaxy (i.e., mass-quenching, e.g., Peng et al. 2010). One interpretation is that environmental quenching operates independently and does not depend strongly on stellar mass (but see also De Lucia et al. 2012; Wetzel et al. 2013; Fillingham et al. 2015). Nevertheless, there is growing evidence that the situation is very different at $z \gtrsim 1$. Recent works have reported a mass dependence in the environmental quenching efficiencies at redshift $z \gtrsim 1$ (Cooper et al. 2010; Balogh et al. 2016; Kawinwanichakij et al. 2017; Fossati et al. 2017; Papovich et al. 2018; Pintos-Castro et al. 2019; van der Burg et al. 2020), which suggests that the effects from both classes are no longer separable. This points to a possible change in the dominant environmental quenching mechanism at high redshift (Balogh et al. 2016).

The studies of environmental quenching efficiency that were mentioned above mostly rely on measuring the stellar mass function of the star-forming and quiescent galaxy population as a function of redshift and environment. The relative fraction of the two populations provides important constraints on galaxy evolution in different environments. Additional, complementary, information about the morphologies or structural properties of the galaxy population is also often considered. Many of the proposed environmental mechanisms have unique implications or predictions on the morphology of the galaxies. The most striking example of them all is the stripped gas tails that are produced by gas removal processes such as ram pressure stripping, which are easily recognizable by their peculiar morphologies in $H\alpha$ imaging and spectroscopy (e.g., Gavazzi et al. 2001; Fumagalli et al. 2014; Yagi et al. 2015; Sheen et al. 2017). Significant efforts have been put into searching for galaxies that exhibit tails reminiscent of a debris trail, known as "jellyfish" galaxies, in groups and clusters at low redshifts (e.g., Poggianti et al. 2016; McPartland et al. 2016; Roberts & Parker 2020). Similarly, galaxies with peculiar morphologies, such as merging pairs, tidal features, and truncated or warped disks, are often treated as the proof of the existence of the corresponding mechanisms (i.e., mergers, harassment, stripping). The structural properties of galaxies have also provided crucial insights into the evolutionary path of the galaxy population. For example, studies of the quiescent galaxy population in clusters and the field at high redshift have shown that they are on average more compact than their local counterparts of the same mass (e.g., Trujillo et al. 2006; Newman et al. 2012; van der Wel et al. 2014a; Chan et al. 2018; Matharu et al. 2019), which suggests that they must have undergone significant evolution in size but only mild growth in mass (but also see Valentinuzzi et al. 2010; Cooper et al. 2012; Lani et al. 2013; Poggianti et al. 2013; Delaye et al. 2014). Repeated minor mergers have been shown to be the primary mechanism that gives rise to the observed size evolution and the inside-out growth of the galaxies (e.g., Naab et al. 2009;

van Dokkum et al. 2010; Shankar et al. 2013; Suess et al. 2019), although the effect of continual arrival of larger quenched galaxies may also play a role (i.e., progenitor bias, e.g., van Dokkum & Franx 2001; Saglia et al. 2010; Carollo et al. 2013; Poggianti et al. 2013; Belli et al. 2015; Matharu et al. 2020).

The projected axis ratio (ellipticity) distribution of the galaxy population has also long been used to study their intrinsic structural properties and shapes (e.g., Sandage et al. 1970; Franx et al. 1991; Tremblay & Merritt 1996). Although individual axis ratios do not carry much information because they are degenerate with the inclination angle, their distribution can be used to infer the intrinsic shape distribution under the assumption of random viewing angles. From the studies that have been published in the last few decades, it is established that the majority of star-forming galaxies in the local Universe are flattened, oblate systems (e.g., Ryden 2004; Padilla & Strauss 2008). Furthermore, van der Wel et al. (2014b) showed that this is also true for the more massive $log(M/M_\odot) > 10.0$ star-forming galaxies at high redshift, up to $z \sim 2.5$. Meanwhile, the projected axis ratio distribution of the early-type population in the local Universe requires a two-component model, which comprises a triaxial set and an oblate set of objects, to well describe its properties (e.g., Tremblay & Merritt 1996; Holden et al. 2012). The exceptions are the massive quiescent galaxies with $log(M/M_\odot) > 11.0$, where they are preferentially round and can be described by a single triaxial population (e.g., van der Wel et al. 2009). Chang et al. (2013b) extended this axis ratio analysis to quiescent galaxies at $1 < z < 2.5$ in the field and found that the fraction of oblate galaxies relative to the total population evolves over redshift. For massive quiescent galaxies with $log(M/M_\odot) > 11.0$, the oblate fraction is almost three times higher at $z > 1$.

This two-component picture is also supported by the observed stellar kinematics of the low-$z$ quiescent galaxies (i.e., the slow rotators and fast rotators, e.g., Emsellem et al. 2011). Recent integral field spectroscopy studies have shown that the intrinsic shape of a galaxy is correlated to the degree of rotational support. For example, Weijmans et al. (2014) found that fast rotators have flattened intrinsic shape distributions that are similar to spiral galaxies, while slow rotators are likely to be mildly triaxial (see also Cortese et al. 2016; Pulsoni et al. 2018). With a larger sample, Foster et al. (2017) showed that galaxies with higher "spin" parameter (Emsellem et al. 2011) have more flattened intrinsic axis ratios and are more likely to be axisymmetric systems.

The projected axis ratio distribution has also been used to study the formation of lenticular galaxies (S0s), which are abundant in local galaxy clusters. For example, Vulcani et al. (2011) studied the axis ratio distributions of a sample of early-type galaxies in intermediate redshift clusters ($z \sim 0.6$) and compared them to those in local clusters. They found that there are fewer flattened objects in the intermediate redshift sample due to a lower fraction of S0 galaxies. Similar studies also found that the S0 fraction drops rapidly with increasing $z$. By $z \sim 0.5$ the fraction of S0 galaxies is found to be <10% (e.g., Fasano et al. 2000; Postman et al. 2005). The high occurrence of S0 in low-$z$ clusters, although not fully understood, is generally believed to be due to environmental effects (e.g., Just et al. 2010; Johnston et al. 2014; Kelkar et al. 2017). Since the mass dependence of environmental quenching efficiencies
emerges at $z \gtrsim 1$, it is interesting to extend the axis ratio distribution studies to even higher redshift.

In this paper, we investigate the axis ratio distributions of the galaxies in 11 clusters of the Gemini Observations of Galaxies in Rich Early ENvironments survey (GOGREEN; Balogh et al. 2017, 2021) at $1.0 < z < 1.4$. This recently completed survey is an imaging and spectroscopic survey that targets 21 known high-redshift overdensities that are representative of the progenitors of the clusters we see today. The deep spectroscopy and imaging of GOGREEN allows us to study the axis ratio distributions of an unprecedentedly large sample of cluster galaxies in this redshift range. The goal of this work is to study the effect of environment on galaxy structures by comparing the axis ratio distributions of cluster galaxies to those in the general field. The field comparison sample is taken from the CANDELS (Grogin et al. 2011; Koekemoer et al. 2011), and 3D-HST Treasury programs (Brammer et al. 2012; Skelton et al. 2014).

This paper is organized as follows. In Section 2, we describe the data set that is used in this work and we present the derivation of structural parameters and other quantities. We present the results of the axis ratio distributions and describe the procedure and results of the axis ratio modeling in Section 3. We then explore the relationship between environmental quenching and morphological transformation and discuss the results in Section 4. In Section 5, we draw our conclusions.

Throughout this paper, we assume the standard flat cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_\Lambda = 0.7$ and $\Omega_m = 0.3$. Magnitudes quoted are in the AB system (Oke & Gunn 1983). The stellar masses in this paper are computed with a Chabrier (2003) initial mass function (IMF).

## 2. Sample and Data

The cluster sample that is used in this work is from the GOGREEN survey (Balogh et al. 2017, 2021). The GOGREEN sample consists of 21 overdensities at $1.0 < z < 1.5$ spanning a wide range of halo masses, including three clusters from the South Pole Telescope (SPT) survey (Brodwin et al. 2010; Foley et al. 2011; Stalder et al. 2013), nine clusters from the Spitzer Adaptation of the Red-sequence Cluster Survey (SpARCS, Muzzin et al. 2009; Wilson et al. 2009), of which five were followed up extensively by the Gemini Cluster Astrophysics Spectroscopic Survey (GCLASS, Muzzin et al. 2012), and nine group candidates selected in the COSMOS and Subaru-XMM Deep Survey (SXDS) fields.

In this study, we focus on 11 GOGREEN clusters at $1.0 < z < 1.4$ that have complete spectroscopic and photometric catalogs at the time of this work.<fnref n="20" /> Eight of the clusters were discovered using the red-sequence or the stellar-bump technique (Muzzin et al. 2009; Wilson et al. 2009; Demarco et al. 2010). The remaining three clusters were discovered via the Sunyaev-Zeldovich effect signature (Bleem et al. 2015). The properties of the clusters are summarized in Table 1.

The main spectroscopic data set of GOGREEN was obtained from a Gemini Large and Long Program (GS LP-1 and GN LP-4; PI Balogh) using the Gemini Multi-Object Spectrographs (GMOS) on Gemini-North and South. This large program

allows us to obtain unbiased spectroscopy of galaxies of all types down to stellar masses of $M_* \gtrsim 10^{10.3} M_\odot$, with the faintest targets having exposure time up to 15 hr. Five of the GOGREEN clusters that are part of GCLASS also have similar GMOS spectroscopy data for the bright galaxies in the clusters. These data have been incorporated into the GOGREEN spectroscopy sample. For the specifics of the targeting selection and data reduction, we refer the reader to the survey and data release papers (Balogh et al. 2017, 2021).

In addition to deep spectroscopy, GOGREEN has obtained deep multi-band imaging (UBVRIzYJK and IRAC 3.6 $\mu$m) for the sample. We begin by describing in detail the derivation of structural properties of galaxies, the key focus of this paper, from HST imaging in Sections 2.1 and 2.2. The deep multi-wavelength photometric data also allows us to characterize galaxy SEDs and derive photometric redshifts, stellar population parameters, and rest-frame colors. A brief summary of the derivation of these properties is given in Sections 2.3, 2.4 and 2.5. We refer the reader to van der Burg et al. (2020) for more details.

### 2.1. HST Observations and Data Reduction

We make use of the near-infrared HST/WFC3 F160W imaging of the GOGREEN clusters to quantify the structural properties of the galaxies. The HST/WFC3 F160W images were obtained in a Cycle 25 program (GO-15294; PI: Wilson), which is dedicated to studying galaxy morphologies. Each cluster was targeted with a 1 × 2 mosaic of WFC3 pointings centered on the cluster, covering a region of $136'' \times 233''$. At the redshift of the GOGREEN clusters, this corresponds to a ~1.1 × 1.9 Mpc rectangular region on the sky. Each pointing has 1-orbit depth. We constrained the ORIENT to within 20° of the GMOS mask orientation to maximize the overlap between the imaging and the GMOS spectroscopy.

The data are reduced and combined using ASTRODRIZZLE (version 2.1.22) (Gonzaga et al. 2012). All of the calibrated frames (`.flt.fits`) were downloaded from the Mikulski Archive for Space Telescopes (MAST) archive and are first examined to check the quality of the cosmic rays and bad pixels identification by the `calwf3` pipeline. We find hot stripes that span across the field of view (FOV) in two of the frames (e.g., due to satellite trails), which are not fully flagged by the pipeline. We mask these regions generously in the data quality array of the `flt` files. In addition, a total of seven frames in five clusters show a smooth background gradient, presumably due to earthshine. To remove the gradient, we follow a similar approach described in Windhorst et al. (2011). Sources on the `flt` image are first masked, and then the gradient is fitted with a fifth-order bivariate spline function and is subtracted from the frame before drizzling.

For the final drizzling, we adopt a pixel scale of 0.0"06 pixel⁻¹, a square kernel, and a `pixfrac` of 0.8. We produce weight maps using both inverse variance map (IVM) and error map (ERR) weighting for different purposes. The IVM weight maps, which contain all background noise sources except Poisson noise of the objects, are used for object detection, while the ERR weight maps are used for structural analysis because the Poisson noise of the objects is included. The final images and the weight maps are included in the first GOGREEN public release (Balogh et al. 2021). The characteristic point-spread function (PSF) of each cluster is constructed by median-stacking of the isolated bright unsaturated stars.

<footnote n="20">²⁰ The one cluster in the 12 GOGREEN clusters that was not included is SpARCS 1033 because deep K-band imaging has not yet been obtained at the time of this work.</footnote>
<table><caption>**Table 1**
Summary of the Gogreen Cluster Sample Used in this Study</caption><thead><tr><th>Name</th><th>R.A.<sub>J2000</sub></th><th>Decl.<sub>J2000</sub></th><th>Redshift</th><th>σ<sub>v</sub><sup>a</sup><br>(kms<sup>-1</sup>)</th><th>M<sub>200</sub><sup>a</sup><br>(10<sup>14</sup>M<sub>☉</sub>)</th><th>R<sub>200</sub><sup>a</sup><br>(Mpc)</th><th>N<sub>mem,HST</sub><sup>b</sup></th></tr></thead><tbody><tr><td>SpARCS J1051+5818</td><td>10:51:11.23</td><td>+58:18:02.7</td><td>1.035</td><td>689 ± 36</td><td>2.51<sup>+0.65</sup><sub>-0.55</sub></td><td>0.88 ± 0.07</td><td>51</td></tr><tr><td>SPT-CL J0546-5345</td><td>05:46:33.67</td><td>-53:45:40.6</td><td>1.067</td><td>977 ± 68</td><td>6.11<sup>+1.52</sup><sub>-1.30</sub></td><td>1.17 ± 0.09</td><td>151</td></tr><tr><td>SPT-CL J2106-5844</td><td>21:06:04.59</td><td>-58:44:27.9</td><td>1.132</td><td>1055 ± 83</td><td>7.65<sup>+2.02</sup><sub>-1.72</sub></td><td>1.23 ± 0.10</td><td>145</td></tr><tr><td>SpARCS J1616+5545</td><td>16:16:41.32</td><td>+55:45:12.4</td><td>1.156</td><td>782 ± 39</td><td>3.29<sup>+0.69</sup><sub>-0.60</sub></td><td>0.92 ± 0.06</td><td>85</td></tr><tr><td>SpARCS J1634+4021</td><td>16:34:37.00</td><td>+40:21:49.3</td><td>1.177</td><td>715 ± 37</td><td>2.66<sup>+0.60</sup><sub>-0.52</sub></td><td>0.85 ± 0.06</td><td>57</td></tr><tr><td>SpARCS J1638+4038</td><td>16:38:51.64</td><td>+40:38:42.9</td><td>1.196</td><td>564 ± 30</td><td>1.52<sup>+0.42</sup><sub>-0.36</sub></td><td>0.71 ± 0.06</td><td>49</td></tr><tr><td>SPT-CL J0205-5829</td><td>02:05:48.19</td><td>-58:28:49.0</td><td>1.320</td><td>678 ± 57</td><td>2.22<sup>+0.89</sup><sub>-0.70</sub></td><td>0.76 ± 0.09</td><td>61</td></tr><tr><td>SpARCS J0219-0531</td><td>02:19:43.56</td><td>-05:31:29.6</td><td>1.325</td><td>810 ± 77</td><td>2.51<sup>+1.33</sup><sub>-0.98</sub></td><td>0.79 ± 0.12</td><td>47</td></tr><tr><td>SpARCS J0035-4312</td><td>00:35:49.68</td><td>-43:12:23.8</td><td>1.335</td><td>840 ± 52</td><td>4.14<sup>+1.00</sup><sub>-0.87</sub></td><td>0.93 ± 0.07</td><td>90</td></tr><tr><td>SpARCS J0335-2929</td><td>03:35:03.56</td><td>-29:28:55.8</td><td>1.368</td><td>542 ± 33</td><td>1.60<sup>+0.65</sup><sub>-0.51</sub></td><td>0.67 ± 0.08</td><td>47</td></tr><tr><td>SpARCS J1034+5818</td><td>10:34:49.47</td><td>+58:18:33.1</td><td>1.385</td><td>250 ± 28</td><td>0.08<sup>+0.03</sup><sub>-0.03</sub></td><td>0.24 ± 0.03</td><td>49</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The cluster mass $M_{200}$ and $R_{200}$ (radius where the mass overdensity is 200 times the critical density at the cluster redshift) are derived from a scaling relation with the velocity dispersions $\sigma_v$. See Biviano et al. (2021) and Old et al. (2020) for details.

<sup>b</sup> $N_{mem,HST}$ is the number of galaxies that are spectroscopically and photometrically selected as cluster members, are within the HST image FOV, and have a good structural fit.

Depending on the cluster, 5–22 stars are used in the stack. The full-width-half-maximum (FWHM) of the PSFs are ~0''17–0''18.

## 2.2. Structural Parameters

We derive structural parameters for all sources in the F160W image of each cluster by fitting them with two-dimensional single Sérsic profiles (Sersic 1968). The parameters are derived using a modified version of GALAPAGOS (based on v.2.3.1) (Barden et al. 2012; Häußler et al. 2013) with GALFITM (v.1.2.1). The five independent parameters of the Sérsic profile —namely the total luminosity ($L_{tot}$), the Sérsic index ($n$), the half-light radius / effective semimajor axis ($R_e$), the axis ratio ($q = b/a$, where $a$ and $b$ are the major and minor axis respectively) and the position angle ($P.A.$), as well as the centroid $(x, y)$ of the source—are left as free parameters. Source detection from the HST image and initial guesses for fitting these parameters were derived by running SEXTRACTOR, which is incorporated in the GALAPAGOS run.

We apply fitting constraints of $0.2 < n < 12$, $0.3 < R_e < 400$ (pix), $0 < \text{mag} < 40$, $0.0001 < q < 1$, and $-180^\circ < P.A. < 180^\circ$. The local sky level for each source is fixed to the value determined by GALAPAGOS, which is derived using an elliptical annulus flux growth method. We derive noise maps (rms noise) from the ERR weight maps output by ASTRODRIZZLE and use them as sigma map input for GALFIT. The noise maps that we generate from ERR weight maps are a more realistic representation of the noise than the internal error estimation in GALFIT because they include pixel-to-pixel exposure time differences originating from image drizzling and dithering patterns in observations, as well as a more accurate estimation of shot noise. The Sérsic model is convolved with the characteristic PSF of each cluster. Nearby objects that are close to the primary source of interest are fitted simultaneously.

We refine the GALAPAGOS configuration parameters, including those that control local sky level estimation and close neighbor treatment, using extensive tests with simulated galaxies. The details and result of these simulations are provided in Appendix B. In brief, we inject a set of 20000 simulated galaxies (20 at a time), with surface brightness

profiles described by a Sérsic profile, to random locations in the sky region of the F160W image and recover their structural parameters with our science setup. Using this set of simulated galaxies, we then compute the biases of our measurements, modify the GALAPAGOS setup, and re-derive the parameters and biases. This process is iterated a few times to get the best configuration parameters that minimize the biases.

The simulation allows us to characterize the biases in our structural parameter measurements. Among the three Sérsic structural parameters that describe the shape of a galaxy ($R_e$, $n$, and $q$), the axis ratio $q$ can typically be measured with the highest accuracy. The axis ratio shows an average bias and dispersion of ~2% and ~12% at F160W = 23 (AB), which corresponds roughly to $\log(M/M_☉) \sim 9.5$, which is the mass limit that we adopted in this work. This is a factor of five (two) better than the average bias (dispersion) we see in Sérsic index $n$. We stress that while Sérsic index is useful for morphological selection, it is not straightforward to compare their distribution between different samples. From our simulations, we find that biases in $n$ depend also heavily on $n$ itself, such that high $n$ values are more uncertain (see also van der Wel et al. 2012, for a description of systematic uncertainties). These systematic errors can have a detrimental effect on the cluster and field comparison, especially for low-mass galaxies, because the expected difference is small at this redshift (see, e.g., Chan et al. 2018; Matharu et al. 2019, for the difference in $n$ for a sample of massive galaxies). Therefore, in this work we focus primarily on the axis ratio distributions.

We also visually inspect outliers that have large sizes for their particular magnitudes or parameters that hit the boundary of the constraints with a procedure similar to Chan et al. (2016). In cases where sources or nearby objects are not correctly deblended, extra Sérsic components are added iteratively if necessary to ensure adjacent sources are well-fitted. Fits that still hit the boundary of our fitting constraints are considered as bad fits and are excluded from subsequent analyses.

## 2.3. Photometric Catalogue

We utilize the $K_s$-band selected photometric catalog that is derived from the multi-band imaging of each cluster. We refer
to van der Burg et al. (2020) for details of the procedure to construct these catalogs. Source detection is performed on the Ks-band image using SExTRACTOR (Bertin & Arnouts 1996). Aperture photometry is measured on the PSF-matched images using circular apertures with a diameter of 2". To preserve the spatial resolution of the ground-based imaging, aperture photometry of Spitzer/IRAC data is measured with a larger aperture of 3" and rescaled, following the approach in van der Burg et al. (2013). The area covered by the catalogs range from ~5' × 5' to ~10' × 10' depending on the cluster. Therefore, the area that is considered for this study is limited by the HST imaging coverage.

## 2.4. Spectroscopic and Photometric Redshifts

Spectroscopic redshifts ($z_{\text{spec}}$) are measured using the Manual and Assisted Redshifting software (MARZ, Hinton et al. 2016), which utilizes a cross-correlation algorithm to match the spectra against a variety of spectral templates. These are supplemented with publicly available $z_{\text{spec}}$ from various surveys. An exhaustive list of surveys can be found in van der Burg et al. (2020). Photometric redshifts ($z_{\text{phot}}$) are derived using EAZY (Brammer et al. 2008) with standard templates. In this work, we use the peak of the posterior probability distribution of the redshift estimated with EAZY as $z_{\text{phot}}$. A correction, in the form of a quadratic function, has been applied to the $z_{\text{phot}}$ to minimize the residual between the measured $z_{\text{spec}}$ and $z_{\text{phot}}$ (van der Burg et al. 2020). In this work, we derive an additional correction to the $z_{\text{phot}}$ of each cluster to better match the $z_{\text{spec}}$ at the cluster redshift. The correction is taken as the median offset between the $z_{\text{phot}}$ and the $z_{\text{spec}}$ of the cluster members (see Section 2.6 for a description of the cluster membership). The magnitude of this correction is generally small, but in some cases can reach up to ~0.09.

## 2.5. Stellar Mass Estimates and Rest-frame Colors

The stellar masses for all of the galaxies are inferred from the multi-band photometry using FAST (Kriek et al. 2009), derived in van der Burg et al. (2020). This includes a rescaling factor that is applied to the input aperture fluxes, such that SED fitting gives the total mass of the galaxy. This factor is taken to be the ratio of Ks-band FLUX_AUTO measurements to the aperture flux from SExTRACTOR (i.e., $F_{\text{auto}}/F_{\text{aper}}$). We use the Bruzual & Charlot (2003) stellar population synthesis models and we assume a Chabrier (2003) IMF, solar metallicity, and the Calzetti et al. (2000) dust law. The star formation history (SFH) is parameterized as an exponentially declining history $SFR \propto e^{-t/\tau}$, where the timescale $\tau$ ranges between 10 Myr and 10 Gyr. We note that stellar masses derived in this way can typically be 0.2 dex lower than those derived from non-parametric star formation histories (Leja et al. 2019; Webb et al. 2020). Nevertheless, as we describe in Section 2.7, the method that we used to derive stellar masses in GOGREEN is largely consistent with our chosen field sample and therefore has the advantage of allowing us to compare the stellar masses directly.

In this work, we utilize the rest-frame UVJ color classification to separate the galaxies into star-forming and quiescent. UVJ classification has become a standard technique in galaxy evolution studies because it can separate “genuine” quiescent galaxies from dusty star-forming galaxies (e.g., Labbé et al. 2005; Williams et al. 2009; Muzzin et al. 2013a). The

rest-frame U-V and V-J colors required for this classification are again derived using EAZY. We adopt the UVJ color classification criteria in Muzzin et al. (2013b). The best redshift estimate of individual galaxies ($z_{\text{phot}}$ for those lacking $z_{\text{spec}}$) is used to measure the rest-frame colors. We also computed the colors by fixing the redshifts of all galaxies in each cluster field to the cluster mean redshift listed in Table 1, and confirm this does not change our conclusion.

## 2.6. Cluster Membership and Sample Selection

We use both the spectroscopic and photometric redshift information of the GOGREEN sample to define cluster membership. For galaxies with spectroscopic redshifts, we define them as cluster members if they are ±2000 km s⁻¹ around the mean redshift of the cluster. This corresponds to $\Delta z_{\text{spec}} \sim 0.015$ at GOGREEN redshifts. Our simple cluster membership selection is slightly different from those used in previous GOGREEN papers (e.g., Old et al. 2020; Balogh et al. 2021; Biviano et al. 2021) but are fully adequate for the scope of the present analysis given that we also include members that are selected photometrically.

For galaxies without spectroscopic information, we use their photometric redshifts to determine membership. Since the main goal of this work is to compare the axis ratio of cluster galaxies with the field, it is important to strike a balance between maximizing the number of cluster members (i.e., sample completeness) and compromising the purity of the sample. To find the optimal $z_{\text{phot}}$ selection criteria, we compare the $z_{\text{spec}}$ of the GOGREEN spectroscopic sample with their $z_{\text{phot}}$ to estimate the completeness and purity as a function of the $z_{\text{phot}}$ selection width. The result of this test is given in Appendix A. This test is possible because the spectroscopic sample is a representative subset of the photometrically selected galaxy population (see Appendix A.2 in van der Burg et al. 2020, for a discussion). Photometric cluster members are defined as those with $\Delta z_{\text{phot}}/(1+z_{\text{phot}}) \le 0.06$ around the mean redshift of the cluster. Our test suggests that this selection gives a completeness of ~85% and a purity of ~80%.

We crossmatch the structural parameter catalog (i.e., detected from the F160W image) with the Ks-band selected photometric catalogs. Adopting the above mentioned cluster membership selection and a stellar mass limit of $\log(M/M_{\odot}) = 9.5$ results in 860 cluster members that are within the HST image FOV. This stellar mass limit corresponds to a ~80% completeness of the catalogs (see van der Burg et al. 2020, for details of the completeness characterization). After excluding bad structural fits (i.e., fits that hit the boundaries of the fitting constraints), our final sample contains 832 cluster members that have robust structural parameters. In the left-hand column of Figure 1, we show the rest-frame U-V and V-J color distribution of the cluster sample, color coded by their axis ratios and Sérsic indices.

## 2.7. Field Comparison Sample

Although a large number of spectroscopically-confirmed field galaxies are available in GOGREEN, the number of field galaxies within the HST image FOV is still too small for a morphology comparison between clusters and the field within GOGREEN. Expanding this field sample with photometric redshifts is not straightforward because the number of galaxies declines sharply with the redshift selection. We end up with
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Rest-frame UVJ diagram of the GOGREEN cluster sample and the field sample used in this study. Top: Data points are color coded by their axis ratio *q*. As expected, both quiescent and star-forming galaxies in clusters and the field exhibit a wide range of axis ratios. Galaxies that have both red *U* − *V* and *V* − *J* colors can be predominantly seen with small *q*, as the fact that they are viewed edge on result in a higher dust extinction value. Bottom: Data points color coded by their Sérsic index *n*. As pointed out in various works, quiescent and star-forming galaxies in both clusters and the field show discernible differences in their *n* distributions. The number of quiescent (*N*<sub>Q</sub>) and star-forming (*N*<sub>SF</sub>) galaxies in each sample are provided next to the UVJ classification line. The cluster BCGs are marked with a star symbol. The cluster sample shows a higher relative abundance of quiescent galaxies compared to star-forming galaxies than the field.</figcaption></figure>

either a small sample or a sample with low purity (see
Appendix A for a discussion).

The field comparison sample that we use in this study is taken from the CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) and 3D-HST Treasury programs (Brammer et al. 2012; Skelton et al. 2014). Among the ~99000 3D-HST grism redshift measurements (Momcheva et al. 2016), ~5200 redshifts are within the range of 0.9 < *z* < 1.5. For structural parameters, we use the F160W-band measurements of all five CANDELS/3D-HST fields (COSMOS, GOODS-N, GOODS-S, EGS, and UDS) derived by van der Wel et al. (2014a). The CANDELS F160W wide imaging has, on average, one and one-third orbit depth, thus being comparable in depth to the GOGREEN imaging. To ensure that our structural parameter measurements are compatible with the van der Wel et al. (2014a) measurements, we apply our methodology described in Section 2.2 to the CANDELS imaging for a sample of galaxies in the redshift range of 0.9 < *z* < 1.5. Overall we find that the median ratio and 1σ uncertainty between our measurements and van der Wel et al. (2014a) are 1.00 ± 0.03 (0% ± 3%) down to the mass limit of this work. The result of the comparison is shown in Appendix C. We also check that using either set of measurements gives a consistent conclusion. Throughout this work, we show the field results using the van der Wel et al. (2014a) measurements.

Stellar masses, photometric redshifts, and rest-frame colors
are estimated using FAST and EAZY from multi-band photo-
metry (Skelton et al. 2014), in a way that is almost identical to
GOGREEN (van der Burg et al. 2020). There are two main

differences. First, Skelton et al. (2014) adopted a minimum
timescale τ of 40 Myr as opposed to 10 Myr. The second subtle
difference is on the definition of the total fluxes, which affects
the stellar mass estimates. On top of rescaling the aperture
fluxes to FLUX_AUTO measurements like in GOGREEN (see
Section 2.5), Skelton et al. (2014) also factored in a correction
to account for the missing flux that falls outside the AUTO
aperture, which is determined from measuring the growth
curves of the F160W PSFs (the detection band of the 3D-HST
catalog). To ensure that the stellar masses are comparable, we
apply a correction to both the 3D-HST and GOGREEN stellar
masses, rescaling the stellar masses to the total F160W fluxes
of the best-fit Sérsic profile of the galaxies. Overall this
correction is small—it only increases the stellar mass by ~0.02
(3D-HST) and ~0.03 (GOGREEN) dex on average, although
in some cases it can exceed 0.1 dex. The fact that GOGREEN
galaxies require a slightly larger correction is accordant with
the additional missing flux correction that is applied in
3D-HST.

We select galaxies that (i) are in the redshift range of
0.9 < z < 1.5, (ii) with a stellar mass of log(M/M<sub>☉</sub>) ≥ 9.5, and
iii) have robust structural parameters as our field sample. The
selection is done using the z<sub>best</sub> catalogs (v4.1.5), which uses
ground-based z<sub>spec</sub> of the galaxies if available, then grism
redshift z<sub>grism</sub>, and finally z<sub>phot</sub> if the other two are not available.
These selection criteria result in a sample of 6471 galaxies. We
verified that applying more sophisticated redshift cuts (e.g.,
1.0 < z < 1.4 for z<sub>spec</sub> and z<sub>grism</sub>, 0.9 < z < 1.5 for z<sub>phot</sub>) does
not affect our conclusion.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Projected axis ratio *q* distribution of star-forming (top) and quiescent (bottom) galaxies as a function of mass in clusters and the field. The cluster BCGs are marked with a star symbol. The right-hand column shows the *q*-mass relation in clusters and the field as a running median in stellar mass bins of 0.2 dex. Only bins with ≥ 10 galaxies in the cluster and field samples are shown. The solid line and black-dashed line correspond to the median axis ratio for clusters and field, respectively. The shaded regions correspond to the standard error of the median *q* (1.253σ/√N). The dotted lines correspond to the running 16th and 84th percentiles of the *q* distributions. The vertical dotted lines correspond to the masses we used to divide the samples into four mass bins in Section 3.2. The median *q* of both star-forming and quiescent galaxies in clusters and the field both show a mass dependence, with high-mass galaxies being rounder. The median relations of the star-forming galaxies in clusters and the field are plotted in the bottom panel as blue solid (cluster) and gray-dashed (field) for comparison. The axis ratio distributions of quiescent galaxies in clusters and the field are distinct, most remarkably at log(*M**/M⊙) ~ 10.0–10.6 (~4.0σ) and log(*M**/M⊙) ~ 11.0–11.3 (~2.4σ). See Section 3.1 for details.</figcaption></figure>

The right-hand column of Figure 1 shows the rest-frame *U* − *V* and *V* − *J* color distribution of the field sample, color coded by their axis ratios and Sérsic indices. We have checked for offsets in the rest-frame *U* − *V* and *V* − *J* colors between the two catalogs by inspecting the color distribution of the 3D-HST galaxies and the GOGREEN field galaxies. We find no evidence of any color offset that is larger than 0.05 mag. To ensure that our results are robust, we move the *UVJ* selection for the clusters in all four directions by 0.05 mag to mimic the effect of potential color offsets and repeat the analyses four times. All of these analyses give consistent results.

## 3. Results

### 3.1. Projected Axis Ratios in Clusters and the Field

Figure 2 shows the axis ratio of the star-forming and quiescent population in clusters and the field as a function of mass. Clusters show a higher relative abundance of quiescent galaxies compared to star-forming galaxies than the field, with an overall quenched fraction of $f_{Q,clus} = (N_Q/N_{SF} + N_Q) = 0.58$ down to our mass limit compared to the field quenched fraction of $f_{Q,field} = 0.16$. This confirms the enhanced quenched fraction in GOGREEN clusters, relative to the field, found by van der Burg et al. (2020). The right-hand panel of Figure 2 shows the *q*-mass relation in clusters and the field as a running median and percentiles in stellar mass bins of 0.2 dex. We can see that the median *q* of both star-forming galaxies and quiescent galaxies in clusters and the field show a mass dependence. For the cluster sample, the median *q* increases from 0.55 ± 0.02 (0.69 ± 0.04) at log(*M**/M⊙) ~ 9.7 to 0.73 ± 0.10 (0.83 ± 0.02) at log(*M**/M⊙) ~ 11.2 for star-forming (quiescent) galaxies.

The median axis ratios of star-forming and quiescent galaxies in both clusters and the field are significantly different,

which suggests that there are fundamental differences in the intrinsic shapes between star-forming and quiescent galaxies. For the field, these differences in median *q* are seen at all masses. Star-forming galaxies in the field show a lower median *q* at a fixed stellar mass than quiescent galaxies, which is consistent with the finding that most star-forming field galaxies are disks at this redshift range (e.g., van der Wel et al. 2014b). Meanwhile, the median *q* of star-forming and quiescent galaxies in clusters show a difference only at low (log(*M**/M⊙) ≤ 10.1) and high (log(*M**/M⊙) ≥ 11.0) masses.

From Figure 2, we can also see that massive quiescent galaxies with log(*M**/M⊙) ≥ 11 in both clusters and the field are not only rounder than their low-mass counterparts but they also have a narrower *q* distribution, which is reflected by their percentiles (dotted lines). Note that this is not an effect merely due to low number statistics because there are 59 and 116 massive cluster and field galaxies, respectively. A similar change is also seen in the axis ratio distributions of local and intermediate-redshift quiescent galaxies (e.g., van der Wel et al. 2009; Holden et al. 2012; Chang et al. 2013a). We will focus more on their distributions in the next sections.

There are some intriguing differences between the axis ratio distributions in clusters and the field. The medians and percentiles of the *q* distributions of star-forming galaxies are largely consistent with each other, although there may be a weak indication that massive star-forming galaxies (log(*M**/M⊙) ≥ 11.0) in clusters show a higher median (*q* ~1.2σ difference). Meanwhile, we note that the distributions of quiescent galaxies in clusters and the field show differences roughly in the following two mass ranges: log(*M**/M⊙) ~ 10.0–10.6 and log(*M**/M⊙) ~ 11.0–11.3, but in the opposite sense. For log(*M**/M⊙) ~ 10.0–10.6, the median *q* are offset to lower values in clusters compared to the field, with a ~4.0σ difference. This effect can be seen in Figure 2 as both the cluster medians and the 16th percentiles extend to lower values
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison of the projected axis ratio *q* distribution of star-forming (top) and quiescent (bottom) galaxies in clusters with the field in different mass bins. In each case the top panels show the *q* histograms, while the bottom panels show the CDFs and the resultant *p*-values of the KS, AD and MW tests. The shaded areas show the 1σ uncertainty of the CDF derived from the bootstrapped sample of the axis ratios. The solid and dashed arrows show the median of the distribution of the cluster and the field sample, respectively. There is no obvious difference between the *q* distribution of star-forming galaxies in clusters and the field. For quiescent galaxies, the *q* distribution in clusters and the field are distinct, most prominently in the 10.1 ≪ log(*M*<sub>*</sub>) < 10.5 mass bin.</figcaption></figure>

than the field. In contrast, at log(M*/M⊙) ~ 11.0–11.3, there is
evidence that the median *q* and the percentiles are offset to higher
values (~2.4σ difference) in clusters compared to the field. We
find similar differences if we limit the cluster sample to only
spectroscopically-confirmed members but with lower significance,
due to the smaller number of galaxies in the sample. In the
following sections, we will explore these differences in more detail
using the axis ratio distributions in different mass ranges and we
will investigate their implications using intrinsic shape reconstruc-
tion techniques.

## 3.2. Reconstructing the Intrinsic Shapes from the Projected Axis Ratio Distributions

In Figure 3 we compare the axis ratio distributions between
clusters and the field in four mass bins: (i) 9.5 ≪ log(M*) < 10.1,
(ii) 10.1 ≪ log(M*) < 10.5, (iii) 10.5 ≪ log(M*) < 10.8, and
(iv) 10.8 ≪ log(M*) < 11.8. The bottom panels show the
cumulative distribution functions (CDFs). To take into account

the measurement uncertainties of *q*, we bootstrap the observed *q*
distribution and at the same time perturb individual *q* with its 1σ
uncertainty. The shaded areas in the bottom panels show the 1σ
uncertainty of the CDF derived from this bootstrapped sample.
The choice of the binning is selected to match the binning used in
Chang et al. (2013b) for the purposes of the *q* distribution
modeling (see Section 3.2.1 for details). Chang et al. (2013b)
adopted three mass bins, with the lowest mass bin being
10.1 ≪ log(M*/M⊙) < 10.5. Since GOGREEN data allow us
to go down to lower masses, here we include an additional mass
bin of 9.5 ≪ log(M*/M⊙) < 10.1²¹. We have repeated the
analysis with a different set of binnings that have more uniform
bin widths and found a consistent conclusion. We apply the
Kolmogorov–Smirnov (KS) test, the Anderson-Darling (AD)

<footnote n="21">²¹ Similarly, the highest mass bin in Chang et al. (2013b) only goes up to log(M*/M⊙) = 11.5. Here we extend it up to log(M*/M⊙) = 11.8. We checked that limiting it to 11.5 gives the same conclusion.</footnote>
test, and the Mann-Whitney U-test (MW)<fnref n="22" /> on the *q* distributions, with the null hypothesis that they come from a common distribution.

Figure 3 confirms the similarities and differences between clusters and the field discussed in Section 3.1. No obvious differences can be seen between star-forming galaxies in cluster and field, which suggests that their intrinsic shapes are likely to be similar. For quiescent galaxies, we see some evidence that the axis ratio distribution between cluster and field are distinct in some mass bins. Cluster galaxies in the 10.1 ≤ log(M*) < 10.5 bin show a flatter distribution with an apparent excess at low *q* compared to the field. All three tests show a small *p* value (p<sub>KS,AD,MW</sub> ≈ 0.00)<fnref n="23" />. There is some weak indication that cluster galaxies in the highest mass bin (10.8 ≤ log(M*) < 11.8) show on average higher *q* than those in the field (p<sub>MW</sub> ≈ 0.05, p<sub>KS,AD</sub> ≲ 0.1). As we have shown in Section 3.1, this is due to the high-mass population (log(M*) > 11). There is no statistically significant difference in the other mass bins.

### 3.2.1. Methodology for Fitting the Observed *q* Distributions

To understand the implication of these differences, we model the projected axis ratio distributions of the *quiescent population* in clusters and the field to reconstruct their intrinsic shapes. We focus only on quiescent galaxies because no difference can be seen for the star-forming population. We adopt the methodology used by previous works (e.g., Holden et al. 2012; Chang et al. 2013b; van der Wel et al. 2014b), assuming that the intrinsic 3D structure of a galaxy can be described by a triaxial ellipsoid. We refer the reader to Section 5 of Chang et al. (2013b) for a description of the relevant equations. The procedure can be briefly described as follows.

A triaxial ellipsoid can be described with three axes (*a*, *b*, *c*), with $a \ge b \ge c$. One can define two intrinsic axis ratios, $\beta = b/a$ and $\gamma = c/a$, the intrinsic ellipticity $E = (1 - \gamma)$ and the triaxiality $T = (1 - \beta^2)/(1 - \gamma^2)$. The triaxial ellipsoid has two axisymmetric cases; the ellipsoid is known as an oblate spheroid if $\beta = 1$ (i.e., $a = b > c$). If $\beta = \gamma$ (i.e., $a > b = c$), then the ellipsoid is known as a prolate spheroid. The goal of the modeling is to find the model galaxy population(s) (i.e., sets of triaxial ellipsoids) that best-reproduces the observed axis ratio distribution. The model population is assumed to have Gaussian distributions of ellipticity and triaxiality. It can, therefore, be described by four parameters ($E, \sigma_E, T, \sigma_T$), where $\sigma_E$ and $\sigma_T$ are the standard deviations of the ellipticity and triaxiality, respectively.

Assuming random viewing angles, we can compute the expected projected axis ratio distribution for such a population. A correction is then applied to include the effects of uncertainties in the *q* measurements (see Rix & Zaritsky 1995, for a description). In practice, the projected axis ratio distribution for a model galaxy population is computed numerically by generating 100,000 galaxies with random viewing angles and input parameters according to the Gaussian distributions. The number of galaxies is chosen so that the resolution of the axis ratio distribution of the model population is sufficient to compare with the observations. This is

essentially the probability distribution function of the projected axis ratio $P(q)$ of the model population given a set of input parameters.

Star-forming galaxies are traditionally modeled with a single model population of triaxial ellipsoids (e.g., van der Wel et al. 2014b)<fnref n="24" />. For quiescent galaxies, it is established that a single model population is unable to reproduce their axis ratio distributions. For example, Holden et al. (2012) modeled the low-redshift quiescent population in SDSS and found that a single-component triaxial model cannot adequately describe the *q* distribution, except for the massive population with log(M*/M☉) > 11.0. They showed that an additional second component, composed of oblate spheroids, is needed to match the observed distributions. Chang et al. (2013b) confirmed that this is also true for quiescent galaxies in the field at higher redshifts (1 < *z* < 2.5). Although it originated from purely empirical needs to reproduce the axis ratio distribution, this two-component ( triaxial + oblate) model is consistent with the dichotomy in local early-type galaxies discovered via stellar kinematics, (i.e., the slow and fast rotators, see Cappellari 2016, for a review).

Following Holden et al. (2012) and Chang et al. (2013b), on top of the single model population we also adopt the two-component model. Since $\beta = 1$ in the oblate model, the triaxiality is always zero, hence the parameters that describe the model are the intrinsic axis ratio $\gamma$ and its standard deviation $\sigma_\gamma$. To be consistent with previous works, we use $b$ and $\sigma_b$ to denote $\gamma$ and $\sigma_\gamma$<fnref n="25" />. Hence, the two-component model can be fully described by seven independent parameters ($E, \sigma_E, T, \sigma_T, b, \sigma_b, f_{ob}$), where the oblate fraction $f_{ob}$ is the fraction of oblate galaxies relative to the total model population.

Nevertheless, high-redshift galaxy samples, including the cluster and field samples used in this work, are often not large enough to constrain all seven parameters simultaneously. Chang et al. (2013b) tackled this by first fitting the model to local quiescent galaxies, and assumed that the same components could be used to describe the axis ratio distributions at high redshift. They fixed the parameters for the triaxial component and only allowed the oblate parameters ($b, \sigma_b, f_{ob}$) to vary to study the redshift evolution of these parameters. They demonstrated that using this approach can reach conclusion that is consistent with other independent analyses. Here we take a similar approach and build on the findings of Chang et al. (2013b). Partly for this reason, we have adopted a similar mass binning as Chang et al. (2013b). We consider three scenarios with different assumptions, summarized below:

1. Case I—Fitting $E, \sigma_E, T, \sigma_T$—We assume the axis ratio distribution can be described by a single-component model, i.e., $f_{ob} = 0$. The single-component model is useful in studying the distributions at the high masses. See Section 3.2.4.

2. Case II—Fitting $f_{ob}, b, \sigma_b$—We assume the values of the remaining four parameters ($E, \sigma_E, T, \sigma_T$) to be the same as the best-fit values in Chang et al. (2013b). The same assumed values are used for cluster and the field, although we find that the conclusion does not depend

<footnote n="22">²² The Mann-Whitney U-statistics tests the hypothesis that the two sample populations are distributed with the same median.</footnote>

<footnote n="23">²³ We have also assessed the significant of this difference using the half-sample method (using only half of the cluster sample) and by jackknifing the cluster sample. Both tests give small *p* values.</footnote>

<footnote n="24">²⁴ Zhang et al. (2019) demonstrated that $R_e$ needs to be taken into account in the modeling, due to the strong correlation between *q* and $R_e$. This correlation is only seen in star-forming galaxies and is not seen in quiescent galaxies.</footnote>

<footnote n="25">²⁵ This definition is first used by Sandage et al. (1970). An oblate galaxy only has two independent axes *a*, *b*, with the intrinsic axis ratio being *b*/*a*. Sandage et al. (1970) assumes *a* = 1 without loss of generality, hence the use of *b*.</footnote>
<table><caption>**Table 2**
Initial Parameters and Setup Used for the Model Fitting</caption><thead><tr><th rowspan="2">Type</th><th rowspan="2">Mass Ranges</th><th colspan="7">Initial Parameters</th><th colspan="4">Fitting Grid Spacing</th></tr><tr><th>f<sub>ob</sub></th><th>b</th><th>σ<sub>b</sub></th><th>T</th><th>σ<sub>T</sub></th><th>E<sup>a</sup></th><th>σ<sub>E</sub></th><th>Δ<sub>T</sub></th><th>Δ<sub>σ<sub>T</sub></sub></th><th>Δ<sub>E</sub></th><th>Δ<sub>σ<sub>E</sub></sub></th></tr></thead><tbody><tr><td colspan="13" style="text-align: left;">Case I—Fitting E, σ<sub>E</sub>, T, σ<sub>T</sub></td></tr><tr><td>Cluster/Field</td><td>9.5 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>0</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.04</td><td>0.02</td><td>0.01</td><td>0.02</td></tr><tr><td>Cluster/Field</td><td>10.1 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>0</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.04</td><td>0.02</td><td>0.01</td><td>0.02</td></tr><tr><td>Cluster/Field</td><td>10.5 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>0</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.04</td><td>0.02</td><td>0.01</td><td>0.02</td></tr><tr><td>Cluster/Field</td><td>10.8 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>0</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.04</td><td>0.02</td><td>0.01</td><td>0.02</td></tr><tr><td colspan="13" style="text-align: left;">Case II—Fitting f<sub>ob</sub>, b, σ<sub>b</sub></td></tr><tr><td>Cluster/Field</td><td>9.5 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>...</td><td>...</td><td>...</td><td>0.48</td><td>0.08</td><td>0.49</td><td>0.12</td><td>0.04</td><td>0.01</td><td>0.01</td><td>...</td></tr><tr><td>Cluster/Field</td><td>10.1 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>...</td><td>...</td><td>...</td><td>0.48</td><td>0.08</td><td>0.49</td><td>0.12</td><td>0.04</td><td>0.01</td><td>0.01</td><td>...</td></tr><tr><td>Cluster/Field</td><td>10.5 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>...</td><td>...</td><td>...</td><td>0.68</td><td>0.08</td><td>0.45</td><td>0.16</td><td>0.04</td><td>0.01</td><td>0.01</td><td>...</td></tr><tr><td>Cluster/Field</td><td>10.8 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>...</td><td>...</td><td>...</td><td>0.64</td><td>0.08</td><td>0.41</td><td>0.19</td><td>0.04</td><td>0.01</td><td>0.01</td><td>...</td></tr><tr><td colspan="13" style="text-align: left;">Case III—Fitting f<sub>ob</sub> only</td></tr><tr><td>Cluster/Field</td><td>9.5 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>...</td><td>0.28</td><td>0.09</td><td>0.48</td><td>0.08</td><td>0.49</td><td>0.12</td><td>0.04</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Cluster/Field</td><td>10.1 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>...</td><td>0.28</td><td>0.09</td><td>0.48</td><td>0.08</td><td>0.49</td><td>0.12</td><td>0.04</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Cluster/Field</td><td>10.5 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>...</td><td>0.28</td><td>0.08</td><td>0.68</td><td>0.08</td><td>0.45</td><td>0.16</td><td>0.04</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Cluster/Field</td><td>10.8 &le; log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>...</td><td>0.29</td><td>0.07</td><td>0.64</td><td>0.08</td><td>0.41</td><td>0.19</td><td>0.04</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

heavily on these assumed values (see Appendix D for a discussion).<fnref n="26" />

3. Case III—Fitting $f_{\text{ob}}$ only—We assume the values of the remaining six parameters to be the same as the best-fit values in Chang et al. (2013b).

The best-fit model population is determined using a maximum likelihood estimation method. To reduce computation time, we first generate a model grid within the parameter space being considered in each case. The spacing of the grid and the assumed values for the remaining parameters can be found in Table 2. We then compute the likelihood for each model population. For each *q* in the observed distribution, we calculate the probability of observing a galaxy with this particular value of *q* according to the probability distribution function $p(q)$ of the model population. The total log-likelihood $\ln(L)$ of the model is then computed by summing the log-probability of all *q* in the observed axis ratio distributions. Then, the best-fit model population is taken as the one with the highest likelihood.

The uncertainty of the fitted parameters is derived by fitting the bootstrapped sample. The $1\sigma$ variation of the best-fit parameters of the bootstrapped sample is taken as the uncertainty. Some examples of the corner plots of the fitting of the bootstrapped sample are shown in Appendix D. The best-fit parameters can be found in Table 3. We provide $p_{\text{KS}}$, $p_{\text{MW}}$, the reduced $\chi^2$ and the $p_{\chi^2}$ values in Table 3 as rough goodness-of-fit indicators. The Akaike information criterion (AIC) value of the best-fit models is also provided in Table 3.

### 3.2.2. Evidence for a Higher Fraction of Oblate Quiescent Galaxies in Clusters

Figure 4 shows the result of the modeling for the four mass bins. The shaded regions show the $1\sigma$ variation of the axis ratio

distribution derived from the bootstrapped sample, which we used to derive the uncertainty of the best-fit parameters.

We find that the observed axis ratio distribution of the field sample in all mass bins can be reasonably described by a single-component triaxial model (Case I). This is consistent with previous findings. Chang et al. (2013b) also found that the single-component model cannot be ruled out with only the use of high-redshift data. Meanwhile, there is some evidence that the single-component model is unable to accurately reproduce the shape of the distribution of the cluster sample in some of the mass bins (see Table 3 for the $p_{\chi^2}$ values). This is intriguing because the axis ratio distribution of local quiescent galaxies in this mass range is known to be not well described by a single component, due to the existence of an oblate component (e.g., Holden et al. 2012).

Assuming the triaxial component as found in Chang et al. (2013b) and fitting the oblate component parameters $f_{\text{ob}}$, $b$, $\sigma_b$ (Case II), we find that the main differences between the cluster and field distribution lie in the fraction of oblate galaxies in the total population. We find tentative evidence that the cluster distribution in the three lower mass bins has a higher oblate fraction than the field, with the $10.1 \le \log(M_*/M_\odot) < 10.5$ mass bin showing the largest difference ($f_{\text{ob,cluster}} = 0.88 \pm 0.14$ versus $f_{\text{ob,field}} = 0.44 \pm 0.17$). The best-fit value of the intrinsic axis ratio $b$ of the cluster sample is low in all four mass bins, with $b = 0.25-0.35$, which is consistent with the value found in Chang et al. (2013b) for their $1 < z < 2.5$ sample ($b \sim 0.29$). The best-fit $b$ for the field models are also consistent with this value, except in the two lower mass bins where both $b$ and $\sigma_b$ are poorly constrained (see Appendix D for the corner plots).

Fixing all parameters to the field values and only allowing $f_{\text{ob}}$ to vary also gives a similar conclusion (Case III). The cluster sample has a much higher oblate fraction of $f_{\text{ob,cluster}} = 0.72 \pm 0.06$ compared to the field $f_{\text{ob,field}} = 0.24 \pm 0.07$ in the mass bin $10.1 \le \log(M_*/M_\odot) < 10.5$. However, we note that the best-fit model of the field in this bin is not a good

<footnote n="26">²⁶ For the lowest mass bin, we use the same assumed values as the $10.1 \le \log(M_*/M_\odot) < 10.5$ bin in Chang et al. (2013b).</footnote>
<table><caption>Table 3
Best-fit Parameters of the Models</caption><thead><tr><th rowspan="2">Type</th><th rowspan="2">Mass Ranges</th><th colspan="6">Best-fit Parameters</th><th colspan="5">GoF</th></tr><tr><th>f<sub>ob</sub></th><th>b</th><th>σ<sub>b</sub></th><th>T</th><th>σ<sub>T</sub></th><th>E<sup>a</sup></th><th>σ<sub>E</sub></th><th>p<sub>KS</sub></th><th>p<sub>MW</sub></th><th>p<sub>χ<sup>2</sup></sub></th><th>χ<sup>2</sup>/ν</th><th>AIC</th></tr></thead><tbody><tr><td colspan="14" style="text-align: center; font-weight: bold;">Case I—Fitting E, σ<sub>E</sub>, T, σ<sub>T</sub></td></tr><tr><td>Cluster</td><td>9.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>...</td><td>...</td><td>...</td><td>0.36<sup>+0.52</sup><sub>-0.00</sub></td><td>0.00<sup>+0.14</sup><sub>-0.00</sub></td><td>0.63<sup>+0.04</sup><sub>-0.17</sub></td><td>0.06<sup>+0.24</sup><sub>-0.06</sub></td><td>0.03</td><td>0.06</td><td>0.02</td><td>2.13</td><td>897.6</td></tr><tr><td>Cluster</td><td>10.1 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>...</td><td>...</td><td>...</td><td>0.88<sup>+0.08</sup><sub>-0.52</sub></td><td>0.00<sup>+0.22</sup><sub>-0.00</sub></td><td>0.53<sup>+0.13</sup><sub>-0.04</sub></td><td>0.22<sup>+0.10</sup><sub>-0.12</sub></td><td>0.63</td><td>0.49</td><td>0.02</td><td>1.99</td><td>1233.27</td></tr><tr><td>Cluster</td><td>10.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>...</td><td>...</td><td>...</td><td>0.36<sup>+0.08</sup><sub>-0.00</sub></td><td>0.10<sup>+0.06</sup><sub>-0.10</sub></td><td>0.72<sup>+0.02</sup><sub>-0.02</sub></td><td>0.02<sup>+0.24</sup><sub>-0.02</sub></td><td>0.12</td><td>0.04</td><td>0.76</td><td>0.66</td><td>978.7</td></tr><tr><td>Cluster</td><td>10.8 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>...</td><td>...</td><td>...</td><td>0.36<sup>+0.04</sup><sub>-0.00</sub></td><td>0.02<sup>+0.02</sup><sub>-0.02</sub></td><td>0.48<sup>+0.05</sup><sub>-0.03</sub></td><td>0.22<sup>+0.04</sup><sub>-0.04</sub></td><td>0.80</td><td>0.34</td><td>0.18</td><td>1.38</td><td>1076.6</td></tr><tr><td>Field</td><td>9.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>...</td><td>...</td><td>...</td><td>0.48<sup>+0.44</sup><sub>-0.12</sub></td><td>0.02<sup>+0.38</sup><sub>-0.02</sub></td><td>0.64<sup>+0.05</sup><sub>-0.17</sub></td><td>0.30<sup>+0.06</sup><sub>-0.06</sub></td><td>0.75</td><td>0.34</td><td>0.89</td><td>0.57</td><td>1757.0</td></tr><tr><td>Field</td><td>10.1 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>...</td><td>...</td><td>...</td><td>0.44<sup>+0.24</sup><sub>-0.08</sub></td><td>0.04<sup>+0.34</sup><sub>-0.04</sub></td><td>0.50<sup>+0.02</sup><sub>-0.04</sub></td><td>0.20<sup>+0.02</sup><sub>-0.02</sub></td><td>0.89</td><td>0.47</td><td>0.98</td><td>0.35</td><td>2433.0</td></tr><tr><td>Field</td><td>10.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>...</td><td>...</td><td>...</td><td>0.92<sup>+0.00</sup><sub>-0.32</sub></td><td>0.00<sup>+0.12</sup><sub>-0.00</sub></td><td>0.51<sup>+0.07</sup><sub>-0.04</sub></td><td>0.24<sup>+0.04</sup><sub>-0.02</sub></td><td>0.32</td><td>0.16</td><td>0.39</td><td>1.06</td><td>2464.3</td></tr><tr><td>Field</td><td>10.8 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>...</td><td>...</td><td>...</td><td>0.48<sup>+0.04</sup><sub>-0.08</sub></td><td>0.02<sup>+0.02</sup><sub>-0.02</sub></td><td>0.53<sup>+0.03</sup><sub>-0.03</sub></td><td>0.20<sup>+0.04</sup><sub>-0.02</sub></td><td>0.31</td><td>0.10</td><td>0.06</td><td>1.67</td><td>2142.4</td></tr><tr><td colspan="14" style="text-align: center; font-weight: bold;">Case II—Fitting f<sub>ob</sub>, b, σ<sub>b</sub></td></tr><tr><td>Cluster</td><td>9.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>0.76<sup>+0.16</sup><sub>-0.16</sub></td><td>0.35<sup>+0.02</sup><sub>-0.02</sub></td><td>0.03<sup>+0.03</sup><sub>-0.03</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.31</td><td>0.20</td><td>0.18</td><td>1.37</td><td>894.4</td></tr><tr><td>Cluster</td><td>10.1 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>0.88<sup>+0.12</sup><sub>-0.20</sub></td><td>0.32<sup>+0.04</sup><sub>-0.03</sub></td><td>0.09<sup>+0.03</sup><sub>-0.03</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.88</td><td>0.41</td><td>0.06</td><td>1.66</td><td>1230.5</td></tr><tr><td>Cluster</td><td>10.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>0.76<sup>+0.16</sup><sub>-0.08</sub></td><td>0.28<sup>+0.02</sup><sub>-0.01</sub></td><td>0.02<sup>+0.02</sup><sub>-0.02</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.78</td><td>0.40</td><td>0.93</td><td>0.46</td><td>971.1</td></tr><tr><td>Cluster</td><td>10.8 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>0.24<sup>+0.08</sup><sub>-0.12</sub></td><td>0.25<sup>+0.06</sup><sub>-0.03</sub></td><td>0.01<sup>+0.03</sup><sub>-0.01</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.22</td><td>0.08</td><td>0.35</td><td>1.11</td><td>1074.9</td></tr><tr><td>Field</td><td>9.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>0.64<sup>+0.16</sup><sub>-0.12</sub></td><td>0.24<sup>+0.08</sup><sub>-0.08</sub></td><td>0.21<sup>+0.08</sup><sub>-0.10</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.85</td><td>0.37</td><td>0.90</td><td>0.57</td><td>1755.0</td></tr><tr><td>Field</td><td>10.1 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>0.44<sup>+0.16</sup><sub>-0.12</sub></td><td>0.45<sup>+0.05</sup><sub>-0.11</sub></td><td>0.22<sup>+0.07</sup><sub>-0.08</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.39</td><td>0.16</td><td>0.82</td><td>0.64</td><td>2434.0</td></tr><tr><td>Field</td><td>10.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>0.64<sup>+0.12</sup><sub>-0.08</sub></td><td>0.29<sup>+0.03</sup><sub>-0.03</sub></td><td>0.08<sup>+0.02</sup><sub>-0.03</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.94</td><td>0.46</td><td>0.64</td><td>0.83</td><td>2458.6</td></tr><tr><td>Field</td><td>10.8 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>0.24<sup>+0.08</sup><sub>-0.08</sub></td><td>0.27<sup>+0.05</sup><sub>-0.07</sub></td><td>0.05<sup>+0.03</sup><sub>-0.05</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.12</td><td>0.25</td><td>0.06</td><td>1.70</td><td>2143.7</td></tr><tr><td colspan="14" style="text-align: center; font-weight: bold;">Case III—Fitting f<sub>ob</sub> only</td></tr><tr><td>Cluster</td><td>9.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>0.64<sup>+0.12</sup><sub>-0.16</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.28</td><td>0.15</td><td>0.38</td><td>1.07</td><td>893.6</td></tr><tr><td>Cluster</td><td>10.1 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>0.72<sup>+0.12</sup><sub>-0.08</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.47</td><td>0.50</td><td>0.07</td><td>1.57</td><td>1228.4</td></tr><tr><td>Cluster</td><td>10.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>0.76<sup>+0.16</sup><sub>-0.08</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.50</td><td>0.34</td><td>0.89</td><td>0.55</td><td>972.0</td></tr><tr><td>Cluster</td><td>10.8 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>0.20<sup>+0.12</sup><sub>-0.08</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.39</td><td>0.17</td><td>0.43</td><td>1.02</td><td>1072.0</td></tr><tr><td>Field</td><td>9.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.1</td><td>0.60<sup>+0.04</sup><sub>-0.12</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.63</td><td>0.25</td><td>0.70</td><td>0.79</td><td>1756.3</td></tr><tr><td>Field</td><td>10.1 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.5</td><td>0.24<sup>+0.08</sup><sub>-0.08</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.02</td><td>0.00</td><td>0.36</td><td>1.09</td><td>2435.4</td></tr><tr><td>Field</td><td>10.5 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 10.8</td><td>0.64<sup>+0.08</sup><sub>-0.08</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.96</td><td>0.46</td><td>0.73</td><td>0.76</td><td>2454.7</td></tr><tr><td>Field</td><td>10.8 ≤ log(M<sub>*</sub>/M<sub>☉</sub>) &lt; 11.8</td><td>0.24<sup>+0.08</sup><sub>-0.08</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.10</td><td>0.20</td><td>0.09</td><td>1.54</td><td>2140.3</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The best-fit intrinsic shape models for the quiescent galaxy population in clusters (top) and the field (bottom) in different mass bins. The solid (dashed) line in each panel corresponds to the observed axis ratio distribution for the clusters (field), while the shaded area represents the $1\sigma$ variation of the distribution derived from the bootstrapped sample. A different binning width than in Figure 3 is used here to highlight the variation. The purple line in each panel corresponds to the best-fit model of case I, where a single-component model is used. The brown line corresponds to the best-fit model of case II, where $f_{\text{ob}}$, $b$, $\sigma_b$ are free parameters. The orange line corresponds to the best-fit model of case III, where only $f_{\text{ob}}$ is being fitted. The best-fit case II and III models for the $10.1 \le \log(M_*) < 10.5$ mass bin both show a higher $f_{\text{ob}}$ in clusters than the field, which suggests that clusters have a higher fraction of disk-like galaxies than the field at this redshift.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The axis ratio distributions for each of the oblate and triaxial components that contribute to the total distribution in the best-fit intrinsic shape models (Case II) for the quiescent galaxy population in clusters (top) and the field (bottom) in different mass bins. The brown line in each panel corresponds to the best-fit model of case II, also shown in Figure 4. The magenta (blue) dashed line corresponds to the model axis ratio distributions of the oblate ( triaxial) component. The shaded magenta area represents the $1\sigma$ variation of the fitted oblate component parameters derived from the bootstrapped sample. For simplicity, only the variation of the oblate component is shown. A different y-scale than in Figure 4 is used to show the components. With a low best-fit *b* value, the oblate component dominates the low-*q* region of the cluster distribution.</figcaption></figure>

representation of the data according to the $p_{\text{KS}}$ and $p_{\text{MW}}$ values, presumably due to the limitation of the assumed models.

Figure 5 shows the axis ratio distributions for each of the oblate and triaxial components separately, which contribute to the total distribution in the best-fit $f_{\text{ob}}$, $b$, $\sigma_b$ models (Case II). In each panel, the brown line corresponds to the same best-fit model in Figure 4, and the shaded magenta area represents the $1\sigma$ variation of the fitted oblate component parameters derived from the bootstrapped sample. We can see that the effect of having a larger $f_{\text{ob}}$ of a low *b* oblate population to the overall distribution. This results in a larger low-*q* contribution relative to the total population and gives rise to a broader *q* distribution that better describes the flatter shape of the cluster distributions, especially in the mass bin $10.1 \le \log(M_*/M_{\odot}) < 10.5$.

Therefore, the difference between the cluster and the field sample that we described in Section 3.1 is consistent with the existence of a larger population of flattened, oblate quiescent galaxies in clusters.

### 3.2.3. The Relationship between Axis Ratio and Sérsic Index

We show that the main difference between the axis ratio distribution of the cluster and the field quiescent population is consistent with the existence of a larger population of flattened, oblate quiescent galaxies in clusters. With such a low *b* value, these oblate quiescent galaxies resemble disk-dominated galaxies as in the majority of the star-forming population. In the literature, the Sérsic index is commonly used to separate
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The Sérsic index of the quiescent galaxies in the cluster sample as a function of axis ratio. Grey points correspond to individual galaxies. Large red circles and error bars correspond to the median *n* of the entire sample and its standard error. Light blue and orange circles correspond to the median *n* in the lowest and the second lowest mass bin, respectively. There is a positive trend between *q* and *n*, as expected. The trend is mass dependent, with low-mass galaxies having on average lower values of *n* for a given *q*.</figcaption></figure>

disk- and bulge-dominated galaxies with the *n* = 2.5 division. Although we do not rely heavily on Sérsic index in this work, here we examine the relation between axis ratio and Sérsic index of the cluster sample as a consistency check.

Figure 6 shows the axis ratio of the quiescent galaxies in the cluster sample as a function of their Sérsic index. We see a strong trend between the two, with low *q* galaxies typically having on average smaller values of *n*. The trend is also mass-dependent, with low-mass galaxies having lower values of *n* for a given *q* bin. The excess population of *q* ~ 0.4 galaxies in clusters have, on average, *n* ~ 2.5, and are consistent with the traditional selection of disky galaxies using Sérsic indices. This strong trend supports our modeling results and interpretation that the difference in the axis ratio distribution originates from a larger population of flattened, disk-like galaxies in clusters compared to the field.

The fact that we see an excess population of disk-like quiescent galaxies in clusters suggests that morphological transformation and quenching in clusters does not operate in the same way as in the field. In Section 4, we discuss this further and explore possible implications together with other quantities.

### 3.2.4. Properties of the Massive Quiescent Galaxies in $10.8 \le \log(M_*/M_\odot) < 11.8$

In Section 3.1 we showed that massive quiescent galaxies $\log(M_*/M_\odot) \ge 11$ in both clusters and the field are not only rounder than their low-mass counterparts but they also have a narrower *q* distribution; i.e., there is a lack of low-*q* galaxies. There is also evidence that the median *q* is offset to higher values in clusters compared to the field. Here we explore the

underlying reason by examining the fitting results of the highest mass bin $10.8 \le \log(M_*/M_\odot) < 11.8$.<fnref n="27" />

There is a stark difference between the best-fit models of the highest mass bin and those of the three lower mass bins, as one can see in Figure 5. To reconstruct the shape of the axis ratio distribution of these massive galaxies, the best-fit case II (and also case III) models have much lower $f_{ob}$ values compared to the three lower mass bins (~2.3 (1.7)$\sigma$ for the cluster (field) sample). In fact, for both the cluster and the field sample, the $f_{ob}$s in the highest mass bin are statistically consistent with 0. This suggests that the contribution of a possible second oblate component is small. Indeed, the single-component triaxial model (case I) can describe the axis ratio distribution of the massive galaxies in both the cluster and field sample well.

In Section 3.2, we note that there is some weak indication that clusters have higher median *q* compared to the field. From the best-fit parameters, the best-fit case I model of the cluster has a trixiality of $T=0.36$ and an ellipticity of $E=0.48$, compared to the field value of $T=0.48$ and $E=0.53$. Nevertheless, they are within $1\sigma$ uncertainty. The best-fit $\sigma_T$ and $\sigma_E$ are consistent with each other in clusters and the field. Therefore, with our sample, there is no evidence of a difference in the intrinsic shape distribution of massive galaxies in clusters and in the field.

### 3.3. Caveats

Dust obscuration, in particular dust lanes or dust gradients within the galaxy, can impact the axis ratio measurements and thus potentially affect the axis ratio distributions. It is also possible that the dust obscuration hides the disk structure of the galaxies, making them appear rounder (van der Wel et al. 2014a). This has a larger effect on star-forming galaxies than quiescent galaxies due to their higher dust content. Previous studies have shown that the dust content of quiescent galaxies, probed via the rest-frame $V-J$ color, is low and does not correlate strongly with their axis ratio (Chang et al. 2013b). Here we check if this is also the case in clusters at this redshift range.

Figure 7 shows the dust extinction $A_V$ of star-forming and quiescent galaxies in the cluster sample, derived via SED fitting with FAST in steps of 0.2, as a function of *q*. We confirm that the quiescent galaxies in our cluster sample have, in general, a low dust content. There is no obvious correlation between the axis ratio and the dust extinction. Meanwhile, the star-forming population have a larger variation in $A_V$ at a certain *q*. Galaxies that have lower *q* show higher dust extinction values on average, which is expected for an inclined / edge-on disk population. A similar correlation can also be seen if the $V-J$ color is used instead. We have also checked that the field sample shows similar correlations as the cluster sample. This suggests that the difference between the cluster and field axis ratios that we see in the quiescent population is unlikely to be driven only by dust, but is due to the difference in their intrinsic shapes.

In addition, one potential caveat is related to the fact that our axis ratio measurements are measured from galaxy surface brightness profiles. Luminosity-weighted structural measurements are not always a reliable measure for the mass

<footnote n="27">²⁷ We note that our mass bin starts from $\log(M_*/M_\odot) = 10.8$. We have also repeated the fits using a mass bin of $11.0 \le \log(M_*/M_\odot) < 11.8$ (59 galaxies) and found consistent results.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The dust extinction $A_V$ of the star-forming (top) and quiescent (bottom) galaxies in the cluster sample as a function of $q$. Orange circles and error bars correspond to the mean $A_V$ and its standard error. The green shading reflects the number of galaxies in a particular bin. The quiescent galaxies in our sample have low dust content. Meanwhile, the dust extinction of the star-forming galaxies shows a negative correlation with $q$.</figcaption></figure>

distribution of a galaxy, due to radial variation in stellar population properties such as age and metallicity. Suess et al. (2019) derived mass-weighted structural parameters for galaxies in the field at $1.0 < z < 2.5$ and found that most galaxies, no matter star-forming or quiescent, show negative color/mass-to-light ratio (M/L) gradients with a strong redshift evolution, resulting in smaller mass-weighted sizes than the luminosity weighted ones. Although their work did not focus on axis ratios, they found that these gradients can account for most of the evolution in the mass-size relation. Similar effects have been observed in clusters at high redshift. Chan et al. (2018) showed that there are strong negative color gradients in quiescent galaxies in three clusters at $z \sim 1.5$ that are consistent with being a combination of age and metallicity gradients. These gradients are even stronger in evolved clusters compared

to the field at a similar redshift. The existence of these negative M/L gradients suggests that the older stellar population in the bulge (which has a larger M/L) could be outshone by younger and bright stars in disks, which can bias the luminosity-weighted axis ratios to lower values.

Our HST data set does not have the necessary imaging bands to derive mass-weighted structural parameters (see, e.g., Chan et al. 2016; Suess et al. 2019, for a discussion of the methodologies and their requirements), therefore we are unable to examine the effects of stellar population gradients on the axis ratio distributions. Instead, we can use the results in Chan et al. (2018) as a reference. In addition to size, they compared both the luminosity-weighted and mass-weighted axis ratio distributions between clusters and the field at $z \sim 1.5$.<fnref n="28" /> They found that the mass-weighted axis ratios are on average smaller, with a broader mass-weighted distribution extending to lower $q$ than the luminosity-weighted $q$ distribution (see Section 6.1.3 in Chan et al. 2018, for more details). This implies that the low $q$ excess and flat $q$ distributions that we see in these clusters are likely to be from genuine oblate structures instead of hidden bulges. Nevertheless, fully ruling out the possibility that the difference we see is caused by stellar population gradients would require a detailed analysis of the mass-weighted properties and color gradients in both clusters and in the field.

Another potential caveat is that the axis ratio distribution modeling relies heavily on the assumption that the galaxy population is observed from random viewing angles. This assumption breaks down when the galaxies are intrinsically aligned with respect to a certain direction. There have been reports that cluster galaxies may be aligned radially toward the center of the cluster in low-redshift clusters (e.g., Huang et al. 2018; Georgiou et al. 2019). We therefore examined the radial alignment signal in our sample and measured the alignment of the galaxies with respect to the center of the cluster. The procedure and result are discussed in Appendix E.

We find that the average radial and tangential alignments for the cluster sample within $1R_{200}$ are consistent with zero. Examining the alignment signal as a function of cluster-centric radius, there is weak evidence that the average radial alignment is positive in the region close to the cluster center ($\sim 1.4\sigma$ for $R < 0.2R_{200}$). We conclude that the potential intrinsic alignments in the sample are unlikely to affect our results.

## 4. Discussion

The goal of this work is to examine the effect of environment on galaxy structural properties. We find that the axis ratio distributions of quiescent galaxies in clusters and the field are distinct. By modeling the axis ratio distribution in different mass bins, we find evidence that quiescent galaxies in clusters have a higher fraction of flattened, oblate galaxies than the field in the intermediate-mass range. The most massive cluster galaxies, with $10.8 \le \log(M_*/M_\odot) < 11.8$, have a low fraction of oblate galaxies, and those in clusters exhibit a lower ellipticity than the field. Here we discuss the implications of these results. We begin with the result of the massive galaxies in Section 4.1. We then discuss the result of the intermediate-mass range in Section 4.2 in the context of a simple toy accretion model. In Section 4.3, we explore the implication of

<footnote n="28">²⁸ Chan et al. (2018) also found evidence that the $q$ distribution in the evolved clusters show an excess of low axis ratio ($q < 0.4$) compared to the field, albeit with low number statistics.</footnote>
our results in the context of the “early mass-quenching” scenario, discussed in van der Burg et al. (2020). In Section 4.4, we combine our results with the measured stellar age of a subset of the population to explore the underlying physical mechanism that drives environmental quenching.

## 4.1. Evolution of the Massive Quiescent Galaxies in Clusters

The fact that massive galaxies, in both clusters and the field, have significantly different axis ratio distributions than their low-mass counterparts has been observed in previous studies at lower redshifts. For example, van der Wel et al. (2009) showed that there is a lack of low-q galaxies ($q<0.6$) in the local massive quiescent population at $log(M_*/M_\odot) \ge 11.0$. Holden et al. (2012) reported a similar transition exists at $z \sim 0.7$. Similar effects has been found in clusters locally and at intermediate redshifts (e.g., Vulcani et al. 2011). We show that this is also true at $1.0 < z < 1.4$.

These results are often interpreted as evidence that massive quiescent galaxies experienced repeated major and minor mergers, which make them appear rounder gradually. In simulations, it is shown that the importance of mergers increases as a function of mass (e.g., Wang & Kauffmann 2008; De Lucia et al. 2010; Qu et al. 2017). This picture is also largely supported by the observed kinematics of these low-z massive quiescent galaxies (i.e., the slow rotators, e.g., Emsellem et al. 2011), their observed number density evolution, and the merger rates (e.g., Man et al. 2016). The offset to higher *q* values in the cluster population may therefore be a consequence of massive galaxies in clusters having experienced more mergers than the field at the epoch of observation.

This can either be due to a) merger rates in clusters are (or were) elevated compared to the field, or b) they are in a more advanced evolutionary stage compared to their field counterparts. Observations have shown that major mergers can play an important role in mass assembly of the brightest cluster galaxies at $z \sim 1$ (e.g., Lidman et al. 2013), although this is less likely for satellite galaxies. Nevertheless, we note that the massive cluster population that we see is likely to be composed of galaxies that were central galaxies for most of their lifetime (e.g., De Lucia et al. 2012), hence the merger could have happened before or when the galaxy was being accreted. Meanwhile, a recent stellar kinematics study of massive galaxies ($log(M_*/M_\odot) \ge 11.0$) at $0.6 < z < 1.0$ found that only those in the densest environments are primarily slow rotators (Cole et al. 2020). They suggest that slow rotators are being built in dense environments first through repeated minor mergers and hence they are more kinematically evolved compared to the field. We note that this relation, which is known as the kinematic morphology-density relation (Cappellari et al. 2011), has been studied extensively in dense environments in the local Universe. Several studies showed that the fraction of massive slow rotators increases with galaxy number density, and the locations of the slow rotators are strongly correlated with peak densities in groups and clusters (e.g., D'Eugenio et al. 2013; Houghton et al. 2013; Graham et al. 2019), although some studies suggesting that this relation is driven by the stellar mass distribution with the environment and not the environment itself (e.g Brough et al. 2017; Veale et al. 2017). Although we are not able to distinguish the effects of major and minor mergers, our result supports the picture that

mergers are a crucial component in the evolution of the massive galaxies in clusters at $1.0 < z < 1.4$.

## 4.2. The Effect of Environmental Quenching on Galaxy Structure

In this section, we aim to combine the results of the axis ratio distribution modeling with the quenched fractions in our samples to quantify the extent of morphological transformation. It is essential to explore the relationship between environmental quenching and morphological transformation at this redshift because the dominant environmental mechanism(s) needs to be able to explain the morphological mix or the morphological signatures of the population. Here we consider a simple accretion model and compute the expected axis ratio distributions under various assumptions. We then compare these distributions with the observed *q* distribution of the quiescent galaxies in clusters.

We start by computing the quenched fraction $f_Q$ in the cluster and field samples for the four mass bins, following the method of van der Burg et al. (2020). Note that the field $f_Q$ value here is different from their work because the adopted field sample comes from CANDELS/3D-HST as opposed to UltraVISTA (Muzzin et al. 2013b) in van der Burg et al. (2020). Over the whole mass range that we considered in this work ($9.5 \le \log(M_*/M_\odot) < 11.8$), the quenched fraction in clusters is more than three times higher than in the field ($f_{Q,clus} = 0.58$ versus $f_{Q,field} = 0.16$).

Assuming that mass quenching occurs in the same way in clusters as in the field, we can compute the fraction of the excess quenched galaxies (i.e., quenched via environmental processes) in the cluster quiescent sample in each mass bin<fnref n="29" />: $f_{EQ} = (f_{Q,clus} - f_{Q,field})/f_{Q,clus}$. The quantity $f_{EQ}$ varies in different mass bins, with the lowest mass bins having the largest value.

We first investigate the effect of having such an excess population on the field quiescent axis ratio distributions in each mass bin. To do this, we compute the expected axis ratio distributions of a galaxy population with a fraction $f_{EQ}$ of “accreted galaxies”. Star-forming galaxies, which are randomly drawn from the field star-forming *q* distributions ($P_{SF}(q)$) of the corresponding mass bin, are added into the field quiescent *q* distributions until the fraction of this accreted population reaches $f_{EQ}$. Using the star-forming *q* distributions from the field population as the parent distribution mimics the effect of having no morphological transformation, in the sense that the accreted galaxies retain the same morphology (axis ratio) as they would have had in the field.

The top row of Figure 8 shows the result of this accretion model in different mass bins. The black line corresponds to the distribution with $f_{EQ}$ of star-forming population mixed in, while the gray lines correspond to the $1\sigma$ variation of the expected distribution derived from bootstrapping, in which the bootstrapped samples contain the same number of galaxies as the cluster quiescent sample. Comparing with the observed *q* distribution in the clusters, we find that the expected distribution of the accretion model matches the overall shape of the cluster distribution in the two intermediate-mass bins

<footnote n="29">²⁹ Note that $f_{EQ}$ describes the fraction of environmentally quenched galaxies in the quiescent galaxy population. This is different from the Quenched Fraction Excess (QFE) in van der Burg et al. (2020), which describes the fraction of galaxies that would have been star-forming in the field but are quenched by the environment.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The expected axis ratio distribution of the quenched galaxies in clusters from a simple accretion model. Top: Comparison of the observed axis ratio distributions in clusters (red line) to the expected axis ratio distribution of the model with a $f_{EQ}$ fraction of star-forming population ($P_{SF}(q)$) (black line) in different mass bins. Grey dotted-dashed lines correspond to the $1\sigma$ variation of the expected distribution. Middle: The cumulative distribution function of the expected axis ratio distribution of models with different accreted star-forming fractions. Bottom: The axis ratio distribution of the models with different accreted star-forming fractions. See Section 4.2 for details.</figcaption></figure>

$10.1 \le \log(M_*/M_\odot) < 10.5$ and $10.5 \le \log(M_*/M_\odot) < 10.8$, with $p_{KS}$ values of $\simeq 0.56$ and 0.39, respectively. However, the model cannot match the shape of the highest and lowest mass bins, with both bins having $p_{KS}$ value of $\simeq 0$.

The middle and bottom rows of Figure 8 illustrate the effect of varying the accreted star-forming fraction to the axis ratio distribution. In general, increasing the fraction of the accreted population increases the abundance of the low *q* galaxies, resulting in a broader *q* distribution. Varying this fraction can also be regarded as changing the amount of morphological transformation of the accreted galaxies. Since we have an independent constraint on the fraction of the environmental quenched population, if an accreted star-forming fraction that is smaller than $f_{EQ}$ matches the data well, then it suggests that part of the accreted population has transformed, such that their distribution matches closer to the field quiescent population than the star-forming ones.

Interestingly, we find that an accreted fraction that is consistent with $f_{EQ}$ best-fits the cluster data in the two intermediate-mass bins, as seen from the cumulative distributions. For $10.1 \le \log(M_*/M_\odot) < 10.5$, the KS test indicates that models with $f_{EQ} \pm 0.2$ are also acceptable representations. Hence, for the full $10.1 \le \log(M_*/M_\odot) < 10.8$ range, our model suggests that the observed axis ratio distributions are statistically consistent with a scenario where no morphological transformation occurs after the galaxy was accreted and quenched environmentally.

For the lowest and highest mass bins, it is intriguing to see that none of the model distributions are a good representation of the data. For the highest mass bin, since the cluster distributions are shifted to higher *q* relative to the field, it is

unsurprising that our model will not work. As we discussed in Section 4.1, mergers are likely to be a crucial component in the evolution of these massive galaxies which changes their axis ratio distributions. However, although the lowest mass bin has a high $f_{EQ}$, our model fails to reproduce the observed distribution. The model may be too simplistic to reproduce the characteristics of the observed distribution, particularly in the low-*q* region ($q < 0.3$). It is also possible that our sample size for the lowest-mass bin is simply too small.

We assumed that the accreted population has the same axis ratio distribution as the star-forming population in the field at the same epoch. However, this might not be true if the star-forming population is accreted earlier. We repeat the analysis using the axis ratio distribution of two higher redshift samples of star-forming galaxies in the field at $1.5 < z < 2.0$ and $2.0 < z < 2.5$ (1.4 and 2.2 Gyr earlier). The two samples have 1189 and 995 galaxies, respectively. We find that the results remain unchanged, primarily due to the fact that their axis ratio distributions are consistent with the distribution at GOGREEN redshifts ($p_{KS} \simeq 0.3$ and 0.2, respectively).

### 4.3. Morphological Transformation in Context of the Early Mass-quenching Scenario

It has been demonstrated by van der Burg et al. (2020) that the shape of the stellar mass function (SMF) in star-forming and quiescent galaxies is indistinguishable between cluster and field at $1.0 < z < 1.4$. This leads to the attractive explanation that galaxies in clusters quench through the same mass-quenching process as those in the field, but at an earlier time and at a higher rate. Under this “early mass-quenching”
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The mass-weighted ages of the quiescent galaxies in clusters (in unit of cosmic time $t_{\text{univ,z}} - t_{\text{mw}}$) as a function of *n* and *q* for the four mass bins. Only galaxies that have both reliable stellar age and structural parameter measurements are shown. There is no clear trend in age with *n* or *q*. Disk-like galaxies (e.g., *n* < 2.5 or *q* < 0.4) show a large variation in age.</figcaption></figure>

scenario, they found that a difference in formation time of $\ge 1$ Gyr can result in a quenched fraction difference that is consistent with the data. This difference in formation time would manifest itself in the age difference between the cluster and field quiescent population (see the discussion in Webb et al. 2020, for more details on the relationship between the two quantities). Although cluster galaxies are on average older than the field (Webb et al. 2020), the observed difference $(0.31_{-0.33}^{+0.51})$ Gyr) is inconsistent with the required difference in formation time, as van der Burg et al. (2020) also pointed out.

It is also unclear how to reconcile our finding of a higher fraction of oblate galaxies in the cluster population with this “early mass quenching” scenario. Observational studies suggest that morphological transformation is a prerequisite for quenching of star formation in central galaxies in the field, presumably as a result of the compaction phase that the galaxies underwent before quenching through internal feedback processes (e.g., Tacchella et al. 2015; Zolotov et al. 2015; Barro et al. 2017). The compaction, which originated from mergers or disk instability, leads to the formation of a bulge while the disk slowly fades due to the declining star formation. The strong association between morphological properties and quiescence (e.g., Lang et al. 2014; Whitaker et al. 2017) is, therefore, a signature of this “mass-quenching” process. If the same process is responsible for the cluster population, then we would naively expect the cluster population to have fewer disk-like oblate galaxies than in the field, given that the cluster population had a head-start. This seems to go against our findings, except for the highest mass bin. This suggests that at least part of the environmentally quenched population originated from a different process than mass-quenching.

We note that other quenching processes that have similar effects to the morphologies are also unlikely to be directly responsible for the “environmentally quenched” oblate population. An example is major mergers<fnref n="30" />, which are expected to lead to the formation of spheroids and quench galaxies through gas funnelling and triggering central starbursts (Hopkins et al. 2009, 2010). Hence, if the excess quenching is due to an elevated merging rate in clusters, then we would also expect to see fewer disk-like galaxies.

### 4.4. Age Variation in the Oblate Quiescent Population

Here we explore the age variation in the cluster quiescent population by combining our results with the measured stellar ages. We utilize the mass-weighted age measurements from Webb et al. (2020), derived using SED fitting of the GOGREEN spectroscopy and photometry. We consider the mass-weighted age in units of cosmic time ($t_{\text{univ,z}} - t_{\text{mw}}$, i.e., the formation time, younger galaxies having a larger/later formation time), taking into account the redshift difference between the clusters in the sample. We refer the readers to Webb et al. (2020) for the methodologies. Limited by the FOV of the HST images, there are 175 quiescent cluster members in the Webb et al. (2020) sample that have both structural parameters and stellar age measurements.

Figure 9 shows the mass-weighted ages of the sample as a function of *n* and *q* in different mass bins. The ages are correlated with mass, with the median formation time and its standard error decreasing from the lowest mass bin ($1.94 \pm 0.23$ Gyr) to the highest mass bin ($1.11 \pm 0.05$ Gyr). We find that the ages do not show a significant trend in *q* alone. While there are formation times that are as late as 2.7 Gyr at low *q*, the median formation time and its standard error at low *q* ($q < 0.4$, $1.31 \pm 0.10$ Gyr) for the whole sample are consistent with those at high *q* ($q \ge 0.4$, $1.21 \pm 0.05$ Gyr). Meanwhile, we find tentative evidence that the median formation time is higher for low *n* galaxies ($n < 2.5$, $1.45 \pm 0.10$) compared to high *n* galaxies ($n \ge 2.5$, $1.15 \pm 0.04$) ($\sim 2\sigma$ difference). We also find similar results by excluding the galaxies in the highest mass bin.

The lack of a strong trend in mass-weighted age with *n* or *q* suggests that disk-like galaxies comprise objects with a mix of ages that are not significantly different from the bulk of the population. Even if we select a strictly “disky” sample with both *n* < 2.5 and *q* < 0.4, the median formation time is $1.42 \pm 0.14$ Gyr with a standard deviation of 0.44 Gyr, which shows that it comprises both young and old galaxies. Hence, not all of the disk-like galaxies in the cluster sample were recently quenched. Instead, some of them are formed and quenched early and remained a disk until the epoch of observation.

The age variation that we see implies that the quenching process that produces the disk excess has been occurring since high redshift. One possibility is that the quenching may happen when or even before the galaxy was accreted into the cluster.

<footnote n="30">³⁰ van der Burg et al. (2020) included a merger-quenching recipe as implemented by Peng et al. (2010) in the early mass-quenching model. They found that the inclusion of merger quenching has no significant effects on the SMF or the quenched fractions.</footnote>
Fossati et al. (2017) studied the fraction of galaxies that were quenched by environmental processes in the five CANDELS/3D-HST fields and showed that satellite galaxies are efficiently environmentally quenched in haloes of all masses at this redshift. Groups in GOGREEN also show higher quenched fraction relative to the field, mostly at the high-mass end (Reeves et al. 2021). Since a cluster grows not only by accreting field central galaxies but also by accreting smaller haloes (McGee et al. 2009), part of the quiescent population in the clusters will be galaxies that have gone through this “pre-processing” stage. Fossati et al. (2017) reported that the inferred quenching time of the satellites is consistent with them being quenched by a gas exhaustion “starvation”-like mechanism, similar to the “over-consumption” model that was proposed by McGee et al. (2014).

In the “over-consumption” model, galaxies are quenched because they essentially run out of fuel to sustain star formation. Outflows from star formation are expected to shorten the quenching time. This effect in turn depends on the star formation rate. No study has yet examined the effect of such a model on galaxy morphologies. However, because galaxies are primarily quenched via gas exhaustion, they will likely retain their disk-like structure and not undergo a drastic morphological change. Simulation studies following the evolution of disk galaxies in group environments also show that their morphological evolution is dependent on the initial inclination of the disk, and that central bulges are not produced or enhanced from interactions with the group environment alone (Villalobos et al. 2012).

Whether or not pre-processing can fully explain the excess quenching that we see in clusters at this redshift range remains to be seen. The effects of pre-processing have been established in clusters at local and intermediate redshifts (e.g., Bianconi et al. 2018; Olave-Rojas et al. 2018; van der Burg et al. 2018). In addition, the observed halo mass dependence of the galaxy ages favors a model with pre-processing at $1 < z < 1.5$ (Reeves et al. 2021). Nevertheless, the effect is expected to be weaker at high redshift and is thought to be negligible in group-size halos above $z \sim 1.5$ (Fossati et al. 2017).

## 4.5. Caveats—More Complex Morphological History?

We have treated morphological transformation as a one-way process throughout this work, in the sense that galaxies only transform from disks to spheroids. While this might be true on a global level, the morphologies of individual galaxies can often have more complex morphological histories that switch back and forth between the two general types. Simulations have shown that morphological transformation is a complex interplay between different processes, including mergers, disk instabilities, and gas accretion (e.g., Brennan et al. 2015); and specifically in the case of clusters, gravitational perturbations, such as tidal shocking at pericenter passages (Joshi et al. 2020).

Given a significant accretion of gas or stars, it is possible for a quiescent spheroid to regrow a disk and subsequently appear as a quiescent disk. De Lucia et al. (2011) studied the rates of disk regrowth in bulge-dominated galaxies in simulations and found that disk regrowth is important for intermediate and low-mass galaxies at high redshifts. Although such a population has not been confirmed at this redshift, there is evidence that some local galaxies have experienced disk regrowth at a certain stage. For example, number density studies by Graham et al. (2015) and de la Rosa et al. (2016) suggested that a population

of the spheroids at high redshift may have regrown a disk and have been hiding in plain sight as the bulges of local spirals and S0 galaxies.

We cannot exclude the possibility that the excess quiescent disk population that we see in clusters is due to disk regrowth, presumably from tidal interactions or minor mergers that predominantly deposit materials in the outer part of the galaxy (e.g., Rodriguez-Gomez et al. 2016). However, the implication of this hypothesis is worth exploring, which requires morphology and age data at different redshifts to track the changes in morphology in both environments over time. One way is to look into the internal gradients of the quiescent population because one would expect the disk-regrown galaxies to show an age difference between their inner structure and the regrown disk. Another way is to study the morphologies of the quiescent population as a function of their local environment (e.g., local surface density or cluster-centric radius). For example, previous studies have found evidence of a relationship between average intrinsic ellipticity of the quiescent cluster population and cluster-centric radius (e.g., D'Eugenio et al. 2015). This analysis would help us to locate where the quiescent disks preferentially reside in the cluster and their relationship with the environment. Unfortunately, given the FOV of the HST data set that we used in this work, we do not have a large enough sample<fnref n="31" /> to perform a similar analysis as a function of cluster-centric radius.

## 5. Summary and Conclusions

In this work, we have presented the axis ratio distributions of 832 galaxies in 11 clusters at $1.0 < z < 1.4$ from the Gemini Observations of Galaxies in Rich Early ENvironments (GOGREEN) Survey. We compare their distributions with a sample of 6471 galaxies in the field, taken from the CANDELS/3D-HST survey, to investigate the effect of the environment on the structural properties of galaxies. Our results can be summarized as follows:

1. Star-forming and quiescent galaxies in clusters show different axis ratio distributions down to the mass limit of $log(M/M_⊙) = 9.5$, similar to the field. The median $q$ of both star-forming galaxies and quiescent galaxies in clusters and the field increases with mass.

2. Massive quiescent galaxies with $log(M/M_⊙) \ge 11$ in both clusters and the field are on average rounder and have a narrower $q$ distribution than their low-mass counterparts.

3. Comparing the axis ratio distribution of star-forming galaxies between cluster and field, we find that they are consistent with each other in all mass bins.

4. The axis ratio distributions of quiescent galaxies in clusters and the field are significantly distinct. For the $10.1 \le log(M/M_⊙) < 10.5$ mass bin, cluster galaxies show a flatter axis ratio distribution with an apparent excess at low $q$ compared to the field.

5. We model the axis ratio distributions of the quiescent population in clusters and the field under different assumptions of their intrinsic shapes, following the methodology of Chang et al. (2013b). We find some evidence that a single-component ( triaxial only) model is unable to reproduce the observed axis ratio distribution of

<footnote n="31">³¹ In the GOGREEN cluster sample, only 27% (47%) of the quiescent (star-forming) galaxies are located in regions with $R \gtrsim 0.5R_{200}$.</footnote>
the cluster galaxies in some mass bins. The axis ratio distribution is instead well-described by a two-component ( triaxial + oblate) model. We find tentative evidence that the cluster distribution in the three lower mass bins has a higher oblate fraction than the field, with the $10.1 \le \log(M/M_{\odot}) < 10.5$ mass bin showing the largest difference ($f_{\text{ob,cluster}} = 0.88 \pm 0.14$ versus $f_{\text{ob,field}} = 0.44 \pm 0.17$).

6. Our modeling shows that both the cluster and field distributions in the highest mass bin are well-described by a single-component model. The contribution of a second oblate component is small, resulting in their distinct distribution shapes.

7. We confirm that two potential sources of bias, the intrinsic shape alignment and dust extinction in cluster galaxies, are unlikely to affect our results. The intrinsic alignment signal in the cluster sample is consistent with zero, and the quiescent galaxies have a low dust content that does not exhibit an axis ratio dependence.

8. We combine the results of the axis ratio distributions, and the quenched fractions of the cluster and the field samples. Using a simple accretion model, we show that the observed axis ratio distributions are statistically consistent with a scenario where no morphological transformation occurs for the environmentally quenched population for the $10.1 \le \log(M/M_{\odot}) < 10.5$ and $10.5 \le \log(M/M_{\odot}) < 10.8$ mass bins. However, the model fails to reproduce the observed axis ratio distributions in the lowest and highest mass bins.

Our results indicate that the environmental quenching mechanism(s) likely produce(s) a population that has a different morphological mix than those resulting from the dominant quenching mechanism in the field. We find that there is also no strong trend in the mass-weighted ages of the quiescent population with *q* or *n*, which suggests that disk-like quiescent galaxies comprise objects with a mix of ages that are not significantly different from the bulk of the population. Our result is consistent with a scenario where the intermediate-mass galaxies are quenched by a starvation-like mechanism(s), such as the “over-consumption” model, that are not expected to drastically change the morphologies. The result of this work suggests that morphology continues to provide important constraints on the underlying physical mechanism that drives environmental quenching at this redshift range.

The authors would also like to thank Arjen van der Wel for useful discussions on the structural parameter measurements of the CANDELS/3D-HST fields.

This work is supported by HST program number GO-15294, and by grant No. 80NSSC17K0019 issued through the NASA Astrophysics Data Analysis Program (ADAP) as well as NSF grants AST-1815475, AST-1517863 and AST-1518257. Support for program number GO-15294 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. B.V. acknowledges funding from the INAF main-stream funding program (PI B. Vulcani) and from the Italian PRIN-Miur 2017 (PI A. Cimatti). J.N. acknowledges support from Universidad Andres Bello internal grant DI-12-19/R.

This paper includes data gathered with the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (NSF. United States), the National Research Council (NRC, Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). This work is based (in part) on observations and archival data obtained with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA; the ESO Telescopes at the La Silla Paranal Observatory under program ID 097.A-0734; the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile; the Canada-France-Hawaii Telescope (CFHT) which is operated by the NRC of Canada, the Institut National des Sciences de l'Université de la Centre National de la Recherche Scientifique of France, and the University of Hawaii; MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA; Subaru Telescope, which is operated by the National Astronomical Observatory of Japan (NAOJ) and the WIRCam, a joint project of CFHT, the Academia Sinica Institute of Astronomy and Astrophysics (ASIAA) in Taiwan, the Korea Astronomy and Space Science Institute (KASI) in Korea, Canada, France.

# Appendix A

## Completeness and Purity of the Cluster Member Selection

Here we expand on the discussion in Section 2.6 about how we chose the selection criteria for the cluster members. The left-hand panel of Figure 10 shows the completeness and purity of the resultant cluster sample as a function of the width of the photometric redshift selection (in units of $\Delta z_{\text{phot}}/(1+z_{\text{phot}})$ around the mean redshift of the clusters). For simplicity, cluster members are defined as those with $|\Delta z_{\text{spec}}/(1+z_{\text{spec}})| \le 0.015$ and above the mass limit of $\log(M/M_{\odot}) \ge 9.5$ in this test. Note that the validity of this test is partially based on the fact that the GOGREEN spectroscopic sample is a representative subset of the photometrically selected galaxy population, for galaxies with stellar mass $\log(M/M_{\odot}) \ge 10.3$. The completeness and purity are a strong function of the $z_{\text{phot}}$ selection width. Increasing this box width results in higher sample completeness at the expense of the sample purity because the sample contains more galaxies that are photometrically selected as cluster members but are interlopers. To optimize both the completeness and purity, we define photometric cluster members as those with $|\Delta z_{\text{phot}}/(1+z_{\text{phot}})| \le 0.06$, which gives a completeness of ~85% and a purity of ~80%.

As we mentioned in Section 2.7, the number of spectroscopically-confirmed GOGREEN field galaxies within the FOV of the HST image is too small for a morphology comparison. While, in theory, we can expand this sample with $z_{\text{phot}}$ like the cluster members, we find that there are no good selection criteria to do so. The right-hand panel of Figure 10 shows the completeness and purity of the “field” sample as a function of the photometric redshift selection width. The genuine field galaxies that can be used for the cluster versus field comparison are defined as those with $|\Delta z_{\text{spec}}/(1+z_{\text{spec}})| > 0.015$ and are within a redshift range of $1.0 < z_{\text{spec}} < 1.4$. In this case, the selection width refers to the width of the redshift region around the mean redshift of the
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Completeness and purity of the sample as a function of the photometric redshift selection width in units of $\Delta z_{phot}/(1+z_{phot})$. Left-hand panel: cluster member selection. Both completeness and purity are a strong function of the selection width. To optimize both the completeness and purity we define photometric cluster members as those with $|\Delta z_{phot}/(1+z_{phot})| \le 0.06$ (gray-dashed line). Right-hand panel: photometric field galaxy selection using GOGREEN data. Since GOGREEN targets massive overdensities, it is difficult to select a clean photometric sample because of contamination from cluster members. We would either end up with a low completeness sample or a sample with low purity. Consequently, we instead choose to draw our field sample from CANDELS and 3D-HST.</figcaption></figure>

cluster that was *avoided* while selecting field galaxies. Increasing the width means that we select galaxies that are further from the cluster redshift, resulting in a higher purity but a less complete field sample. Hence, the trend in completeness and purity are opposite to those shown in the left-hand panel. Overall, the completeness (and therefore size of the resultant sample) drops sharply with the selection width. Even when imposing a wider photometric redshift cut ($0.9 < z_{phot} < 1.5$), we either end up with a low completeness sample (i.e., small sample) or a sample with low purity. For example, using the $\Delta z_{phot}/(1+z_{phot}) \le 0.06$ cut as done for the cluster member selection will result in a low completeness of ~50% and a purity of ~60%. We therefore use galaxies in the CANDELS/ 3D-HST survey as our field comparison sample.

# Appendix B

## Details of the Simulation and the Biases of the Axis Ratio Measurements

In this section we describe the setup and the result of the simulation we used to characterize the biases of our measurements.

We use a similar method as described in Chan et al. (2016, 2018). In this work, we generate a set of 20,000 simulated galaxies, which are uniformly distributed within a magnitude range of $19.0 \le F160W \le 25.0$. Each galaxy has its surface brightness profiles described by a random Sérsic profile. The input structural parameters used to generate these profiles are taken from the parameter distribution of galaxies in the field at a similar redshift range as the GOGREEN clusters (van der Wel et al. 2014a) to mimic real galaxies. These galaxies are then convolved with the appropriate PSF and injected, 20 at a time, to random locations in the sky region of the F160W image. To avoid direct overlap with existing sources on the image, the segmentation maps from SextracTOR are used as a reference of the sky region. We then run these 1000 images through our GALAPAGOS setup to recover the structural parameters of the simulated galaxies. The

simulations were run on the images of five different clusters with a range of richness. We verified that the biases and uncertainties we obtain do not depend on the cluster used.

The simulations are used to refine the GALAPAGOS configuration parameters. Using the bias between the input structural parameters of the simulated galaxies and the recovered galaxies, we then modify the GALAPAGOS configuration parameters and rerun the simulation. Through iterating this process a few times, we optimize our configuration setup by minimizing the biases in the recovered parameters.

Using the simulated galaxies, we also characterize the potential biases for the derived structural parameters. The biases are derived as a function of input magnitudes, *n*, and $\log(R_e)$. We find that the biases are not only a strong function of magnitude (i.e., S/N) but also depend on *n* and $R_e$ (albeit more weakly). Galaxies with higher *n* at a given magnitude show higher biases and larger uncertainties (see also Häussler et al. 2007; van der Wel et al. 2012, for discussions of these second-order effects). Overall, the biases are only significant at faint magnitudes. For example, the Sérsic index *n*, the parameter that is hardest to constrain shows an average bias and dispersion of ~13% and ~24% at F160W = 23 (AB), which corresponds roughly to $\log(M/M_{\odot}) \sim 9.5$. Meanwhile, the effective radius $R_e$ shows an average bias and dispersion of ~6% and ~27% at the same magnitude. The axis ratio *q* shows an average bias and dispersion of ~2% and ~12% at the same magnitude. Figure 11 shows the fractional difference between the recovered and input *q* as a function of F160W magnitude and surface brightness of the simulated galaxies. Among the three Sérsic structural parameters, the axis ratio *q* can typically be measured with the highest accuracy. These bias relations can be used to bias correct the structural parameter measurement for individual galaxies. Nevertheless, since the biases for *q* is almost negligible, we have not applied the bias corrections that we derived from the simulated galaxies on the cluster sample.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Fractional differences between recovered and input structural parameters as function of the input F160W mag (left-hand panel) and surface brightness of the simulated galaxies (right-hand panel). The blue line and the error bars correspond to the median and $1\sigma$ dispersion in different bins (0.5 mag/mag arcsec$^{-2}$ bin width). Green-shaded 2D histogram shows the number density distribution of the simulated galaxies.</figcaption></figure>

## Appendix C

### Axis Ratio Comparison between van der Wel et al. (2014a) and this Work

The axis ratio measurements of the field sample used in this work are taken from the structural parameter catalog of van der Wel et al. (2014a). Although the van der Wel et al. (2014a) measurements are derived with an overall consistent method using GALAPAGOS, different setups and treatment of the images may induce a bias in the cluster and field comparison. To ensure that our measurements are compatible with van der Wel et al. (2014a), we apply our methodology described in Section 2.2 to galaxies in the redshift range of $0.9 < z < 1.5$ in the CANDELS imaging.

Figure 12 shows the result of the comparison as a function of mass. Overall, our derived axis ratios are very consistent with those measured by van der Wel et al. (2014a). There are a total of 6225 galaxies that have a good structural fit in both our measurements and those by van der Wel et al. (2014a) down to $\log(M_*/M_\odot) \sim 9.5$, the mass limit of this work. The entire sample has a median $q$ ratio (ours / van der Wel et al. 2014a) and $1\sigma$ of $1.00 \pm 0.03$ (0% bias and 3% scatter). We also split the sample into UVJ star-forming and quiescent galaxies to check if the ratio differences depend on galaxy types. The blue and red histograms in Figure 12 show the distributions of the $q$ ratios. We find that both galaxy types show consistent ratio distributions.

The bias between the axis ratios derived using our methodologies and those by van der Wel et al. (2014a) is small. It is even smaller than the average bias of $q$ in the cluster sample compared to simulated galaxies. Therefore, this is not a

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of axis ratio derived using the method used in this work with the van der Wel et al. (2014a) measurements as a function of mass. Blue points correspond to star-forming galaxies and red points correspond to the quiescent galaxies. The median ratio and $1\sigma$ between our measurements and van der Wel et al. (2014a) are $1.00 \pm 0.03$ down to $\log(M_*/M_\odot) \sim 9.5$. Both galaxy types show similar ratio distributions, as shown in the histogram on the right-hand side.</figcaption></figure>

major source of uncertainty in our results. Meanwhile, we note that the biases of $n$ and $R_e$ between our measurements and van der Wel et al. (2014a) are non-negligible. Since we have not compared $n$ and $R_e$ between the clusters and the field in this work, we defer the discussion on these biases to an upcoming paper on the mass-size relations of the two samples (J. C. C. Chan et al. 2021 in preparation).
**Appendix D**

**Results of the Axis Ratio Distribution Modeling**

Here we expand on Section 3.2 and give more details on the
fitting results. Figures 13 and 14 show the corner plots of the
model fitting of the bootstrapped sample for the Case II models,

where $f_{\text{ob}}$, $b$, $\sigma_b$ are left as free parameters. The best-fit parameter values, and the median and $1\sigma$ values of the parameters derived from the bootstrapped sample are provided in the histograms. We can see that the $b$ and $\sigma_b$ parameters for the field sample in the two lower mass bins are poorly constrained.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Corner plots of the Case II model fitting of the bootstrapped cluster and field sample in the four mass bins. The color of the green shade in each panel represents the distribution of the bootstrap results. The gray-dashed lines correspond to the best-fit parameter values derived from the original sample. The vertical lines in the histogram correspond to the 16th, 50th, and 84th percentile of the distribution.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Corner plots of the Case II model fitting of the bootstrapped cluster and field sample in the four mass bins (cont'd).</figcaption></figure>

For the case II models we have assumed the values of the four parameters ($E$, $\sigma_E$, $T$, $\sigma_T$) to be the same as the best-fit values in Chang et al. (2013b). However, we find that our results are not very sensitive to these parameters. For example, we tested that using the initial parameters for the $10.5 \le \log(M/M_\odot) < 10.8$ bin to fit the $10.1 \le \log(M/M_\odot) < 10.5$ bin gives consistent results. The best-fit parameters for the bootstrapped cluster sample are $f_{\text{ob}} = 0.88 \pm 0.15$, $b = 0.33 \pm 0.04$, $\sigma_b = 0.09 \pm 0.03$, which is completely consistent with the result using the “correct” initial parameters.

## Appendix E
Intrinsic Alignments in Cluster Galaxies

The axis ratio distribution modeling used in this work relies on the assumption that the galaxy population is observed from random viewing angles, which breaks down when the galaxies are not randomly oriented (i.e., aligned with respect to a certain direction). The alignment of, or between, galaxies within a halo has been a subject of heated interest (See Joachimi et al. 2015, for a review). There are numerous types of intra-halo galaxy
alignment. In this section, we focus on the alignment of the shape of the satellite galaxies, which is the only one that will impact the observed axis ratio distribution.

In numerical simulations, it is established that major axes of satellite galaxies are preferentially oriented toward the center of mass of the halo or the central galaxy (e.g., Faltenbacher et al. 2008; Knebe et al. 2020; Tenneti et al. 2021). The strength of this radial alignment is expected to be strongest at the cluster core and to decrease with cluster-centric radius. Observationally, whether this alignment exists is still a matter of debate. Sifón et al. (2015) measured the alignments of galaxies in 90 clusters at $0.05 < z < 0.55$ and detected no alignments out to $3R_{200}$. Meanwhile, more recent works by Huang et al. (2018) and Georgiou et al. (2019) reported a radial alignment signal that is stronger in satellites with smaller distance to the cluster/group center, although the measured signal is a few times smaller than predicted in the simulations. The alignment studies are limited to low-redshift clusters. Here we would like to examine the radial alignment signal in our sample.

We follow the procedure outlined in Sifón et al. (2015) to measure the alignment of the galaxies with respect to the center of the cluster. In this work, the center of each cluster is taken to be the location of the BCG. To quantify the alignment signal, we adopt the commonly-used ellipticity components $\epsilon_+$ and $\epsilon_-$, which are defined as:

$$ \epsilon_+ = \epsilon_1 \cos 2\theta + \epsilon_2 \sin 2\theta, \\ \epsilon_- = \epsilon_1 \sin 2\theta - \epsilon_2 \cos 2\theta \qquad (E1) $$

where $\epsilon_1$ and $\epsilon_2$ are the galaxy ellipticities in the Cartesian frame and $\theta$ is the azimuthal angle of the individual galaxy with respect to the BCG of the cluster. The ellipticity $\epsilon_1$ measures the ellipticity in the R.A. and decl. directions and $\epsilon_2$ in diagonal directions. The two quantities are related to the axis ratio $q$ through the following:

$$ \epsilon_1 = \frac{1-q}{1+q} \cos 2\phi, \\ \epsilon_2 = \frac{1-q}{1+q} \sin 2\phi \qquad (E2) $$

where the angle $\phi$ is the position angle of the major axis of the galaxy. The two Equations (E1) show that $\epsilon_+$ and $\epsilon_-$ are rotated to the frame with one axis pointing toward the radial direction from the cluster center. A positive (negative) $\epsilon_+$ therefore indicates a radial (tangential) alignment of the galaxies toward the center of the cluster. Meanwhile, $\epsilon_-$ measures the shape alignment at $\pm 45^\circ$ from the radial direction. This component is commonly used as a check for systematic effects, as $\epsilon_-$ should be consistent with zero due to symmetry.

Figure 15 shows the average alignment of our cluster sample ($\langle\epsilon_+\rangle$ and $\langle\epsilon_-\rangle$) as a function of cluster-centric radius. Given the field of view of the HST images we can measure the alignment only out to $1R_{200}$. We find that the average radial/tangential alignment for the entire cluster sample within $1R_{200}$ is consistent with zero, with $\langle\epsilon_+\rangle = -0.0011 \pm 0.0079$. The cross component is also consistent with zero, with $\langle\epsilon_-\rangle = 0.0018 \pm 0.0077$. Examining the alignment signal as a function of radius, there is weak evidence that the average radial alignment is positive in the region close to the cluster center ($\sim 1.4\sigma$ for $R < 0.2R_{200}$).

<figure><img src="image_18.png" /><figcaption>Figure 15. Average alignments ($\langle\epsilon_+\rangle$ and $\langle\epsilon_-\rangle$) of the galaxies in GOGREEN clusters. The circles correspond to $\epsilon_+$, while squares correspond to $\epsilon_-$. The error bars correspond to the standard error of the mean. Both alignment signals are consistent with zero.</figcaption></figure>

Since the alignment is expected to be strongest in regions close to the center, we repeat our axis ratio distribution analysis excluding galaxies in the $R < 0.2R_{200}$ region and find that it does not affect our conclusion. The exclusion only results in an increase in uncertainty due to lower number statistics. The potential intrinsic alignments in the sample are therefore unlikely to affect our results.

## ORCID iDs

Jeffrey C. C. Chan https://orcid.org/0000-0001-6251-3125
Gregory Rudnick https://orcid.org/0000-0001-5851-1856
Adam Muzzin https://orcid.org/0000-0002-9330-9108
Kristi A. Webb https://orcid.org/0000-0002-8610-0672
Andrea Biviano https://orcid.org/0000-0002-0857-0732
M. C. Cooper https://orcid.org/0000-0003-1371-6019
Gabriella De Lucia https://orcid.org/0000-0002-6220-9104
Ricardo Demarco https://orcid.org/0000-0003-3921-2177
Ben Forrest https://orcid.org/0000-0001-6003-0541
Pascale Jablonka https://orcid.org/0000-0002-9655-1063
Chris Lidman https://orcid.org/0000-0003-1731-0497
Sean L. McGee https://orcid.org/0000-0003-3255-3139
Irene Pintos-Castro https://orcid.org/0000-0002-9133-4457
Bianca Poggianti https://orcid.org/0000-0001-8751-8360
Andrew M. M. Reeves https://orcid.org/0000-0003-2618-6408
Benedetta Vulcani https://orcid.org/0000-0003-0980-1499
Howard K. C. Yee https://orcid.org/0000-0003-4935-2720
Dennis Zaritsky https://orcid.org/0000-0002-5177-727X

## References

* Baldry, I. K., Balogh, M. L., Bower, R. G., et al. 2006, MNRAS, 373, 469

* Balogh, M. L., Baldry, I. K., Nichol, R., et al. 2004, ApJL, 615, L101

* Balogh, M. L., Gilbank, D. G., Muzzin, A., et al. 2017, MNRAS, 470, 4168

* Balogh, M. L., McGee, S. L., Mok, A., et al. 2016, MNRAS, 456, 4364

* Balogh, M. L., Morris, S. L., Yee, H. K. C., Carlberg, R. G., & Ellingson, E. 1997, ApJL, 488, L75

* Balogh, M. L., Navarro, J. F., & Morris, S. L. 2000, ApJ, 540, 113

* Balogh, M. L., van der Burg, R. F. J., Muzzin, A., et al. 2021, MNRAS, 500, 358

* Barden, M., Häußler, B., Peng, C. Y., McIntosh, D. H., & Guo, Y. 2012, MNRAS, 422, 449

* Barro, G., Faber, S. M., Koo, D. C., et al. 2017, ApJ, 840, 47

* Bassett, R., Papovich, C., Lotz, J. M., et al. 2013, ApJ, 770, 58
Belli, S., Newman, A. B., & Ellis, R. S. 2015, ApJ, 799, 206

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bianconi, M., Smith, G. P., Haines, C. P., et al. 2018, MNRAS, 473, L79

Biviano, A., van der Burg, R. F. J., Balogh, M. L., et al. 2021, A&A, 650, 105

Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

Boselli, A., & Gavazzi, G. 2006, PASP, 118, 517

Boselli, A., & Gavazzi, G. 2014, A&ARv, 22, 74

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Brammer, G. B., van Dokkum, P. G., et al. 2012, ApJS, 200, 13

Brennan, R., Pandya, V., Somerville, R. S., et al. 2015, MNRAS, 451, 2933

Brodwin, M., Ruel, J., Ade, P. A. R., et al. 2010, ApJ, 721, 90

Brough, S., van der Sande, J., Owers, M. S., et al. 2017, ApJ, 844, 59

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Cappellari, M. 2016, ARA&A, 54, 597

Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 416, 1680

Carollo, C. M., Bschorr, T. J., Renzini, A., et al. 2013, ApJ, 773, 112

Chabrier, G. 2003, PASP, 115, 763

Chan, J. C. C., Beifiorsi, A., Mendel, J. T., et al. 2016, MNRAS, 458, 3181

Chan, J. C. C., Beifiorsi, A., Saglia, R. P., et al. 2018, ApJ, 856, 8

Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013a, ApJ, 762, 83

Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013b, ApJ, 773, 149

Cole, J., Bezanson, R., van der Wel, A., et al. 2020, ApJL, 890, L25

Cooper, M. C., Coil, A. L., Gerke, B. F., et al. 2010, MNRAS, 409, 337

Cooper, M. C., Griffith, R. L., Newman, J. A., et al. 2012, MNRAS, 419, 3018

Cortese, L., Fogarty, L. M. R., Bekki, K., et al. 2016, MNRAS, 463, 170

de la Rosa, I. G., La Barbera, F., Ferreras, I., et al. 2016, MNRAS, 457, 1916

De Lucia, G., Boylan-Kolchin, M., Benson, A. J., Fontanot, F., & Monaco, P. 2010, MNRAS, 406, 1533

De Lucia, G., Fontanot, F., Wilman, D., & Monaco, P. 2011, MNRAS, 414, 1439

De Lucia, G., Weinmann, S., Poggianti, B. M., Aragón-Salamanca, A., & Zaritsky, D. 2012, MNRAS, 423, 1277

Delaye, L., Huertas-Company, M., Mei, S., et al. 2014, MNRAS, 441, 203

Demarco, R., Wilson, G., Muzzin, A., et al. 2010, ApJ, 711, 1185

D'Eugenio, F., Houghton, R. C. W., Davies, R. L., & Dalla Bontà, E. 2013, MNRAS, 429, 1258

D'Eugenio, F., Houghton, R. C. W., Davies, R. L., & Dalla Bontà, E. 2015, MNRAS, 451, 827

Dressler, A. 1980, ApJ, 236, 351

Emsellem, E., Cappellari, M., Krajnović, D., et al. 2011, MNRAS, 414, 888

Faltenbacher, A., Jing, Y. P., Li, C., et al. 2008, ApJ, 675, 146

Fasano, G., Poggianti, B. M., Couch, W. J., et al. 2000, ApJ, 542, 673

Fillingham, S. P., Cooper, M. C., Wheeler, C., et al. 2015, MNRAS, 454, 2039

Foley, R. J., Andersson, K., Bazin, G., et al. 2011, ApJ, 731, 86

Fossati, M., Wilman, D. J., Mendel, J. T., et al. 2017, ApJ, 835, 153

Foster, C., van der Sande, J., D'Eugenio, F., et al. 2017, MNRAS, 472, 966

Franx, M., Illingworth, G., & de Zeeuw, T. 1991, ApJ, 383, 112

Fumagalli, M., Fossati, M., Hau, G. K. T., et al. 2014, MNRAS, 445, 4335

Gavazzi, G., Boselli, A., Mayer, L., et al. 2001, ApJL, 563, L23

Georgiou, C., Chisari, N. E., Fortuna, M. C., et al. 2019, A&A, 628, A31

Gonzaga, S., Hack, W., Fruchter, A., & Mack, J. 2012, The DrizzlePac Handbook (Baltimore, M.D.: STScI)

Graham, A. W., Dullo, B. T., & Savorgnan, G. A. D. 2015, ApJ, 804, 32

Graham, M. T., Cappellari, M., Bershady, M. A., & Drory, N. 2019, arXiv:1910.05139

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Gunn, J. E., & Gott, J. R., III 1972, ApJ, 176, 1

Häußler, B., Bamford, S. P., Vika, M., et al. 2013, MNRAS, 430, 330

Häussler, B., McIntosh, D. H., Barden, M., et al. 2007, ApJS, 172, 615

Hinton, S. R., Davis, T. M., Lidman, C., Glazebrook, K., & Lewis, G. F. 2016, A&C, 15, 61

Holden, B. P., Illingworth, G. D., Franx, M., et al. 2007, ApJ, 670, 190

Holden, B. P., van der Wel, A., Rix, H.-W., & Franx, M. 2012, ApJ, 749, 96

Hopkins, P. F., Cox, T. J., Younger, J. D., & Hernquist, L. 2009, ApJ, 691, 1168

Hopkins, P. F., Croton, D., Bundy, K., et al. 2010, ApJ, 724, 915

Houghton, R. C. W., Davies, R. L., D'Eugenio, F., et al. 2013, MNRAS, 436, 19

Huang, H.-J., Mandelbaum, R., Freeman, P. E., et al. 2018, MNRAS, 474, 4772

Joachimi, B., Cacciato, M., Kitching, T. D., et al. 2015, SSRv, 193, 1

Johnston, E. J., Aragón-Salamanca, A., & Merrifield, M. R. 2014, MNRAS, 441, 333

Joshi, G. D., Pillepich, A., Nelson, D., et al. 2020, MNRAS, 496, 2673

Just, D. W., Zaritsky, D., Sand, D. J., Desai, V., & Rudnick, G. 2010, ApJ, 711, 192

Kawinwanichakij, L., Papovich, C., Quadri, R. F., et al. 2017, ApJ, 847, 134

Kelkar, K., Gray, M. E., Aragón-Salamanca, A., et al. 2017, MNRAS, 469, 4551

Knebe, A., Gámez-Marín, M., Pearce, F. R., et al. 2020, MNRAS, 495, 3002

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Labbé, I., Huang, J., Franx, M., et al. 2005, ApJL, 624, L81

Lang, P., Wuyts, S., Somerville, R. S., et al. 2014, ApJ, 788, 11

Lani, C., Almaini, O., Hartley, W. G., et al. 2013, MNRAS, 435, 207

Larson, R. B., Tinsley, B. M., & Caldwell, C. N. 1980, ApJ, 237, 692

Leja, J., Johnson, B. D., Conroy, C., et al. 2019, ApJ, 877, 140

Lidman, C., Iacobuta, G., Bauer, A. E., et al. 2013, MNRAS, 433, 825

Man, A. W. S., Zirm, A. W., & Toft, S. 2016, ApJ, 830, 89

Matharu, J., Muzzin, A., Brammer, G. B., et al. 2019, MNRAS, 484, 595

Matharu, J., Muzzin, A., Brammer, G. B., et al. 2020, MNRAS, 493, 6011

McGee, S. L., Balogh, M. L., Bower, R. G., Font, A. S., & McCarthy, I. G. 2009, MNRAS, 400, 937

McGee, S. L., Bower, R. G., & Balogh, M. L. 2014, MNRAS, 442, L105

McPartland, C., Ebeling, H., Roediger, E., & Blumenthal, K. 2016, MNRAS, 455, 2994

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Moore, B., Lake, G., & Katz, N. 1998, ApJ, 495, 139

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013a, ApJ, 777, 18

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013b, ApJS, 206, 8

Muzzin, A., Wilson, G., Yee, H. K. C., et al. 2009, ApJ, 698, 1934

Muzzin, A., Wilson, G., Yee, H. K. C., et al. 2012, ApJ, 746, 188

Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJL, 699, L178

Newman, A. B., Ellis, R. S., Bundy, K., & Treu, T. 2012, ApJ, 746, 162

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Olave-Rojas, D., Cerulo, P., Demarco, R., et al. 2018, MNRAS, 479, 2328

Old, L. J., Balogh, M. L., van der Burg, R. F. J., et al. 2020, MNRAS, 493, 5987

Padilla, N. D., & Strauss, M. A. 2008, MNRAS, 388, 1321

Papovich, C., Kawinwanichakij, L., Quadri, R. F., et al. 2018, ApJ, 854, 30

Peng, Y.-J., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Pintos-Castro, I., Yee, H. K. C., Muzzin, A., Old, L., & Wilson, G. 2019, ApJ, 876, 40

Poggianti, B. M., Calvi, R., Bindoni, D., et al. 2013, ApJ, 762, 77

Poggianti, B. M., Fasano, G., Omizzolo, A., et al. 2013, AJ, 151, 78

Postman, M., Franx, M., Cross, N. J. G., et al. 2005, ApJ, 623, 721

Pulsoni, C., Gerhard, O., Arnaboldi, M., et al. 2018, A&A, 618, A94

Qu, Y., Helly, J. C., Bower, R. G., et al. 2017, MNRAS, 464, 1659

Reeves, A. M. M., Balogh, M. L., van der Burg, R. F. J., et al. 2021, MNRAS, 506, 3364

Rix, H.-W., & Zaritsky, D. 1995, ApJ, 447, 82

Roberts, I. D., & Parker, L. C. 2020, MNRAS, 495, 554

Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371

Ryden, B. S. 2004, ApJ, 601, 214

Saglia, R. P., Sánchez-Blázquez, P., Bender, R., et al. 2010, A&A, 524, A6

Sandage, A., Freeman, K. C., & Stokes, N. R. 1970, ApJ, 160, 831

Sersic, J. L. 1968, Atlas de Galaxias Australes

Shankar, F., Marulli, F., Bernardi, M., et al. 2013, MNRAS, 428, 109

Sheen, Y.-K., Smith, R., Affé, Y., et al. 2017, ApJ, 763, L7

Sifón, C., Hoekstra, H., Cacciato, M., et al. 2015, A&A, 575, A48

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Stalder, B., Ruel, J., Suhada, R., et al. 2013, ApJ, 763, 93

Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2019, ApJ, 877, 103

Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Sci, 348, 314

Tenneti, A., Kitching, T. D., Joachimi, B., & Di Matteo, T. 2021, MNRAS, 501, 5859

Tremblay, B., & Merritt, D. 1996, AJ, 111, 2243

Trujillo, I., Förster Schreiber, N. M., Rudnick, G., et al. 2006, ApJ, 650, 18

Valentinuzzi, T., Poggianti, B. M., Saglia, R. P., et al. 2010, ApJL, 721, L19

van der Burg, R. F. J., McGee, S., Aussel, H., et al. 2018, A&A, 618, A140

van der Burg, R. F. J., Muzzin, A., Hoekstra, H., et al. 2013, A&A, 557, A15

van der Burg, R. F. J., Rudnick, G., Balogh, M. L., et al. 2020, A&A, 638, A112

van der Wel, A., Bell, E. F., Häussler, B., et al. 2012, ApJS, 203, 24

van der Wel, A., Chang, Y.-Y., Bell, E. F., et al. 2014b, ApJL, 792, L6

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014a, ApJ, 788, 28

van der Wel, A., Rix, H.-W., Holden, B. P., Bell, E. F., & Robaina, A. R. 2009, ApJL, 706, L120
van Dokkum, P. G., & Franx, M. 2001, ApJ, 553, 90

van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018

Veale, M., Ma, C.-P., Greene, J. E., et al. 2017, MNRAS, 471, 1428

Villalobos, Á., De Lucia, G., Borgani, S., & Murante, G. 2012, MNRAS, 424, 2401

Vulcani, B., Poggianti, B. M., Dressler, A., et al. 2011, MNRAS, 413, 921

Wang, L., & Kauffmann, G. 2008, MNRAS, 391, 785

Webb, K., Balogh, M. L., Leja, J., et al. 2020, MNRAS, 498, 5317

Weijmans, A.-M., de Zeeuw, P. T., Emsellem, E., et al. 2014, MNRAS, 444, 3340

Wetzel, A. R., Tinker, J. L., & Conroy, C. 2012, MNRAS, 424, 232

Wetzel, A. R., Tinker, J. L., Conroy, C., & van den Bosch, F. C. 2013, MNRAS, 432, 336

Whitaker, K. E., Bezanson, R., van Dokkum, P. G., et al. 2017, ApJ, 838, 19

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Wilson, G., Muzzin, A., Yee, H. K. C., et al. 2009, ApJ, 698, 1943

Windhorst, R. A., Cohen, S. H., Hathi, N. P., et al. 2011, ApJS, 193, 27

Yagi, M., Gu, L., Koyama, Y., et al. 2015, AJ, 149, 36

Zhang, H., Primack, J. R., Faber, S. M., et al. 2019, MNRAS, 484, 5170

Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327