Redshift and stellar mass dependence of intrinsic shapes
of disc-dominated galaxies from COSMOS observations
below z = 1.0

K. Hoffmann, C. Laigle, N. E. Chisari, P. Tallada-Crespí, R. Teyssier, Y.

Dubois, J. Devriendt

► To cite this version:

K. Hoffmann, C. Laigle, N. E. Chisari, P. Tallada-Crespí, R. Teyssier, et al. Redshift and stellar mass dependence of intrinsic shapes of disc-dominated galaxies from COSMOS observations below z = 1.0. Monthly Notices of the Royal Astronomical Society, 2022, 10.1093/mnras/stac1988. insu-03748237

HAL Id: insu-03748237

https://insu.hal.science/insu-03748237

Submitted on 9 Aug 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Redshift and stellar mass dependence of intrinsic shapes of disc-dominated galaxies from COSMOS observations below $z = 1.0$

K. Hoffmann<sup>D, 1,2★</sup> C. Laigle,<sup>3</sup> N. E. Chisari<sup>D, 4</sup> P. Tallada-Crespí,<sup>5,6</sup> R. Teyssier<sup>D, 2</sup> Y. Dubois<sup>3</sup>
and J. Devriendt<sup>7</sup>

Accepted 2022 July 8. Received 2022 June 17; in original form 2020 December 14

## ABSTRACT

The high abundance of disc galaxies without a large central bulge challenges predictions of current hydrodynamic simulations of galaxy formation. We aim to shed light on the formation of these objects by studying the redshift and mass dependence of their intrinsic 3D shape distributions in the COSMOS galaxy survey below redshift $z = 1.0$. This distribution is inferred from the observed distribution of 2D shapes, using a reconstruction method which we test using hydrodynamic simulations. Our tests reveal a moderate bias for the inferred average disc circularity and relative thickness, but a large bias on the dispersion of these quantities. Applying the reconstruction method on COSMOS data, we find variations of the average disc circularity and relative thickness with redshift of around ~1 per cent and ~10 per cent, respectively, which is comparable to the error estimates on these quantities. The average relative disc thickness shows a significant mass dependence which can be accounted for by the scaling of disc radius with galaxy mass. We conclude that our data provides no evidence for a strong dependence of the average circularity and absolute thickness of disc-dominated galaxies on redshift and mass that is significant with respect to the statistical uncertainties in our analysis. These findings are expected in the absence of disruptive merging or feedback events that would affect galaxy shapes. They hence support a scenario where present-day discs form early ($z > 1.0$) and subsequently undergo a tranquil evolution in isolation. However, more data and a better understanding of systematics are needed to reaffirm our results.

**Key words:** galaxies: evolution – galaxies: spiral – galaxies: disc – galaxies: statistics – methods: statistical.

## 1 INTRODUCTION

The structure of late-type galaxies is comprised of different components with distinct kinematic and morphological characteristics. The defining component is the thin stellar disc with its spiral arm over densities (e.g. Hubble 1926, 1936; Vaucouleurs 1959a, b; Freeman 1970; van der Kruit & Searle 1982). Observations at low redshifts and in our own Milky Way revealed that this thin disc is often enclosed by a thick stellar disc as well as a stellar halo with relatively low densities (e.g. Carollo et al. 2010; Trujillo & Bakos 2013; Martínez-Lombilla & Knapen 2019, among others). Near the Galactic Centre, the spiral arms of most discs transition either into a bar or into a bulge. The contribution of these components to the overall mass and morphology of a given galaxy is determined by the galaxies' formation history (e.g. Binney & Tremaine 2008; Vogelsberger et al. 2020). Therefore, studying how the shapes of late-type galaxies (hereafter referred to as disc galaxies) are distributed at different

epochs of the Universe can provide significant insights into the processes which dominate galaxy formation, as detailed below.

During the recent years, it was pointed out that the observed abundance of a particular type of objects, namely disc-dominated galaxies with only a small or no central bulge, challenges our current understanding of galaxy formation in the ΛCDM model (e.g. Kormendy et al. 2010). As a matter of fact, bulges are expected to form in galaxies for several reasons. van den Bosch (1998) suggested that bulges can form ‘inside out’ from the low-angular momentum components of the initial gas overdensities in which the galaxies are born (see also Kepner 1999). In another popular scenario, bulges result from the redistribution of angular momentum within the disc by central bars, spiral arms, or bending instabilities which distort the stellar orbits (Kormendy & Kennicutt 2004; Debattista et al. 2006). This so-called ‘secular evolution’ leads typically to the formation of a small ‘pseudo’ bulge, but may also promote the emergence of a larger ‘classical’ bulge when occurring at high redshifts (Elmegreen, Bournaud & Elmegreen 2008; Bournaud 2016). A third and important scenario for bulge formation is the accretion of satellite galaxies, globular clusters or clumpy streams

<footnote n="1"><sup>1</sup>Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="2"><sup>2</sup>Institute for Computational Science, University of Zurich, Winterthurerstr. 190, CH-8057 Zürich, Switzerland</footnote>

<footnote n="3"><sup>3</sup>CNRS and UPMC Univ. Paris 06, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Institute for Theoretical Physics, Utrecht University, Princetonplein 5, NL-3584 CC Utrecht, The Netherlands</footnote>

<footnote n="5"><sup>5</sup>Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, E-28040 Madrid, Spain</footnote>

<footnote n="6"><sup>6</sup>Port d'Informació Científica (PIC), Campus UAB, C. Albareda s/n, E-08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="7"><sup>7</sup>Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford, OX1 3RH, UK</footnote>

<footnote>*E-mail: kai.d.hoffmann@gmail.com</footnote>
of cold gas, which can move low-angular momentum material to the centre of the disc (e.g. through violent disc instabilities) without destroying it (e.g. Walker, Mihos & Hernquist 1996; Dekel, Sari & Ceverino 2009; Hopkins et al. 2009; Dubois et al. 2012; Kretschmer, Agertz & Teyssier 2020). In relatively rare occasions, even major mergers can lead to the formation of a bulge-dominated discs, as discs can survive or reform if one of the progenitors is gas-rich (Toomre 1977; Negroponte & White 1983; Hernquist 1992; Naab & Burkert 2003; Hopkins et al. 2008, 2009; Jackson et al. 2020), especially in the case of prograde mergers (Martin et al. 2018). Since mergers and cold streams are believed to play a significant role in the formation and growth of disc galaxies (e.g. Kauffmann, White & Guiderdoni 1993; Baugh, Cole & Frenk 1996; Aguerre, Balcells & Peletier 2001; Eliche-Moral et al. 2006) it is difficult to explain how a significant fraction of high mass discs (i.e. $M_* \gtrsim 10^{10} M_\odot$) could grow without forming a large bulge (Kautsch et al. 2006; Kautsch 2009; Weinzirl et al. 2009; Buta et al. 2015).

One possible solution to this problem are feedback processes, driven by starbursts, supernovae, or active Galactic nuclei, which can suppress the formation of a bulge by removing low-angular momentum material from the disc, in particular after gas-rich major mergers (e.g. Governato et al. 2010; Brook et al. 2011, 2012; Hopkins et al. 2012; Übler et al. 2014; Dubois et al. 2016; Grand et al. 2017). However, simulations suggest that this suppression is only effective for discs with stellar masses of the Milky Way ($\simeq 6 \times 10^{10} M_\odot$) or below, but not at higher masses (Brooks & Christensen 2016). An additional challenge for formation scenarios including major mergers are the low-density environments in which bulgeless discs typically reside (e.g. Grossi et al. 2018). Kormendy et al. (2010) argued that this finding speaks for a gentle, rather than a violent mass accretion history (see also Kormendy 2016; Jackson et al. 2020). Recently, Peebles (2020) therefore discussed an alternative solution to the problem, according to which small-scale non-Gaussianities in the initial conditions of a warm or mixed dark matter universe could lead to a large fraction of bulge-less discs. However, this approach remains to be explored with simulations of galaxy formation.

With this work, we aim to shed light on the formation of bulge-less discs. We therefore use a novel approach for discriminating the first scenario in which feedback processes suppress bulge formation after mergers from alternative, less violent mass accretion scenarios by comparing the distributions of shapes of observed disc-dominated galaxies at different redshift and stellar mass ranges. This approach relies on the hypothesis that disruptive events, such as mergers and strong feedback should dramatically affect the morphology of the remaining Galactic disc, besides the absence of a bulge. Such changes can occur in the form of an increase in disc thickness which results from vertical heating by feedback and merging events (e.g. Quinn, Hernquist & Fullagar 1993; Grand et al. 2016) or tidal debris of mergers (Abadi et al. 2003). The latter can furthermore lead to large sub-structures which decrease the disc circularity, such as accreted satellites or spiral arms that form due to tidal interactions during merging events (Springel & Hernquist 2005; Robertson et al. 2006; Peschken, Łokas & Athanassoula 2020). We therefore expect mergers and strong feedback to cause a significant mass and redshift dependence of the discs thickness and circularity.

On the contrary, if the bulge-less discs underwent an early and regular accretion of mass, for instance through infalling smooth streams of cold gas at high redshifts followed by a calm evolution without major mergers, their morphologies should exhibit a much weaker or no dependency on mass and redshift. It has been discussed in the literature that a secular redshift evolution of the disc thickness can be expected from vertical heating by bars, spiral arms, stellar

clumps, and giant molecular clouds (e.g. Bournaud, Elmegreen & Martig 2009; Saha, Tseng & Taam 2010; Aumer, Binney & Schönrich 2016). However, Grand et al. (2016) find that the effect of such internal perturbations is weak compared to those caused by merging events. Moreover, Park et al. (2021) do not find a significant redshift dependence of the disc thickness in an ensemble of simulated disc galaxies that underwent a quiescent growth without significant mergers since $z = 1.0$.

A challenging aspect of studying observed galaxy morphologies is to interpret the projected two-dimensional (2D) galaxy shape distributions in terms of 3D models of galaxy formation. In this work we address this challenge by reconstructing the distribution of 3D galaxy shapes from the observed distribution of 2D shapes. We thereby assume a simple ellipsoidal model for the 3D light distribution within a given galaxy. The 3D shape of the ellipsoids is thereby fully characterized by two of the three possible ratios between the major, intermediate and minor axes ($A_{3D}$, $B_{3D}$, and $C_{3D}$, respectively),

$$q_{3D} \equiv \frac{B_{3D}}{A_{3D}}, r_{3D} \equiv \frac{C_{3D}}{B_{3D}}, s_{3D} \equiv \frac{C_{3D}}{A_{3D}}. \qquad (1)$$

For disc galaxies, the $q_{3D}$ parameter can be regarded as a measure for the circularity, while $r_{3D}$ and $s_{3D}$ both quantify the relative disc thickness. The 2D galaxy shapes of the projected ellipsoid are described by the axial ratio

$$q_{2D} \equiv \frac{B_{2D}}{A_{2D}}. \qquad (2)$$

The main advantage of such an ellipsoidal approximation is that the distribution of 2D axis ratios can be predicted from a given model for the 3D axial ratio distribution at low-computational cost. This allows for the reconstruction of the 3D axes ratio distribution by tuning the corresponding model parameters such that the predicted distribution of 2D axis ratios matches observations. This reconstruction technique dates back to Hubble (1926) who derived first constraints on the intrinsic 3D axis ratios of galaxies by modelling them as oblate ellipsoids. Over the last century it has been shown that, despite its simplicity, this methodology reproduces the observed 2D axial ratio distribution of late- as well as of early-type galaxies, assuming oblate or prolate ellipsoidal models (e.g. Sandage, Freeman & Stokes 1970; Binney 1978; Noerdlinger 1979). The agreement with observed axial ratio distributions was further improved by modelling galaxies as triaxial ellipsoids, which allowed for more detailed interpretations of the observations (Benacchio & Galletta 1980; Binney & de Vaucouleurs 1981; Lambas, Maddox & Loveday 1992). These early studies where continued using larger samples to study the relation between intrinsic galaxy shapes and properties, such as size, luminosity, and colour in the local universe, observed by the Sloan Digital Sky Survey (SDSS, e.g. Ryden 2004; Vincent & Ryden 2005; Padilla & Strauss 2008; Rodríguez & Padilla 2013). The evolution of intrinsic shapes with redshift has been studied with the same approach in the galaxy surveys SDSS, 3D-HST, GOODS, COSMOS, and CANDELS (Yuma et al. 2011; Holden et al. 2012; Yuma, Ohta & Yabe 2012; Chang et al. 2013; van der Wel et al. 2014; Takeuchi et al. 2015; Satoh, Kajisawa & Himoto 2019; Zhang et al. 2019). Physical interpretations drawn from these reconstructed axial ratio distributions rely on the validity of the ellipsoidal model for the 3D galaxy shapes as well as on the accuracy of the model for the 3D axial ratio distribution. In addition, observational sources of systematics on the observed 2D axial ratio distribution, for example those induced by dust extinction within the galaxy, need to be taken into account in the analysis.
The objective of our analysis is therefore two fold. Our main goal is to constrain potential explanations for the absence of large bulges in disc galaxies based on the reconstructed 3D shape distribution of disc-dominated galaxies as outlined earlier in this section. Our investigation is based on data from the COSMOS galaxy survey (Scoville et al. 2007), which provides excellent space-based imaging of galaxies over a wide range of redshifts and stellar masses, and is therefore ideal for our analysis.

However, the assumptions on which the reconstruction of the 3D shape distribution is based are highly simplistic which may induce uncertainties in our as well as in previous analyses. Our second goal is therefore to test several of these assumptions and to assess the overall performance of the reconstruction method. For that purpose we employ for the first time two state-of-the-art cosmological hydrodynamic simulations of galaxy formation, Horizon AGN and Illustris TNG. These simulations provide 3D as well as projected 2D galaxy morphologies which allows for an examination of the reconstruction method under controlled conditions.

This paper is structured as follows. In Section 2, we present the galaxy catalogues from the COSMOS survey and the hydrodynamic simulations and explain our sample selection. Section 3 provides details on the shape reconstruction method together with validations of the model assumptions and accuracy tests. The method is then applied on the COSMOS data in Section 4. A summary of our results can be found together with our conclusions in Section 5.

# 2 DATA

## 2.1 COSMOS observations

The 2 deg<fnref n="2" /> COSMOS field (Scoville et al. 2007) has been observed extensively by different ground and space based telescopes, including Hubble, Spitzer, VISTA, CFHT, and Subaru. The joint analysis of these observations led to different catalogues providing estimates of galaxy properties such as stellar masses, star formation rates and morphological characteristics over a large range in redshift and luminosity. Our analysis is based on three of these catalogues, which are described below.

### 2.1.1 Photometry

The public COSMOS2015 catalogue<fnref n="1" /> (Laigle et al. 2016) comprises photometry in 30 bands, covering ultra-violet to mid-infrared wavelengths. In our analysis, we use redshift, stellar mass, and specific star formation rate (sSFR)<fnref n="2" /> estimates provided in this catalogue, which were derived for each galaxy by fitting templates of spectral energy distributions (SEDs) to the photometric data (Ilbert et al. 2006). We discard objects which are i) residing in regions flagged as ‘bad’ (mostly because they are close to stars or to the edge of the field), ii) saturated, and iii) not classified as galaxies by requiring the corresponding catalogue flags to be [flag_hjmcc, flag_peter, type] = [0,0,0]. We further impose cuts at the limiting AB magnitudes in the near-infrared Ks-band of 24.0 and 24.7 in the deep and ultra ultra-deep fields, respectively. These magnitudes are defined within a fixed 3'' diameter aperture ($K_s M A G_A P E R 3$) and the limits correspond to a $3\sigma$ detection. After applying these cuts, the catalogue contains

252 527 galaxies. For this sample, the standard deviation of the relative differences in redshift with respect to the zCOSMOS-bright spectroscopic control sample (Lilly et al. 2007) is $\sigma_{\Delta z/(1+z_s)} = 0.007$ (catastrophic failure fraction $\eta = 0.5$ per cent). Note that this error is an optimistic estimate, since it was inferred by comparison to spectroscopic redshifts of bright galaxies and is higher for dimmer objects. However, due to the bright magnitude cuts used in this work (see Section 2.1.5), we expect this inaccuracy to be a realistic estimate for the galaxies in our samples. The accuracy on the stellar mass and star formation rate estimates is expected to be ~0.1 dex and ~0.2–0.6 dex, respectively at $z < 1.5$ (Laigle et al. 2019). These values may be lower for our samples as we focus on bright objects only. Nevertheless, Laigle et al. (2019) find evidence that the strong scatter (and bimodality) in the SFR estimates is mainly driven by the inaccurate modelling of dust extinction within the galaxies at the SED-fitting stage (see their Figure B3). We also find strong evidence of dust extinction in the galaxies from our sample (e.g. Fig. D1). We therefore study the sSFR only for discs with high apparent axis ratios, which we consider to be inclined towards a face-on orientation at which the impact of dust extinction (and dust extinction modelling) is expected to be minimal.

### 2.1.2 Shapes

We use galaxy shape estimates from the public Advanced Camera for Surveys General Catalog (ACS-GC,<fnref n="3" /> Griffith et al. 2012). This catalogue is based on Hubble Space Telescope (HST) imaging in the optical red $I_{AB}$ broad-band filter F814W. The absence of atmospheric distortions allows for an excellent image resolution, which is mainly limited by the width of the HST point spread function (PSF) of 0.085'' in the F814W filter and the pixel scale of 0.03''. Sources were detected using the GALAPAGOS software (Häußler et al. 2011), which runs SExTRACTOR (Bertin & Arnouts 1996) and GALFIT (Peng et al. 2002) in two subsequent steps. Galaxy shapes are described by the two-dimensional major over minor axis ratios $q_{2D}$, which are derived by GALFIT from fits of a single Sérsic model to each objects image. The surface brightness in this model is given by

$$ \Sigma(r) = \Sigma_e \exp \left\{ -\kappa_n \left[ \left( \frac{r}{r_e} \right)^{1/n} - 1 \right] \right\}, \quad (3) $$

