IllustrisTNG in the HSC-SSP: No Shortage of Thin Disk Galaxies in TNG50

The HSC-SSP (Aihara et al., 2018) is an optical imaging survey that uses the HSC (Miyazaki et al., 2012, 2018) of the $8.2$-meter Subaru Telescope. HSC-SSP has released three public data releases (PDRs; Aihara et al., 2018, 2019, 2022) and covers about $1000\,\mathrm{deg^{2}}$ in $grizy$ bands in the WIDE layer. One can reliably extract the surface brightness profiles of galaxies down to $28.5\,\mathrm{mag\,arcsec^{-2}}$ in $r$ band and $i$ band using PDR2 WIDE images (Huang et al., 2018; Li et al., 2022). The best spatial resolution is in the $i$-band, with a median seeing of $0\farcs 61$ for the PDR3 WIDE layer. The HSC has a pixel scale of $0\farcs 168$.

The HSC-SSP photometric catalog¹¹1The photometric catalogs are described in https://hsc-release.mtk.nao.ac.jp/schema/. has incorporated publicly available spectroscopic redshifts (spec-$z$) from various sources such as the Galaxy And Mass Assembly (GAMA; Driver et al., 2009, 2011) DR2 (Liske et al., 2015) and the Sloan Digital Sky Survey (SDSS; York et al., 2000) DR16 (Ahumada et al., 2020). A homogenized flag in the photometric catalog indicates which spec-$z$ has been adopted out of the various available sources, opting for the spec-$z$ with the smallest error (Aihara et al., 2018). We select all primary objects with known spec-$z$ by querying the photometric catalog with $\mathtt{isprimary=True}$ and $\mathtt{specz\_redshift}>0$. Galaxies are selected from the photometric catalog based on $\mathtt{extendedness=1}$. Furthermore, we select full-depth galaxy images of superior quality based on the following flags: $\mathtt{g\_inputcount\_value\geq 4}$, and $\mathtt{(i\&y)\_inputcount\_value\geq 5}$, NOT pixelflags_edge, NOT pixelflags_offimage, NOT mask_brightstar_halo, NOT mask_brightstar_blooming. We cross-match the HSC-SSP photometric catalog with the stellar mass and size measurements from GAMA (Taylor et al., 2011) and SDSS (Ahn et al., 2014). The stellar mass for a galaxy whose spec-$z$ from either SDSS or GAMA, is obtained from the measurements of the corresponding source. Taylor et al. (2011) estimate the stellar masses using the Chabrier (2003) initial mass function and Bruzual & Charlot (2003) stellar population synthesis models, while Ahn et al. (2014) use the Kroupa (2001) initial mass function and adopts the Flexible Stellar Population Synthesis code (FSPS; Conroy et al., 2009). We confirm a systematic offset between the two measurements of stellar mass, as reported in Kado-Fong et al. (2020a, b). The SDSS stellar masses are systematically higher than the GAMA ones, with a mean offset of $0.17$ dex and a scatter of $0.26$ dex estimated from the overlapped galaxies. To align the measurements and to project consistency with the generation of the synthetic images (Section 2.2), we subtract the offset of $0.17$ dex from the SDSS stellar masses. To compare with the TNG50 simulation, we select our main sample of 4657 galaxies within a redshift range of $0.005\leq z\leq 0.05$ and with a stellar mass range of $9\leq\log(M_{\star}/M_{\sun})\leq 11.5$. Furthermore, to gain a comprehensive view of the axis ratio distribution including dwarf galaxies, we also include galaxies with smaller stellar masses of $8.5\leq\log(M_{\star}/M_{\sun})<9$.

In this work, we use the $g$, $i$, and $y$-band images from the HSC-SSP PDR3 WIDE layer. We extract cutout images with a size of $15\,r_{90}\times 15\,r_{90}$, where $r_{90}$ is the radius within which the integrated flux constitutes $90\%$ of the total galaxy flux. The value of $r_{90}$ is obtained from the same source as the spec-$z$ (either SDSS or GAMA). We also retrieve point spread function (PSF) images at the center of each individual galaxy. We use the command-line SQL tool²²2https://hsc-gitlab.mtk.nao.ac.jp/ssp-software/data-access-tools/-/tree/master/pdr3 to download all image products.

