The High Fraction of Thin Disk Galaxies Continues to Challenge ΛCDM Cosmology

In order to compare simulations with observations, we have to determine q_sky for a galaxy with λ_i , where i = 1–3 (Section 3). Only the ratios of the λ_i are relevant for our analysis, but it will be helpful to think of them as actual lengths. We approximate that the galaxy is an ellipsoid and find what this looks like when viewed by a distant observer. We work in Cartesian coordinates in a reference frame centered on the galaxy and aligned with the eigenvectors of its inertia tensor. Thus, the "edge" of the galaxy is given by

$$\begin{eqnarray}&&\displaystyle \frac{{x}^{2}}{{{\lambda }_{1}}^{2}}\,+\,\displaystyle \frac{{y}^{2}}{{{\lambda }_{2}}^{2}}\,+\,\displaystyle \frac{{z}^{2}}{{{\lambda }_{3}}^{2}}\,=\,1.\end{eqnarray} \tag{ 1 }$$

Suppose a distant observer is located toward the direction $\widehat{{\boldsymbol{n}}}$, where we use hats to denote unit vectors. Our approach is to find the extent of the image along the direction ${\widehat{{\boldsymbol{n}}}}_{3}$, which lies entirely within the sky plane (i.e., $\widehat{{\boldsymbol{n}}}\cdot {\widehat{{\boldsymbol{n}}}}_{3}=0$). Our main goal is to find the position vector r of the point corresponding to the edge of the galaxy image along the direction ${\widehat{{\boldsymbol{n}}}}_{3}$. For this purpose, it will be useful to define the unit vector within the sky plane in the orthogonal direction, which we call ${\widehat{{\boldsymbol{n}}}}_{2}$ (orthogonal to both $\widehat{{\boldsymbol{n}}}$ and ${\widehat{{\boldsymbol{n}}}}_{3}$).

We know that ${\boldsymbol{r}}\cdot {\widehat{{\boldsymbol{n}}}}_{2}=0$, or else the point would not appear to be in the direction ${\widehat{{\boldsymbol{n}}}}_{3}$ from the galaxy's center. To get an additional constraint, we note that the tangent plane to the galaxy at the point r must not contain ${\widehat{{\boldsymbol{n}}}}_{3}$, or else it would be possible to move along the galaxy's boundary and reach a larger apparent separation along ${\widehat{{\boldsymbol{n}}}}_{3}$. The plane normal must therefore be $\pm {\widehat{{\boldsymbol{n}}}}_{3}$, but for our analysis, it is sufficient to know that the plane normal is orthogonal to $\widehat{{\boldsymbol{n}}}$.