where $\Sigma_e$ is the surface brightness at the effective radius $r_e$ and the parameter $\kappa_n$ is chosen such that $r_e$ encloses half of the total flux. The Sérsic index $n$ quantifies the concentration of the surface brightness profile. The 2D axial ratio $q_{2D}$ (provided as BA_GALFIT HI in the catalogue) enters equation (3) via $r = \sqrt{x^2 + (y/q_{2D})^2}$, where $x$ and $y$ are the coordinates on the major and minor axis, respectively. The modelled surface brightness profiles are convolved with the ACS PSF before being compared to the reference observation during the fit. The morphological parameters from GALFIT hence describe the intrinsic 2D galaxy shapes and do not require further PSF correction. We select objects from the catalogue which were classified by SExTRACTOR as galaxies (CLASS_STAR HI < 0.1) and with good fits to the GALFIT model (FLAG_GALFIT HI = 0 and CHI2NU HI < 2). From the remaining sample we reject 65 objects which have axis ratios equal to zero or larger than unity or effective radii from GALFIT above 750 ACS pixels (22.5''). After these cuts the final catalogue contains 128 365 objects.

<footnote n="1">¹https://www.eso.org/qi/</footnote>

<footnote n="2">²We are aware that this sSFR estimates is not very robust, but we checked that it is sufficient for the purpose of our analysis.</footnote>

<footnote n="3">³vizier.u-strasbg.fr/viz-bin/VizieR-3?-source = J/ApJS/200/9/acs-gc</footnote>
### 2.1.3 Morphological classification

For the selection of disc-dominated galaxies we use the Zurich Structure & Morphology Catalog<fnref n="4" /> which is derived from the same HST imaging data as the ACS-GC. It provides a morphological classification for each galaxy, derived with the Zurich Estimator of Structural Type and is referred to as ZEST catalogue in the following. The ZEST classification is based on a principal component analysis of five non-parametric diagnostics: asymmetry, concentration, Gini coefficient, 2nd-order moment of the brightest 20 per cent of galaxy pixels (M<sub>20</sub>) and ellipticity (see Scarlata et al. 2006, 2007; Sargent et al. 2007, who also provide validations of their classification). The catalogue contains galaxies in the COSMOS field brighter than I<sub>AB</sub> = 24. We select objects, which i) are classified as galaxies ([acs mu class, stellarity] = [1,0]), ii) do not reside in automatically or manually masked regions ([acs mask, acs masked] = [0,1]) and iii) are not flagged as unusable or spurious ([acs clean, junk-flag] = [1,0]). After applying these conditions the remaining sample contains 108 800 galaxies. The morphological classification is considered to be unreliable for galaxies with half-light radii smaller than twice the size of the ACS F814W PSF (i.e. 0.17", that is ≳5.6 times larger than the ACS pixel size), which we take into account in our sample selection (Section 2.1.5). We validate the morphological classification for the disc-dominated galaxies used in our analysis in Section 2.1.5 and in Appendix A.

### 2.1.4 Matched catalogue

Matching objects in these three COSMOS catalogues is not straightforward due to variations in common properties such as positions and magnitudes, which can lead to spurious mismatches. These variations can originate from atmospheric distortions in the ground-based COSMOS2015 data, differences in the employed image analysis software (i.e. SExTRACTOR and GALFIT) and its configuration, in the employed telescopes, cameras and filters, as well as different quality cuts applied before matching. In order to minimize the chance for mismatches, galaxies are matched based on angular positions as well as on magnitudes. We start by matching objects in the COSMOS2015 and ACS-GC catalogue in three steps. 1) We select pairs of galaxies as candidate matches if their angular separation is <0.6", which is slightly below the typical seeing of the ground-based telescopes contributing to the COSMOS survey. 2) We discard candidate matches with a difference in brightness of more than 1.0 magnitude. 3) We finally select the matches as those with the smallest difference in angular positions. The COSMOS2015 magnitudes used in step 2) are measured in the Subaru *i*+ band within a fixed 3" aperture (referred to as ip\_mag\_aper3 in the catalogue). They are compared to the SExtractor magnitudes MAG\_BEST\_HI from ACS-GC which are defined as magnitudes measured within a flexible elliptical aperture (referred to as MAG\_AUTO) or corrected isophotal magnitudes if contaminating sources are located in the vicinity. The average wavelengths weighted by transmission in the Subaru *i*+ and ACS F814W filters are 7683.88 Å and 8073.43 Å, respectively, while the width and shape of their transmission curves differ significantly. These differences are accounted for by the relatively large tolerance in magnitude. We identify 98 604 objects in the matched COSMOS2015 and ACS-GC catalogue. From these

matched objects 50 (95) per cent have less than 0.07 (0.22)" and 0.12 (0.42) magnitude differences in their matched angular positions and luminosities, respectively which is well below the chosen tolerances described above and indicates that the match is robust. Subsequently we matched the joint COSMOS2015 and ACS-GC catalogue with the ZEST catalogue using the same three-step method with the same tolerances for magnitudes and angular positions. The ACS-GC MAG\_BEST\_HI are now compared to the SExtractor ACS\_MAG\_AUTO provided in ZEST.

The final matched catalogue contains 707 080 objects from which 50 (95) per cent have less than 0.03 (0.2)" and 0.03 (0.16) magnitudes differences in their matched angular positions and luminosities, respectively. These smaller differences compared to the first matching between COSMOS2015 and ACS-GC can be attributed to the fact that the positions and magnitudes used for the second matching are all derived with SExtractor from the same HST ACS imaging data. A second reason is the cut at I<sub>AB</sub> = 24 in ZEST, which excludes the dimmest objects with highest uncertainties on position and luminosity estimates. This latter cut also explains the strong drop in the number of objects in the final catalogue.

### 2.1.5 Volume limited main sample of disc-dominated galaxies

Our study is focused on disc-dominated late-type galaxies in the matched catalogue, which we identify as those with ZEST parameters type = 2 and bulg = 2 or 3. Among them, we select a volume limited main sample adopting cuts in photometric redshift, absolute Subaru *i*+ magnitude and comoving effective radius. The selection is displayed in Fig. 2 and described below. Examples of the galaxies in our main sample are shown in Figs 1 and A1.

**Redshift cuts** The redshift range of our main sample is set to 0.2 < z < 1.0 (marked by vertical red dashed lines in Fig. 2). The upper limit is a compromise between a deep selection in redshift and a sufficiently high number of galaxies in the volume limited sample. The lower limit is defined by the end of the distribution, below which very few objects are found due to the small volume of the light cone. The redshifts on which the cuts are applied are the median of the likelihood distribution from the SED template fits, which are referred to as photoz in the COSMOS2015 catalogue.

**Magnitude cuts** Galaxies are further selected to have apparent isophotal AB magnitudes in the Subaru *i*+ band (referred to as ip\_MAG\_ISO in the COSMOS2015 catalogue and hereafter as *m*<sub>i</sub>) brighter than *m*<sub>i</sub><sup>max</sup> = 24 above which the ZEST morphological classification becomes unreliable (Scarlata et al. 2007). The apparent magnitude cut introduces a redshift dependent selection by the absolute magnitude which hampers the comparison of galaxy populations at different redshifts (see bottom panel of Fig. 2). We therefore require the absolute restframe Subaru *i*+ magnitudes (hereafter referred to as *M*<sub>i</sub>) to be brighter than *M*<sub>i</sub><sup>max</sup> = −21.5, ensuring that all objects in our redshift range are sufficiently bright to be unaffected by the apparent magnitude cut.

Note that the cut in absolute magnitude is chosen to be brighter than a naive cut at $M_i^{\max} = m_i^{\max} - DM(z_{\max} = 1.0) = 20.1$ (indicated as blue dotted line in the bottom panel of Fig. 2) to mitigate effects of dust extinction on the observed axial ratio distribution. In the top panel of Fig. 3, we show that the apparent magnitude of discs is fainter for objects with low-axis ratios. This effect can be expected from an increased dust extinction in objects that are inclined towards an edge-on orientation (e.g. Graham & Worley 2008). As a consequence, high redshift discs (at $z \simeq 1.0$) that have low-axis ratios ($q_{2D} \le 0.5$) and where selected by $M_i^{\max} = -20.1$ can fall below the apparent

<footnote n="4">⁴irsa.ipac.caltech.edu/data/COSMOS/tables/morphology/cosmos_morph.zurich.1.0.tbl</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Random examples of ACS images of late-type galaxies in our volume limited COSMOS sample in the redshift range 0.2 < z < 0.4. The galaxies are selected to be disc-dominated according to the ZEST morphological classification scheme (Section 2.1.3). Galaxies in the top and bottom panels are further classified by ZEST as being face-on and edge-on oriented, respectively. The three columns on the left (right) show galaxies in our low (high) mass sample, with stellar masses below (above) $M_{\star}^{\text{cut}} = 10^{10.35} M_{\odot}$. The 2D axis ratios $q_{2D}$, provided in the ACS-GC catalogue, were obtained from single Sérsic profile fits (Section 2.1.2). The image sizes are adjusted for each galaxy.</figcaption></figure>

magnitude cut of $m_l^{\text{max}} = 24$, marked as horizontal dashed line in Fig. 3. Dust extinction can therefore introduce a bias in the observed axial ratio distribution towards apparently rounder (i.e. face-on) galaxies, as we discuss in more detail in Appendix D. In the bottom panel of Fig. 3, we demonstrate that selecting galaxies with absolute magnitudes brighter than $M_l = -21.5$ ensures that the apparent magnitudes are brighter than $m_l^{\text{max}} = 24$, even if the galaxies are seen edge-on ($q_{2D} \ll 1$) at the highest redshifts considered in this work.

**Size cuts** Our final selection addresses the sizes of galaxies. We require objects in our sample to have angular radii $R_{\theta}$ which are at least twice as large as the width of the PSF in the HST F814W imaging (i.e. $R_{\theta}^{\text{min}} = 0.17''$) to minimize the effect of inaccuracies in the PSF correction on shape measurements as well as PSF effects on the morphological classifications (Scarlata et al. 2007; Griffith et al. 2012). Simply applying a cut on $R_{\theta}$ would lead to objects with smaller physical sizes entering the main sample at lower redshifts. This effect could introduce an apparent evolution of the galaxy shapes with redshift since the apparent physical sizes and shapes are correlated (see top panels of Fig. 4 and discussion below). We therefore apply a cut on the transverse comoving radii $R_{\perp}(z) = R_{\theta}D_{A}(z)$, which ensures that the observed angular radii of the galaxies in our sample are always larger than $R_{\theta}^{\text{min}}$, even for the most distant objects at $z_{\text{max}} = 1.0$. Here $D_{A}(z)$ is the angular diameter distance at the source redshift $z$. Assuming a flat $\Lambda$CDM universe we obtain the condition $R_{\perp} > R_{\theta}^{\text{min}} D_{A}(z_{\text{max}}) = 0.64$ kpc. $D_{A}$ and DM are computed using cosmological parameters from the Planck Collaboration (2018), i.e. $(\Omega_{M}, \Omega_{b}, H_{0}, \sigma_{8}, n_{s}) = (0.31, 0.049, 67.66, 0.81, 0.9665)$.

We show $R_{\perp}$ corresponding to one and two times the angular width of the HST PSF ($0.085''$) at a given redshift as dotted and dash-dotted blue lines respectively in the top panel of Fig. 2 together with the cut on $R_{\perp}$ which defines the main sample (red dashed horizontal line). The angular radii $R_{\theta}$ are the PSF corrected effective radii from the GALFIT single Sérsic model fits, which quantify the galaxy size along the projected major axis (referred to as RE_GALFIT HI in the ACS-GC).

For opaque discs with zero thickness the size of the projected major axis would be equal to the three-dimensional major axis, independent of its inclination with respect to the plane of projection. In that case a cut by projected size would not introduce a cut by shape or inclination with respect to the observer. However, real

disc galaxies have a non-zero thickness and are not opaque but to a certain degree transparent. Their observed surface brightness profile is hence affected by projection effects, which depend on the inclination angle with respect to the observer. As a consequence, the diameter of the observed 2D isophotes is larger for edge-on than for face on discs (e.g. Holmberg 1946; de Vaucouleurs 1959; Heidmann, Heidmann & de Vaucouleurs 1972). In the top panels of Fig. 4, we show for galaxies brighter than our absolute magnitude cut that $R_{\perp}$ increases significantly with decreasing apparent axis ratios, indicating the expected dependence of $R_{\perp}$ on the disc inclination. However, for our bright sample and within the redshift range we consider, the lower limit on $R_{\perp}$ is sufficiently small to have only a negligible impact on the observed distribution of axis ratios as shown in the bottom panels of Fig. 4. We therefore do not expect a relevant impact of the size cut on the disc inclinations of our sample, which could otherwise bias the axial ratio distribution towards edge-on objects (e.g. Huizinga & van Albada 1992). All cuts defining the volume limited sample are summarized in Table 1. After applying these cuts, the final volume limited main sample contains 3739 disc-dominated galaxies in total. The redshift distributions of these galaxies is shown in Fig. 5.

**Type and colour distributions** In Fig. 6, we display the relative abundance of the disc-dominated galaxies in our main sample with respect to the total galaxy population in that sample in three redshift bins. These abundances are compared to those of galaxies classified as discs with large bulges (ZEST parameter type = 2 and bulg = 0 or 1), as elliptical (type = 1) and irregular (type = 3). The figure confirms reports from the literature that a large fraction of discs galaxies has no significant bulge. The fraction of discs remains roughly constant at $\le 80$ per cent, with a slight increase with redshift. It is further interesting to note that at $z \simeq 1.0$, the majority of discs has no large bulge, while the opposite is the case at $z \simeq 0.4$. The decline of the fraction of bulgeless discs with decreasing redshift lines up with the findings from Sachdeva (2013) based on Chandra Deep Field observations. Fig. 6 further shows that the rate at which this fraction declines is similar to the rate at which the fraction of discs with large bulges increases. This finding is consistent with the scenario in which discs with large bulges form out of bulgeless discs, potentially during mergers. Note, the absolute values of the different fractions may be specific to our
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Selection of the volume limited main sample from disc-dominated galaxies in the matched COSMOS catalogue by photometric redshift *z*, transverse comoving radii $R_\perp$, and apparent and absolute *i* band magnitudes $m_i$, and $M_i$, respectively. The cuts on each variable are marked by red dashed lines, enclosing the selected sample in the red area. The blue dotted and dash-dotted lines in the top panel show the comoving radii which correspond to one and two times the angular size of the HST PSF, respectively. The blue dash-dotted line in the bottom panel shows the limit on $M_i$, given by the apparent magnitude cut $m_i^{\text{max}}$ and the distance modulus DM. The horizontal blue dotted line in the same panel shows a naive cut on $M_i$ (see Section 2.1.5 for details). The over-abundances at $z \sim 0.35$ and $z \sim 0.7$ are known large-scale structures (see e.g Fig. 4 in Guzzo et al. 2007; Gozaliasl et al. 2019), while some sharp vertical spikes in the redshift distribution might be artefacts of the photometric redshift estimation.</figcaption></figure>

sample selection. In particular, the lower limit on the size separates out many elliptical galaxies, which tend to be more compact than discs, increasing the fraction of discs in our sample.

For further validation of our data set, we compare the colour–colour diagrams of the disc-dominated and elliptical galaxies in our main sample in Fig. 7. We find that the discs and the ellipticals reside preferentially in the blue star forming and red quenched sequence respectively, which are well separated by the colour–colour cuts from Laigle et al. (2016). This result can be expected in general from the well known correlation between galaxy morphology and colour (e.g. Larson, Tinsley & Caldwell 1980; Strateva et al. 2001; Baldry et al. 2004; Martig et al. 2009). We find that the bulge dominated discs in

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Apparent 2D axis ratios $q_{2D}$ of disc-dominated galaxies in COSMOS versus their apparent *i*-band magnitude $m_i$ at the highest redshifts considered in this work. Horizontal bars show the average $m_i$ in bins of $q_{2D}$ to elucidate the correlation between both observables. Face-on discs with $q_{2D} \simeq 1$ appear brighter than edge-on discs, which indicates a dependence of dust extinction on the disc inclination relative to the observer. Top panel shows the volume limited sample selected as shown in Fig. 2, but with a more conventional cut on the absolute *i* band magnitude of $M_i < -20.1$. The $m_i$ distribution of this sample is cut off at the apparent magnitude limit at $m_i = 24$, which biases the $q_{2D}$ distribution, as illustrated in Fig. D1. We mitigated this bias with a more conservative cut of $M_i < -21.5$. The resulting sample, shown in the bottom panel, is only weakly affected by the cut on $m_i$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top: Comoving transverse effective radii $R_\perp$ of disc-dominated galaxies in our volume limited COSMOS sample verses the apparent 2D axial ratio $q_{2D}$. Red bars mark the average $R_\perp$ in bins of $q_{2D}$ and show that the apparent size depends on the apparent axial ratio, presumably as an effect of the disc inclination. The horizontal dashed lines show minimum value for $R_\perp$ used as cut in our sample selection. Bottom: Distribution of $2D$ for disc galaxies before and after the size cut (black dashed and red solid lines, respectively). This figure demonstrates that our size cut does not bias the $q_{2D}$ distribution of our sample.</figcaption></figure>
<table><caption>**Table 1.** Selection cuts for the volume limited COSMOS main sample, shown in Fig. 2.</caption><thead><tr><th>galaxy property</th><th>constraint</th></tr></thead><tbody><tr><td>photometric redshift</td><td>0.2 &lt; z &lt; 1.0</td></tr><tr><td>apparent 3" aperture AB Subaru i+ magnitude</td><td>m<sub>i</sub> &lt; 24</td></tr><tr><td>absolute rest-frame Subaru i+ magnitude</td><td>M<sub>i</sub> &lt; -21.5</td></tr><tr><td>transverse comoving effective radius</td><td>R<sub>⊥</sub> &gt; 0.64 kpc</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Selection of redshift sub-samples from our volume limited main sample of disc-dominated galaxies in COSMOS. The photometric redshift distribution is split into three redshift bins, which are enclosed by vertical red dashed lines at $z = 0.2, 0.7, 0.9$, and $1.0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Fractions of galaxy types in our volume limited COSMOS sample in the three redshift bins shown in Fig. 5. Blue circles, red triangles, and green crosses show results for galaxies classified as discs, ellipticals, and irregulars, respectively. Results for discs with and without a large bulge are connected by dash-dotted and dashed lines, respectively.</figcaption></figure>