The TNG project³³3http://www.illustris-project.org/ is a suite of magneto-hydrodynamic cosmological simulations of the formation and evolution of various elements, including gas, stars, dark matter, supermassive black holes, and magnetic fields from $z=127$ to $z=0$ (Weinberger et al., 2017; Pillepich et al., 2018; Nelson et al., 2019a, b; Pillepich et al., 2019). The TNG suite consists of three simulations, respectively referred to as TNG50, TNG100, and TNG300. They run at $V_{\mathrm{box}}$ volumes of $51.7^{3}$, $106.5^{3}$, and $302.6^{3}$ cMpc³, respectively. Two commonly used resolution indicators are (1) the Plummer-equivalent gravitational softening length for collisionless particles, $\epsilon_{\star,\mathrm{dm}}=(0.29,0.74,1.48)$ ckpc, and (2) the average size of star-forming gas cells, $\bar{r}_{\mathrm{cell,sf}}=(0.138,0.355,0.715)$ ckpc.

We use the synthetic HSC-SSP WIDE images generated by Bottrell et al. (2024) from the TNG50 simulation. Our work focuses on simulated galaxies in snapshot 99 ($z=0$) mock observed at $z=0.05$. The photo-realistic but noise-free (idealized) synthetic images in HSC $grizy$ bands were first generated using SKIRT Monte Carlo radiative transfer (Camps & Baes, 2020, SKIRT version 9). Using a bespoke version of RealSim⁴⁴4https://github.com/cbottrell/RealSim for the HSC-SSP (Bottrell et al., 2019), the idealized images made by SKIRT were then convolved with reconstructed HSC-SSP seeing kernels and statistically injected into final-depth HSC-SSP cutouts so that these synthetic images have consistent properties with observed HSC-SSP WIDE layer images of real galaxies. The TNG50 snapshot 99 includes $3092$ galaxies with $9\leq\log(M_{\star}/M_{\sun})\leq 11.5$, and their synthetic HSC-SSP galaxy images are generated for each TNG50 galaxy at four different viewing angles, resulting in a final sample of $12368$ synthetic images. We refer readers to Bottrell et al. (2024) for a detailed description of the procedure.

Figure 1 compares the stellar mass distributions of the observed HSC-SSP sample ($0.005\leq z\leq 0.05$) and the TNG50 sample (snapshot 99 mock observed at $z=0.05$).

We use the sep package⁵⁵5https://sep.readthedocs.io/en/v1.1.x (Barbary et al., 2016), a Python implementation of SourceExtractor, to detect and deblend objects using the cold + hot approach (e.g., Rix et al., 2004). We run sep with two sets of parameters. The hot mode detects bright objects and ensures a complete segmentation of the target galaxy by using a $5\sigma$ detection threshold and $\mathtt{deblend\_cont}=0.1$. The cold mode detects faint galaxies and stars by using a $1\sigma$ detection threshold and $\mathtt{deblend\_cont}=0.001$. We first obtain the position and size information of the target from the hot mode. Within $3$ times the Kron radius of the target, we remove all the objects detected by the cold mode, except for those that cover more than $80\%$ of the area of the target galaxy in the hot mode. Then, we remove the target from both the cold mode and the hot mode, and combine the remaining objects to create masks that we use for further analysis. We dilate the object masks based on the distance of each object from the galaxy center. For objects located more than $6\,r_{50}$ away from the target center, where $r_{50}$ is the radius within which $50\%$ of the total flux is enclosed, we apply a Gaussian filter with a sigma of $3$ times the median FWHM of PSFs in each band. For those objects located within $6\,r_{50}$ of the center, we use an elliptical aperture mask whose semi-major (minor) axis lengths are the second order moments of the object’s flux along the major (minor) axis given by sep. We combine the mask with the bad pixel mask (either ‘NO_DATA’ or ‘SAT’ mask planes) provided by HSC-SSP PDR3 to form the final masks.

The above analysis is performed for both real and synthetic HSC-SSP images. Moreover, the final masks for the TNG50 synthetic images also include the detected source footprints as part of the original HSC-SSP image products, to err on the safe side despite that they should mostly overlap with the sep masks. This extra step does not apply to the real images because the footprints include the galaxy of interest and are not deblended.