These two constraints are sufficient to determine the direction of r . Its magnitude is found through Equation (1). Thus, the extent of the image along ${\widehat{{\boldsymbol{n}}}}_{3}$ is d, which we find through the following procedure involving the intermediate vectors q and v :

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{i}\,\equiv \,\displaystyle \frac{{\widehat{{\boldsymbol{n}}}}_{i}}{{{\lambda }_{i}}^{2}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\,{\boldsymbol{v}}\,\equiv \,{\widehat{{\boldsymbol{n}}}}_{2}\times {\boldsymbol{q}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\,d\,=\,\displaystyle \frac{{\boldsymbol{v}}\cdot {\widehat{{\boldsymbol{n}}}}_{3}}{\sqrt{\sum _{i=1}^{3}{\left(\tfrac{{{\boldsymbol{v}}}_{i}}{{\lambda }_{i}}\right)}^{2}}}.\end{eqnarray} \tag{ 4 }$$

We repeat this for a low-resolution grid of ${\widehat{{\boldsymbol{n}}}}_{3}$, whose direction we parameterize using the so-called position angle. Starting from the ${\widehat{{\boldsymbol{n}}}}_{3}$ in this grid which gives the lowest d, we apply the gradient descent algorithm (Fletcher & Powell 1963) to find the minimum value of d to high precision. We then rotate ${\widehat{{\boldsymbol{n}}}}_{3}$ through a right angle, and start a gradient ascent stage to search for the maximum extent of the image. The ratio of these d values is the sky-projected aspect ratio q_sky ≤ 1, which forms the heart of our analysis.

To build up the q_sky distribution, we repeat this procedure for a 2D grid of viewing angles $\widehat{{\boldsymbol{n}}}$, with each result weighted according to the solid angle it represents. The observer is thus assumed to be in a random direction relative to the galaxy.

The Galaxy And Mass Assembly survey (GAMA; Driver et al. 2009, 2011) is a multiwavelength photometric and spectroscopic survey. An overview of the survey regions and their corresponding magnitude limits is given in Table 1 of Baldry et al. (2018). Here, we download the stellar masses, redshifts, and ellipticities by submitting an sql query to the GAMA DR3 database (Baldry et al. 2018). ⁹ In detail, we extract the galfit (Peng et al. 2002, 2010) r-band ellipticity (1 − b/a) from the SérsicPhotometry (v09) - SérsicCatSDSS catalog (Kelvin et al. 2012). We use the stellar mass labeled as "logmstar" from the StellarMasses (v20) catalog (Taylor et al. 2011), which is derived from matched aperture photometry in the r band, missing therewith flux beyond the AUTO aperture. ¹⁰ Consequently, the stellar masses M_*,AUTO within the photometric aperture ("logmstar") are corrected by

$$\begin{eqnarray*}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)={\mathrm{log}}_{10}\left(\displaystyle \frac{{M}_{* ,\mathrm{AUTO}}}{{M}_{\odot }}\right)+{\mathrm{log}}_{10}\left(\displaystyle \frac{{f}_{{\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}}}{{f}_{\mathrm{AUTO}}}\right),\end{eqnarray*}$$

where ${f}_{{\rm{S}}\acute{{\rm{e}}}\mathrm{rsic}}/{f}_{\mathrm{AUTO}}$ is the so-called "fluxscale" parameter, i.e., the linear ratio between the total r-band flux from a Sérsic profile fit cut at 10 effective radii and the r-band AUTO aperture flux (see also, e.g., Taylor et al. 2011; Kelvin et al. 2014; Lange et al. 2016; Vázquez-Mata et al. 2020). The stellar masses are calculated by assuming concordance cosmology with the cosmological parameters being Ω_m,0 = 0.3, Ω_Λ,0 = 0.7, and h=0.7 (Taylor et al. 2011).

For our analysis, we only consider galaxies with a fluxscale correction of 0.5 < fluxscale < 2 and a spectral energy distribution (SED) fit with a posterior predictive P-value PPP > 0, which removes failed SED fits. In addition, we exclude objects with a heliocentric redshift z < 0.005 to remove stars. ¹¹ Also requiring M_* > 10¹⁰ M_⊙ and z < 0.1 yields a final sample of 5304 galaxies that pass the above quality cuts.

The SAMI Galaxy Survey Data Release 3 (DR3) includes 3068 galaxies in the redshift range 0.004 < z < 0.095. The ellipticity distribution of a subsample of the SAMI Galaxy Survey consisting of 826 galaxies is shown in the third panel of Figure 6 in Oh et al. (2020).

In order to increase the sample size, we analyze the aspect ratio distribution of galaxies listed in the input catalogs of the three equatorial regions of the GAMA survey (i.e., G09, G12, and G15; Bryant et al. 2015) and the eight cluster regions (i.e., APMCC0917, A168, A4038, EDCC442, A3880, A2399, A119, and A85; Owers et al. 2017) for the SAMI DR3 Galaxy Survey. For this, we download the stellar masses, redshifts, and projected ellipticities by using an sql/Astronomical Data Query Language (adql) query to the Data Central database servers ¹² of the InputCatGAMADR3 and InputCatClustersDR3 catalogs (for more details, see also Table 1 of Croom et al. 2021). The ellipticities are derived from Sérsic fits in the r band (Bryant et al. 2015; Owers et al. 2019).

Requiring M_* > 10¹⁰ M_⊙ gives a sample of 4252 galaxies, of which 3238 are from GAMA and 1014 are from the cluster regions. As we show in Section 5, the aspect ratio distributions of the GAMA survey and the here described sample disagree only at the 0.013σ confidence level. Although the ellipticity values of the same galaxy listed in the InputCatGAMADR3 and the GAMA survey are slightly different, these catalogs are not fully independent of each other. Because of these reasons and since the GAMA survey contains more galaxies (5304), we do not use the input catalogs of the GAMA and cluster regions for SAMI DR3 in our statistical comparison with simulations (Section 5).

The SDSS is a flux-limited galaxy survey including galaxies with a Petrosian r magnitude brighter than 17.77 (SDSS Collaboration 2000). Here, we download the stellar masses, redshifts, and photometric parameters from its DR16 ¹³ (Ahumada et al. 2020) using an sql query to the SDSS database server. ¹⁴

We selected galaxies listed in the PhotoPrimary catalog, which also contains the aspect ratios of galaxies derived from an exponential and a de Vaucouleurs profile (de Vaucouleurs 1948) for the SDSS r, i, u, z, and g filters. Throughout this analysis, we use the aspect ratio parameters based on the SDSS r-band magnitude. In order to distinguish between spiral and elliptical galaxies, we access the fracDeV parameter, the weight of the de Vaucouleurs component in a linear combination of the exponential and de Vaucouleurs profiles (for a more detailed description, see Section 3.1 of Abazajian et al. 2004). Following Padilla & Strauss (2008), we use the aspect ratio derived from the exponential fit if fracDeV < 0.8 and from the de Vaucouleurs fit if fracDeV ≥ 0.8.

The here used stellar masses from the galSpecExtra table correspond to the median of the estimated logarithmic stellar mass probability density function using model photometry (Kauffmann et al. 2003; Salim et al. 2007). The redshifts are extracted from the SpecObj table. ¹⁵

Selecting galaxies with z < 0.1 and M_* > 10¹⁰ M_⊙ gives a sample of 232,315 galaxies. The so-obtained and here analyzed q_sky distribution is consistent with Padilla & Strauss (2008), as shown in Section 5. They obtained a volume-limited sample by weighting each galaxy in the SDSS DR6 (Adelman-McCarthy et al. 2008) by $1/{V}_{\max }$, where ${V}_{\max }$ is the volume corresponding to the maximum distance at which we can observe a galaxy with its absolute magnitude (see Section 3.1 of Padilla & Strauss 2008). In total, their sample contains 303,290 "spirals" (fracDeV < 0.8) and 282,203 "ellipticals" (fracDeV ≥ 0.8). Here, we extract the $1/{V}_{\max }$ weighted q_sky distributions of the spiral and elliptical samples from their Figure 1 (bottom panels), and combine these in order to obtain the q_sky distribution of the total galaxy sample.

Zoom In Zoom Out Reset image size

As a consistency check, we also investigate the q_sky distribution of galaxies within the Local Volume (LV), a sphere of radius 11 Mpc centered on the Sun. The stellar masses are calculated from K-band luminosities using a mass-to-light ratio of 0.6 (e.g., McGaugh & Schombert 2014). We select 52 galaxies with $10.0\lt {\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\leqslant 11.65$ from the Updated Nearby Galaxy Catalog (Karachentsev et al. 2013), a renewed version of The Catalog of Neighboring Galaxies ¹⁶ (Karachentsev et al. 2004). These galaxies cover a wide range of q_sky (0.13–0.94), as shown in the right panel of Figure 1. Due to the small sample size and high resulting Poisson uncertainties, we do not calculate its tension with the simulated q_sky distributions.

These are nearby galaxies with well-established morphological types, allowing direct well-resolved images to inform us of which galaxies are thin disks. Therefore, we show the distribution of the morphological T-types of our LV galaxy sample in the left panel of Figure 1. The T-types of these galaxies are listed in the Updated Nearby Galaxy Catalog and discussed and analyzed in detail in Karachentsev et al. (2018). Of the 52 galaxies, 50% are early-type spirals (Sa, Sab, Sb, Sbc, Sc; morphological type T = 1–5), 23% are late-type spirals (Scd, Sd, Sdm, Sm; T = 6–8), and the remaining 27% are E, S0, dSph, or S0a galaxies (T < 1). Thus, the morphological classification scheme used by Karachentsev et al. (2013, 2018) is slightly different to the de Vaucouleurs system.

It is not in general possible to disentangle the contributions to q_sky from inclination and intrinsic thickness. The inclination i between disk and sky planes listed in the Updated Nearby Galaxy Catalog is derived with Equation (5) of Karachentsev et al. (2013):

$$\begin{eqnarray}&&{\sin }^{2}i\,=\,\displaystyle \frac{1-{q}_{\mathrm{sky}}^{2}}{1-{\widetilde{q}}_{\mathrm{int}}^{2}},\end{eqnarray} \tag{ 5 }$$

where the assumed intrinsic aspect ratio ${\widetilde{q}}_{\mathrm{int}}$ is assigned a value depending on the morphological T-type (T) of the considered galaxy as given by their Equation (6): ¹⁷

$$\begin{eqnarray}{\mathrm{log}}_{10}{\widetilde{q}}_{\mathrm{int}}\,=\,\left\{\begin{array}{ll}-0.43\,-0.053T & \quad \left(T\leqslant 8\right),\\ -0.38 & \left(T=9,10\right),\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where according to Karachentsev et al. (2018) T = 9 and T = 10 refer to Im/BCD and Ir galaxies, respectively. Thus, the so-obtained intrinsic aspect ratios are a best guess based on the observed morphology. As an example, the M31 galaxy with T = 3 has with this approach ${\widetilde{q}}_{\mathrm{int}}=0.26$ (q_sky = 0.33).

Applying Equation (6) to the 52 selected LV galaxies shows that their intrinsic aspect ratio distribution covers the range ${\widetilde{q}}_{\mathrm{int}}=0.07\mbox{--}0.55$ with a global peak at ${\widetilde{q}}_{\mathrm{int}}\approx 0.23$, as presented in the right panel of Figure 1. Converting q_sky to q_int depends on the relation between ${\widetilde{q}}_{\mathrm{int}}$ and the morphological T-type, so the intrinsic aspect ratio distribution of our LV sample is mainly considered for illustrative purposes. However, it helps to demonstrate that the LV galaxies seem to be mostly thin disks. According to this, if the LV were to be representative of the universe, then 81% ± 13% (42/52) of all galaxies heavier than M_* = 10¹⁰ M_⊙ are thin disk galaxies with ${\widetilde{q}}_{\mathrm{int}}\lt 0.4$. For completeness, we quantify the tension between the here obtained intrinsic aspect ratio distribution of LV galaxies and the ΛCDM models in Section 6.1.

In order to avoid any bias on the shapes of galaxies introduced by model assumptions when converting q_sky to q_int, we only consider the former when statistically comparing observations to simulations in our main analysis. As explained in Section 3.1, this requires us to generate the q_sky distribution given the intrinsic shapes of simulated galaxies. We can then directly compare the resulting distribution with the observed q_sky distribution (Section 6.2).

Figure 3 shows the q_int distribution of galaxies in central and noncentral subhalos with M_* > 10¹⁰ M_⊙ in the TNG50-1, TNG100-1, Illustris-1, EAGLE50, and EAGLE100 simulations. These all show a unimodal distribution with different peak positions. We show later that increasing the resolution broadens the peak and shifts it to smaller q_int (Section 7.1). Galaxies in the EAGLE and TNG50-1 runs are typically thinner than in the Illustris-1 and TNG100-1 runs. The TNG50-1 run contains 133 galaxies with q_int < 0.3 and M_* > 10¹⁰ M_⊙, with the thinnest galaxy of this sample having q_int ≈ 0.19. Clearly, the resolution and the adopted temperature floor of about 10⁴ K (10 kK; see Trayford et al. 2017 for EAGLE and Nelson et al. 2019 for TNG) allow for the formation of massive thin disk galaxies. This is consistent with the work of Sellwood et al. (2019), who managed to maintain a thin disk in a hydrodynamical CDM-based simulation of M33 with a slightly higher temperature floor of 12 kK.

Zoom In Zoom Out Reset image size

The EAGLE runs have a very similar q_int distribution—the peak position is similar to TNG50, and EAGLE forms galaxies as thin as q_int = 0.18 (the thinnest galaxies in each run are q_int = 0.20 for EAGLE25, q_int = 0.19 for EAGLE50, and q_int = 0.18 for EAGLE100). This is in agreement with TNG50-1 (right panel of Figure 3). By using different resolution realizations of the TNG50 and EAGLE simulations, we test the numerical convergence of the aspect ratio distribution in Section 7.1. While there is no evidence that the TNG50 simulation has numerically converged, the aspect ratio distributions of different EAGLE runs closely agree with each other, possibly because they have the same resolution. The similarity between EAGLE runs is consistent with Lagos et al. (2018), who showed that higher-resolution realizations of the EAGLE project yield galaxies with a similar ellipticity distribution (see their Figures A2 and A3). Moreover, the q_int distributions of EAGLE and TNG50-1 are very similar despite differences in resolution and other details of the models.

In the LV, 81% ± 13% of galaxies with M_* > 10¹⁰ M_⊙ have ${\widetilde{q}}_{\mathrm{int}}\lt 0.4$ (Section 4.4). In contrast, the fractions of galaxies with such masses that have q_int < 0.4 in the highest-resolution simulations analyzed here (TNG50-1, EAGLE25) are only 39% ± 2% and 46% ± 8%, respectively. The stated Poisson uncertainties are estimated as $\sigma =\sqrt{{N}_{\lt 0.4}+1}/{N}_{\mathrm{tot}}$, where N_<0.4 and N_tot are the number of galaxies with q_int < 0.4 and the total number of galaxies, respectively. The overall intrinsic aspect ratio distribution of the LV galaxies (blue dashed line in the left panel of Figure 1) is in 5.42σ tension with the TNG50-1 run. We emphasize again that the intrinsic aspect ratios of LV galaxies are based on the morphological T-type (Equation (6)). Therefore, we provide in the following section another test of the ΛCDM framework in which we compare the sky-projected aspect ratio distributions of the ΛCDM simulations with the GAMA survey and SDSS.

In order to compare these simulation results with observations from the GAMA survey and SDSS (Section 4), the simulated galaxies are projected onto the sky from a grid of viewing angles to find the distribution of q_sky (Section 3.1). The resulting q_sky distribution of each simulation is compared with the GAMA survey and SDSS in Figure 4. This reveals that the simulations significantly underproduce the fraction of galaxies with a low q_sky. In particular, only about 10% of simulated galaxies with $10.0\lt {\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\leqslant 11.65$ have q_sky < 0.4, while these galaxies make up about 22%–25% of locally observed galaxies (Table 2).

Zoom In Zoom Out Reset image size

The tension between the simulated and observed q_sky distributions is quantified using a standard χ² statistic (Section 5). All simulations significantly disagree with local observations. The smallest tension is 12.52σ, which arises from comparing the GAMA survey with TNG50-1. The smallest tension for the SDSS is 14.82σ, corresponding to a comparison with EAGLE50 (we do not consider results from the EAGLE25 run as it only has 81 galaxies within the analyzed M_* range, but see Section 7.1). The total χ² values and the corresponding levels of tension are listed in Table 3 for different simulations.

The fact that zoomed-in ΛCDM simulations can produce thin disk galaxies (e.g., Wetzel et al. 2016) suggests that further increasing the numerical resolution may solve the discrepancy reported in Section 6.2. For example, Ludlow et al. (2021) recently argued that the coarse-grained implementation of dark matter halos causes an artificial heating of the stellar particles, which increases the vertical velocity dispersion of the galaxies. This yields thicker disks than better resolved galaxies, which could be the reason for the here reported tension. According to their Table 3, the vertical velocity dispersion σ_z of TNG50-1 galaxies embedded in a dark matter halo of M₂₀₀ < 7.9 × 10¹¹ M_⊙ is numerically increased by Δσ_z > 0.1 V₂₀₀, where V₂₀₀ is the virial velocity of the halo, and M₂₀₀ is its virial mass.

In this section, we analyze the effect of numerical resolution on the distribution of galaxy shapes using the TNG50 simulation. The flagship of the TNG project is called TNG50-1, which has the highest resolution among the here analyzed simulation runs. In addition, the TNG50 simulation suite includes the runs TNG50-2, TNG50-3, and TNG50-4, where a higher suffixed number indicates a lower-resolution realization as summarized in Table 1. The particle mass differs by a factor of eight between any TNG50 run and the next higher-resolution realization.

The left panel of Figure 5 shows the intrinsic aspect ratio distribution for different resolution realizations of the TNG50 simulation. The formation of thin galaxies and the peak position of the q_int distribution strongly depend on the resolution: the mode shifts from q_int = 0.88 for TNG50-4 to q_int = 0.33 for TNG50-1.

Zoom In Zoom Out Reset image size

The convergence of different simulation runs is statistically quantified in Table 4. The q_int distributions of the TNG50-1 and TNG50-2 runs disagree with each other at 7.60σ. Although this is still a high tension, it is much lower compared to that between TNG50-2(-3) and TNG50-3(-4), which differ at the 13.75σ (14.17σ) confidence level. Based on the q_sky distribution, the TNG50-1 and TNG50-2 runs disagree at only 0.83σ.

Table 4. Testing the Numerical Convergence of Different Simulation Runs (Column 1) for Subhalos with M_* > 10¹⁰ M_⊙ by Showing the Total χ² (Equation (10)) between Their Distributions of q_int (Column 2) and q_sky (Column 3), along with the Corresponding Level of Tension (Bracketed Numbers)

Simulations	${\chi }_{\mathrm{int}}^{2}$ between Them (Tension)	${\chi }_{\mathrm{sky}}^{2}$ between Them (Tension)
TNG50-3 versus TNG50-4	268.26 (14.17σ)	91.80 (6.69σ)
TNG50-2 versus TNG50-3	255.73 (13.75σ)	105.10 (7.45σ)
TNG50-1 versus TNG50-2	107.85 (7.60σ)	20.82 (0.83σ)
TNG50-1 versus EAGLE100	87.83 (6.35σ)	18.35 (0.58σ)
TNG50-1 versus EAGLE50	56.78 (4.24σ)	6.04 (0.0014σ)
TNG50-1 versus EAGLE25	19.70 (0.71σ)	2.98 (4.83 × 10−6 σ)
EAGLE100 versus EAGLE25	19.23 (0.66σ)	0.50 (2.80 × 10−13 σ)
EAGLE50 versus EAGLE25	12.87 (0.15σ)	0.95 (1.26 × 10−10 σ)

Note. We use 20 bins in q_int and q_sky, giving 20 degrees of freedom. EAGLE50 and EAGLE100 use the same resolution, so the similarity between their results (χ² not shown) is not a strong test of numerical convergence. Note that the higher-resolution realization EAGLE25 only has 81 galaxies with M_* > 10¹⁰ M_⊙, which reduces the total χ² because of the higher Poisson uncertainties. The intrinsic aspect ratio distributions underlying these comparisons are plotted in Figure 5.

Download table as:  ASCIITypeset image

The right panel of Figure 5 compares the q_int distribution of TNG50-1 with EAGLE100, EAGLE50, and the higher-resolution realization EAGLE25 (Recal-L0025N0752). The similarity of the EAGLE results (Table 4) implies that the EAGLE simulations have numerically converged, so the deficit of intrinsically thin galaxies cannot be explained by resolution effects—as also concluded by Lagos et al. (2018) and apparent in their Figures A2 and A3. Although the q_int distribution of the EAGLE25 run is very similar to EAGLE50 or EAGLE100, the tension between the sky-projected aspect ratio distribution of EAGLE25 and GAMA DR3 (SDSS) is only 1.44σ (0.37σ). This is because EAGLE25 only has 81 galaxies in the stellar mass range $10.0\lt {\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\leqslant 11.65$, which results in large Poisson uncertainties that decrease therewith the total χ².

Importantly, the q_sky distributions of galaxies in TNG50-1 and EAGLE50 (EAGLE100) are consistent with each other at 1.4 × 10⁻³ σ (0.58σ) confidence (Table 4). Thus, we demonstrate for the first time that ΛCDM simulations with SPH-based (EAGLE) and adaptive grid (Illustris/TNG) methods produce the same sky-projected aspect ratio distribution. The baryonic algorithms are also quite different between these simulations.

In order to quantify the thickening of simulated galaxies due to limited numerical resolution (e.g., Ludlow et al. 2021), we apply as an ansatz a parametric correction to the intrinsic shape distribution. The short axis λ₁ of each galaxy is scaled down by the factor x such that its q_int is corrected by

$$\begin{eqnarray}&&{q}_{\mathrm{int},\ \mathrm{corr}}\,=\,{{xq}}_{\mathrm{int}},\quad x\,=\,\alpha +\left(1-\alpha \right){q}_{\mathrm{int}}.\end{eqnarray} \tag{ 12 }$$

α is a variable ranging from 0–1 in steps of Δα = 0.01. Applying this correction to the q_int of simulated galaxies keeps a galaxy with a high q_int round, but makes a thin disk galaxy even thinner. α = 1 does not change the q_int of a galaxy, while α = 0 yields a corrected ${q}_{\mathrm{int},\ \mathrm{corr}}={q}_{\mathrm{int}}^{2}$, therewith causing the most substantial thinning of the galaxy population in our parameterization.

The tension between the observed and rescaled simulated q_sky distributions is shown in Figure 6 for different resolution realizations of the TNG50 simulation. As expected from our previous analysis, the tension systematically decreases with higher resolution. The tension between TNG50-1 and the GAMA survey (SDSS) reaches the 5σ confidence level if galaxies in the TNG50-1 run are corrected by α = 0.668 (α = 0.738). The tension between TNG50-1 and GAMA (SDSS) is minimized for a correction factor of α = 0.32 (α = 0.47), the tension in this case being only 4.6 × 10⁻³ σ (1.2 × 10⁻² σ). This demonstrates that our simple parametric correction with a single parameter α (Equation (12)) can make the simulations agree very well with observations.

Zoom In Zoom Out Reset image size

We can use this procedure to quantify how improving the resolution affects the simulated aspect ratio distribution, which then allows us to quantify if the simulations are numerically converged. To do this, we first find which α value must be applied to a TNG50 run to minimize the tension with the uncorrected q_int distribution of the next higher-resolution realization with an eight-times-higher mass resolution. We find that the tensions between the q_int distributions of the rescaled TNG50-4 and TNG50-3, rescaled TNG50-3 and TNG50-2, and rescaled TNG50-2 and TNG50-1 runs become minimal for α values of 0.03 (0.85σ), 0.44 (5.9 × 10⁻³ σ), and 0.69 (0.80σ) applied to the TNG50-4, TNG50-3, and TNG50-2 runs, respectively. The true q_int distribution of each TNG50 run and the so-corrected distribution of the next lower-resolution TNG50 run is compared in the left panel of Figure 7, demonstrating that our parametric correction to, e.g., TNG50-2 galaxy shapes reproduces quite well the actual q_int distribution of TNG50-1. The required α value increases with higher resolution, but we cannot confirm that the TNG50 simulations have numerically converged, as that would be achieved if α = 1.

Zoom In Zoom Out Reset image size

In a second step, we extrapolate the α values by plotting $\mathrm{ln}\left(1-\alpha \right)$ against the logarithmic dark matter particle mass of the corresponding TNG50 run (right panel of Figure 7). The α value required for galaxies in the TNG50-1 run to mimic the result of an eight-times-lower dark matter particle mass simulation is estimated using a parabolic and a linear fit in the ${\mathrm{log}}_{10}\left({m}_{\mathrm{dm}}/{M}_{\odot }\right)$ versus $\mathrm{ln}\left(1-\alpha \right)$ diagram of Figure 7. This predicts α = 0.824 (0.834) for a linear (parabolic) fit. Being more conservative, we apply the linearly extrapolated result that α = 0.824 to the TNG50-1 run, which yields an 8.68σ (8.71σ) tension with GAMA (SDSS). Thus, an eight-times higher-resolution realization of the TNG50-1 run would almost certainly not reduce the here reported tension below the 5σ confidence level.

In the third and final step, we estimate the tension with observations by extrapolating the data further to five more refinement levels than in the TNG50-1 run. Therefore, we extract the α values not only for an eight-times-lower m_dm but also for an 8²×, 8³×, 8⁴×, and 8⁵× lower m_dm than in the TNG50-1 run. These so-obtained α values for the five different resolution levels are listed in Table 5 for the parabolic and linear extrapolations. The α values are then successively applied to estimate the intrinsic aspect ratio of each subhalo in an 8⁵ × higher-resolution realization than TNG50-1. In other words, the aspect ratio distribution for an 8²× higher-resolution run than TNG50-1 is obtained by correcting the aspect ratios of subhalos in the eight-times higher-resolution run by applying Equation (12) with α = 0.915 (parabolic fit) or α = 0.900 (linear fit). This procedure is repeated until we reach an 8⁵× higher refinement level, which on the last step requires applying α = 0.991 (parabolic fit) or α = 0.982 (linear fit) to the q_int distribution of the 8⁴× higher-resolution realization than TNG50-1. Since α is now close to unity, this almost reaches full convergence in q_int. Further improvements to the resolution can be expected to have sub-percent level effects on the intrinsic aspect ratios.

The q_sky distributions for an 8× and 8⁵× higher dark matter mass resolution than the TNG50-1 run derived from a parabolic (linear) extrapolation are presented in the left (right) panel of Figure 8. Raising the resolution of the dark matter mass in TNG50-1 by a factor of 8⁵ yields an estimated tension of 5.58σ (parabolic fit) and 4.36σ (linear fit) with GAMA. In contrast, the discrepancy is significantly alleviated with respect to the older SDSS data set—we estimate a tension of 3.27σ with a parabolic fit and only a 1.61σ tension with a linear fit. The difference between the two observational samples is mainly because the fraction of galaxies with q_sky ≲ 0.3 is higher in the GAMA survey compared to SDSS (Figure 8). The total χ² values and corresponding levels of tension with the GAMA and SDSS surveys for different resolution realizations and extrapolation methods are summarized in Table 5.

Zoom In Zoom Out Reset image size

Since the parabolic fit considers more information than the linear fit and exactly matches the data, we consider it as the nominal case in our main analysis. The linear fit is still shown in order to illustrate the effect of uncertainties in the extrapolation procedure. Therefore, our estimate is that increasing the mass resolution of TNG50-1 by 8⁵× still leaves a significant tension of 5.58σ with GAMA, exceeding the 5σ plausibility threshold. The SDSS data set is then in 3.27σ tension, which still points to an underestimated fraction of very thin galaxies in the ΛCDM simulations (see Figure 8).

Finally, we note that our approach is very conservative with respect to the ΛCDM framework. First of all, the extrapolation to five more refinement levels than TNG50-1 is based only on the four resolution realizations of the TNG50 simulation. For example, the TNG50-4 simulation has m_dm = 2.3 × 10⁸ M_⊙, which is a factor of 8⁸ = 1.7 × 10⁷ more massive than dark matter particles in an 8⁵× higher-resolution realization than the TNG50-1 run. The parabolic and linear fits have been derived from simulations with 4.5 × 10⁵ (TNG50 − 1) ≤ m_dm/M_⊙ ≤ 2.3 × 10⁸ (TNG50 − 4). This introduces additional uncertainties because the shape of the relation between ${\mathrm{log}}_{10}\left({m}_{\mathrm{dm}}/{M}_{\odot }\right)$ and $\mathrm{ln}\left(1-\alpha \right)$ may deviate significantly from the derived polynomial fits at lower m_dm. Second, we performed the extrapolation to very high-resolution realizations without assuming that numerical convergence is reached between two different refinement levels. It is not impossible that numerical convergence (α = 1) is already attained before reaching the 8⁵× higher-resolution level. Consequently, higher-resolution runs of self-consistent ΛCDM simulations are still required to test when numerical convergence is reached and if the tension with observations is then alleviated.

In contrast to our findings, Bottrell et al. (2017b) found a significant deficit of bulge-dominated galaxies based on the photometric bulge/total ratios of subhalos in the Illustris-1 simulation compared to SDSS. In particular, Bottrell et al. (2017a) photometrically derived ${\left(B/T\right)}_{\mathrm{phot}}$ ratios of subhalos with M_* > 10¹⁰ M_⊙ at z = 0 in the Illustris-1 simulation by performing a bulge-disk decomposition of the surface brightness profile with fixed Sérsic indices n_b = 4 for the bulge and n_d = 1 for the disk. By applying the same decomposition analysis to observed SDSS galaxies, they found a significant deficit of bulge-dominated galaxies with M_*/M_⊙ = 10¹⁰–10¹¹ in the Illustris-1 simulation (see Figures 4 and 6 in Bottrell et al. 2017b).

If the Illustris-1 simulation indeed lacks bulge-dominated galaxies, this deficit should also be evident in the q_int parameter (Section 6.1). To test if ${\left(B/T\right)}_{\mathrm{phot}}$ reflects the shapes of simulated galaxies as quantified by q_int, we show in Figure 9 its distribution for four simulated galaxy subsamples of Bottrell et al. (2017a) with different photometric morphologies, i.e., ${\left(B/T\right)}_{\mathrm{phot}}$: 0, (0,0.2], (0.2,0.5], and (0.5,1]. ¹⁸ Throughout this work, we use the ${\left(B/T\right)}_{\mathrm{phot}}$ ratio derived for the r band from camera angle 0. Remarkably, all four galaxy subsamples have a very similar aspect ratio distribution, with the peak between q_int = 0.6–0.8. The distributions have very few galaxies with q_int < 0.5: the proportions for the above-mentioned ${\left(B/T\right)}_{\mathrm{phot}}$ bins are only 3.4%, 7.1%, 3.6%, and 0.26%, respectively, implying that the sample of Bottrell et al. (2017a) is actually bulge-dominated.

Zoom In Zoom Out Reset image size

In order to understand this mismatch between ${\left(B/T\right)}_{\mathrm{phot}}$ and q_int, we begin by presenting color images ¹⁹ generated from the Friends-of-Friends (FoF) halo finder. Subhalos with low $0\lt {\left(B/T\right)}_{\mathrm{phot}}\lt 0.2$ and q_int < 0.4 indeed have spiral structures and a disk (Figure 10). The κ_rot morphological parameter (Equation (2) in Rodriguez-Gomez et al. 2015) indicates that these subhalos are rotation-dominated (κ_rot > 0.5).

Zoom In Zoom Out Reset image size

The mismatch between the morphological parameters becomes evident in Figure 11, which shows the significant fraction of subhalos (13.1%) with low $0\lt {\left(B/T\right)}_{\mathrm{phot}}\lt 0.2$ but high q_int > 0.8 for their stellar component. ²⁰ These are featureless-looking dispersion-dominated (κ_rot < 0.5) galaxies without a pronounced disk or spiral structures. Based on these images and given that the intrinsic aspect ratio is a very robust parameter to quantify the actual shape of a galaxy, we conclude that the 2D parametric surface brightness decomposition applied by Bottrell et al. (2017a) is inadequate if applied to simulated galaxies.

Zoom In Zoom Out Reset image size

In addition to a bulge-disk decomposition, Bottrell et al. (2017a) also applied a pure Sérsic model to the galaxies in subhalos by varying the Sérsic index between 0.5 and 8. Figure 12 shows a significant correlation between ${\left(B/T\right)}_{\mathrm{phot}}$ and the Sérsic index, highlighting the inevitable confusion between an elliptical galaxy with n_s ≈ 1 and a thin exponential disk if only considering the surface brightness profile. Choosing a different camera angle leads to the same results. In combination with Figure 9, this shows that neither n_s nor ${\left(B/T\right)}_{\mathrm{phot}}$ correlates with q_int.

Zoom In Zoom Out Reset image size

This might be due to the resolution of the Illustris-1 simulation not being sufficient to apply a photometric decomposition. Importantly, we have shown in this contribution that the observed and simulated galaxy morphology distributions differ significantly (Table 3). The vast majority of observed galaxies are spirals (e.g., Loveday 1996; Delgado-Serrano et al. 2010), while the ΛCDM simulations form a far-too-large fraction of bulge-dominated galaxies (Figures 3 and 4). Thus, assuming that all galaxies with low ${\left(B/T\right)}_{\mathrm{phot}}$ are intrinsically thin is quite accurate in the real universe where spirals are quite common, but not in a ΛCDM universe where they are rare (Section 6.2). In the Illustris-1 simulation, there are so few disk galaxies that galaxies with low ${\left(B/T\right)}_{\mathrm{phot}}$ are nearly always ellipticals with an exponential-like surface brightness profile (Sérsic index close to 1).

These problems with the ${\left(B/T\right)}_{\mathrm{phot}}$ parameter are avoided by using the sky-projected aspect ratio (Figure 4), which is closely linked to the intrinsic aspect ratio and thus yields a more robust measurement of the galaxy morphology. Since q_sky is observable, it allows for a much more direct test of the model, provided the sample size is sufficient to statistically sample over projection effects. This is true in our case because the TNG50-1, GAMA DR3, and SDSS DR16 samples contain 882, 5304, and 232,128 galaxies, respectively, in the range $10.0\lt {\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\leqslant 11.65$ used for our comparisons.

The ΛCDM theory strictly implies a hierarchical merger-driven build-up of the galaxy population. Galaxy mergers are driven by dynamical friction on the extended dark matter halos of interacting galaxies (e.g., Kroupa 2015). Mergers grow the bulge component and thicken the stellar disk. By studying galaxy–galaxy interactions in N-body/hydrodynamical simulations, Hwang et al. (2021) showed that the disk angular momentum of the late-type galaxy decreases by about 15%–20% after a prograde collision. Thus, galaxies with a quiescent merger history are expected to have lower bulge fractions than galaxies that have undergone a major merger (see also, e.g., Bournaud et al. 2005; D'Onghia et al. 2006). Since most observed galaxies are late types (e.g., Delgado-Serrano et al. 2010) and ≈50% of these have no classical bulge (Graham & Worley 2008; Kormendy et al. 2010), we might expect that galactic mergers are less frequent in the universe than in ΛCDM simulations (Disney et al. 2008; Stewart et al. 2008; Fakhouri et al. 2010; Kroupa 2015; Wu & Kroupa 2015). This may underlie the tension between the observed and simulated sky-projected aspect ratio distributions (Figure 4).

We therefore investigate if galaxies with a quiescent merger history in the ΛCDM framework are typically much thinner. We focus on the TNG50-1 run, as it has the highest resolution (Table 1) and yields the lowest tension in the here analyzed simulations (Table 3). The merger trees of the Illustris and TNG projects can be downloaded from their webpages, ²¹ with a detailed description available in Rodriguez-Gomez et al. (2015). We quantify the merger history of a galaxy by considering the total mass ratio μ ≤ 1 of the two progenitors and the maximum lookback time ${t}_{\max }$ during which we consider mergers. If there are multiple mergers in this timeframe, we use the most major merger, i.e., that with the highest μ. This μ value is used to construct two galaxy samples that differ according to whether the galaxy had at least one major merger with μ ≥ 1/12 within a lookback time of ${t}_{\max }=12\,\mathrm{Gyr}$ (z = 3.7). We restrict ourselves to galaxies with M_* > 10¹⁰ M_⊙ and M_dm/M_* > 1 at z = 0. The constraint on the M_dm/M_* ratio is applied to exclude dark matter–poor galaxies, which cannot be traced back accurately. This excludes 26 galaxies from the original sample. About 89% of all subhalos at z = 0 have undergone at least one major merger, which is broadly consistent with other ΛCDM simulations (Stewart et al.2008; Fakhouri et al. 2010). The remaining 11% of all subhalos had at most only minor merger(s) with μ ≤ 1/12 during the past 12 Gyr.

As expected from zoomed-in simulations of galaxies in underdense environments, intrinsically thin galaxies are more frequent in the sample of galaxies with a quiescent merger history compared to that with an active merger history (see the left panel of Figure 13). In particular, 39% ± 2% of galaxies with an active merger history have q_int < 0.4, but this rises to 50% ± 7% for galaxies with a quiescent merger history. The aspect ratio of the thinnest such galaxy (q_int = 0.22) is similar to that in the sample with an active merger history (q_int = 0.19). Additionally, 9.5% ± 1.1% of galaxies with an active merger history have q_sky < 0.4, but this rises to 12.5% ± 3.7% for the galaxies with a quiescent merger history (right panel of Figure 13).

Zoom In Zoom Out Reset image size

Thus, the lack of intrinsically thin galaxies in ΛCDM could be partly due to major mergers. More likely, it is due to a combination of major and minor mergers, which are unavoidable in ΛCDM, as galaxies grow their mass through mergers. We emphasize that minor mergers are expected to be much less frequent in an alternative framework where galaxies lack extended dark matter halos, as occurs in MOND (Milgrom 1983). In addition, secular processes like disk-halo angular momentum exchange could also drive the formation of significant bars and bulges in ΛCDM (Athanassoula 2002; Sellwood et al. 2019), but perhaps not in the real universe (Banik et al. 2020; Roshan et al. 2021a) where the fraction of bars differs substantially from ΛCDM expectations (Reddish et al. 2021). There is also a highly significant discrepancy between the pattern speeds of bars in observations and in ΛCDM simulations, where the results seem to have converged with respect to the ratio of bar length to corotation radius (Roshan et al. 2021b). The tension is caused by the fact that bars are expected to be slowed down by dynamical friction with the dark matter halo, but observed bars are fast. This is another indication against dynamical friction from massive dark matter halos.

In the previous section, we discussed that the disagreement between the observed and ΛCDM simulated galaxy shapes could be partly due to the frequency of mergers being too high in this framework. Another possibility is that the tension is caused by the feedback description used in the simulations. In particular, Lagos et al. (2018) discussed the link between loss/gain of angular momentum and dry/wet mergers. Using the EAGLE and HYDRANGEA (Bahé et al. 2017; Barnes et al. 2017) hydrodynamical simulations, they showed that dry mergers typically decrease the stellar spin parameter while wet mergers increase it. Therefore, galaxies that further accrete cold gas are able to reform their disks. This process is sensitive to the implemented feedback description.

As shown in Sections 6.1 and 7.1, the TNG50-1 and EAGLE simulations produce very similar aspect ratio distributions despite relying on different sub-grid models. This is a strong indication that the tension is likely not caused by the implemented sub-grid feedback models. Even so, improved models might yet alleviate or resolve the tension. For example, the feedback recipe in both the TNG and EAGLE projects could be too strong to allow the (re)formation of disks. In this case, the shape of the galaxies would be set by the merger rate independently of whether the mergers are dry or wet. However, a weaker feedback description could potentially increase the efficiency of disk formation. We note that strong feedback is required in ΛCDM simulations for them to explain why the Newtonian dynamical mass of a galaxy or galaxy group often greatly exceeds 6.4× its baryonic mass (e.g., Müller et al. 2021), even though ΛCDM needs this to be the cosmic ratio between baryonic and total mass (Planck Collaboration VI 2020). A strong feedback prescription is also needed to solve the missing satellite galaxy problem (e.g., Brooks et al.2013).

The difficulties faced by ΛCDM with regards to the high fraction of thin disk galaxies motivate us to consider MOND one of the main alternative frameworks that is generally considered to perform better on galaxy scales (for reviews, see, e.g., Famaey & McGaugh 2012; Kroupa 2015; Banik & Zhao 2022). For illustrative purposes, we show in the left panel of Figure 14 the present-day q_int distribution of six disk galaxies with 10¹⁰ < M_*/M_⊙ ≤ 9.56 × 10¹⁰ formed in hydrodynamical MOND simulations of collapsing gas clouds conducted by Wittenburg et al. (2020) using the phantom of ramses code (Lüghausen et al. 2015), an adaptation of the N-body and hydrodynamics solver ramses (Teyssier 2002) for MOND gravity (for a user guide, see Nagesh et al. 2021). These non-cosmological simulations have a box size of 960 kpc per side, a minimum stellar mass of M_* ≈ 3 ×10⁴ M_⊙, and depending on the run a maximum grid cell resolution of 117.19 pc, 234.38 pc, or 468.75 pc. The temperature floor for the gas is set to 10 kK (for further information on the initial conditions and numerical parameters of the individual simulations, see Table 1 of Wittenburg et al. 2020). The Milgromian disk galaxies are systematically thinner than in the here analyzed ΛCDM simulations, as evidenced by their ${q}_{\mathrm{int}}=0.18\mbox{--}0.30\,\left(0.31\mbox{--}0.54\right)$ for r < 30 kpc (r < 2 r_0.5,*), with the global peak at ${q}_{\mathrm{int}}\approx 0.2\,\left(0.3\right)$. While these isolated MOND results cannot yet be directly compared with self-consistent cosmological ΛCDM simulations that allow for galaxy interactions and mergers, we note that neglecting mergers may be a good approximation in MOND as mergers are expected to be rare due to the absence of dynamical friction with the dark matter halo (Renaud et al. 2016).

Zoom In Zoom Out Reset image size

Interestingly, the right panel of Figure 14 shows that the sky-projected aspect ratio distribution of these MONDian disk galaxies is very consistent with that of SDSS spiral galaxies (fracDeV <0.8; see Section 4.3). This also demonstrates that the distinction of spiral from elliptical galaxies based on the linear combination of the exponential and de Vaucouleurs models (Abazajian et al. 2004) succeeds in SDSS, in contrast to the bulge-disk decomposition applied to the Illustris-1 simulation (possible reasons were discussed in Section 7.2). In the future, whether a self-consistent cosmological MOND simulation would be able to reproduce the observed fraction of spiral and elliptical galaxies needs to be explicitly shown. Such MOND simulations are not available at the moment, but are underway in the Bonn-Prague group based on the promising cosmological MOND framework detailed in Haslbauer et al. (2020). Given the problems in reproducing the observed distribution of galactic morphologies with ΛCDM simulations, it is often argued that these simulations depend sensitively on the implemented feedback model and that the problem is likely to be alleviated once the correct feedback implementation has been found. Hydrodynamical MOND simulations of galaxy formation out of post–Big Bang gas clouds, on the other hand, naturally lead to realistic galaxies similar to the observed ones with the available feedback implementations (Wittenburg et al. 2020, Eappen et al. 2022, submitted).