our sample populate both, the quenched as well as the star forming sequences (not shown in the figure for clarity). This result lines up with the large bulges found in discs on the red quenched sequence, for instance in data from COSMOS and the Sloan Digital Sky Survey (see Bundy et al. 2010; Guo et al. 2020). Overall, the expected correlation between colour and morphology shown in Fig. 7 confirms *a posteriori* that the morphological classification from ZEST and the photometric properties from the COSMOS2015 are consistent with each other. This test shows that the matching between both catalogues is sufficiently reliable for deriving physical interpretations from their combined morphological and photometric information.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Colour-colour diagram, based on estimates of the absolute rest-frame magnitudes in the *r*, *j*, and *NUV* filters from the COSMOS2015 catalogue for galaxies which pass our main sample selection on size, magnitude, and redshift. Red and blue dots show objects which are classified in the ZEST catalogue as early-type and disc-dominated late-type, respectively. The black dashed line is taken from Laigle et al. (2016) and separates the quenched and the star forming populations.</figcaption></figure>

### 2.1.6 Stellar mass – redshift samples

In order to investigate the redshift and stellar mass dependence of the galaxy shapes, we select three redshift sub-samples and two stellar mass sub-samples from the COSMOS main sample as shown in Fig. 5 and in the top panel of Fig. 8, respectively. For more detailed investigations we further select six sub-samples, which are selected by both redshift and stellar mass, as shown in the lower panels of Fig. 8. The width of the redshift bins are chosen such that each sample contains a sufficiently large number of objects required for our statistical analysis. The stellar mass cut at $M_{*}^{\text{cut}} = 10^{10.35}M_{\odot}$ corresponds to the median stellar mass of the main sample and lies close to the median stellar mass of the redshift sub-samples. The cuts in mass and redshift are summarized in Table 5 together with the number of galaxies in each sub-sample.

We show in the left-hand column of Fig. 9, the distribution of the observed shapes from the ACS-GC catalogue, quantified by the axial ratio $q_{2D}$, in the three redshift bins. We find that at all redshifts, the distributions of the low and high mass samples, are skewed towards high and low-2D axis ratios, respectively. The difference between the mass samples increases with decreasing redshifts, which is mainly driven by an increasing fraction of high mass galaxies with low-2D axis ratios ($q_{2D} \simeq 0.3$). The axial ratio distribution of low-mass galaxies on the other hand shows no strong change with redshift. In fact, we argue in Section 4 that the redshift dependence of both samples is not significant when considering shot-noise errors on the binned distributions.

Nevertheless, a weak trend in redshift can also be seen in the transverse comoving radii $R_{\perp}$ in the central column of Fig. 9. These radii are shown for discs inclined towards a face-on orientation with apparent axis ratios $q_{2D} > 0.5$ to reduce the impact of inclination on the observed size (see Fig. 4). As for the axis ratios, we find that at high redshift, the size distribution for low and high mass samples is similar. At lower redshift, the deviation between both distributions is
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Selection of low- and high-stellar mass sub-samples from the COSMOS main sample, defined over the entire redshift range, and from the three redshift sub-samples shown in Fig. 5 (top and lower panels, respectively). The stellar mass distributions are split at $M_{\star}^{\text{cut}} = 10^{10.35}$ (red dashed lines), which corresponds to the median stellar mass of the main sample and lies close to the median stellar masses of the redshift sub-samples (black dotted lines).</figcaption></figure>

slightly larger, mainly due to an increase of the sizes in the high mass disc sample. This trend can be seen in the average radii, displayed as dotted and dash-dotted vertical lines for the high and low mass sample, respectively.

A non-geometric galaxy property which follows the same behaviour is the specific star formation rate (sSFR) from the COSMOS2015 catalogue, whose distribution is displayed in the right-hand column of Fig. 9. The selection here is again restricted to galaxies with $q_{2D} > 0.5$ to mitigate a potential bias in the sSFR estimates induced by inaccuracies of dust attenuation model employed in the SED fitting, as reported by Laigle et al. (2019). We see that the sSFR distribution is similar for the low and high mass samples at high redshifts. At low redshifts, both distributions differ significantly, which is mainly driven by a strong decrease of the sSFR for high mass discs. This result lines up with those of Grossi et al. (2018), who find that the sSFR of disc dominated galaxies in COSMOS field decreases with increasing stellar mass and decreasing with redshift. It is worth noting that the mass and redshift dependence is weaker for the SFR (not shown here) than for the sSFR. This weaker dependence is consistent with the overall weak redshift dependence of star formation in nearby galaxies, suggested by findings of Kroupa et al. (2020) and may result from a partial compensation by the mass

differences between the sub-samples. Overall, our findings indicate a correlation of the galaxy shapes, sizes, and specific star formation rates with the stellar masses, which appear to increase as galaxy formation proceeds. We will discuss this result in more detail together with the reconstructed 3D axial ratio distribution in Section 5.

We stress here that our discussion of differences between the $R_{\perp}$ and sSFR distributions from different mass and redshift samples discussed above is based on our visual impressions. Clarifying whether or not these differences are significant would require a statistical analysis, which is beyond the scope of this work.

## 2.2 disc galaxies in hydro-dynamic simulations

We use the state-of-the-art hydrodynamic simulations of galaxy formation Horizon-AGN<fnref n="5" /> and Illustris TNG100⁶ (hereafter referred to as HAGN and TNG100, respectively) to test the methods and assumptions on which our analysis of the COSMOS data is based. These tests are performed at redshift $z = 1.0$, were the available snapshots of these simulations are closest to the maximum of the redshift distribution of our volume limited COSMOS sample. Both simulations assume a $\Lambda$CDM cosmology with recently constrained parameters, cover similar cosmological volumes and include sub-grid mechanisms to model the formation of galaxies at a similar degree of complexity, as detailed below. However, their simulation techniques differ significantly, which allows for testing the robustness of our conclusions. The main characteristics of these simulations are compared to each other in Table 2. Note that these simulations where run at relatively low resolutions in order to cover large volumes, which has a noticeable impact on the galaxy morphologies as discussed in Section 3.1.

### 2.2.1 Horizon-AGN

The hydrodynamic simulation HAGN was produced with the grid-based adaptive-mesh-refinement code RAMSES (Teyssier 2002), using cosmological parameters compatible with the constraints from WMAP-7 (Komatsu et al. 2011). The simulation includes the key processes relevant for galaxy formation: cooling, heating, and chemical enrichment of gas, the formation and evolution of stars and black holes as well as feedback from stellar winds, supernovae and Active Galactic Nuclei (AGNs) (Dubois et al. 2014). Galaxies were identified in the distribution of stellar particles as groups with more than 50 members using the AdaptaHOP finder (Aubert, Pichon & Colombi 2004). A key achievement of HAGN with respect to previous generations of cosmological hydrodynamic simulations is the broad diversity of galaxy morphologies, with realistic fractions of discy and elliptical objects (Dubois et al. 2016). A reasonable agreement with observed luminosity functions and colour distributions has been shown by Kaviraj et al. (2017), while moderate deviations from observed angular clustering have been found by Hatfield et al. (2019).

For our analysis, we select disc galaxies, which we define via the ratio $r_v = V_{\text{rot}}/\sigma_v$ between the stellar rotation $V_{\text{rot}}$, defined as the mean tangential velocity of stellar particles with respect the galaxies’ spin axis and the velocity dispersion $\sigma_v$ (see Chisari et al. (2015) for details). High values of $r_v$ indicate relatively high rotational velocities compared to those radial or parallel to the spin axis. For our analysis we select the 20 per cent of galaxies with the highest values of $r_v$,

<footnote n="5">⁵www.horizon-simulation.org</footnote>

<footnote n="6">⁶www.tng-project.org</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Probability distributions of apparent 2D axis ratios, transverse comoving radii and specific star formation rates (left, central, and right-hand panels, respectively) for disc-dominated galaxies in our volume limited COSMOS main sample with stellar masses below and above $10^{10.35}M_{\odot}$ (blue and red histograms, respectively). Vertical panels show results in three redshift bins with limits indicated on the right. The comoving radii and specific star formation rates are shown for discs with apparent axis ratios above $q_{2D} > 0.5$, to minimize bias induced by projection and dust extinction (see Fig. 4 and D1). Vertical black dashed lines mark the minimum radius of the volume limited sample. Vertical solid and dotted lines indicate the mean radii for each mass sample over the full redshift range and in each redshift bin, respectively.</figcaption></figure>

<table><caption>**Table 2.** Main properties of the hydro-dynamic simulations Illustris TNG 100 and Horizon-AGN.</caption><thead><tr><th></th><th></th><th>TNG100</th><th>HAGN</th></tr></thead><tbody><tr><td>Ω<sub>Λ</sub></td><td>–</td><td>0.6911</td><td>0.728</td></tr><tr><td>Ω<sub>m</sub></td><td>–</td><td>0.3089</td><td>0.272</td></tr><tr><td>Ω<sub>b</sub></td><td>–</td><td>0.0486</td><td>0.045</td></tr><tr><td>H<sub>0</sub></td><td>[s<sup>-1</sup>km]</td><td>67.74</td><td>70.4</td></tr><tr><td>σ<sub>8</sub></td><td>–</td><td>0.8159</td><td>0.81</td></tr><tr><td>n<sub>s</sub></td><td>–</td><td>0.9667</td><td>0.967</td></tr><tr><td>L<sub>box</sub></td><td>[h<sup>-1</sup>Mpc]</td><td>75</td><td>100</td></tr><tr><td>m<sub>*</sub></td><td>[M<sub>☉</sub>]</td><td>–</td><td>2 × 10<sup>6</sup></td></tr><tr><td>m<sub>baryon</sub></td><td>[M<sub>☉</sub>]</td><td>1.4 × 10<sup>6</sup></td><td>–</td></tr><tr><td>m<sub>dm</sub></td><td>[M<sub>☉</sub>]</td><td>7.5 × 10<sup>6</sup></td><td>8 × 10<sup>7</sup></td></tr></tbody></table>

which corresponds to a cut at $r_v > 1.06$ at redshift $z = 1.0$ (Fig. 10). We attribute this low value of $r_v$ to the relatively low resolution of the simulation and demonstrate in Section 3.1 that the selected galaxies are indeed discy, although with ‘puffed-up’ shapes, characterized by relatively high $C_{3D}/B_{3D}$ axis ratios. We further limit the selection to galaxies with more than 500 particles, which corresponds to a stellar mass cut at $M_{*} > 10^{9}M_{\odot}$, to ensure reliable measurements of the morphological properties (see also Section 2.2.3). The final sample contains 14198 disc galaxies.

### 2.2.2 Illustris TNG100

The TNG100 (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018b; Springel et al. 2018) is a magneto-hydrodynamic simulation, produced with the moving-mesh code AREPO (Springel 2010), which was run with cosmological parameters from the Planck Collaboration (2016). It includes the same key processes for modelling galaxy formation as HAGN,

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Left: Distribution of the ratio between stellar rotation $V_{rot}$ and the velocity dispersion $\sigma_v$ for galaxies in the HAGN simulation. Right: Distribution of the fraction of disc particles in each galaxy of the TNG100 simulations. The populations of galaxies which we classify as discs are marked in red at the tails of the distributions and make up 20 per cent of the entire sample. Results are shown for $z = 1.0$.</figcaption></figure>

although with significantly different implementations which are described in Pillepich et al. (2018a). Galaxies are identified in dark matter subhaloes with non-zero stellar components. These subhaloes are detected in friends-of-friends groups by the SUBFIND algorithm (Davis et al. 1985; Springel et al. 2001; Dolag et al. 2009). While it has been calibrated to match the galaxy mass function at $z = 0$, the simulation has been shown to predict main characteristics of observed galaxy populations reasonably well, such as morphological diversity (Rodriguez-Gomez et al. 2019), stellar mass functions at higher redshift (Pillepich et al. 2018b), the colour bimodality (Nelson et al. 2018) as well as colour-dependent two-point clustering statistics (Springel et al. 2018). The basic properties of the simulation are compared to those of the HAGN simulation in Table 2.

To identify discs in the TNG100 simulation, we use each galaxy’s fraction of stellar mass in the disc component with respect to the total stellar mass. For this purpose stellar disc particles have been defined
via their circularity parameter $\epsilon \equiv J_z/J(E)$, where $J_z$ is the specific angular momentum around a selected z-axis and $J(E)$ is the maximum specific angular momentum possible at the specific binding energy $E$ of a given stellar particle.

The z-axis is thereby given by the principle angular momentum vector of the star forming gas, or the stellar particles, if there is no star forming gas in the system (Vogelsberger et al. 2014; Genel et al. 2015). Particles with $\epsilon > 0.7$ are considered to belong to the disc component (Abadi et al. 2003). The circularity parameter is provided on the Illustris data base for all subhalos with $M_* > 3.4 \times 10^8 M_\odot$ within the half-mass diameter $2R_{1/2,*}$ and at least 100 stellar particles. We classify TNG100 galaxies as discs if their fraction of disc particles is above 0.35, which corresponds to the upper 20 per cent of the distribution at $z = 1.0$, as shown in Fig. 10. As for HAGN, we attribute this low fraction of disc particles used for the cut to the limited resolution of the simulation and demonstrate in Section 3.1 that it is a reasonable choice for selecting discy objects. Note that we choose the disc particle fraction for the morphological classification in TNG100 since this quantity is provided for the $z = 1.0$ snapshot on the public TNG data base in contrast to $r_v$. Besides this practical motivation, using two different morphological classification schemes in both simulations allows for drawing more robust conclusions regarding the validation of our analysis methods. The final sample contains 5 674 disc galaxies.

### 2.2.3 Axis ratio measurements

The galaxies' axis ratios provided for the HAGN and TNG100 simulations are computed from the moment of inertia of their stellar mass distributions

$$I_{i,j} = \frac{1}{M_*} \sum_{n}^{N_*} m_*^n r_i^n r_j^n, \qquad (4)$$

where $N_*$ is the number of stellar particles in the galaxy, $m_*^n$ is the stellar mass of the $n^{th}$ particle, $M_* = \sum_{n}^{N_*} m_*^n$, and $r_i^n$ are the components of the particle position vectors, defined with respect to the centre of mass (Chisari et al. 2015; Genel et al. 2015).

The moment of inertia can be defined via the particle positions in 3D as well as in 2D. In the latter case, the positions are projected along one coordinate axis of the simulation, assuming an observer at infinity. In the 3D case, the square roots of the three absolute eigenvalues $\lambda_1 \ge \lambda_2 \ge \lambda_3$ provide a measure of the 3D major, intermediate and minor axis lengths, respectively, i.e. $(A_{2D}, B_{2D}, C_{2D}) = \sqrt{(\lambda_1, \lambda_2, \lambda_3)}$. Accordingly, the 2D major and minor axis lengths are given by $(A_{2D}, B_{2D}) = \sqrt{(\lambda_1, \lambda_2)}$, respectively. The axis ratios are defined according to equation (1) and (2). We expect that the bias on such axial ratio measurements, which results from the discreteness of the particle distributions (e.g. Joachimi et al. 2013; Hoffmann et al. 2014; Chisari et al. 2015), is negligible, due to the lower mass limits imposed on our sample (see Section 2.2.1 and 2.2.2). The 2D axis ratios of the projected stellar mass distributions in HAGN are used to test our method for reconstructing the 3D axial ratio distribution. For TNG100 such measurements are currently not publicly available. However, we use 2D axis ratios measured from second-order moments in synthetic images from TNG100 galaxies from Rodriguez-Gomez et al. (2019). These images have been produced using the SKIRT radiative transfer code<fnref n="7" /> (Baes et al. 2011; Camps & Baes 2015) by the Sloan Digital

Sky Survey in the g and i broad-band filters. The axis ratios are obtained from the second-order moments of the pixels in these images that belong to the source. Since these measurements are not PSF-corrected, we do not use them for testing the shape reconstruction method. Instead, we use them to obtain a rough estimation of how strongly the observed galaxy shapes depend on the wavelength range in which they are measured in order to interpret the COSMOS observations in Appendix E.

## 3 RECONSTRUCTING 3D GALAXY SHAPE DISTRIBUTIONS FROM 2D OBSERVATIONS

We now present the method used for constraining the distribution of 3D galaxy axis ratios from the observed distribution of 2D axis ratios in COSMOS and test this method using galaxy shapes from the TNG100 and HAGN simulations. The method is based on the assumption that each galaxy can be represented by an absorption-free, self-similar, coaxial ellipsoidal stellar system, to which we refer to as 3D ellipsoid in the following. The shape of its 3D luminosity density is fully described by the two axis ratios $q_{3D}$ and $r_{3D}$, given by equation (1). We expect such a one-component model to be a simplistic, but useful description for the objects in our observed sample: disc-dominated galaxies, discarding bulge-dominated, and irregular objects (see Section 2.1.3). However, sub-structures such as spiral arms could bias our results (see Fig. 1). The isodensity contours of the projected luminosity density (i.e. the isophotes) of such model galaxies are self-similar, coaxial ellipses (Stark 1977), whose 2D axis ratios can be obtained from the 3D axis ratios analytically, as detailed in Appendix B. This allows for an efficient prediction of the 2D axial ratio distribution, $P(q_{2D})$, for an ensemble of randomly oriented model galaxies with a given distribution of 3D axis ratios, $P(q_{3D}, r_{3D})$. This distribution can hence be constrained by comparing the corresponding $P(q_{2D})$ prediction to observations. This approach relies on a physically meaningful model for $P(q_{3D}, r_{3D})$, which we motivate in Section 3.1. The free model parameters are obtained from the COSMOS data using Bayesian interference as detailed in Section 3.2 and 3.3. In the latter section we also study the impact of inaccuracies of the employed $P(q_{3D}, r_{3D})$ model on the inferred 3D axial ratio distribution.