AutoProf⁶⁶6https://autoprof.readthedocs.io/en/latest is an automatic nonparametric tool for fitting isophote and extracting surface brightness profiles of extended sources (Stone et al., 2021). The code has several modular steps that can be customized by the user to form a desired pipeline of data processing. For our application, the full feature of AutoProf is more than enough. We only use the initialization step to measure the global axis ratio ($q$) and position angle (PA) of the galaxies. We supply the masks prepared in Section 3.1 and PSF FWHMs obtained via 2D Gaussian fits as input for AutoProf. Moreover, we need quality flags produced by the checking fit step to help discard poor fits. Below, we provide a brief overview of these two steps. For further details, please refer to Sections 2.4 and 2.7 of Stone et al. (2021).

The initialization is done in two steps. First, fast Fourier transform (FFT) is performed over fluxes along circular paths concentric with the galaxy, and the phases are determined from the second Fourier coefficient $\mathcal{F}_{2}$. The global PA is the average direction of the phases for the five outermost circular paths whose average fluxes are still above the background level. Second, with the global PA being fixed, the global axis ratio is obtained by minimizing $\left|\mathcal{F}_{2}/\left(\tilde{\mu}+\sigma_{b}\right)\right|$, where $\tilde{\mu}$ is the median flux, $\sigma_{b}$ is the background noise, and $\tilde{\mu}+\sigma_{b}$ is approximately equal to $\mathcal{F}_{0}$.

The checking fit step performs four checks on the quality of the fits by examining the flux distribution along the best-fit isophotes: (1) and (2) the magnitudes of the first and the second Fourier mode, respectively; (3) the variability of fluxes; and (4) the variability of fluxes compared with that along the initialized isophotes. The first check failing due to significant power in the first Fourier mode often indicates a failed center determination or lopsidedness in galaxy flux distribution, and the second suggests that the optimization was not achieved or the presence of elongated structures (e.g., bars). The third check compares the interquartile range of the fluxes to their median and fails when a minimum fraction of the isophotes have such high variabilities. The fourth check is similar to the third, only that the variability is contrasted against that obtained from the initialized isophotes. The thresholds for triggering the failures were chosen by Stone et al. (2021) to catch visually identified cases.

Upon visual inspection of a random sample of 500 galaxies, we find that the fits that have passed at least three checks are of good quality. Hence, galaxies that failed more than one check are removed from the sample, which constitutes a small fraction of our sample (e.g., $\sim 3\%$ for HSC-SSP $i$-band observations). Additionally, we implement a resolution cut to ensure sufficient spatial resolution for the initialized isophotes where $q$ was measured. Specifically, we required a semi-minor axis length of the initialized isophote larger than $1.5$ times the PSF FWHM, i.e., $R\times q\geq 1.5\times\mathrm{FWHM}$, where $R$ is the semi-major axis length of the initialized isophote. The resolution cut has minimal impact on the sample, even in the case of the $y$-band images, which exhibit the worst seeing statistics (see Appendix A for details).

In this section, we describe how we measure the concentration index from the (synthetic) HSC-SSP images and how the spheroidal-to-total ratio (S/T) is derived from stellar kinematics in the simulation.

We adopt the curve-of-growth method to derive the concentration index for both the HSC-SSP galaxies and the TNG50 galaxies. We use AutoProf to extract 1D, azimuthally averaged surface brightness profiles of the galaxies. Subsequently, based on the growth curve of their flux truncated at the limiting surface brightness, we derive $r_{20}$ and $r_{80}$, where $r_{20}$ ($r_{80}$) represents the radius within which the integrated flux constitutes $20\%$ ($80\%$) of the total galaxy flux. We define the concentration index as $C_{82}=5\times\log{(r_{80}/r_{20})}$. We exclude galaxies with $r_{20}<1.5\times{\rm FWHM}$ to ensure sufficient spatial resolution for accurate measurements. The index characterizes the degree of flux concentration of a galaxy and is often used to separate bulge-dominated and disk-dominated galaxies (e.g., Conselice, 2003, 2014).

S/T of a galaxy quantifies the fraction of its stellar mass in the spheroidal (kinematically hot) components. We follow the definition of Tacchella et al. (2019) that stars of spheroidal components have $\epsilon\leq 0.7$, where $\epsilon$ is the circularity of stellar particles. Genel et al. (2015) provided the fractional mass of stars with $\epsilon>0.7$ computed within ten times the stellar half-mass radii of TNG galaxies as CircAbove07Frac⁷⁷7https://www.tng-project.org/data/docs/specifications/#sec5c. So, we simply derive $S/T=1-\mathtt{CircAbove07Frac}$.