### 3.1 Model for the 3D axial ratio distribution

The reconstruction of 3D axial ratio distributions relies on a physically meaningful model for those distributions. Several of such models have been proposed in the literature. Sandage et al. (1970) found that the $q_{2D}$ distribution of spiral galaxies can be fitted reasonably well with a simple oblate disc model according to which $q_{3D} = 1$ for all objects while $s_{3D}$ is normal distributed around $\langle s_{3D} \rangle \simeq 0.25$. Later studies based on larger samples found that the $q_{2D}$ distribution of spirals is better fitted using slightly triaxial disc models which describe the absence of perfectly circular face-on spirals in observations (e.g. Binney & de Vaucouleurs 1981; Fasano et al. 1993). Using normal distributions for $q_{3D}$ and $s_{3D}$ Lambas et al. (1992) obtained good fits to the observations from the APM Bright Galaxy Survey. The good performance of this model has been confirmed by recent results from Satoh et al. (2019) based on COSMOS data. Ryden (2004) found that a normal distribution for $s_{3D}$ and a log-normal distribution for the ellipticity $\epsilon_{q_{3D}} \equiv 1 - q_{3D}$ delivers good fits to the $q_{2D}$ distribution of spirals in SDSS observations, which was confirmed by later studies (e.g. Padilla & Strauss 2008; Rodríguez & Padilla 2013). An alternative model based on normal distributions of the ellipticity $\epsilon_{s_{3D}} \equiv 1 - s_{3D}$

<footnote n="7">⁷https://skirt.ugent.be</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Distributions of 3D axis ratios for disc galaxies in hydrodynamic simulations at z = 1.0, obtained from the simple moment of inertia. Solid and dotted lines show fits to a normal and a skew normal distribution, respectively.</figcaption></figure>

and the triaxiality $T = (1 - q_{3D}^2)/(1 - s_{3D}^2)$ proved to match the $q_{2D}$ distribution from SDSS, 3D-HST, COSMOS, and CANDELS observations (Chang et al. 2013; van der Wel et al. 2014). These different models were motivated by the fact that their prediction for the 2D axis ratios (or ellipticities) fit observations reasonably well, which indicates that observations do not provide a strong constrain on the functional form of the 3D axes ratio distribution.

In our analysis, we use an approach that overcomes this limitation by employing a new model for the 3D axial ratio distribution which we validate using the HAGN and TNG100 simulations, described in Section 2.2. Since cosmological hydrodynamical simulations reproduce observed galaxy morphologies to varying degree of success and with a known dependence on simulation techniques, sub-grid models, and resolution (Snyder et al. 2015; Dubois et al. 2016; Correa et al. 2017; Park et al. 2019; Rodriguez-Gomez et al. 2019; Tacchella et al. 2019), we only use them as a tool to validate our reconstruction method.

In Fig. 11, we show the marginalized probability distributions of $q_{3D}$ and $r_{3D}$ for galaxies classified as discs in the HAGN and TNG100 simulations (see Section 2.2 for details). We find very similar distributions for both simulations, which is interesting to note given a) the significant differences in the physical models employed in the simulations and b) the different parameters used for the morphological classification. The latter appears to be reasonable, as the discs tend to be oblate with $r_{3D} < q_{3D} \lesssim 1$. Note that the absence of thin discs with $r_{3D} \ll 1$ can be attributed to the effective pressurized equation of state employed in these simulations and their relatively low resolution. Thin discs are present for instance in the TNG50 run (see Pillepich et al. 2019) or the New Horizon Simulation (Park et al. 2019), which are higher resolution, lower volume versions of the runs analysed here.

We fit the probability distributions with normal and skew normal distributions, given by $G(x) \propto \exp\{-\delta_x^2\}$, and $G_s(x) = G(x)(1 + \text{erf}\{\gamma\delta_x\})$, respectively, with $\delta_x = (x - x_0)/(\sqrt{2}\sigma)$.<fnref n="8" /> Both functions are truncated outside of the interval (0,1] to ensure that $\infty > A_{3D}$

$\ge B_{3D} \ge C_{3D}$. We show in Fig. 11 that the fit of the skew normal distribution matches the axis ratios slightly better than the normal distribution, in particular $q_{3D}$. However, we decide to neglect the skewness in our modelling, since the improvement in the 3D axial ratio fits is relatively small. Furthermore we show in Appendix C that a skewness on $q_{3D}$ has a very minor impact on the corresponding distribution of $q_{2D}$ and could hence be constrained very poorly by observations. In the same appendix, we also show that $q_{3D}$ and $r_{3D}$ are uncorrelated in HAGN as well as in TNG100. Based on these findings we approximate the joint 3D axial ratio distribution with the product $G(q_{3D})G(r_{3D})$, i.e.

$$ \tilde{P}(q_{3D}, r_{3D}) = \exp\left\{-\frac{1}{2}\left[\left(\frac{q_{3D} - q_0}{\sigma_q}\right)^2 + \left(\frac{r_{3D} - r_0}{\sigma_r}\right)^2\right]\right\}, \quad (5) $$

where $q_0, \sigma_q, r_0, \sigma_r$ are the free model parameters. The normalized truncated distribution is then given by

$$ P = \begin{cases} \tilde{P}_{3D}/N & \text{if } q_{3D}, r_{3D} \in (0, 1] \\ 0 & \text{else} \end{cases} \quad (6) $$

$$ \text{with } N = \int_{0}^{1} \int_{0}^{1} \tilde{P}_{3D}(q_{3D}, r_{3D}) dr_{3D} dq_{3D}. $$

## 3.2 Error estimation and parameter inference

We obtain constraints on our model parameter vector $\theta = (q_0, \sigma_q, r_0, \sigma_r)$ from the observed data using Bayesian inference. We assume a multivariate normal distribution of the likelihood for observing the data vector $\mathbf{d}$ given the parameters $\theta$,

$$ \ln \mathcal{L}(\mathbf{d}|\theta) = -\frac{1}{2} \chi^2(\mathbf{d}|\theta) + \text{const.} \quad (7) $$

with

$$ \chi^2(\mathbf{d}|\theta) = [\mathbf{d} - \mathbf{m}(\theta)]^T C^{-1} [\mathbf{d} - \mathbf{m}(\theta)]. \quad (8) $$

The data vector is the observed distribution $P(q_{2D})$, measured in 25 bins of equal width within the $q_{2D}$ interval (0,1] and $m(\theta)$ is the corresponding prediction from our $P(q_{3D}, r_{3D})$ model for the parameters $\theta$. The posterior distribution of the parameters $\theta$ given the data $\mathbf{d}$ is given by Bayes' theorem as

$$ P(\theta|\mathbf{d}) \propto \mathcal{L}(\mathbf{d}|\theta)\Pi(\theta), \quad (9) $$

where $\Pi(\theta)$ is the prior, which we set to be flat in the intervals (0,1] and [0,1] for the parameters ($q_0, r_0$) and ($\sigma_q, \sigma_r$), respectively. We thereby include parameter combinations in our sampling that describe axial ratio distributions not only for discs ($r_0 < < q_0$), but for any kind of ellipsoids, e.g. $r_{3D} \gtrsim q_{3D}$.

The covariance matrix $C$ is assumed to be given by $C_{ij} = \delta_{ij} \sigma_i^2$, with Poisson shot noise variance $\sigma_i^2 \propto N_i$, where $N_i$ are the counts of galaxies in a given bin $i$. The model prediction $m(\theta)$ is obtained by drawing a sample of 3D axial ratio pairs ($q_{3D}, r_{3D}$) from the $P(q_{3D}, r_{3D})$ distribution for a given set of parameters $\theta$. The corresponding $q_{2D}$ axis ratios are then computed for a random orientation, as detailed in Appendix B. The prediction for $P(q_{2D})$ is measured from the resulting $q_{2D}$ sample in the same bins as the observed data. For generating the predictions we choose a sample size of $10^5$, which is a compromise between a fast computation, needed for efficiently estimating the posterior, and having errors on the prediction which are negligible compared to those on the data vector. The latter condition is satisfied, as the number of galaxies in our six COSMOS sub-samples is more than two magnitudes smaller than the samples used for generating the predictions (see Table 5). We estimate $P(\theta|\mathbf{d})$ by sampling the parameter space with the Markov chain-Monte Carlo

<footnote n="8">⁸The fits are obtained by minimizing the $\chi^2$ deviation between model and data, assuming shot noise errors on the latter.</footnote>
(MCMC) method using the code `emcee`<fnref n="9" /> (Foreman-Mackey et al. 2013). For each posterior we run 32 independent chains with at last 1000 steps each. The best fit parameters are obtained from the position of the maxima of the marginalized posterior distribution.

### 3.3 Testing the reconstruction method

The method for reconstructing the distribution of 3D axis ratios $P(q_{3D}, r_{3D})$ from the distribution of 2D axis ratios $P(q_{2D})$ is now tested using the HAGN simulation. We begin by validating two assumptions on which this method is based. Those are 1) that the galaxies' 3D stellar mass isodensities are coaxial, self-similar 3D ellipsoids, whose 3D axis ratios can be related analytically to the 2D axis ratios of the projected 2D stellar mass isodensities as described in Appendix B, and 2) that our Gaussian model for the 3D axial ratio distributions from equation (6) is sufficiently accurate to provide good predictions for the 2D axial ratio distribution. Subsequently, we test how well the distribution of 3D axis ratios in the simulation can be recovered from the corresponding distribution of 2D axis ratios.

To test assumption 1) we compare in Fig. 12, the $q_{2D}$ distribution in the HAGN simulation, measured directly from the projected distributions of each galaxy's stellar particles as described in Section 2.2.3 (top panel), to the distribution of 2D axis ratios, obtained analytically from each galaxy's 3D axes ratios assuming ellipsoidal 3D stellar mass isodensities (central panel). We find that the distributions differ significantly from each other, which indicates that the assumption of ellipsoidal stellar mass isodensities can only serve as a very rough approximation. This conclusion can already be expected from the images of face-on late-type galaxies in our COSMOS sample, shown in Fig. 1. However, certain characteristics are present in both distributions, namely the cut-offs at $q_{2D} \simeq 0.4$ and $q_{2D} \simeq 0.9$, as well as a skewness in the distributions towards low-axis ratios. To test assumption 2) we show in the bottom panel of Fig. 12 the $q_{2D}$ distribution predicted for an ensemble of $10^6$ 3D axis ratios, drawn from our $P(q_{3D}, r_{3D})$ model fit to the HAGN simulation (see Fig. 11 and C2). We find that this predicted distribution is similar to the results obtained from the analytical projection of 3D axis ratios (shown in the central panel of Fig. 12), which indicates that our Gaussian model for $P(q_{3D}, r_{3D})$ from equation (6) is an appropriate approximation for our analysis.

We study the performance of the reconstruction method, starting with a self consistency test (a) in which we fit the $P(q_{3D}, r_{3D})$ model to match its own $P(q_{2D})$ prediction, using parameters from the HAGN fits shown in Fig. 12. We find that the fit, shown as green line in the bottom panel of Fig. 12, is in good agreement with the reference measurements. The MCMC estimate of the posterior distribution is displayed in Fig. 13 as light and dark green contours, indicating the 68 per cent and 95 per cent confidence intervals, respectively. We find that the parameter constraints are in good agreement with the input parameters of the $P(q_{3D}, r_{3D})$ model, shown as black dashed lines in Fig. 13. This result demonstrates that the 3D axial ratio distribution can be reconstructed from the corresponding distribution of 2D axis ratios for a sample of idealized disc galaxies, which satisfy the model assumptions outlined above. Note that this is not necessarily the case for elliptical galaxies, which are not subject of this work. In Appendix F, we show that the input parameters are also recovered for samples that are as small as the observational samples used on this work, but with larger confidence intervals.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Probability distribution of 2D axis ratios. Histograms in the top and central panels show results from the HAGN simulation, measured from the projected stellar mass distributions and analytic projections derived from 3D axis ratios, assuming coaxial, self-similar ellipsoidal stellar mass isodensities, respectively. The histogram in the bottom panel shows the 2D axial ratio distribution for a Gaussian model of the 3D axis ratio distribution with parameters from the 3D fits, shown in Fig. 11. In each panel, we show fits of the model prediction as solid lines. The corresponding parameter constraints and best-fitting values are compared in Fig. 13 and Table 3. Poisson shot-noise estimates of the standard deviation are shown as error bars on the histograms. Note that errors are smaller than the line width in the bottom panel.</figcaption></figure>

We perform the same test under more realistic conditions by fitting the $q_{2D}$ distribution obtained by the projection of idealized 3D ellipsoids with axis ratios measured from HAGN galaxies (test b). The fit, shown as red line in the central panel of Fig. 12, is in an good agreement with the corresponding measurements. We attribute the remaining deviations to inaccuracies of the $P(q_{3D}, r_{3D})$ model and not to the assumption of ellipsoidal stellar mass isodensities used for obtaining the 2D axis ratios, since this assumption is employed for generating both, the prediction as well as the reference measurements. The parameter constraints, shown as red contours in Fig. 13, are significantly biased with respect to the reference parameters from the $P(q_{3D}, r_{3D})$ model fits to the 3D axial ratio distribution in HAGN. However, the relative deviations of the best fits (defined as maxima of the marginalized probabilities) are at the per cent level for $q_0, r_0$, and $\sigma_r$ (see Table 3). The $\sigma_q$ parameter shows the strongest deviation from the reference values of roughly 100 per cent.

Our most realistic performance test of the reconstruction method consists in fitting the model to the distribution of 2D axis ratios measured from each galaxy's projected stellar mass distribution (test c). We find that the fit, shown as blue line in the top panel of Fig. 12, strongly deviates from the corresponding measurements. This finding can be expected from the previously discussed deviation

<footnote n="9">⁹emcee.readthedocs.io</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Posteriors of the parameters of the Gaussian model for the 3D axial ratio distribution, given by equation (5). Results are derived from the 2D axial ratio distributions of a) an ensemble of projected ellipsoids with 3D axis ratios that follow the same Gaussian model (model input values are marked as black dashed lines), b) projected 3D ellipsoids with 3D axial ratio distributions of disc galaxies in the HAGN simulation, c) projected stellar mass distributions of disc galaxies in HAGN. Light and dark areas mark 95 and 68 per cent confidence levels, respectively.</figcaption></figure>

between the $q_{2D}$ distribution obtained from the galaxies' projected stellar masses and the one obtained analytically from the 3D axial ratio measurements assuming idealized 3D ellipsoidal isodensities since the latter assumption is employed in the modelling. However, some main characteristics, such as the cut-off of the distribution, the position of the local maxima and the skewness towards low-axis ratios, are captured by the fit. We find that the parameter constraints, shown as blue contours in Fig. 13, are biased more strongly with respect to the reference values from the direct $P(q_{3D}, r_{3D})$ model fits compared to the previous test case b) (see also Table 3).

A potential shortcoming of these tests is that we defined discs in HAGN as the upper 20 per cent of galaxies with the highest values of $r_v = V_{rot}/\sigma_v$, which corresponds to a selection cut at $r_v > 1.06$ (see Section 2.2.1). We thereby might include discs with large bulges for which a one-component ellipsoidal model is inaccurate. In order to investigate how robust our accuracy tests are towards variations of the disc selection criteria, we repeat our analysis using galaxies with $r_v > 1.17$ and $r_v > 1.27$. This corresponds to selecting the upper 10 per cent and 5 per cent of galaxies with the highest $r_v$.

<table><caption>**Table 4.** Relative deviations between the parameters of the Gaussian $P(q_{3D}, r_{3D})$ model, reconstructed from the 2D axial ratio distribution of the projected stellar masses in HAGN (test case c) in Section 3.3) with respect to the reference parameters obtained from direct fits to the 3D axial ratio distribution. Results are shown for three samples, consisting of the 20%, 10%, and 5% of galaxies at $z = 1.0$ with the highest values of $V_{rot}/\sigma_v$, respectively. Each samples minimum value of $V_{rot}/\sigma_v$ is given in the left column.</caption><thead><tr><th>min(V<sub>rot</sub>/σ<sub>v</sub>)</th><th>Δq<sub>0</sub>[%]</th><th>Δσ<sub>q</sub>[%]</th><th>Δr<sub>0</sub>[%]</th><th>Δσ<sub>r</sub>[%]</th></tr></thead><tbody><tr><td>1.06</td><td>3.0</td><td>-64.1</td><td>-8.0</td><td>-14.6</td></tr><tr><td>1.17</td><td>3.1</td><td>-74.2</td><td>-8.1</td><td>-16.6</td></tr><tr><td>1.27</td><td>3.6</td><td>-51.2</td><td>-8.3</td><td>-24.3</td></tr></tbody></table>

values, respectively, and should reduce the fraction of discs with large bulges at the price of higher noise on the axial ratio distributions due to the smaller sample size. Our results from these different samples show no significant improvement in the accuracy of the reconstructed parameters when increasing the minimum values of $r_v$. This finding speaks against the possibility that discs with large bulges are the main reason for the inaccuracies of the reconstruction method. For test c), that is based on the projected 2D stellar mass distributions, we obtain overall relative deviations between the reconstructed model parameters and those derived from fits to the 3D stellar mass distributions of around ~3 per cent, ~60 per cent, ~8 per cent, and ~20 per cent for $q_0$, $\sigma_q$, $r_0$, and $\sigma_r$, respectively (Table 4).

Refocusing on our sample of discs in HAGN with $r_v > 1.06$, we compare in Fig. 14, the reconstructed distributions for $q_{3D}$ and $r_{3D}$ directly with the HAGN measurements. This figure illustrates that the location of the maxima of the reconstructed $q_{3D}$ and $r_{3D}$ distributions, given by the model parameters $q_0$ and $r_0$, are close to the positions of the maxima in the measurements. The width of the reconstructed $q_{2D}$ distribution, quantified by the $\sigma_q$ parameter, differs strongly from the HAGN measurements and is therefore not considered in the final discussion of our results.

We emphasize that the conclusion drawn in this section could change when conducting these tests under more realistic conditions. An ideal test would be based on synthetic images including effects of dust extinction, lensing, PSF convolution and pixelization, as well as image noise. These images would then need to be analysed using the same software as used for the observations. Such a realistic test is beyond the scope of this work, but would be an interesting objective for future investigations.

## 4 3D AXIAL RATIO DISTRIBUTIONS IN COSMOS

We apply the method for reconstructing 3D axial ratio distributions to the different samples from the COSMOS data, that were selected

<table><caption>**Table 3.** Best-fitting parameters of the model for the 3D axial ratio distribution, derived from disc galaxies in the HAGN simulation at $z = 1.0$. The two top rows show results for the Gaussian and the skew Gaussian model from direct fits to the marginalized distributions of $q_{3D}$ and $r_{3D}$ (see Section 3.1), shown in Fig. 11. The two bottom rows show the parameters of the Gaussian model inferred from the distribution of 2D axis ratios using the reconstruction method. The 2D axis ratios used for the reconstruction are obtained either analytically from the galaxies' 3D shapes, assuming that the stellar mass distributions are perfect 3D ellipsoids, or directly from shape measurements of the projected stellar mass distribution (labelled as projected ellipsoids and projected stellar mass respectively, see Fig. 12).</caption><thead><tr><th></th><th>model</th><th>q<sub>0</sub></th><th>σ<sub>q</sub></th><th>γ<sub>q</sub></th><th>r<sub>0</sub></th><th>σ<sub>r</sub></th><th>γ<sub>r</sub></th></tr></thead><tbody><tr><td rowspan="2">fits to 3D axial ratio distribution:</td><td>skew Gauss</td><td>0.977</td><td>0.081</td><td>-5.009</td><td>0.426</td><td>0.112</td><td>2.016</td></tr><tr><td>Gauss</td><td>0.921</td><td>0.049</td><td>—</td><td>0.504</td><td>0.073</td><td>—</td></tr><tr><td rowspan="2">reconstruction from 2D axial ratio distribution:</td><td>projected 3D ellipsoids</td><td>Gauss</td><td>0.934</td><td>0.093</td><td>—</td><td>0.515</td><td>0.071</td><td>—</td></tr><tr><td>projected stellar mass</td><td>Gauss</td><td>0.949</td><td>0.017</td><td>—</td><td>0.464</td><td>0.062</td><td>—</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Testing the reconstruction of the 3D axial ratio distribution of disc galaxies in the HAGN simulation at $z = 1.0$. Grey histograms show measurements in the simulation with $1\sigma$ shot-noise errors displayed as black bars for each bin. Black solid lines show the fits to a truncated normal distribution, which are also shown as solid lines in the top panel of Fig. 11. Red dashed lines show the distribution reconstructed from fits to the 2D axial ratio distribution which was derived from the 3D axis ratios using the ellipsoidal galaxy model (see central panel of Fig. 12). Dash-dotted blue lines show the reconstruction based on fits to the 2D axis ratios measured from the projected stellar mass (see top panel of Fig. 12). Thin lines show predictions for 500 random sampling points of the posterior distributions and reflect the uncertainties expected on the prediction.</figcaption></figure>

by stellar mass and redshift as described in Section 2.1.6. The fits to the observed distribution of 2D axis ratios from which we infer the parameters of our model for the 3D axial ratio distribution are shown in Fig. 15. For the samples that are selected by both, stellar mass and redshift, we find that the fits match the measurements well with $\chi^2$ values per degree of freedom around unity (see Table 5), capturing the skewness towards low and high values of $q_{2D}$ for the high and low mass samples, respectively (as shown in Fig. 9).

The fits to the samples that are selected only by stellar mass or only by redshift, shown in the top row and left-hand column of Fig 15, respectively, are less accurate as indicated by the higher $\chi^2$ values per degree of freedom. This finding can be expected on

one hand from the larger sample size that decreases the errors on the measurements and hence increases the significance of deviations between measurements and model with respect to the errors. On the other hand the underlying model assumption that the 3D axis ratios, $q_{3D}$ and $r_{3D}$, follow a Gaussian distribution may be less accurate for samples that are more broadly defined. This could in particular be the case if sub-samples have very distinct $q_{2D}$ distributions, like the low and high mass sub-samples.

However, we note that the model fits the observed $q_{2D}$ distributions overall better than corresponding distributions from the projected stellar mass in the HAGN simulation, shown in the top panel of Fig. 12.

This finding might result from the fact that we selected all disc galaxies in HAGN, including those with large bulges for which a one-component ellipsoidal shape model is clearly an inadequate approximation. Additional reasons may be differences between the 2D shapes of the projected stellar mass and the projected luminosity densities or shortcomings of the simulation, caused for instance by resolution effects as discussed in Section 2.2.

In Fig. 16, we show the MCMC estimates for the posterior probability distribution of the $P(q_{3D}, r_{3D})$ model parameters, $q_0, r_0, \sigma_q$, and $\sigma_r$, quantifying the average disc circularity, the average relative disc thickness and the corresponding dispersions, respectively. The best fit parameters, defined as the positions of the maxima of the posteriors are summarized in Table 5. In this table, we also provide results for the parameter $s_0 \equiv q_0 r_0$, whose posterior distribution we obtain by multiplying the coordinates of the sampling points of the posterior distributions of $q_0$ and $r_0$. It specifies the position of the maximum in the $P(q_{3D}, s_{3D})$ distribution and quantifies the disc thickness relative to the major axis, facilitating a comparison to the disc circularity parameter $q_0$ and the estimation of the absolute disc thickness (see Section 4.2.2). The posteriors shown in Fig. 16 and the corresponding model parameters in Table 5 were derived from $q_{2D}$ distributions measured in 25 bins. When repeating the same analysis using 20 bins, we find only marginal changes of the posteriors and no significant change in the corresponding model parameters. However, larger fluctuations may occur for larger variations of the binning. In Appendix F, we show that the posteriors can be strongly affected by noise in the data. Varying the binning would change the noise and therefore the posteriors. However, a rigorous investigation of how the binning impacts the posteriors is beyond the scope of our analysis.

## 4.1 Redshift dependence

We start examining the redshift dependence of the galaxy shape distribution by comparing the fits to the different redshift sub-samples with those for the entire redshift range (shown as red solid and blue dashed lines in Fig. 15, respectively). We find that within a given mass range, these different fits are close to each other, which is a first indication that the redshift evolution of the axial ratio distribution in our samples is weak.

This visual impression lines up with the relatively small variation of the best fit parameters for the average disc circularity and the relative disc thickness ($q_0$ and $r_0$) with redshift. For the full mass redshift sub-samples, we find the values of $q_0$ and $r_0$ to vary by ~1 per cent and ~10 per cent around the results for the full redshift range, respectively (see Table 5). For the higher and low-mass sub-samples these variation increase up to ~10 per cent and ~25 per cent for $q_0$ and $r_0$, respectively. The redshift variations of the corresponding dispersions, $\sigma_q$ and $\sigma_r$ are much higher with up to ~50 per cent. However, the amplitudes of these different redshift variations are
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Probability distribution of apparent axis ratios $q_{2D}$ of disc-dominated galaxies from our COSMOS samples in different redshift and stellar mass ranges, shown as black histograms. Error bars show shot noise estimates of the standard deviation in each bin. Red lines are model fits for each redshift sub-sample in a given stellar mass range. Blue lines are fits for the stellar mass samples defined over the full redshift range (shown in the top panels). The model parameters derived from these fits are summarized in Table 5. The corresponding posterior distributions are shown in Fig. 16.</figcaption></figure>

overall consistent with the ~68 per cent confidence intervals from the marginalized posteriors.

In order to obtain a more detailed insight into how the parameter variations with redshift compare to the uncertainties, we show the posterior probability distributions of the $P(q_{3D}, r_{3D})$ model

parameters from the different sub-samples in Fig. 16. The 68 per cent confidence levels of the posteriors for the full mass samples in the different redshift bins, shown in the top left of Fig. 16, are mostly overlapping. A weak indication for a redshift dependence in the data might be given by the fact that the 68 per cent confidence levels of the
<table><caption>**Table 5.** Parameters of the 3D axial ratio distribution model (equations (5) and (6)), inferred from fits to the 2D axial ratio distribution from the COSMOS survey, shown in Fig. 15. Results are given for different stellar mass – redshift samples, which are selected as indicated in the two left columns. For each model parameter we provide the lower and upper limit of the 68% confidence level, obtained from the marginalized posteriors. The right column shows the $\chi^2$ goodness of fit per degree of freedom.</caption><thead><tr><th>stellar mass range</th><th>redshift range</th><th>N<sub>gal</sub></th><th>q<sub>0</sub></th><th>σ<sub>q</sub></th><th>r<sub>0</sub></th><th>σ<sub>r</sub></th><th>s<sub>0</sub></th><th>χ<sup>2</sup>/d.o.f.</th></tr></thead><tbody><tr><td rowspan="4">full mass range</td><td>0.2 &lt; z &lt; 1.0</td><td>3749</td><td>0.88<sup>0.92</sup><sub>0.87</sub></td><td>0.08<sup>0.12</sup><sub>0.07</sub></td><td>0.29<sup>0.30</sup><sub>0.27</sub></td><td>0.09<sup>0.12</sup><sub>0.09</sub></td><td>0.25<sup>0.27</sup><sub>0.24</sub></td><td>2.20</td></tr><tr><td>0.2 &lt; z &lt; 0.7</td><td>971</td><td>0.89<sup>0.91</sup><sub>0.87</sub></td><td>0.04<sup>0.10</sup><sub>0.03</sub></td><td>0.29<sup>0.31</sup><sub>0.28</sub></td><td>0.10<sup>0.13</sup><sub>0.09</sub></td><td>0.26<sup>0.28</sup><sub>0.25</sub></td><td>1.39</td></tr><tr><td>0.7 &lt; z &lt; 0.9</td><td>1551</td><td>0.89<sup>0.92</sup><sub>0.87</sub></td><td>0.09<sup>0.12</sup><sub>0.07</sub></td><td>0.27<sup>0.28</sup><sub>0.26</sub></td><td>0.10<sup>0.12</sup><sub>0.09</sub></td><td>0.24<sup>0.25</sup><sub>0.23</sub></td><td>1.18</td></tr><tr><td>0.9 &lt; z &lt; 1.0</td><td>1217</td><td>0.87<sup>0.93</sup><sub>0.86</sub></td><td>0.10<sup>0.14</sup><sub>0.07</sub></td><td>0.32<sup>0.33</sup><sub>0.30</sub></td><td>0.12<sup>0.14</sup><sub>0.11</sub></td><td>0.28<sup>0.29</sup><sub>0.26</sub></td><td>1.01</td></tr><tr><td rowspan="4">M<sub>*</sub> &lt; 10<sup>10.35</sup> M<sub>☉</sub></td><td>0.2 &lt; z &lt; 1.0</td><td>1913</td><td>0.87<sup>0.90</sup><sub>0.86</sub></td><td>0.07<sup>0.10</sup><sub>0.06</sub></td><td>0.37<sup>0.39</sup><sub>0.36</sub></td><td>0.21<sup>0.23</sup><sub>0.20</sub></td><td>0.33<sup>0.34</sup><sub>0.31</sub></td><td>1.28</td></tr><tr><td>0.2 &lt; z &lt; 0.7</td><td>443</td><td>0.89<sup>0.91</sup><sub>0.82</sub></td><td>0.05<sup>0.22</sup><sub>0.04</sub></td><td>0.43<sup>0.63</sup><sub>0.41</sub></td><td>0.21<sup>0.29</sup><sub>0.18</sub></td><td>0.38<sup>0.53</sup><sub>0.36</sub></td><td>0.58</td></tr><tr><td>0.7 &lt; z &lt; 0.9</td><td>789</td><td>0.86<sup>0.88</sup><sub>0.84</sub></td><td>0.07<sup>0.11</sup><sub>0.06</sub></td><td>0.28<sup>0.31</sup><sub>0.11</sub></td><td>0.35<sup>0.52</sup><sub>0.30</sub></td><td>0.23<sup>0.27</sup><sub>0.10</sub></td><td>0.98</td></tr><tr><td>0.9 &lt; z &lt; 1.0</td><td>681</td><td>0.83<sup>0.86</sup><sub>0.80</sub></td><td>0.06<sup>0.21</sup><sub>0.05</sub></td><td>0.35<sup>0.46</sup><sub>0.33</sub></td><td>0.20<sup>0.29</sup><sub>0.18</sub></td><td>0.30<sup>0.40</sup><sub>0.28</sub></td><td>1.13</td></tr><tr><td rowspan="4">M<sub>*</sub> &gt; 10<sup>10.35</sup> M<sub>☉</sub></td><td>0.2 &lt; z &lt; 1.0</td><td>1836</td><td>0.90<sup>0.97</sup><sub>0.89</sub></td><td>0.15<sup>0.18</sup><sub>0.12</sub></td><td>0.26<sup>0.27</sup><sub>0.25</sub></td><td>0.06<sup>0.06</sup><sub>0.05</sub></td><td>0.24<sup>0.25</sup><sub>0.23</sub></td><td>1.39</td></tr><tr><td>0.2 &lt; z &lt; 0.7</td><td>532</td><td>0.88<sup>0.94</sup><sub>0.86</sub></td><td>0.07<sup>0.15</sup><sub>0.06</sub></td><td>0.24<sup>0.26</sup><sub>0.24</sub></td><td>0.04<sup>0.06</sup><sub>0.03</sub></td><td>0.21<sup>0.23</sup><sub>0.21</sub></td><td>0.79</td></tr><tr><td>0.7 &lt; z &lt; 0.9</td><td>768</td><td>0.98<sup>0.98</sup><sub>0.88</sub></td><td>0.25<sup>0.29</sup><sub>0.20</sub></td><td>0.27<sup>0.29</sup><sub>0.26</sub></td><td>0.06<sup>0.07</sup><sub>0.04</sub></td><td>0.26<sup>0.28</sup><sub>0.24</sub></td><td>1.18</td></tr><tr><td>0.9 &lt; z &lt; 1.0</td><td>536</td><td>0.92<sup>0.97</sup><sub>0.88</sub></td><td>0.12<sup>0.15</sup><sub>0.08</sub></td><td>0.28<sup>0.30</sup><sub>0.27</sub></td><td>0.09<sup>0.10</sup><sub>0.07</sub></td><td>0.26<sup>0.28</sup><sub>0.25</sub></td><td>0.87</td></tr></tbody></table>

intermediate and the high redshift samples are slightly disconnected in the $r_0$ direction. Nonetheless, the 95 per cent confidence levels still overlap clearly in that case. More significant differences between the posteriors from different redshift bins are present for the low and high mass sub-samples, as shown in the bottom left and right of Fig. 16, respectively. However, for the low-mass sample the 68 per cent confidence levels still overlap plainly. For the high mass sample the 68 per cent confidence level of the central redshift bin is slightly disjoint from those of the other bins, although the 95 per cent confidence levels from all redshifts still overlap. Overall our Bayesian analysis does not provide convincing evidence for a redshift dependence of the $P(q_{3D}, r_{3D})$ model parameters that is significant with respect to the uncertainties. It is thereby important to note that the differences in the shapes of the posteriors are not an indication for such a dependence since they can be expected from the sampling noise in the $q_{2D}$ distributions, as we demonstrate in Appendix F.

It is further important to note that the insignificant redshift dependence of the parameters with respect to the uncertainties does not imply that the shapes of disc-dominated galaxies do not evolve since $z = 1.0$. A potential evolution could just be too weak to be detected reliably with the Bayesian analysis presented here. In particular, the large uncertainties on the dispersions of the circularity and relative thickness may obscure even a relatively strong redshift dependence of these parameters. We therefore employ the two-sample Kolmogorov-Smirnov test (e.g. Hodges 1958, hereafter referred to as K-S test) as an independent tool for testing the null hypothesis that the axial ratio distributions from two different redshift samples are drawn randomly from the same redshift independent distribution. We apply this test directly on the observed $q_{2D}$ distributions and are therefore independent of any model assumptions. In Table 6, we show the p-values of the K-S test, which is the probability of obtaining a K-S test result that is as large or larger than the one obtained for the two samples under the null hypothesis. A typical, but arbitrary rejection criteria is $p < 0.05$. For the redshift samples defined over the full mass range the K-S test does not reject the null hypothesis. This result lines up with the strong overlap of the model parameter contours for these samples, which we discussed above. The situation is less conclusive for the p-values from the low- and high-mass sub-samples. For the low-mass sub-sample the K-S test suggests, that the $q_{2D}$ distributions from the lowest and highest redshift bins are drawn

from different underlying distributions, since the p-value of 0.0273 is below 0.05. This finding is interesting, given that the parameter contours from these two redshift bins are strongly overlapping (see bottom left-hand panel of Fig. 16). Given this apparent disagreement and that fact that i) our arbitrarily chosen rejection criteria just slightly surpassed for this redshift bin combination and ii) the null hypothesis has not been rejected for the other redshift bins combinations, we conclude that the low-mass sub-sample provides no strong indication for a substantial redshift evolution of the galaxy shapes.

For the high mass sub-sample the K-S test suggests, that the $q_{2D}$ values from the intermediate and the highest redshift bin are drawn from different underlying distributions. In this case, the rejection criteria is clearly fulfilled by the p-value of 0.0053. This finding lines up with the weaker overlap of the corresponding parameter contours, shown in the bottom right of Fig. 16, and could be an indication for a redshift dependence of the shapes of high mass disc-dominated galaxies. However, it could also result from systematic effects, such as potential misclassifications of galaxy morphologies in the ZEST catalogue. Another explanation for a redshift dependence might be cosmic variance, as galaxy clusters in the COSMOS field (see Fig. 2) could change the mean values and dispersions of the disc thickness and circularity in a given redshift bin. It is thereby interesting to note that the $q_{2D}$ distributions from the lowest and the highest redshift sample are consistent with each other, according to both, the K-S test as well the overlap of the parameter contours for these two samples shown in Fig. 16, which speaks against a physical evolution of the shapes. However, the same systematic effects may also change the $q_{2D}$ distributions in such a way that a detection of a redshift evolution is prevented. Samples with more objects that are distributed over larger areas and morphological classifications based on different techniques would be needed to better understand the impact of misclassification or cosmic variance on our results.