We aim to conduct a fair comparison of axis ratios between the HSC-SSP galaxies and the TNG50 galaxies. We apply the same method defined in Section 3.2 to both the synthetic HSC-SSP images of TNG50 galaxies and the real HSC-SSP images, in order to obtain comparable measurements of galaxy shapes. The comparison is performed within the stellar mass range of $10\leq\log(M_{\star}/M_{\sun})\leq 11.5$. The upper limit of the stellar mass range was set to ensure nonzero presence of galaxies in the highest stellar mass bins. We quantify the difference between the observation and the simulation using chi-square ($\chi^{2}$) statistics and the associated statistical significance. The formulation of $\chi^{2}$ that describes the difference between the axis ratio distributions of the observed and simulated galaxies is given below.

where $N_{\mathrm{sim},i}$ denotes the number of simulated galaxies in the $i$th axis ratio bin and $N^{\prime}_{\mathrm{obs},i}$ the weighted number of observed galaxies therein; $\sigma_{i}$ are their respective Poisson uncertainties. Following the modified Neyman’s $\chi^{2}$, when counts of galaxies are zero, a minimum counts of one is enforced to avoid the situation of $\sigma_{i}=0$. $N_{\mathrm{obs},i,j}$ is the number of observed galaxies in both the $i$th axis ratio bin and the $j$th stellar mass bin. We will clarify the definition of their weights $w_{\mathrm{obs},j}$ when describing the details of the stellar-mass weighting scheme. $N_{{\rm bins},q}=20$ is the number of axis ratio bins with a fixed width of 0.05 and $N_{{\rm bins},M_{\star}}=10$ is the number of stellar mass bins with a fixed width of 0.15 dex.

It is worth noting that we construct the $\chi^{2}$ test statistic based on the difference in the number of galaxies between observation and simulation, rather than using the ratio of the number of galaxies in each axis ratio bin to the total number, as done in a previous study by H22. This approach is chosen because, in the former case, we can approximate the Poisson counts by a Gaussian distribution, allowing the sum of the test statistics to follow a $\chi^{2}$ distribution. In contrast, when using ratios, the total number of galaxies itself is a Poisson random variable. This complicates the distribution of the ratios, making it difficult to approximate their sum with a $\chi^{2}$ distribution.

The different stellar mass distributions of the observed and the simulated galaxies should be taken into account when comparing their axis ratios (see Figure 1), as galaxy shape is known to depend on stellar mass. We follow the same stellar mass-weighting method in H22 to correct the axis ratio distribution of the observed galaxies before comparing it with that of the simulated galaxies. We calculate a weight for each observed galaxy as the number ratio of simulated galaxies to observed galaxies in the corresponding stellar mass bin. Specifically, the weight of an observed galaxy in the $j$th stellar mass bin is defined as $w_{\mathrm{\rm obs},j}=N_{\mathrm{sim},j}/N_{\mathrm{obs},j}$. In the $i$th axis ratio bin, $N^{\prime}_{\mathrm{obs},i}=\sum_{j}{w_{\mathrm{obs},j}N_{\mathrm{obs},i,j}}$ represents the total weight of galaxies therein, a summation of their $w_{{\rm obs},j}$. So, $N^{\prime}_{\mathrm{obs},i}$ is the adjusted number of the observed HSC-SSP galaxies in the axis ratio bin after the stellar-mass weighting, as if the observed galaxies obey the same stellar mass distribution of the TNG50 galaxies. The total weight of all the observed galaxies by definition equals the total number of simulated galaxies.

We then compute the $p$-value for the test statistic $\chi^{2}$ as $1-$scipy.stats.chi2.cdf($\chi^{2}$, dof=$N_{{\rm bins},q}-1$). In addition, we translate the $p$-value to $n\sigma$ significance as if the test statistic is a Gaussian variable by scipy.stats.norm.ppf($1-p/2$), in order to help comparison with previous studies. The statistics for all three bands are summarized in Table 1. We also include the measurements from H22 for comparison, even if their definition of $\chi^{2}$ is different from ours and the upper limit of our stellar mass range is slightly smaller by 0.15 dex. We find that using their definition of $\chi^{2}$ results in a minor difference of at most $\sim 10\%$ reduction in our $\chi^{2}$ values. The small fraction of galaxies in the stellar mass range of $11.5\leq\log(M_{\star}/M_{\sun})\leq 11.65$ has negligible impacts on the statistics. To sum up, the stark difference between their extreme $\chi^{2}$ ($\gtrsim 200$) and ours cannot be accounted for by the difference in the implementation details.