## 4.2 Mass dependence

We now proceed with studying the dependence of the galaxy shapes on stellar mass. Given that the evidence for a redshift evolution of the shapes is overall weak, we will thereby focus on the mass samples that are defined over the entire redshift range ($0.2 < z < 1.0$) to maximize the statistical power of the samples.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Posteriors of the 3D axial ratio model parameters from equation (5), derived from the 2D axial ratio distributions in COSMOS (Fig. 15). Light and dark areas mark 95 and 68 per cent confidence levels, respectively. Results for i) redshift samples defined over the full stellar mass range (top left), ii) stellar mass samples, defined over the full redshift range (top right), iii, iv) redshift samples with stellar masses below and above $10^{10.35} M_{\odot}$ (bottom left and right, respectively). The marginalized posterior distributions are shown on the diagonal panels of each sub-figure.</figcaption></figure>

<table><caption>**Table 6.** *p*-values of the two sample Kolmogorov–Smirnov test. The test is applied on pairs of 2D axial ratio distributions measured in three redshift bins (Fig. 15). The redshift bins $z_0$, $z_1$, $z_2$ correspond to the ranges [0.2, 0.7], [0.7, 0.9], [0.9, 1.0], respectively.</caption><thead><tr><th></th><th>full mass range</th><th><i>M</i><sub>*</sub> &lt; 10<sup>10.35</sup>1</th><th><i>M</i><sub>*</sub> &gt; 10<sup>10.35</sup></th></tr></thead><tbody><tr><td><i>z</i>0–<i>z</i>1</td><td>0.3796</td><td>0.2177</td><td>0.2496</td></tr><tr><td><i>z</i>1–<i>z</i>2</td><td>0.1684</td><td>0.1648</td><td>0.0053</td></tr><tr><td><i>z</i>2–<i>z</i>0</td><td>0.8569</td><td>0.0273</td><td>0.1221</td></tr></tbody></table>

### 4.2.1 Intrinsic disc circularity

The intrinsic circularity of the discs in our samples is quantified by the axial ratio $q_{3D}$. The maximum and the width of the $q_{3D}$ distribution are described by the parameters $q_0$ and $\sigma_q$ of our $P(q_{3D}, r_{3D})$ model.

In the top right of Fig. 16, we compare the constraints on these model parameters from the low and high-mass samples in the entire redshift range. We find that the high-mass sample prefers slightly higher values of $q_0$ and $\sigma_q$ compared to the low-mass samples. However, the corresponding confidence intervals from both mass samples are overlapping significantly, which indicates that our data provides no evidence for a mass dependence of the disc circularity.

We find $q_0$ values of around 0.9 for both samples (see Table 5), which is consistent with estimates reported for discs at lower redshifts in previous studies (Sandage et al. 1970; Fasano et al. 1993; Ryden 2004, 2006; Rodríguez & Padilla 2013) and describes the deficit of circular face-on discs with $q_{2D} = 1.0$ in the observed data (see Fig. 15). Different reasons for this deficit have been discussed in the literature (e.g. Bertola, Vietri & Zeilinger 1991; Huizinga & van Albada 1992; Rix & Zaritsky 1995; Bernstein & Jarvis 2002; Joachimi et al. 2013). On the one hand, it is argued that deviations
from perfect circularity could result from observational effects, such as noisy isophotes, artifacts like cosmic rays in the galaxy images, or simply the fact that the images are pixelized. However, we expect the impact of such effects on our measurements to be minor, since i) we selected objects with ‘good’ fits to the Sérsic profile (see Section 2.1.2), which should remove objects whose images are heavily distorted by artefacts from the sample and ii) the ACS pixel scale of 0.03" is well below our lowest limit for the effective angular radius of 0.17" at z = 1.0 (see Fig. 2). However, since the deviations from perfect circularity are predicted to be relatively small, even minor systematics might be relevant. On the other hand, one could expect that the disc galaxies are intrinsically not perfectly circular due to patchy star formation activity and sub-structures, such as spiral arms (see Fig. 1) or Galactic warps (see e.g. Binney 1992; Gómez et al. 2017, for the latter). This expectation lines up with the non-circularity of discs in the HAGN and TNG100 simulations, which we see in Fig. 11 and 12 as the lack of galaxies with $q_{3D}$ and $q_{2D}$ close to unity.

Note that we omit entering a detailed discussion of our $\sigma_q$ constraints, since the test of our reconstruction method in Section 3.3 and Appendix F suggested that these constraints may be strongly biased.

### 4.2.2 Intrinsic disc thickness

The relative disc thickness is quantified by the minor to intermediate axial ratio $r_{3D}$. The position of the maximum and width of the $r_{3D}$ distribution are described in our $P(q_{3D}, r_{3D})$ model by the parameters $r_0$ and $\sigma_r$. In the top right of Fig. 16, we show that $r_0$ exhibits a significant dependence on stellar mass. It is predicted to be around 0.3 with lower and higher values for the high and low-mass sample, respectively (see Table 5). The corresponding constraints on the relative thickness $s_0 \equiv q_0 r_0$ (quantifying the position of the maximum of the $s_{3D} \equiv C_{3D}/A_{3D}$ distribution) are 0.33$^{0.34}_{0.31}$ and 0.24$^{0.25}_{0.23}$ for the low-mass and high-mass sample, respectively, where the upper and lower values here are the limits of the 68 per cent confidence intervals. These values are consistent with results reported for discs at low redshifts which were obtained using similar reconstruction techniques as used in this work (Fasano et al. 1993; Ryden 2004, 2006; Rodríguez & Padilla 2013) or from direct measurements based on edge-on oriented discs (Mosenkov et al. 2015; Reshetnikov et al. 2016; Mosenkov et al. 2020). The dispersion of the relative disc thickness, $\sigma_r$, shows a significant mass dependence as well, as it takes lower and higher values for the high and low-mass sample, respectively. This mass dependence of the thickness dispersion explains why the cuff-off on the left side of the $q_{2D}$ distribution in Fig. 15 is sharper for high than for low-mass discs, while we expect also this latter parameter to be significantly biased by ~20 per cent (Section 3.3).

In Fig. 17, we show the reconstructed marginalized distributions of the 3D axis ratios $q_{3D}$ and $r_{3D}$ for the low- and high-mass sample over the full redshift range, as predicted from the best-fitting parameters given in Table 5. The uncertainties on the prediction is illustrated as the density of thin lines which are obtained from randomly selected sampling points of the posterior in the $P(q_{3D}, r_{3D})$ model parameters space. The corresponding values $q_0$ and $r_0$ are shown as vertical lines. They are enclosed by shaded areas which indicate the 68 per cent confidence intervals for these parameters, given in Table 5. Note that these uncertainties are similar to the systematic bias of ~3 per cent and ~8 per cent for $q_0$ and $r_0$, respectively which we found in our systematic tests (Section 3.3).

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Distributions of galaxy 3D axis ratios, reconstructed from the distribution of 2D galaxy axial ratio in the COSMOS survey. Results are shown for the low and high mass sample, defined within 0.2 < z < 1.0. Vertical lines mark the positions of the maxima, described by the best-fitting model parameters $q_0$ and $r_0$ from Table 5. Vertical shaded areas indicate the corresponding 68 per cent confidence intervals. Thin lines show predictions for 500 random sampling points of the posterior distribution.</figcaption></figure>

When comparing the distributions of both mass samples, we see that the values of $r_{3D}$ are well below those of $q_{3D}$ as expected for disc galaxies. This indicates that our results carry physically meaningful information, despite the expected biases. Fig. 17 further illustrates that i) high-mass discs tend to be thinner with respect to their diameter than low-mass discs, ii) high-mass discs tend to be slightly more circular than low-mass discs, iii) the dispersion of the circularity is larger for high mass than for low-mass discs and iv) the dispersion of the relative thickness is larger for low mass than for high-mass discs. However, the mass dependence of the peak of the circularity distribution, quantified by $q_0$, is statistically not significant, since the 68 per cent confidence intervals for $q_0$ overlap (see also discussion in Section 4.2.1). The lower relative thickness for high-mass discs could indicate that these galaxies tend to be more relaxed than low-mass discs, which are more prone to perturbations by feedback, merging and tidal interactions. This interpretation is supported by the higher dispersion of the relative thickness of low-mass discs and lines up with the overall higher star formation rates for low-mass discs, which we see in the right-hand panels of Fig. 9.
However, another reason to expect a lower relative thickness of massive discs with respect to their radius is that the radius is on average larger for discs in the high-mass sample, as indicated by the mass dependence of the effective radii shown in the central panels of Fig. 9. Such an increase of the radius with mass could lead to a decrease of the relative thickness, even if the size of the absolute thickness is mass independent.

In order to test, if this scenario is consistent with our data we derive an estimate of the expectation value for the absolute thickness, $\langle C_{3D} \rangle$, for the low- and high-mass sample via our constraints on the peaks of the $s_{3D} \equiv C_{3D}/A_{3D}$ axial ratio distributions, quantified by $s_0$. Note that we use $s_{3D}$ instead of $r_{3D} \equiv C_{3D}/B_{3D}$, since we can approximate $A_{3D}$ directly from the observed 2D radii. We estimate the expectation value for $C_{3D}$ using the approximation $\langle C_{3D} \rangle \simeq \langle A_{3D} \rangle \langle C_{3D}/A_{3D} \rangle = \langle A_{3D} \rangle \langle s_{3D} \rangle \simeq \langle A_{3D} \rangle s_0$. This approximation is based on two assumptions: 1) We assume $\langle s_{3D} \rangle \equiv \langle C_{3D}/A_{3D} \rangle \simeq \langle C_{3D} \rangle / \langle A_{3D} \rangle$. This approximation corresponds to the zero-order term in a perturbative expansion of $f(A, C) = C/A$ around the $\langle A \rangle$ and $\langle C \rangle$. We validated this relation using our sample of disc galaxies in the HAGN simulation and found it to be accurate at the sub-per cent level. 2) We assume that the expectation of the $s_{3D}$ distribution is located at its maximum, i.e. $\langle s_{3D} \rangle \simeq s_0$. This implies that the $s_{3D}$ distribution is symmetric, which is commonly assumed to be the case for observed discs (e.g. Ryden 2004). Again, we also verified that this is a good approximation using the disc galaxies in our HAGN sample.

We obtain the posterior distribution of $s_0$ from the posteriors of $q_0$ and $r_0$, using $s_0 = q_0 r_0$ as detailed in the beginning of this section. It is shown in the top panel of Fig. 18 for the low-mass and high-mass disc samples, defined over the entire redshift range. The distributions differ significantly from each other, as it can be expected from the constraints on $r_0$, shown in Fig. 16. We proceed by estimating the posterior for the expectation value of the minor axis size from the posterior of $s_0$, using the aforementioned relation $\langle C_{3D} \rangle \simeq \langle A_{3D} \rangle s_0$, which we apply on each MCMC sampling point. For this last step, we need an estimate for the average comoving size of the major axis $A_{3D}$. This size is approximated with the mean effective comoving radius $R_\perp$ of the disc galaxies in our sample, using objects with $q_{2D} > 0.5$. This value is shown in the central panels of Fig. 9 as vertical red and blue solid lines for the high and low-mass sample, respectively. The cut on $q_{2D}$ was thereby chosen to mitigate the impact the projection effects on the effective radius as discussed in Section 2.1.5.

The estimated posterior of $\langle C_{3D} \rangle$ is shown for the samples of high-mass and low-mass discs in the bottom panel of Fig. 18. The two posteriors overlap each other almost completely, which means that our approximation of the average 3D disc thickness shows no significant mass dependence. These posterior distributions are compared to the estimates of the mean major axis sizes $\langle A_{3D} \rangle$ which are shown as vertical lines in the same panel. The strong mass dependence, which we find for $\langle A_{3D} \rangle$ suggests that the mass dependence in the relative 3D thickness $s_{3D}$ is driven by the mass dependence of the major axis sizes.

# 5 SUMMARY AND CONCLUSIONS

We studied the 3D shapes of disc-dominated late-type galaxies from the COSMOS survey in different mass and redshift ranges, tackling the question of how these galaxies could grow without forming a large central bulge. We approximated the 3D light distribution of these galaxies as 3D ellipsoids described by the two ratios $q_{3D} \equiv B_{3D}/A_{3D}$ and $r_{3D} \equiv C_{3D}/B_{3D}$, which quantify the circularity and relative thickness of the discs, respectively. We inferred the distribution of these 3D axis ratios from the observed distribution

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Top: Marginalized posterior of the parameter $s_0 \equiv q_0 r_0$ that serves as an approximation for the average 3D axial ratio $\langle s_{3D} \rangle = \langle C_{3D}/A_{3D} \rangle$. Bottom: Dash-dotted vertical lines mark the average major axes length $\langle A_{3D} \rangle$ that is approximated as the average comoving effective radius $R_\perp$ of discs with $q_{2D} > 0.5$. Dashed lines show the marginalized posterior of the average minor axes length $\langle C_{3D} \rangle$ which is approximated as $\langle C_{3D} \rangle \simeq s_0 \times \langle A_{3D} \rangle$. The different approximations are discussed in Section 4.2.2. Results for the low and high mass sample, defined within $0.2 < z < 1.0$, are shown as blue and red lines, respectively.</figcaption></figure>

of 2D axis ratios, using a reconstruction method based on the assumption that the distribution of $q_{3D}$ and $r_{3D}$ is well approximated by a two-dimensional Gaussian. This Gaussian is characterized by the average disc circularity $q_0$, the average disc thickness $r_0$ and the corresponding dispersions $\sigma_q$ and $\sigma_r$.

Variations of this method have been widely used in the literature, but their accuracy as well as the assumptions they employ remained to be tested. We began our analysis by performing such tests for the first time using two state-of-the-art hydrodynamic simulations of galaxy formation in cosmological volumes, Horizon-AGN and Illustris TNG100. We demonstrated that the 3D axial ratio distributions of disc galaxies in these simulations are adequately described by a Gaussian model. Reconstructing the 3D axial ratio distribution from the distribution of 2D axis ratios, we found that the inferred parameters of the Gaussian model are biased with respect to those derived directly from fits to the 3D distributions. For our most realistic test, based on 2D galaxy shapes from projected stellar mass distributions, we find this bias to be moderate for the parameters $q_0$ and $r_0$ (~3 per cent and ~8 per cent, respectively) but strong for the corresponding dispersions $\sigma_q$ and $\sigma_r$ (~60 per cent and ~20 per cent, respectively, see Table 4). We concluded that the bias is mainly driven by an inaccuracy of the ellipsoidal model for the stellar mass distribution, which we use for relating 3D to 2D galaxy shapes analytically during the reconstruction (see Section 3.3). The strong simplification implied when approximating late-type galaxies as 3D ellipsoids becomes obvious in the COSMOS images of such objects, shown in Fig. 1. Nevertheless, our bias estimates derived
from the HAGN simulation may be overly pessimistic, since our test
is based on all disc galaxies found in this simulation, including those
which have a bulge for which a one-component ellipsoidal galaxy
model is expected to be inadequate. More realistic tests should be
based on synthetic galaxy images, take into account the effect of
dust extinction and be analysed in the same way as the reference
observational data.

After having tested the shape reconstruction method, we applied it
on a volume limited sample of disc-dominated galaxies in COSMOS
which is limited in redshift (0.2 < z < 1.0), absolute magnitude
(Mᵢ < -21.5), and transverse comoving size (R⊥ > 0.64 kpc). We
demonstrated that a conservative cut on Mᵢ is required in order to
minimize the impact of dust extinction on the observed 2D axial
ratio distribution. Otherwise, this effect can lead to an apparent
redshift evolution of the distribution (Fig. D1) and consequently to
incorrect physical interpretations of the observations. We found that
the ellipsoidal galaxy shape model in conjunction with the Gaussian
model for the 3D axial ratio distribution provides good fits to the
distribution of the observed 2D axis ratios (Fig. 15).

The constraints on the parameters *q*₀ and *r*₀ inferred from these fits
show a variation with redshift around ~1 per cent and ~10 per cent,
respectively when considering the full mass range (Table 5). This
finding indicates that the redshift evolution of the shapes is relatively
weak. Splitting the sample into a high and a low mass sub-sample
leads to larger variations of ~10 per cent and ~25 per cent for *q*₀ and
*r*₀, respectively. The parameters *σ*_q and *σ*_q vary up to ~50 per cent
across the different redshift samples. However, overall the variations
with redshift lie within or close to the ~68 per cent confidence
intervals on these parameters. Studying the joint posteriors on the
different shape model parameters, we find no clear indication for a
redshift evolution that is significant with respect to the uncertainties
as the 95 per cent confidence intervals of these parameters from
different redshift samples overlap (Fig. 16). We crosschecked this
result using the Kolmogorov-Smirnov (K-S) statistic to test the null
hypothesis that the 2D axial ratio distributions from two redshift
bins are drawn from the same underlying redshift independent
distribution. Studying all redshift bin combinations for different
mass samples, we find that the K-S test is mostly consistent with
the results from the Bayesian analysis. As an exception, the K-S test
clearly rejects the null hypothesis that the 2D axis ratios of the high-
mass samples in our intermediate and high-redshift bins are drawn
from the same distributions. This finding appears to be inconsistent
with the agreements between the 2D axial ratio distributions from the
high-mass samples in the low and high-redshift bins and could be an
indication for cosmic variance or misclassified objects contaminating
the high-mass sample at intermediate redshifts. In any case, the fact
that the Bayesian analysis revealed no significant redshift dependence
indicates that any dependence detected by the K-S test is relatively
weak compared to the errors on the data.

This weak dependence on redshift is contrasted by a relatively
strong dependence on mass (Fig. 16). In particular, the parameter
$r_0$ is significantly higher for discs in our sample with stellar masses
below $10^{10.35} M_{\odot}$ compared to discs with higher masses, even within
the expected uncertainties of the reconstruction method (e.g. Fig. 17).
This finding indicates that the relative disc thickness decreases with
mass. However, the absolute disc thickness, estimated from the
relative thickness and a proxy for the disc diameter, shows no mass
dependence (Fig. 18), which suggests that the relative thickness of
low-mass discs is higher mainly because of their smaller diameters
(e.g. Fig. 9).