Figure 2 shows the axis ratio distributions for the TNG50 galaxies (red solid) and the HSC-SSP galaxies after stellar-mass weighting (blue dashed). We find that the stellar-mass weighting has minor impacts on the axis ratio distribution of the HSC-SSP galaxies, because they have similar stellar mass distribution to that of the TNG50 galaxies (see right panel of Figure 1). The excellent agreement of the axis ratio distributions between the observed and simulated galaxies contradicts the findings of H22 that TNG50 produces too few thin galaxies.

In addition, we adopt from H22 the parameter that quantifies the fraction of galaxies with axis ratios smaller than 0.4, $f_{q\leq 0.4}=N_{q\leq 0.4}/N_{\mathrm{tot}}$, where $N_{q\leq 0.4}$ is the number of galaxies with $q\leq 0.4$ and $N_{\mathrm{tot}}$ is the total. For the observed HSC-SSP galaxies, the number counts have been adjusted in the aforementioned stellar-mass weighting, and therefore their total number equals that of the TNG50 galaxies. The parameter is not a direct measure of the fraction of intrinsically thin galaxies but is a function of the distribution of intrinsic thickness, in a way that more intrinsically thin galaxies yield larger $f_{q\leq 0.4}$. The uncertainty of $f_{q\leq 0.4}$ for the TNG50 galaxies is calculated as the Gaussian approximation for binomial proportion $\sqrt{f_{q\leq 0.4}(1-f_{q\leq 0.4})/N_{\mathrm{tot}}}$. For the stellar-mass weighted HSC-SSP sample, because the stellar-mass weighting has minor impacts on the $q$ distribution, we can approximate the uncertainty of $f_{q\leq 0.4}$ using the same formula as for the TNG50 galaxies.

Figure 3 shows the measurements of $f_{q\leq 0.4}$ for various samples of observed and simulated galaxies, after applying the stellar mass weighting to the observational samples. As expected, the $f_{q\leq 0.4}$ of the TNG50 galaxies ($\sim 24\%$) aligns closely with those of the HSC-SSP galaxies ($\sim 25\%$). Moreover, the $f_{q\leq 0.4}$ of the HSC-SSP galaxies ($\sim 25\%$) more or less agrees with the observational results collected by H22 from other imaging surveys ($\sim 22\%$ for SDSS and $\sim 24\%$ for GAMA). Interestingly, they are all significantly larger than the $f_{q\leq 0.4}$ reported for simulated galaxies ($\sim 10\%$) in H22, which implies that the disagreement they found was rooted in their analysis of the simulations. We will discuss more about this in Section 5.1.

Building on the remarkable consistency in the overall axis ratio distributions between the observations and TNG50, this section will delve deeper into exploring the dependence of axis ratios on stellar mass over a broader stellar mass range. The $q\textendash\log M_{\star}$ diagram shown in Figure 4 juxtaposes the galaxies in HSC-SSP and TNG50 with stellar masses of $9.0\leq\log(M_{\star}/M_{\sun})\leq 11.2$. We truncate the sample at $\log(M_{\star}/M_{\sun})=11.2$ to ensure enough galaxies at the highest stellar mass for decent sampling of their axis ratio distribution. To simulate the uncertainties in real-world stellar mass estimates, we also randomly perturb the stellar masses of the TNG50 galaxies as suggested by Rodriguez-Gomez et al. (2019) and show the results in the rightmost column of Figure 4. To illustrate the dependence of the axis ratio distributions on stellar mass, we show the 1st, 16th, 50th, 84th, and 99th percentiles of the axis ratios in each stellar mass bin of $0.25$ dex width. The percentiles of the HSC-SSP and TNG50 galaxies show broad agreement with each other. However, we note that there are noticeable differences at lower stellar masses ($\log(M_{\star}/M_{\sun})\lesssim 9.5$) and in redder bands ($y$ band). We perform similar comparisons as those in Section 4.1 but in much smaller stellar mass bins and find that the difference between the axis ratio distributions of the HSC-SSP galaxies and the TNG50 galaxies are generally small ($\leq 2\sigma$) except for the $y$-band axis ratios in the stellar mass range $\log(M_{\star}/M_{\sun})<9.5$, which surges to $>8\sigma$ (Table 2). Moreover, the statistical significance of the differences at low stellar mass increases monotonically toward longer wavelengths that better probe stellar structures of galaxies, although the increase from $g$ to $i$ band is marginal compared to that from $i$ to $y$ band. The differences in the $y$-band axis ratios are still significant even after smoothing the trend with perturbing the TNG50 stellar masses.