In summary, we found indications that the distribution of the 3D
axis ratios of disc-dominated galaxies in our sample is not (or if,

then just weakly) redshift dependent in the range 0.2 < z < 1.0 and
the absolute disc thickness does not depend significantly on mass.
These findings speak against mass accretion by major mergers and
a subsequent suppression of bulge formation by strong feedback
after z ≲ 1.0, since one could expect that such disruptive events
would decrease the discs circularity and increase their thickness due
to vertical heating (e.g. Quinn et al. 1993; Grand et al. 2016; Park
et al. 2021). The decrease of the specific star formation rates of
massive discs with redshift, shown in Fig. 9, suggests that feedback
occurs in these galaxies, suppressing star formation by removing
cold gas, but may not be sufficiently strong to significantly affect the
observed shapes. The absence of major mergers in disc-dominated
galaxies since z ≲ 1 is further supported by the fact that the mean
comoving sizes, displayed in Fig. 9, show no redshift evolution for
low-mass galaxies and a very weak increase with redshift for high-
mass galaxies. We conclude that disc-dominated galaxies accreted
most of their mass before z = 1.0 and lived preferentially in isolation
ever since. This picture lines up with the results from Grossi et al.
(2018), who find that the star formation rates of disc-dominated
galaxies show no significant dependence on the density of their
environment, which indicates the absence of major interactions (see
also Sachdeva & Saha 2016). A tranquil evolution is further supported
by the weak redshift dependence of the transverse comoving size
distribution (Fig. 17), which has also been reported for disc galaxies
in the GOODS fields by Ravindranath et al. (2004).

However, we emphasize that the validity of our results is limited
for several reasons. The small sizes of our samples lead to large
statistical errors, which may obscure a weak evolution of the galaxy
shapes with redshift and could hence affect our conclusions. This
limitation can be overcome in future analysis by using measurements
of galaxy shapes in larger volumes from upcoming weak lensing
surveys like *Euclid* or the Legacy Survey of Space and Time on
the Vera Rubin Observatory. Such high precision measurements will
require an improved accuracy of the shape reconstruction method,
based for instance on corrections which can be developed using
realistic mock images.

Our results may further be affected by inaccuracies of the
assumptions on which the reconstruction method is based. One
inaccurate assumption is the aforementioned ellipsoidal model for
galaxy morphologies, which differs strongly from the complex
shapes of real disc galaxies. Another potential inaccuracy may result
from the Gaussian model, that we assume for the 3D axial ratio
distribution. Our choice for this model was based on 3D axial ratio
distributions measured in two hydro-dynamic simulations. However,
we also demonstrated that the 2D axial ratio distributions in these
simulations differ significantly from the COSMOS observations due
to the lack of thin discs with $q_{2D} \ll 1$, which we attribute to resolution
limits and approximations in the simulation method (see Section 2.2).
More realistic simulations could reveal that our Gaussian model
is inadequate, which could bias the outcome of the reconstruction
method.

Potential errors in our analysis may further result from a miss-
classification of galaxies in the ZEST catalogue. Although this
catalogue has been carefully calibrated and passed rigorous vali-
dations, it may contain a fraction of objects that are miss-classified
as disc-dominated late-type galaxies and contaminate our samples.
By inspecting the colour-colour distribution of the objects in our
catalogue, we find no indication of a significant fraction of miss-
classified objects. However, visual inspection of randomly selected
galaxies revealed the presence of some objects, in particular at high
redshifts and low masses, whose morphological class is ambiguous
to us (Appendix A). In order to assess the uncertainty from potential
miss-classifications on our results, one could compare the mor-
phological classification in ZEST against alternative classification
methods that are for instance based on machine learning techniques.

More detailed interpretations of the redshift evolution of observed galaxy shapes may also require a characterization of two observational systematic effects which we assumed to be insignificant with respect to the errors on our measurements. A first potential systematic could result from the fact that galaxy shapes are observed in the same filter at different redshifts. High-redshift galaxies are therefore seen at bluer rest frame wavelengths than galaxies at lower redshifts. This may affect the observed shapes for two reasons. The first is that the age of stellar populations is not uniformly distributed, as star formation takes place in distinct regions such as spiral arms (e.g. Martin & Kennicutt 2001). The second reason is that extinction and reddening by dust in the source galaxy has a stronger impact on observations in blue than in red rest-frame wavelengths and may hence distort the observed shapes more strongly at high than at lower redshifts. This effect is further complicated by the fact that the distribution and overall density of dust evolves with time. In fact, we find strong indications of dust extinction in our data (Fig. D1). However, studying shapes measured from synthetic galaxy images at z = 0.0, we find no strong dependence of the 2D axial ratio distribution on the filter band. This finding suggests that observing shapes in the same filter at different redshift has only a mild effect on the observed shape distribution (see Appendix D). A second systematic effect can be expected from gravitational lensing, while this effect is typically weak with a contribution of less than 1 per cent to the observed ellipticity (e.g. Kirk et al. 2015).

An important outcome of our analysis is that the model for the 3D galaxy shapes and their distribution provides good fits to the observed distribution of 2D shapes for disc-dominated galaxies, despite the expected inaccuracies. This model can hence be employed to generate mock catalogues of galaxy shapes in large cosmological dark matter-only simulations in which stellar mass distributions of individual galaxies are not provided. In order to construct such mocks, Joachimi et al. (2013) approximated disc galaxies by flat opaque cylinders and showed that the resulting 2D ellipticity distribution differs significantly from COSMOS observations. An improvement on this aspect is crucial for building precision mock galaxy catalogues with intrinsic alignments for the preparation of future weak lensing surveys. Disc galaxies are expected to dominate the lensing source samples at high redshifts in which intrinsic alignments contaminate the gravitational shear induced by the large scale structure at lower redshifts (e.g. Fig. 6). Building and testing such improved mocks is subject of our ongoing work.

AVAILABILITY OF DATA

The galaxy catalogues from COSMOS (COSMOS2015, ACS-GC, ZEST) and Illustris TNG used in this work are publicly available at the links provided in Section 2. The matched COSMOS catalogue described in that section is available upon request to K. Hoffmann (kai.d.hoffmann@gmail.com). The Horizon AGN catalogue is available upon reasonable request to N. E. Chisari (n.e.chisari@uu.nl).

ACKNOWLEDGEMENTS

KH acknowledges support by the Swiss National Science Foundation (Grant No. 173716, 198674) and the Forschungskredit grant of the University of Zurich (Projekt K-76106-01-01). CL acknowledges support from the Programme National Cosmology and Galaxies (PNCG) of CNRS/INSU with INP and IN2P3. This work is part of

the Delta ITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). The research of JD is supported by the Beecroft Trust and STFC. Our analysis is based on data products from observations made with ESO Telescopes at the LaSilla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. Our analysis was performed using JUPYTER NOTEBOOK (Kluyver et al. 2016), with the Python packages COLOSSUS (Diemer 2018), EMCEE (Foreman-Mackey et al. 2013), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), PANDAS (Wes McKinney 2010), PYGTC (Bocquet & Carter 2016), SCIPY (Virtanen et al. 2020).

REFERENCES

Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ, 597, 21

Aguerri J. A. L., Balcells M., Peletier R. F., 2001, A&A, 367, 428

Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

Aumer M., Binney J., Schönrich R., 2016, MNRAS, 459, 3326

Baes M., Verstappen J., De Looze I., Fritz J., Saftly W., Vidal Pérez E., Stalevski M., Valcke S., 2011, ApJS, 196, 22

Baldry I. K., Glazebrook K., Brinkmann J., Ivezić Ž., Lupton R. H., Nichol R. C., Szalay A. S., 2004, ApJ, 600, 681

Baugh C. M., Cole S., Frenk C. S., 1996, MNRAS, 283, 1361

Benacchio L., Galletta G., 1980, MNRAS, 193, 885

Bernstein G. M., Jarvis M., 2002, AJ, 123, 583

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

Bertola L., Vietri M., Zeilinger W. W., 1991, ApJ, 374, L13

Binney J., 1978, MNRAS, 183, 501

Binney J., 1985, MNRAS, 212, 767

Binney J., 1992, ARA&A, 30, 51

Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn., Princeton University Press, Princeton, NJ

Binney J., de Vaucouleurs G., 1981, MNRAS, 194, 679

Bocquet S., Carter F. W., 2016, J. Open Source Soft., 1, 46

Bournaud F., 2016, Astrophysics and Space Science Library, 418, 355

Bournaud F., Elmegreen B. G., Martig M., 2009, ApJ, 707, L1

Brook C. B. et al., 2011, MNRAS, 415, 1051

Brook C. B., Stinson G., Gibson B. K., Roškar R., Wadsley J., Quinn T., 2012, MNRAS, 419, 771

Brooks A., Christensen C., 2016, Astrophysics and Space Science Library, 418, 317

Bundy K. et al., 2010, ApJ, 719, 1969

Buta R. J. et al., 2015, ApJS, 217, 32

Camps P., Baes M., 2015, Astronomy and Computing, 9, 20

Carollo D. et al., 2010, ApJ, 712, 692

Chang Y.-Y. et al., 2013, ApJ, 773, 149

Cheng T.-Y., Huertas-Company M., Conselice C. J., Aragón-Salamanca A., Robertson B. E., Ramachandra N., 2021, MNRAS, 503, 4446

Chisari N. et al., 2015, MNRAS, 454, 2736

Correa C. A., Schaye J., Clauwens B., Bower R. G., Crain R. A., Schaller M., Theuns T., Thob A. C. R., 2017, MNRAS, 472, L45

Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371

de Vaucouleurs G., 1959, AJ, 64, 397

Debattista V. P., Mayer L., Carollo C. M., Moore B., Wadsley J., Quinn T., 2006, ApJ, 645, 209

Dekel A., Sari R., Ceverino D., 2009, ApJ, 703, 785

Diemer B., 2018, ApJS, 239, 35

Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497

Dubois Y., Pichon C., Haehnelt M., Kimm T., Slyz A., Devriendt J., Pogosyan D., 2012, MNRAS, 423, 3616

Dubois Y. et al., 2014, MNRAS, 444, 1453

Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948
Eliche-Moral M. C., Balcells M., Aguerri J. A. L., González-García A. C., 2006, A&A, 457, 91

Elmegreen B. G., Bournaud F., Elmegreen D. M., 2008, ApJ, 688, 67

Fasano G., Amico P., Bertola F., Vio R., Zeilinger W. W., 1993, MNRAS, 262, 109

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Freeman K. C., 1970, ApJ, 160, 811

Genel S., Fall S. M., Hernquist L., Vogelsberger M., Snyder G. F., Rodriguez-Gomez V., Sijacki D., Springel V., 2015, ApJ, 804, L40

Georgiou C. et al., 2019, A&A, 622, A90

Gómez F. A., White S. D. M., Grand R. J. J., Marinacci F., Springel V., Pakmor R., 2017, MNRAS, 465, 3446

Governato F. et al., 2010, Nature, 463, 203

Gozaliasl G. et al., 2019, MNRAS, 483, 3545

Graham A. W., Worley C. C., 2008, MNRAS, 388, 1708

Grand R. J. J., Springel V., Gómez F. A., Marinacci F., Pakmor R., Campbell D. J. R., Jenkins A., 2016, MNRAS, 459, 199

Grand R. J. J. et al., 2017, MNRAS, 467, 179

Griffith R. L. et al., 2012, ApJS, 200, 9

Grossi M., Fernandes C. A. C., Sobral D., Afonso J., Telles E., Bizzocchi L., Paulino-Afonso A., Matute I., 2018, MNRAS, 475, 735

Guo R., Hao C.-N., Xia X., Shi Y., Chen Y., Li S., Gu Q., 2020, ApJ, 897, 162

Guzzo L. et al., 2007, ApJS, 172, 254

Harris C. R. et al., 2020, Nature, 585, 357

Hatfield P. W., Laigle C., Jarvis M. J., Devriendt J., Davidzon I., Ilbert O., Pichon C., Dubois Y., 2019, MNRAS, 490, 5043

Häußler B., Barden M., Bamford S. P., Rojas A., 2011, in Evans I. N., Accomazzi A., Mink D. J., Rots A. H., eds, ASP Conf. Ser., Vol. 442, Astronomical Data Analysis Software and Systems XX. Astron. Soc. Pac., San Francisco, p. 155

Heidmann J., Heidmann N., de Vaucouleurs G., 1972, Mem.RAS, 75, 121

Hernquist L., 1992, ApJ, 400, 460

Hodges J. L., 1958, Ark. Mat., 3, 469

Hoffmann K. et al., 2014, MNRAS, 442, 1197

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

Holmberg E., 1946, Meddelanden fran Lunds Astronomiska Observatorium Serie II, 117, 3

Hopkins P. F., Hernquist L., Cox T. J., Dutta S. N., Rothberg B., 2008, ApJ, 679, 156

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

Hopkins P. F., Kereš D., Murray N., Quataert E., Hernquist L., 2012, MNRAS, 427, 968

Hubble E. P., 1926, ApJ, 64, 321

Hubble E. P., 1936, Realm of the Nebulae. Yale University Press, New Haven, Connecticut

Huizinga J. E., van Albada T. S., 1992, MNRAS, 254, 677

Hunter J. D., 2007, Comput Sci Eng., 9, 90

Ilbert O. et al., 2006, A&A, 457, 841

Jackson R. A., Martin G., Kaviraj S., Laigle C., Devriendt J. E. G., Dubois Y., Pichon C., 2020, MNRAS, 494, 5568

Joachimi B., Semboloni E., Bett P. E., Hartlap J., Hilbert S., Hoekstra H., Schneider P., Schrabback T., 2013, MNRAS, 431, 477

Kauffmann G., White S. D. M., Guiderdoni B., 1993, MNRAS, 264, 201

Kautsch S. J., 2009, PASP, 121, 1297

Kautsch S. J., Grebel E. K., Barazza F. D., Gallagher J. S. I., 2006, A&A, 445, 765

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kepner J. V., 1999, ApJ, 520, 59

Kirk D. et al., 2015, Space Sci. Rev., 193, 139

Kluyver T. et al., 2016, in Loizides F., Schmidt B., eds, Proc. 20th Elpub, Positioning and Power in Academic Publishing: Players, Agents and Agendas. IOS Press BV, Netherlands, p. 87

Komatsu E. et al., 2011, ApJS, 192, 18

Kormendy J., 2016, Astrophysics and Space Science Library, 418, 431

Kormendy J., Kennicutt Robert C. J., 2004, ARA&A, 42, 603

Kormendy J., Drory N., Bender R., Cornell M. E., 2010, ApJ, 723, 54

Kretschmer M., Agertz O., Teyssier R., 2020, MNRAS, 497, 4346

Kroupa P., Haslbauer M., Banik I., Nagesh S. T., Pflamm-Altenburg J., 2020, MNRAS, 497, 37

Laigle C. et al., 2016, ApJS, 224, 24

Laigle C. et al., 2019, MNRAS, 486, 5104

Lambas D. G., Maddox S. J., Loveday J., 1992, MNRAS, 258, 404

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

Lilly S. J. et al., 2007, ApJS, 172, 70

Marinacci F. et al., 2018, MNRAS, 480, 5113

Martig M., Bournaud F., Teyssier R., Dekel A., 2009, ApJ, 707, 250

Martin C. L., Kennicutt Robert C. J., 2001, ApJ, 555, 301

Martin G., Kaviraj S., Devriendt J. E. G., Dubois Y., Pichon C., 2018, MNRAS, 480, 2266

Martin G., Kaviraj S., Hocking A., Read S. C., Geach J. E., 2020, MNRAS, 491, 1408

Martínez-Lombilla C., Knapen J. H., 2019, A&A, 629, A12

Mosenkov A. V., Sotnikova N. Y., Reshetnikov V. P., Bizyaev D. V., Kautsch S. J., 2015, MNRAS, 451, 2376

Mosenkov A. et al., 2020, MNRAS, 494, 1751

Naab T., Burkert A., 2003, ApJ, 597, 893

Naiman J. P. et al., 2018, MNRAS, 477, 1206

Negroponte J., White S. D. M., 1983, MNRAS, 205, 1009

Nelson D. et al., 2018, MNRAS, 475, 624

Noerdlinger P. D., 1979, ApJ, 234, 802

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

Park M.-J. et al., 2019, ApJ, 883, 25

Park M. J. et al., 2021, ApJS, 254, 2

Peebles P. J. E., 2020, MNRAS, 498, 4386

Peng C. Y., Ho L. C., Impey C. D., Rix H.-W., 2002, AJ, 124, 266

Peschken N., Łokas E. L., Athanassoula E., 2020, MNRAS, 493, 1375

Pillepich A. et al., 2018a, MNRAS, 473, 4077

Pillepich A. et al., 2018b, MNRAS, 475, 648

Pillepich A. et al., 2019, MNRAS, 490, 3196

Planck Collaboration, 2016, A&A, 594, A13

Planck Collaboration, 2018, A&A, 641, 67

Quinn P. J., Hernquist L., Fullagar D. P., 1993, ApJ, 403, 74

Ravindranath S. et al., 2004, ApJ, 604, L9

Reshetnikov V. P., Mosenkov A. V., Moiseev A. V., Kotov S. S., Savchenko S. S., 2016, MNRAS, 461, 4233

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

Robertson B., Bullock J. S., Cox T. J., Di Matteo T., Hernquist L., Springel V., Yoshida N., 2006, ApJ, 645, 986

Rodríguez S., Padilla N. D., 2013, MNRAS, 434, 2153

Rodríguez-Gomez V. et al., 2019, MNRAS, 483, 4140

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

Ryden B. S., 2006, ApJ, 641, 773

Sachdeva S., 2013, MNRAS, 435, 1186

Sachdeva S., Saha K., 2016, ApJ, 820, L4

Saha K., Tseng Y.-H., Taam R. E., 2010, ApJ, 721, 1878

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

Sargent M. T. et al., 2007, ApJS, 172, 434

Satoh Y. K., Kajisawa M., Himoto K. G., 2019, ApJ, 885, 35

Scarlata C. et al., 2006, preprint(astro-ph/0611644)

Scarlata C. et al., 2007, ApJS, 172, 494

Scoville N. et al., 2007, ApJS, 172, 1

Snyder G. F. et al., 2015, MNRAS, 454, 1886

Springel V., 2010, MNRAS, 401, 791

Springel V., Hernquist L., 2005, ApJ, 622, L9

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726

Springel V. et al., 2018, MNRAS, 475, 676

Stark A. A., 1977, ApJ, 213, 368

Strateva I. et al., 2001, AJ, 122, 1861

Tacchella S. et al., 2019, MNRAS, 487, 5416

Takeuchi T. M., Ohta K., Yuma S., Yabe K., 2015, ApJ, 801, 2
Teyssier R., 2002, A&A, 385, 337

Toomre A., 1977, in Tinsley B. M., Larson Richard B., Gehret D. C., eds, Evolution of Galaxies and Stellar Populations, Proceedings of a Conference at Yale University. Yale University Observatory, New Haven, Conn, p. 401

Trujillo I., Bakos J., 2013, MNRAS, 431, 1121

Übler H., Naab T., Oser L., Aumer M., Sales L. V., White S. D. M., 2014, MNRAS, 443, 2092

van den Bosch F. C., 1998, ApJ, 507, 601

van der Kruit P. C., Searle L., 1982, A&A, 110, 61

van der Wel A. et al., 2014, ApJ, 792, L6

Vaucouleurs G., 1959a, Handbuch der Physik, 11, 275

Vaucouleurs G., 1959b, Handbuch der Physik, 11, 311

Vincent R. A., Ryden B. S., 2005, ApJ, 623, 137

Virtanen P. et al., 2020, Nature Methods, 17, 261

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42

Walker I. R., Mihos J. C., Hernquist L., 1996, ApJ, 460, 121

Weinzirl T., Jogee S., Khochfar S., Burkert A., Kormendy J., 2009, ApJ, 696, 411

Wes McKinney, 2010, in van der S., Millman W. J., eds, Proceedings of the 9th Python in Science Conference. p. 56

Yuma S., Ohta K., Yabe K., Kajisawa M., Ichikawa T., 2011, ApJ, 736, 92

Yuma S., Ohta K., Yabe K., 2012, ApJ, 761, 19

Zhang H. et al., 2019, MNRAS, 484, 5170

## APPENDIX A: TESTING THE CLASSIFICATION OF DISC-DOMINATED LATE-TYPE GALAXIES IN THE ZEST

We extend our visual inspection from Fig. 1, now including objects of any orientation and any redshift from our volume limited samples to further validate if these objects have been classified correctly as disc-dominated late-type galaxies.

In Fig. A1, we show 16 randomly selected examples of galaxies for each of our 6 stellar mass-redshift sub-samples. We find visually that the majority of objects appears to be classified correctly. However, we find several objects with irregular shapes, in particular in the high redshift bins, which could potentially be misclassified.

As an additional validation of the morphological classification we show in Fig. A2 the distributions of the Sérsic index *n* from equation (3), as provided in the ZEST catalogue. Results are shown for objects which are classified as late-type disc-dominated, late-type bulge-dominated and early-type galaxies (see Section 2.1.3 for details on the classification). We find that the Sérsic indices of late-type disc-dominated galaxies are distributed in the range $0 \le n \le 2$, while the peak of the distribution is located at $n \simeq 0.5$, which is consistent with values that are typically found for disc-dominated galaxies. The Sérsic indices of early-type galaxies are distributed between $2 \le n \le 7$ and are hence clearly separated from those of late-type disc-dominated galaxies. The Sérsic indices of late-type, bulge-dominated galaxies lie in between the two other distributions with significant overlaps in both directions. Note that such overlaps can be expected as the Sérsic index alone is not sufficient for a detailed morphological classification.

In summary, we conclude that the low Sérsic indices for late-type disc-dominated galaxies and the larger Sérsic indices of late-type bulge-dominated and early-type galaxies are consistent with the morphological ZEST classification. However, the image inspection revealed some poorly resolved objects at high redshifts, which are potentially misclassified as late-type disc-dominated galaxies and could hence contaminate the samples on which our investigation is based. Addressing this problem will require higher resolution imaging surveys of high redshift galaxies, which will become achievable with future space based telescopes, such as the James Webb Space Telescope. The reliability of the ZEST morphological classification could further be tested with a comparison against alternative classification methods already with current data (e.g. Martin et al. 2020; Cheng et al. 2021).
<figure><img src="image_20.png" /><figcaption>**Figure A1.** Object classified as disc-dominated late-type galaxies in the ZEST catalogue of the COSMOS field. The figure shows 16 randomly selected galaxies for each of our 6 redshift-stellar mass samples.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure A2.** Sérsic indices *n* of galaxies in our volume limited main sample with different morphological classifications provided in the ZEST catalogue. The different panels show results from the different redshift bins studied in our analysis. The vertical line is drawn at *n* = 1.</figcaption></figure>

## APPENDIX B: FROM 3D TO 2D AXIS RATIOS

We obtain the 2D axis ratios of a projected 3D ellipsoidal stellar system following Joachimi et al. (2013), under the assumption that the 3D system is absorption-free, self-similar and coaxial. A 3D ellipsoid is thereby expressed in a coordinate system defined by the two orthogonal unit vectors {$\hat{e}_u, \hat{e}_v$} which span the projection plane and the unit vector $\hat{e}_\parallel$, which is orthogonal to {$\hat{e}_u, \hat{e}_v$}, pointing along the observer's line of sight. In this new coordinate system, the principal axes are given by $\tilde{S}_\mu = \{\langle\hat{e}_u S_\mu\rangle, \langle\hat{e}_v S_\mu\rangle, \langle\hat{e}_\parallel S_\mu\rangle\}^\tau \equiv \{\tilde{S}_{u,\mu}, \tilde{S}_{v,\mu}, \tilde{S}_{\parallel,\mu}\}^\tau$ with $S_\mu \in \{A_{3D}, B_{3D}, C_{3D}\}$. The projected 2D ellipse is given by all points *r* in the projection plane which fulfil $r^\tau W^{-1}r = 1$, where

$$W^{-1} = \sum_{\mu=1}^{3} \frac{\tilde{S}_{\perp,\mu} \tilde{S}_{\perp,\mu}^{\tau}}{\tilde{S}_{\mu}^{2}} - \frac{kk^{\tau}}{\alpha^{2}}, \qquad (B1)$$

with

$$k = \sum_{\mu=1}^{3} \frac{\tilde{S}_{\parallel,\mu} \tilde{S}_{\perp,\mu}}{\tilde{S}_{\mu}^{2}} \quad \text{and} \quad \alpha^{2} = \sum_{\mu=1}^{3} \left( \frac{\tilde{S}_{\parallel,\mu}}{\tilde{S}_{\mu}} \right)^{2} \qquad (B2)$$

and $\tilde{S}_{\perp,\mu} \equiv \{\tilde{S}_{u,\mu}, \tilde{S}_{v,\mu}\}^\tau$ is the principal axes component in the projection plane. The 2D ellipticity vector of the projected ellipsoid

is then given by

$$\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \end{pmatrix} = \frac{1}{\mathcal{N}} \begin{pmatrix} W_{11} - W_{22} \\ 2 W_{12} \end{pmatrix} \qquad (B3)$$

with $\mathcal{N} \equiv W_{11} + W_{22} + \sqrt{\det \tilde{W}}$. The absolute value of the ellipticity, $\epsilon = \sqrt{\epsilon_1^2 + \epsilon_2^2}$, is related to the 2D axial ratio $q_{2D} \equiv B_{2D}/A_{2D}$ of the projected ellipsoid as

$$q_{2D} = \frac{1 - \epsilon}{1 + \epsilon}, \qquad (B4)$$

where $A_{2D}$ and $B_{2D}$ are the principle axes of the 2D ellipse. An alternative approach for obtaining $q_{2D}$ has been derived by Binney (1985) (see also Stark 1977; Benacchio & Galletta 1980). However, we find this latter calculation to be computationally less efficient.

## APPENDIX C: VALIDATING THE MODEL FOR THE 3D AXIAL RATIO DISTRIBUTION

We use the shapes of disc galaxies, measured in the HAGN and TNG100 simulations to further validate the approximations that are implied by using the Gaussian model for the 3D axial ratio distribution, $P(q_{3D}, r_{3D})$, introduced in Section 3.1.

A possible improvement of that model could be the inclusion of a skewness in the $q_{3D}$ dimension, as demonstrated in Fig. 11. In the top panel of Fig. C1, we show the joint distribution of $q_{3D}$ and $r_{3D}$

<figure><img src="image_22.png" /><figcaption>**Figure C1.** Top: model for the 3D axial ratio distribution with and without skewness in the $q_{3D}$ dimension (blue and red contours, respectively). Bottom: Corresponding predictions for the 2D axis ratios distribution. The skewness as a very small affect on the 2D axis ratios.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure C2.** Distributions of the axis ratios $q_{3D} \equiv C_{3D}/A_{3D}$ and $r_{3D} \equiv C_{3D}/A_{3D}$ of disc galaxies in the HAGN and TNG100 simulations at $z = 1.0$. These distributions are compared to the $P(q_{3D}, r_{3D})$ model from equation (5), which was fitted to the marginalized distributions of $q_{3D}$ and $r_{3D}$ as shown in Fig. 11.</figcaption></figure>

for the truncated Gaussian model with and without skewness, using parameters which are typical for the disc galaxies in our analysis. The corresponding distribution of the apparent 2D axis ratios are shown in the bottom panel of the same figure. We find that the effect of the skewness is small compared to the errors on our measurements (see Fig. 15). As a consequence the skewness in the $q_{3D}$ distribution cannot be constrained from our observed data and is therefore not included in our $P(q_{3D}, r_{3D})$ model.

A second possible improvement of our model could be the inclusion of a correlation between $q_{3D}$ and $r_{3D}$ which could be

described by a covariance matrix in equation (5). In order to assess the necessity for such a model extension, we inspect the kernel density estimates of the $P(q_{3D}, r_{3D})$ distribution from HAGN and TNG100 in Fig. C2. We find no evidence for a significant correlation between $q_{3D}$ and $r_{3D}$ in both simulations. For a visual impression of the accuracy of our Gaussian model we compare it to the measurements in Fig. C2, using the model parameters from the fits to the marginalized 3D axial ratio distributions, shown in Fig. 11. The model describes the simulation data reasonably well, while the relatively weak deviations from the measurements appear to result from the neglected skew rather from a neglected covariance.

## APPENDIX D: BIAS IN THE 2D AXIAL RATIO DISTRIBUTIONS FROM APPARENT MAGNITUDE CUTS

We study the effect of dust extinction on the colour and apparent magnitude of the disc-dominated galaxies from our matched catalogue with $M_i < -21.5$ in Fig. D1. In this figure, we display the apparent *i*-band magnitude and the *i* − *j* colour index against the apparent 2D axial ratio $q_{2D}$ in three redshift bins.

We see at all redshifts that galaxies with small apparent axis ratios are significantly redder (i.e. higher colour index) and dimmer than those with apparent axis ratios close to unity. These effects can be expected from the extinction by dust in the interstellar medium of the source galaxy, as the pathway of light through the dust of the source towards the observer is longer for discs which are seen edge-on than for face-on objects (i.e. $q_{2D} \ll$ and $q_{2D} \simeq 1$, respectively). As a consequence, dust extinction can shift discs with inclined orientations below the apparent magnitude limit (marked as solid red horizontal line in Fig. D1), in particular at high redshifts (e.g. Binney & de Vaucouleurs 1981; Huizinga & van Albada 1992). This effect can introduce a redshift-dependent bias in the observed distribution of axis ratios towards apparently round (face-on) discs, which could be mistaken for an evolution in the intrinsic shape distribution, if not taken properly into account. We demonstrate this effect by applying a cut at $m_i = 23$, shown as dashed blue horizontal line in Fig. D1. The apparent evolution of the axial ratio distribution introduced by this cut can be seen in the bottom panels of Fig. D1.

An additional problem caused by this selection effect is that orientations of galaxies in samples affected by the apparent magnitude cut cannot be expected to be randomly distributed with respect to the observer, which violates a basic requirement for the 3D shape reconstruction method employed in this work. In order to mitigate these biases, we select objects with absolute *i*-band magnitudes below −21.5, which appear not to be affected by the apparent magnitude cut of $m_i < 24$ for $z < 1.0$, used for selecting our volume limited sample.
<figure><img src="image_24.png" /></figure>

**Figure D1.** Apparent 2D axis ratios $q_{2D}$ of disc-dominated galaxies in COSMOS with absolute *i*-band magnitudes brighter than $M_i = -21.5$ versus the rest-frame colour index and the apparent *i*-band magnitude (top and central panels respectively). The solid red lines in the central panels indicate the apparent magnitude cut of our volume limited sample at $m_i = 24$. The dashed blue lines mark a cut at $m_i = 23$. The corresponding $q_{2D}$ distributions, shown in the bottom panels, display how the apparent magnitude cut can bias the $q_{2D}$ distribution, in particular at high redshifts.

## APPENDIX E: DEPENDENCE OF GALAXY SHAPES ON COLOUR

The galaxy shape measurements used in this work are based on ACS imaging in the F814W filter. This filter corresponds to a different rest-frame wavelength range at each source redshift, which can affect the observed shape of a given galaxy if its colour is not uniformly distributed (for instance, due to extinction and reddening by dust or patchy star formation). We study the impact of this systematic effect on the distribution of apparent axis ratios measured from second-

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

**Figure E1.** Filters of the SDSS *g*- and *i*-bands and the ACS F814W filter used for the imaging in the COSMOS survey (dashed and solid lines, respectively). The F814W filter covers the near-infrared wavelengths at $z = 0.0$, which corresponds to red and green wavelengths at the lowest and highest redshift bin used in our analysis (with mean redshifts of $z = 0.56$ and $z = 95$, respectively).

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

**Figure E2.** Distribution of apparent 2D axis ratios, measured from synthetic images of disc galaxies in the Illustris TNG100 simulation at $z = 0.0$ in the SDSS *g*- and *i*-bands.

order moments in synthetic images of disc galaxies from the TNG100 simulations at $z = 0.0$, produced by Rodriguez-Gomez et al. (2019, see section 2.2.3). These measurements were performed in the SDSS *i*- and *g*-bands, which correspond roughly to the ACS F814 band at $z = 0.0$ and $0.5 < z < 1.0$, respectively, as illustrated in Fig. E1. In Fig. E2, we show that the change of the apparent axial ratio distribution is weak, compared to the shot-noise errors which we expect for our COSMOS samples (see Fig. 15). We find the same result when using PSF corrected axis ratios obtained by Rodriguez-Gomez et al. (2019) from Sérsic profile fits to the same synthetic images. These findings line up with those from Georgiou et al. (2019), who report a minor difference between galaxy ellipticities measured in different filters of the KiDS survey. They are further consistent with from Ryden (2006), who shows that axis ratios measured in the $K_s$- and *B*-band are strongly scattered, but not biased. We therefore do not expect our results to be significantly affected by the fact that
the ACS F814W filter probes different rest-frame wavelengths at different redshifts.

APPENDIX F: DEPENDENCE OF PARAMETER CONTOURS ON NOISE IN THE 2D AXIAL RATIO DISTRIBUTION

We aim to test in this section, if the strong variations of the posteriors from different redshifts samples, shown in the central and bottom panels of Fig. 13, can be expected if the axis ratios of galaxies in these samples are drawn randomly from the same redshift independent axial ratio distribution.

We tackle this question with a numerical experiment which is based on two sets of $10^6$ artificial ellipsoids whose 3D axial ratio distributions we generate with our Gaussian model from equation (5) and (6). The parameters used to generate these two samples are the best-fitting values for the high and low-mass samples in the entire redshift range of the COSMOS main sample, given in Table 5. These low and high-mass parameters are referred to as parameter set 1

<figure><img src="image_27.png" /><figcaption>**Figure F1.** 2D axial ratio distributions, generated from our model for the 3D axial ratio distribution for two different sets of parameters. The top panels show results for the parent samples, consisting of $10^6$ objects. The bottom panels show results from three random sub-samples (a, b, c) of the parent sample, consisting of 500 objects each. Red lines show fits from 3D axial ratio distribution model to each sub-sample.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure F2.** Posteriors of the parameters of the Gaussian model for the 3D axial ratio distribution, derived from the 2D axis ratios of projected ellipsoids, whose 3D axis ratios follow the same Gaussian model (Fig. F1). Results are shown for the three random sub-samples (a,b,c). Light and dark areas mark 95 per cent and 68 per cent confidence levels respectively. The parameter values used for generating the parent samples are marked as horizontal and vertical lines. Top and bottom panels show results for the model parameter sets 1 and 2, respectively.</figcaption></figure>

and 2, respectively. The distributions of 2D axis ratios from the corresponding samples are derived assuming random orientations and are shown in the top panels of Fig. F1. From each of the two sets we draw three sub-sets with 500 randomly selected objects, which corresponds roughly to the size of our observational samples. The 2D axial ratio distributions of these sub-sets are shown in the three lower panels of Fig. F1. We proceed by performing the same MCMC parameter inference as for the observational samples, described in Section 3.2. The best-fitting model prediction for each sub-set is shown as red line in Fig. F1. The corresponding posterior distributions are shown in Fig. F2. As for the observational data we
find strong variations across results for the different sub-sets. This
finding can be expected, since the posterior is computed directly
from the data. Variations in the data hence translate into variations in
the posterior. The true values of the input model parameters, shown
in Fig. F2 by horizontal and vertical lines, lie within the 95 per cent
confidence intervals, but often outside of the 68 per cent intervals.

This test shows that strong variations of the posteriors can be expected
for samples that are drawn from the same underlying distribution if
the sample size is similar to the sizes of our observational redshift-
stellar mass samples.

This paper has been typeset from a T\_EX/L\_TEX file prepared by the author.