Overall, the $q_{0}$⁸⁸8We use the 1st percentile of the $q$ distribution as a robust approximation of the minimum axis ratio. shows a subtle (except for in the $y$ band) U-shape dependence on stellar mass for the TNG50 galaxies, which is reminiscent of the one reported by Sánchez-Janssen et al. (2010). But, the vertex masses of $\log(M_{\star}/M_{\sun})\approx 10.4$ is different from $\log(M_{\star}/M_{\sun})\approx 9.3$ of Sánchez-Janssen et al. (2010). We also find that the $q_{0}$ of the HSC-SSP galaxies has no discernible dependence on stellar masses in all three bands, and therefore there is no local minimum of $q_{0}$ at $\log(M_{\star}/M_{\sun})\approx 9$, even with the extended stellar mass range. We suspect that the difference in the approaches to measure axis ratio may be the cause, because Sánchez-Janssen et al. (2010) used axis ratios measured at a constant surface brightness of 25 mag arcsec^-2. However, even if we replace our $q$ measurements with those measured at the same surface brightness level as Sánchez-Janssen et al. (2010), we are still not able to reproduce the U-shape trend in the HSC-SSP sample. The difference in sample selection may also contribute to the difference in results, as the sample of Sánchez-Janssen et al. (2010) is closer and spans a smaller distance range (25–80 Mpc) than ours. The origin of the difference between our results and those of Sánchez-Janssen et al. (2010) is unclear, and further investigation is beyond the scope of this study.

In Section 4.1, our self-consistent comparisons of HSC-SSP and TNG50 galaxies revealed remarkable agreement between their shapes. Our results contrast strongly with those of H22, which reported statistically significant tensions between observed and simulated galaxies—citing a deficit of thin simulated galaxies as a key characteristic of this tension. In this section, we will discuss plausible causes for the discrepancies between our results and theirs.

Most obviously, their comparison used heterogeneous observational measurements of galaxy axis ratios collected from literature and applied a qualitatively different method to measure axis ratios from the simulations. On the contrary, we applied the same method to data of the same format and quality (real and synthetic HSC-SSP images) for the observed sample and the TNG50 galaxies, thanks to the forward modeling of the simulations by Bottrell et al. (2024). The approach incorporates various complications that impact the measurements of galaxy shapes, such as stellar population gradients, dust attenuation and scattering, and instrumental effects (PSF smearing and noise), which would otherwise be ignored in measurements based on idealized stellar mass distributions. We note that the observational measurements collected by H22 agree with both the real HSC-SSP measurements and the mock HSC-SSP measurements of the TNG50 galaxies obtained in this study (Figure 3), albeit the adopted techniques are different. Only the axis ratios of the simulated galaxies obtained by H22 consistently stand out as outliers. This leads us to wonder why the method that H22 applied to the simulations yielded such discrepant results.

H22 first extracted the 3D axis lengths of the simulated galaxies defined as $\sqrt{\sum_{j}m_{j}r_{j,i}^{2}}/\sqrt{\sum_{j}m_{j}}$ from supplementary catalogs of the corresponding simulations (e.g., Genel et al., 2015, for TNG), where $m_{j}$ is the stellar mass of the $j$th particle and $r_{j,i}$ is its $i$th axis coordinate. They then projected an ellipsoid extended by such axis lengths onto a sky plane and built a distribution of projected axis ratios under the assumption of random lines of sight. This approach may not be appropriate for deriving comparable axis ratios to those obtained from the images. Firstly, the assumption of a single ellipsoid for their shape is oversimplified for disk galaxies. It blends all structural components together, therefore the spheroidal components, such as stellar bulges and halos, inevitably contribute to the 3D axis ratio calculations and then to the projected axis ratios at all inclinations. The fact that the 3D axis ratios are defined as the second moment of stellar mass distribution may exaggerate the contamination from the spheroidal components when compared with galaxy shapes measured from optical images, because (1) stellar bulges and halos are often fainter than disks at fixed stellar mass density due to their higher optical mass-to-light ratios and (2) they contribute substantially to the intrinsic flatness $q_{0}$—especially the stellar halo particles far away from the disk plane. In the EAGLE simulation, de Graaff et al. (2022) found that galaxy axis ratio distributions measured from projected stellar mass maps are consistent with those from their $r$-band images, which suggests that the former may not be the primary cause and the latter plays a more dominant role.

Our method of fitting an isophote in the outskirt of the projected galaxy images strives to avoid the center. Also, it does not reach the surface brightness level of stellar halos (Merritt et al., 2016). This approach provides reliable measurements of the projected shape of disks, except for edge-on galaxies where bulges inevitably come into play. It makes sense that the axis ratios obtained via parametric fitting (i.e., fitting 2D parametric surface brightness models to images) are broadly consistent with ours in terms of $f_{q\leq 0.4}$ (Figure 3), because both our method and parametric fitting are essentially isophote-based methods. Still, differences in the details remain. Parametric fitting assumes model isophotes to be concentric and to have the same shape, so the best-fit axis ratio is the average axis ratio of the overall flux distribution and would be biased by central components that have a different shape such as bulges. Our method, in contrast, is more robust against such bias.

In summary, the significant disagreement in galaxy axis ratios between the cosmological simulations and the observations reported by H22 can be attributed to the fundamentally different methods applied to obtain the measurements from both sides. Crafting synthetic images that are compatible with observational data enables the application of consistent methods to extract the same observables for comparison with observations. Our study proves such an exercise pivotal by presenting robust and distinct results from those obtained via a different approach.

In Section 4.2, we noted that the TNG50 sample has axis ratio distributions skewed toward larger $q$ compared with those of the HSC-SSP galaxies in the low-mass regime, while the differences are much less significant at higher stellar masses. This implies that the low-mass galaxies in TNG50 are intrinsically rounder/thicker than their observed counterparts. Spurious heating, the gravitational scattering between stellar and dark matter particles due to finite resolution, is known to increase stellar velocity dispersions and the 3D thicknesses of galaxies in cosmological simulations (Ludlow et al., 2021; Wilkinson et al., 2023). Spurious heating should have the most profound effects on low-mass galaxies comprising smaller numbers of stellar and dark matter particles (Ludlow et al., 2023). This can explain the rounder shapes of the low-mass TNG50 galaxies than the observations.

It is worth noting that the generation of the TNG50 synthetic images involves post-processing of the dust content in the simulated galaxies. Because TNG does not explicitly track dust, an empirical model motivated by dust-to-gas ratios in nearby galaxies (Rémy-Ruyer et al., 2014) was used to infer the dust mass from the gas mass in each gas cell based on its metallicity and temperature (Popping et al., 2022; Bottrell et al., 2024). Rodriguez-Gomez et al. (2019) showed that whether modeling dust radiative transfer or not modeling dust at all has a negligible effect on observed median ellipticity derived from synthetic images of simulated galaxies (their Figure 7). Therefore, we argue that the choice of dust post-processing model cannot have such a significant effect to account for the difference between the observed and simulated median $q$ of the low-mass galaxies (Figure 4).

In addition, as highlighted by D. Xu et al. (2024, in preparation), axis ratio measurements for edge-on galaxies are affected by kinematically hot structures such as stellar bulges and other spheroidal components. Consequently, variations in structural compositions of galaxies may play an important role in shaping the dependence of axis ratios on stellar masses. The U-shape dependence of $q_{0}$ on stellar mass is reminiscent of that of S/T on stellar mass in TNG, first reported by Tacchella et al. (2019) in TNG100 and reproduced in Figure 5 for TNG50. A similar trend was also found in the EAGLE simulation (Clauwens et al., 2018). The median S/T reaches its minimum at $10<\log(M_{\star}/M_{\sun})<10.5$ in both TNG and EAGLE. The coincidence of the vertex stellar masses of the $q_{0}$–$\log M_{\star}$ and S/T–$\log M_{\star}$ relations and the stellar mass where early-type galaxies start to dominate in TNG (Huertas-Company et al., 2019) leads us to investigate whether the stellar mass dependence of galaxy structural compositions also contributes to the curvature in the $q_{0}$–$\log M_{\star}$ relation of TNG50.

In Figure 6, we compare the $C_{82}$–$\log M_{\star}$ relation of the observed galaxies and the TNG50 galaxies. The $C_{82}$ is commonly recognized as an economical substitute for the Sérsic (1968) index and bulge-to-total ratio (e.g., Gadotti, 2009). We find broad agreement between the two as well as some interesting differences. The median relations of the TNG50 galaxies exhibit a minimum at a similar stellar mass as the aforementioned U-shape relations, while the median relations of the HSC-SSP galaxies are monotonic. Because it is unclear how the spheroidal components quantitatively contribute to the $q_{0}$ measurements, the absolute difference of $C_{82}$ at fixed stellar masses between the observations and TNG50 is irrelevant for interpreting the shape of the $q_{0}$–$\log M_{\star}$ relations. What matters is the relative difference of the trends with stellar mass between them.

Toward lower stellar masses, the TNG50 relations show a less significant decrease in $C_{82}$ compared to the observational trends in all three bands. Toward higher stellar masses, the TNG50 galaxies exhibit a more prominent increase in $C_{82}$ than the observed ones in all three bands. These differences persist even after smoothing the trends by randomly perturbing the stellar masses of the simulated galaxies. This is also true for the comparison of concentrations and bulge statistics between TNG100 galaxies and real galaxies observed in the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; (Chambers et al., 2016) $i$-band images at $z\approx 0.05$ (Rodriguez-Gomez et al., 2019). Assuming that more prominent spheroidal components at a given stellar mass lead to larger $q_{0}$ estimates of galaxies in the same sample, the difference in the $q_{0}$–$\log M_{\star}$ relations is compatible with the more curved TNG50 $q_{0}$–$\log M_{\star}$ relations than the HSC-SSP ones. It also explains why the U-shape is most noticeable in the reddest $y$ band, where the spheroidal components should be most prominent in flux (as seen in the elevation of $C_{82}$ from the $g$ to $y$ band; Figure 6), because they are mostly quenched structures in massive galaxies. However, we argue that this cannot be the dominant factor for the low-mass thicker galaxies in TNG50, otherwise we should observe a similar degree of mismatch between the observations and TNG50 at the high-mass end and should expect to observe thicker low-mass TNG50 galaxies in the $g$ and $i$ bands as well.

In conclusion, we suggest that the numerical effects of spurious heating are the main cause for the low-mass thicker galaxies in TNG50. The larger concentration of the TNG50 galaxies may contribute to the difference in axis ratios, but we argue that this should only have a minor effect.

The difference in concentration index between the observations and TNG50 is intriguing. As mentioned in Section 5.2, our results are consistent with the trends reported by Rodriguez-Gomez et al. (2019), despite their different definition of total fluxes from ours. However, Guzmán-Ortega et al. (2023) reported different results that TNG50 galaxies have systematically larger average concentration than galaxies observed in the Kilo-Degree Survey (KiDS; de Jong et al., 2013) in the range of $8.5\leq\log(M_{\star}/M_{\sun})\leq 11$ in a similar redshift range of $z<0.05$. Both studies performed similar research using synthetic images and reported broad agreement between the observations and the simulations as our work, and yet high-order differences persist. Difference in image qualities such as survey depth and resolution may contribute to the high-order difference in the trends.

The difference between observations and simulations at the low-mass end can be attributed to poorer numerical resolution. In addition, Rodriguez-Gomez et al. (2019) showed that different dust models yielded $\sim 0.1$ dex variations in the mock observed $i$-band concentration index of simulated galaxies at $\log(M_{\star}/M_{\sun})\approx 9.5$ (their Figure 6). Understanding the difference at the high-mass end requires more dedicated effort. Dust modeling in generating the synthetic images has a negligible effect since massive galaxies are mostly gas-poor (Rodriguez-Gomez et al., 2019). Guzmán-Ortega et al. (2023) found slightly higher merger rates in TNG50 than in KiDS galaxies using both random forest and asymmetry-based classifications, which seems a plausible explanation for the higher concentrations in TNG50 galaxies. However, they also showed that merger and non-merger TNG50 galaxies have similar $r$-band concentrations. Moreover, classification of mergers suffers from incompleteness and impurity issues (e.g., Bickley et al., 2021; Bottrell et al., 2022). The uncertainties in stellar mass estimates of the observed galaxies could also contribute to the difference in the $C_{82}$–$\log M_{\star}$ relations (Rodriguez-Gomez et al., 2019). We strive to emulate such uncertainties by randomly perturbing the stellar masses of the simulated galaxies and find that the difference between TNG50 and the observations remains significant despite the smoother trends in TNG50 (Figure 6). However, we caution that this approach does not account for the inclination dependent uncertainties in stellar mass estimates.