EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 644 (December 2020) A&A, 644 (2020) A56 Full HTML
Open Access
Issue
A&A
Volume 644, December 2020
Article Number A56
Number of page(s) 17
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201936439
Published online 01 December 2020

© T. Devergne et al. 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

According to the standard theory of galaxy formation, the dissipative infall of gas in the gravitational potential wells of dark-matter (DM) haloes forms discs (Fall & Efstathiou 1980); elliptical galaxies are formed by mergers (Toomre & Toomre 1972). Semianalytic models (SAMs) of galaxy formation build on this theory and describe a galaxy as the sum of two components: a disc and a bulge. Observers perform a similar decomposition when they fit galaxies with the sum of an exponential and a Sérsic (1963) profile to compute quantitative morphologies (Simard et al. 2011; Meert et al. 2015, 2016; Dimauro et al. 2018).

This simplification brushes over the complexity and diversity of galactic morphologies, for example, the distinction between normal and barred spirals (Hubble 1926). Gadotti (2009), Weinzirl et al. (2009), Salo et al. (2015), and Erwin (2018, 2019) have considered more detailed models that decompose galaxies into three components: a disc, a bulge, and a bar if present. Despite the clear merit of such decompositions, this approach is possible only for relatively small samples of a few thousand nearby galaxies, such as the Spitzer Survey of Stellar Structure in Galaxies (S4G; Sheth et al. 2010). It would be much more difficult to perform the same analysis on large samples, such as the Sloan Digital Sky Survey (SDSS), or on high-redshift data with even poorer spatial resolution1.

Explaining the broad statistical properties of galaxies in large surveys is the main purpose of SAMs. If the goal is a detailed study of the morphological structure of galaxies, then hydrodynamic simulations are a much better tool. If the observations that we aim to explain cannot distinguish between different types of bulges, then it makes sense to compute the bulge-to-total mass ratio B/T in such a way that any stellar surface-density excess above an exponential fit is assigned to the bulge component, independently of its origin, structure, and kinematics. That is not to say that all bulges are the same.

Many spiral bulges ressemble miniature ellipticals, especially those in galaxies with stellar mass M⋆ ≳ 1011 M⊙ (Fisher & Drory 2011). They are called classical bulges. This similarity suggests a common formation mechanism. In the earliest SAMs (for example, Kauffmann et al. 1993), all bulges were formed through mergers. Some of them never accreted any gas. We call them elliptical galaxies. Others regrew a disc and became spiral galaxies (Baugh et al. 1996).

In this picture, most spirals should be bulgeless because mergers make a negligible contribution to the mass growth of galaxies with M⋆ <  1011 M⊙ (Cattaneo et al. 2011). Bulgeless galaxies are observed (Kormendy et al. 2010), but they do not constitute the majority of the spiral population unless we include dwarf galaxies. Fisher & Drory (2011) find that the fraction of galaxies in which a bulge is detected increases smoothly from ∼20% at M⋆ = 109 M⊙ to ∼100% at M⋆ = 1010.7 M⊙.

Most of the bulges in spirals with M⋆ <  1011 M⊙ are pseudobulges with different kinematics than elliptical galaxies and do not follow the fundamental plane (Kormendy 1982; Kormendy et al. 1993; Kormendy & Kennicutt 2004). To explain these systems, SAMs began to incorporate a second formation mechanism: disc instabilities (Cole et al. 2000; Hatton et al. 2003; Shen et al. 2003).

Self-gravitating thin discs are dynamically unstable (Hohl 1971; Kalnajs 1972). The observation of disc galaxies despite such instability has been one of the historical arguments for DM (Ostriker & Peebles 1973). Morphological features, such as spiral arms, bars, peanut-shaped boxy pseudobulges, rings, and ovals, show that haloes do not completely stabilise discs, however.

Combes & Sanders (1981) used N-body simulations to study the stability of a truncated Toomre (1963) disc and demonstrated that, if the mass of the disc was larger than the mass of the DM within its maximum radius, a persistent bar developed quickly and, after some time, took a more or less pronounced peanut shape when seen edge-on. Observations of peanut-shaped pseudobulges confirm that they are connected with bars and owe their origin to them (Kormendy & Kennicutt 2004). The stronger conclusion that peanut-shaped pseudobulges are nothing more nor less than bars seen edge-on (Combes & Sanders 1981; Combes et al. 1990; Pfenniger & Friedli 1991; Berentzen et al. 1998; Athanassoula & Misiriotis 2002; Athanassoula 2003) is less straightforward from an observational standpoint, but observations of galaxies such as NGC 7582, where the bar is very flat and three times longer than the pseudobulge (Quillen et al. 1997), are consistent with a picture in which the peanut is the vertical extension of a longer, flatter bar (Athanassoula 2005; Wegg et al. 2015).

The literature above demonstrates that classical bulges and pseudobulges are different entities. GALICS 2.0 (Cattaneo et al. 2017) has been the first (and to date only) SAM to treat them as separate components. More detailed investigations aimed at separating pseudobulges from bars or at distinguishing different types of pseudobulges2 are beyond the scope of SAMs3.

Efstathiou et al. (1982, ELN) extended the analysis by Combes & Sanders (1981) to the more realistic case of an exponential disc and found a condition for the circular velocity Vc at 2.2 exponential scale-lengths, where the rotation curve of a self-gravitating exponential disc peaks (Freeman 1970). A thin exponential stellar disc embedded in a static spherical DM halo becomes unstable and develops a bar if

(1)

where Md is the disc mass, Rd is the exponential scale-length and ϵ = 1.1. Christodoulou et al. (1995) used analytic arguments to conclude that a similar criterion with ϵ = 0.9 should apply to gaseous discs.

Since Mo et al. (1998) and van den Bosch (1998, 2000), the ELN criterion (Eq. (1)) has provided the standard description of disc instabilities that all current SAMs adopt (GALACTICUS: Benson 2012; GALICS 2.0: Cattaneo et al. 2017; MORGANA: Lo Faro et al. 2009; SAG: Gargiulo et al. 2015; SANTACRUZ: Porter et al. 2014; YSAM: Lee & Yi 2013; GALFORM: Gonzalez-Perez et al. 2014; LGALAXIES: Henriques et al. 2015; SAGE: Croton et al. 2016).

More realistic simulations (Athanassoula 2008) found that, even in cases where the criterion predicts stability, a bar can still form if resonances destabilise the disc by transferring angular momentum to the halo (ELN’s assumption of a static halo prevents this possibility in their simulations). Moreover, in cases where the ELN criterion predicts instability, the disc can still be stabilised by random motions, which ELN did not consider because of the assumption of thin discs.

Several reasons explain why SAMs have kept using the ELN criterion despite these criticisms:

– The goal of SAMs is to explain the global properties of galaxies (such as the trend of B/T with M⋆) in a cosmological context. A detailed description of galactic dynamics is beyond their scope.

– SAMs separate the formation of haloes from the evolution of baryons within haloes. In reality, baryons affect the radial distribution of subhaloes within groups (Libeskind et al. 2010) as well as the density profiles (Pontzen & Governato 2012; Macciò et al. 2012; Teyssier et al. 2013; Di Cintio et al. 2014; Tollet et al. 2016), spins, and shapes (Bryan et al. 2013) of DM haloes, and the clustering of galaxies on scales as large as ∼1 Mpc (van Daalen et al. 2014). Assuming a static halo is a simplification, but it is consistent with the semi-analytic approximation.

– From the perspective of SAMs, the ELN criterion is first and foremost a criterion for the ratio of the disc mass Md to the halo mass Mvir. This ratio is largely determined by feedback processes (Dekel & Silk 1986; Silk & Rees 1998; Brook et al. 2012; Anglés-Alcázar et al. 2017; Tollet et al. 2019) that cannot be modelled accurately. Disc sizes are based on the crude assumption (Kimm et al. 2011; Stewart et al. 2013; Jiang et al. 2019) that angular momentum is conserved (Mo et al. 1998) and they, too, affect the circular velocity in Eq. (1). The errors from these uncertainties are likely to be more significant than those from the modelling of disc instabilities themselves.

– A thick-disc component with higher velocity dispersion makes discs more stable (Toomre 1964) and there is observational evidence that thick discs are ubiquitous (Yoachim & Dalcanton 2006; Comerón et al. 2011), but what is their origin? If it is secular heating of the thin disc (Villumsen 1985; Villalobos & Helmi 2008; Schönrich & Binney 2009; Steinmetz 2012), then it is logical that simulators should not put in their inititial conditions the effects of processes they want to simulate. If it is interactions with satellites (Quinn et al. 1993; Abadi et al. 2003), then it is legitimate to neglect them in a model for the evolution of isolated galaxies that have not experienced any merger. A third possibility is that thick discs are relics from gas-rich, turbulent, clumpy discs at high redshift and that thin discs formed later (Brook et al. 2004; Bournaud et al. 2009). The problem is that there is no SAM to compute disc scale-heights (but see Efstathiou 2000 for an attempt in this direction). In absence of physical arguments for one scale-height or another, the simplest assumption (thin discs) is the most reasonable.

The main problem is another. The ELN criterion tells us whether a disc is likely to become unstable but not the mass of the bar or pseudobulge that is likely to form because of that instability. To compute this mass, SAMs must incorporate additional assumptions. Cole et al. (2000; also see Gonzalez-Perez et al. 2014 and Gargiulo et al. 2015) considered an extreme model in which discs evolve into bulges whenever the instability criterion in Eq. (1) is satisfied. Hatton et al. (2003) and Shen et al. (2003) proposed a more conservative model in which matter is transferred from the disc to the bulge so that Md decreases until the disc becomes stable again. In discs where Vc(2.2Rd) falls just slightly short of the critical value required for stability, the results obtained with the two methods vary wildly. A galaxy for which B/T = 1 in the first model may have B/T = 0.1 in the second one. Therefore, even if the ELN criterion were fully reliable, its implications for galactic morphologies would still be considerably uncertain.

In this article, we present a series of simulations in which we follow the evolution of a thin exponential stellar disc embedded in a static spherical Navarro et al. (1997, NFW) halo. Our computational set-up is very similar to the one by ELN, although we explore a larger space of parameters and have better resolution. Unsurprisingly, our simulations confirm ELN’s previous findings, but that is not the purpose of our research. Our question is: as different SAMs make very different predictions for B/T even if they are all based on the ELN criterion, which approach (if any) agrees better with the values of B/T measured in the simulations used to establish the ELN criterion?

As any approximations in the simulations are passed on to SAMs through the ELN criterion, the agreement of a SAM with the simulations is no guarantee of its being a correct description of disc instabilities, but, at least, it proves self-consistency (the SAM faithfully reproduces any biases of the simulations). On the contrary, lack of agreement indicates that the SAM is adding biases of its own on top of those already present in the simulations. This is why it is high time for a sanity check and, possibly, new prescriptions that may improve those used to calculate B/T in SAMs.

The structure of the article is as follows. In Sect. 2, we describe our initial conditions, which are entirely specified by three parameters (the ratio rd = Rd/Rvir of the disc scale-length Rd to the virial radius Rvir, the ratio md = Md/Mvir of the disc mass Md to the total mass Mvir within the virial radius, and the concentration c of the DM halo), the explored parameter space and the computational strategy (we use the adaptive-mesh-refinement [AMR] code RAMSES; Teyssier 2002). In Sect. 3, we present our findings for the dependence of B/T on rd, md, and c. In Sects. 4 and 5, we compare our results with previous models and the observed morphology–luminosity relation, respectively. In Sect. 6 we explore the effects of incorporating our findings into the GALICS 2.0 SAM and their implications for galactic morphologies. Section 7 discusses our results and summarises the conclusions of the article.

2. Computational set-up

2.1. Initial conditions

We make our problem dimensionless by expressing all lengths in units of the virial radius Rvir, all masses in units of the virial mass Mvir and all speeds in units of the virial velocity . Owing to the axial symmetry of our initial configuration, we adopt cylindrical coordinates (r, z, ϕ), where z is the direction of the disc’s rotation axis. Here and throughout this article, upper and lower-case letters refer to dimensional and dimensionless quantities, respectively. As Mvir is the total mass within Rvir, md = Md/Mvir and 1 − md are the dimensionless masses of the disc and the DM halo, respectively.

The dimensionless mass of the DM within a sphere of radius r is given by:

(2)

where f(x) = ln(1 + x) − x/(1 + x) and c is the concentration of the NFW profile.

We assume that the disc is exponential and isothermal in the vertical direction (for example, Villumsen 1985; Efstathiou 2000). These assumptions give the density distribution:

(3)

where h is disc’s vertical scale-length (all quantities in Eq. (3) are adimensional). As our goal is to study the stability of thin discs (discussion in Sect. 4), we run all our simulations for h/rd = 0.044. This value is small but not unrealistic. Bland-Hawthorn & Gerhard (2016) find that the Milky Way has a thin-disc vertical scale-height of 220–450 pc and an exponential scale-length of ∼2.5 kpc. Hence, for the Milky Way, h/rd is in the range 0.088–0.18. Figure 1 shows the isodensity contours that contain 40%, 50%, 60 %, 70%, 80%, and 90% of the disc mass for our initial configuration.

thumbnail Fig. 1.

Initial conditions and refinement regions. The black, blue, violet, green, yellow, and red curves show the isodensity contours that contain 90%, 80%, 70%, 60%, 50%, and 40% of the disc mass at t = 0, respectively. The black, blue, violet, green, yellow, and red dashed lines show the cylinders within which the cell size equals 1/8, 1/16, 1/32, 1/64, 1/128, and 1/256 of the disc exponential scale length, respectively.

Equation (3) is used to generate random coordinates for a million stellar particles within the optical radius ropt = 3.2rd that contains 83% of the mass of the disc (Fujii et al. 2011 demonstrated that numerical heating through close encounters becomes negligible when the disc is resolved with at least a million particles). All stellar particles have equal mass. This article is on stellar discs, but our discs contain a small amount of gas to pave the way for a second article on the formation of bulges in discs with gas. We use a gas fraction of 2% because in massive low-redshift galaxies the gas fraction varies from close to zero to a few percent and is rarely larger than 5–10% (for example, Combes et al. 2013), but we also ran four simulations without gas to check the impact that even a small gas fraction could have on our results. The gas in our simulations is isothermal at 104 K and is not allowed to form stars. The total mass of the stellar particles is 0.98 × 0.83 md in dimensionless virial units.

The circular velocity vc(r) is the speed that a star must have to be on a circular orbit of radius r. Its value is the sum in quadrature of the contributions from the disc and the halo:

(4)

where vd(r) is computed as in Freeman (1970) and

(5)

Stars have velocity dispersion:

(6)

determined from Eq. (3) through the requirement that our initial condition should be in equilibrium (albeit unstable). Hence, their velocities:

(7)

will be the sum of an ordered rotational component (oriented as ) and a random deviate from a Maxwellian distribution with velocity dispersion σ. The assumption of an isotropic velocity dispersion is motivated by simplicity, but it is not unreasonable, since the radial and vertical velocity dispersions in the Solar Neighbourhood are (35 ± 5) km s−1 and (25 ± 5) km s−1, respectively (Bland-Hawthorn & Gerhard 2016).

The rotation speed vrot(r) = ⟨vϕ(r)⟩ equals the circular velocity vc(r) only for σ = 0. For σ >  0, vrot and vc are linked by the condition:

(8)

which gives:

(9)

since ⟨Δvϕ⟩ = 0 and .

Our rotation speeds and thus our particle velocities are computed using Eq. (9) everywhere except in a small central region where we set vrot = 0 because vc <  σ. This central region has r ≪ rd by construction, since the normalisation of σ(r) is determined by the disc scale-height (Eq. (6)) and we have assumed that h ≪ rd. However, the fact that σ is computed considering the vertical equilibrium only and that at r ≲ 0.07 rd implies that the central region will expand a little bit when the initial conditions are allowed to relax.

The global stability of our initial conditions is explored through the simulations presented in this article (Sect. 3). Their local stability is a different problem. The discussion of how local disc instabilities can affect the final B/T of globally unstable discs is postponed to Sect. 4.

2.2. Parameter space

Having fixed the disc scale height h = 0.044 rd and the gas fraction fgas = 0.02, our initial conditions are entirely determined by three parameters: rd, md, and c.

The dimensionless radius of a thin exponential disc of dimensionless mass md embedded in an NFW halo with concentration c,

(10)

(Cattaneo et al. 2017) is completely determined by the disc’s spin parameter

(11)

which is identical to the halo’s spin parameter if we follow (Mo et al. 1998) and we assume that the disc and the halo have the same specific angular momentum (for Jd/Md = Jvir/Mvir, Eq. (11) corresponds to the Bullock et al. 2001 definition of the spin parameter).

Equation (10) implies that we can use the spin distribution of DM haloes from cosmological N-body simulations to derive plausible values for rd. This is true even if conservation of angular momentum is a crude approximation (Kimm et al. 2011; Stewart et al. 2013; Jiang et al. 2019) because the model by Mo et al. (1998) reproduces the correct mass – relation for disc galaxies (for example, Cattaneo et al. 2017).

For a flat rotation curve with vc = 1, Eq. (10) gives rd = λ/2. We do not make this approximation. We substitute Eqs. (4) into (10) and solve Eq. (10) numerically for each combination of λ, md, and c that we wish to consider. However, we find that the approximation rd ≃ λ/2 is accurate at the 20% level.

Muñoz-Cuartas et al. (2011) found that λ follows a log-normal distribution with and σln λ = 0.57. Burkert et al. (2016) found a similar distribution with and σln λ = 0.46. So did Cattaneo et al. (2017) with and σln λ = 0.57 (we use the same symbol λ for the spin of the halo and the spin of the disc because we have assumed they are identical). From these figures, we conclude that λ = 0.05 corresponds to a typical value and that ∼80% of all haloes lie in the interval 0.025 <  λ <  0.1. Haloes with λ <  0.0125 comprise less than 2% of all systems. Based on these considerations and the fact that rd is the parameter to which our results are most sensitive, we have explored six values of λ (0.011, 0.018, 0.025, 0.035, 0.05, 0.1).

We sample the parameter space by λ rather than rd because we started this research to find a prescription that may be useful to assign morphologies to galaxies in both semianalytic and halo-occupation-distribution models (for example, Behroozi et al. 2013, 2019; Moster et al. 2013, 2018; Tollet et al. 2017). If the objective is to find a recipe to populate haloes with galaxies, then λ and c are the quantities that we directly measure in N-body simulations, and we should like to find B/T as a function of them. However, the halo spin has no direct effect on our simulations because our haloes are spherical and static. Hence, it is the dependence of B/T on rd that we probe in reality.

In contrast to λ, which appears to be a true random quantity, md strongly correlates with Mvir. Considerable observational evidence shows that the stellar-to-halo mass ratio increases with Mvir over the range of halo masses where spiral morphologies are prevalent (Papastergis et al. 2012; Leauthaud et al. 2012; Reyes et al. 2012; Behroozi et al. 2013; Moster et al. 2013; Wojtak & Mamon 2013; Cattaneo et al. 2017; Tollet et al. 2017). In a halo with Mvir = 1012 M⊙, the typical stellar mass of the central galaxy is M⋆ = 4 × 1010 M⊙ (md = 0.04 if we assume that the galaxy started as a pure disc). In a halo with Mvir = 1011 M⊙, the typical stellar mass is two orders of magnitudes lower and md is lower by about a factor of ten. Hence, we explore four values of md (0.005, 0.01, 0.02, 0.04) that cover the typical range of M⋆/Mvir from dwarf to Milky-Way-type galaxies.

Concentration has a much weaker dependence on halo mass. Dutton & Macciò (2014) find with M12 = Mvir/1012 M⊙ at redshift zero (also see Muñoz-Cuartas et al. 2011). Thus, we expect the mean concentration to vary from c = 13 for a dwarf galaxy to c = 10 for a Milky-Way-type one. The real concentration range is larger because the scatter is significant. We thus consider three concentration values (5, 10, 15) such that the interval 5 <  c <  15 contains ∼80% of the haloes with Mvir >  1011 M⊙.

Six values for λ, four values for md and three values for c make 72 combinations in total. We have simulated only 34 of these 72 combinations because we have explored the cases λ = 0.018 and λ = 0.035 only for c = 10 and because during our work it became obvious that certain sets of parameters would not form a bulge (Sect. 3). Hence, it would have been pointless to study them.

In addition to these 34 simulations, we have run another 14 simulations to address specific issues, such as: (1) relaxation effects due to the assumption of a razor-thin disc, (2) the assumption of a static halo, and (3) the impact of a small gas fraction. The parameter values for these simulations are listed in Tables 1–3, respectively. The simulations run for this study are 34 + 14 = 48 in total.

Table 1.

Simulations with relaxed initial conditions: parameter values and difference in B/T with respect to the unrelaxed simulations.

Table 2.

Simulations with a live halo: parameter values and difference in B/T with respect to those with a static halo.

Table 3.

Simulations without gas: parameter values and difference in B/T with respect to those with fgas = 0.02.

2.3. Refinement strategy

We evolved our initial conditions with the AMR code RAMSES (Teyssier 2002) until the galaxies converged to a stable configuration. This usually occurs within 2 Gyr.

As we use a Eulerian code, the most important numerical aspect of our work is the choice of the grid on which we integrate the equations of motion for the stellar fluid. We centre our discs on a cubical grid with 1283 cells and side length 32 rd. Hence, the resolution on the scale of the entire computational volume is l = 0.25 rd. The resolution is increased within six nested cylinders by a factor of two each time (Fig. 1). The radii r/rd ≃ 3.8,  3.2,  2.7,  2.3,  1.9, and 1.7 of the six cylinders enclose the isodensity contours that contain 80, 70, 60, 50, 40, and 30% of the disc mass and correspond to l/rd = 1/8,  1/16,  1/32,  1/64,  1/128, and 1/256, respectively.

Height equals radius in all cylinders except the innermost one, which has a height-to-radius ratio of 1 : 5 to ensure than the vertical structure of the thin gaseous disc is well resolved even though the gas fraction in our simulations is so small that it has no dynamical effects. The second innermost cylinder (the one marked in yellow in Fig. 1) has a semi-height (0.95 rd) large enough to ensure that the disc remains well resolved when it buckles up, since the scale-length (and therefore the scale-height) of a pseudobulge is usually several times smaller than rd (Gadotti 2009).

3. Results

Some of our discs our stable. They do not show any bar or pseudobulge even after 2 Gyr. In unstable discs, a bar forms rapidly, but after 2 Gyr the growth of B/T becomes very slow and, in many cases, almost inexistent (see below for details of how we measure B/T). Figures 2–4 show 27 simulated galaxies from our main sample of 34 galaxies in total (we have not shown the simulations λ = 0.018 or λ = 0.035). The galaxies are shown face-on at the first time t >  2 Gyr when the growth of B/T has stabilised. In most galaxies, B/T has converged at t <  2.5 Gyr.

thumbnail Fig. 2.

Face-on view of the simulated galaxies at the first time t when B/T has stopped growing (2 Gyr <  t <  3 Gyr). The image brightness is proportional to the logarithm of stellar surface density. This figure show the simulations with c = 5, sorted by md = Md/Mvir (lines) and halo spin λ (columns). On each image, we have shown the corresponding value of rd = Rd/Rvir. Bulge-to-total stellar mass ratios B/T have been computed by decomposing the stellar surface-density distribution into an exponential and a Sérsic (1963) profile. Galaxies that could be fitted with a single exponential profile have B/T = 0. When the fit required a bulge component, its Sérsic index n has been shown.

thumbnail Fig. 3.

Same as Fig. 2 with c = 10 instead of c = 5. To avoid overcrowding the figure, we have not shown the simulations with λ = 0.018 and λ = 0.35.

thumbnail Fig. 4.

Same as Fig. 2 with c = 15 instead of c = 5. The column λ = 0.1 is now missing because we know that simulations with λ = 0.1 will not form any (pseudo)bulge.

Figures 2–4 correspond to c = 5,  10,  15, respectively. They portray the variety of morphologies that we can reproduce, from Scd to SBa types, simply by changing the values of our parameters. Since our galaxies are isolated and the gas fraction is very small, bar instabilities are the only physical mechanisms through which we expect that a bulge should form in our simulations (see the discussion of this article). In agreement with this expectation, 10 out of 17 galaxies with B/T >  0.1 display morphologies that are clearly barred. We have also looked at our galaxies edge-on. In most bulges, an X-shaped peanut is clearly visible (Fig. 5).

thumbnail Fig. 5.

Galaxy in the simulation with c = 10, md = 0.04, λ = 0.025, viewed edge-on. An X-shaped pseudobulge is clearly visible. This galaxy has been chosen as an example and is by no means atypical. The face-on image (Fig. 3) shows that this galaxy has a bar, but not a prominent one.

The galaxy with md = 0.01 and λ = 0.011 on Fig. 2 is the most conspicuous example of a prominent bulge (B/T = 0.24) without any evidence for a bar. Such cases are rare. We have looked at this galaxy edge on. The bulge looks boxy, although the peanut shape is much less obvious than in Fig. 5.

The stellar surface-density profile Σ*(r) of each galaxy has been fitted with the sum of an exponential and a Sérsic (1963) profile (Fig. 6a). This decomposition has been used to compute a B/T ratio for each galaxy and a Sérsic index n for the pseudobulges of all galaxies with B/T >  0. We assume B/T = 0 when a visual inspection shows the surface-density profile is consistent with a single exponential function (Fig. 6b).

thumbnail Fig. 6.

Stellar surface-density for two of our simulated galaxies (black symbols). Galaxy a has been fitted with the sum of an exponential profile (blue line) and a Sérsic profile, with n = 1 in this particular case (red line). Its bulge-to-total mass ratio is B/T = 0.22. The black curve shows the sum of the two components. Galaxy b is consistent with a single exponential profile (blue line).

The decrease in Σ* in the central region (Fig. 6a) is an artefact due to how we set up the initial conditions for the velocity field. Its cause is the relaxation effect discussed in Sect. 2.1 after Eq. (9) coupled to the lower resolution in the central region. The radius r ∼ 0.1 rd within which Σ* begins to decline is so small that it does not affect our results, except when B/T is very low, in which case it limits our accuracy in fitting the surface-density profile of the bulge component.

All galaxies with B/T >  0.1 have n ≲ 1.2, as expected for pseudobulges. In systems with B/T ≤ 0.1, the bulge is small and the Sérsic index may not be robust, also because of the effect mentioned above. The galaxy with md = 0.04 and λ = 0.05 on Fig. 4 exemplifies this situation. Bulges with n >  2 are usually classical bulges and are seldom encountered in galaxies with B/T ≲ 0.1 (Fisher & Drory 2008). Hence, our simulations should not form any bulges with n = 4, let alone one with B/T = 0.1. A visual inspection of the bulge/disc decomposition shows a poor fit to the stellar surface density in the central region. The Sérsic index n = 4 measured for this galaxy cannot be considered meaningful.

For given c and λ, higher values of md correspond to earlier-type morphologies in the Hubble sequence and to higher B/T (Figs. 2–4). We can also increase B/T by reducing λ at constant c and md. At constant λ and md, B/T usually decreases with c. Hence, if a simulation does not form a pseudobulge, simulations with lower md, higher λ or higher c will not form one either (supposing the other two parameters were kept fixed). Thanks to this finding, which we have verified in a few cases, we could nearly half the simulations required for this work (there is no need to simulate galaxies for which we know in advance that we shall find B/T = 0).

The qualitative findings above have a simple interpretation, which will become apparent once we have examined the ELN criterion in greater detail.

For an exponential disc:

(12)

where md(2.2rd) = 0.65md is the disc mass within 2.2rd (Freeman 1970). The first equality in Eq.(12) implies that Eq. (1) can be rewritten as:

(13)

By introducing the new parameter α = 0.31/ϵ, Eq. (13) can be written in the alternative form:

(14)

where α = 0.28 for ϵ = 1.1 (ELN) and α = 0.26 for ϵ = 1.2 (Christodoulou et al. 1995)4.

The second equality in Eq. (12) shows that Eq. (14) is equivalent to a criterion on the disc fraction within 2.2rd because is the total gravitational acceleration at r = 2.2rd in dimensionless units and is the disc’s contribution. Therefore,

(15)

is the disc’s fractional contribution to the gravitational acceleration at r = 2.2rd and is related to disc’s mass fraction within 2.2rd by the equation:

(16)

Van den Bosch (1998) was the first to use the ELN criterion to predict disc stability in a SAM, but he applied Eq. (14) at r = 3rd rather than at r = 2.2rd for the reason that the rotation curve vc(r) peaks at larger radii when the halo’s contribution is considered. We have verified that Eq. (16) keeps holding at r = 3.2rd (the factor 1.3 does not change to the first decimal digit).

Equations (12) and (14) explain why B/T increases with md and why it decreases with λ and c (Figs. 2–4). The higher md, the higher the disc fraction within a given radius. The larger the disc, the more DM it will contain. The more concentrated the halo is, the more DM there will be within the disc. Our next step is to go beyond these semi-qualitative considerations and to explore the relation between B/T and fd in quantitative detail.

Figure 7 shows B/T versus fd at r = 2.2rd (where the rotation curve of an exponential disc peaks) and at the optical radius ropt = 3.2rd, although we have investigated the relation at other radii down to r = 0.5rd. The vertical blue lines correspond to:

(17)

thumbnail Fig. 7.

Relation between bulge-to-total mass ratio B/T and initial disc fraction fd(r) for r = 2.2rd (left) and r = 3.2rd (right). Triangles pointing up, circles, and triangles pointing down correspond to simulations with c = 5, c = 10, and c = 15, respectively. The colours correspond to different disc-to-virial mass ratios: md = 0.005 (green), md = 0.01 (black), md = 0.02 (blue), and md = 0.04 (red). Symbols with B/T = 0.01 correspond to bulgeless galaxies. They have been assigned B/T = 0.01 merely to be able to show them on a logarithmic diagram. The thick solid black lines are log-log linear least-squares fits to the symbols with B/T >  0.01. They correspond to Eqs. (18) and (19). The vertical blue lines show the critical fd that corresponds to the α = 0.28 (equivalent to assuming ϵ = 1.1 in Eq. (1)) The coloured curves compare our results to previous models (Cole et al. 2000; Hatton et al. 2003; Shen et al. 2003; Cacciato et al. 2012).

for α ≃ 0.28.

The qualitative picture is the same at r = 2.2rd and r = 3.2rd. At , all galaxies are pure discs (the galaxies with B/T = 0.01 in Fig. 7 are bulgeless galaxies; we have assigned them B/T = 0.01 merely to be able to show them on a logarithmic plot). At , all galaxies develop a pseudobulge and display a tight correlation between B/T and fd. At intermediate fd, galaxies with B/T = 0 and B/T >  0 coexist. Therefore, the ELN criterion may fail to discriminate between stable and unstable discs (Athanassoula 2008; Fujii et al. 2018). Nevertheless, the notion of a critical fd that separates the two populations remains fundamentally valid.

If we restrict our attention to galaxies with B/T >  0 and apply a linear least-squares fit to the relation between log(B/T) and log fd, we find:

(18)

and

(19)

(thick black solid lines in Fig. 7a and b, respectively). Equation (18) implies B/T = 0.04 for . With Eq. (19), B/T = 0.04 is attained for the slightly lower value . If we use B/T = 0.04 as the minimum bulge-to-total stellar mass ratio below which galaxies are effectively bulgeless, then these figures imply that the critical fd is lower at 3.2rd than it is at 2.2rd. This makes sense because fd(r) decreases at large radii in both stable and unstable discs (the density of the baryons in the disc decreases faster than the density of the DM). Hence, it is logical that the critical fd that separates them should decrease with radius, too.

The standard deviation of log(B/T) from the linear least-squares fit of B/T versus fd(r) (that is, the root-mean-square residuals) decreases slowly but systematically from 0.12 dex at r = 0.5rd to 0.11 dex at r = 3.2rd. Hence, the ELN criterion is fairly insensitive to the radius at which it is applied, but a criterion at r = 3.2rd appears to be preferable.

If the condition is applied at r = 3.2rd, then all discs with are correctly identified as unstable (Fig. 7b). For these galaxies, Eq. (19) gives B/T with an accuracy of ∼30%. In contrast, we find four galaxies with where the ELN criterion failed to predict instability. However, all these galaxies have B/T <  0.06. To assume that galaxies with B/T ≲ 0.06 are bulgeless is not such a great approximation. Even observers may not be able to detect bulges that small, particularly for galaxies at distances ≳100 Mpc.

4. Comparison with previous models

The blue curves and the red curves in Fig. 7 compare our results to the two most common assumptions encontered in SAMs. The blue curves correspond to the model by Cole et al. (2000): as soon as a disc becomes unstable, it collapses into a bulge. The difference from the black lines is plainly obvious. The red solid curves correspond to the model by Hatton et al. (2003) and Shen et al. (2003): as soon as a disc becomes unstable, matter is transferred from the disc to the bulge until the disc is marginally stable.

We compute the analytic form of the red curves by assuming that initially all the stars are in the disc, so that the stellar masses within 2.2 rd and 3.2 rd are 0.65 m* and 0.83 m*, respectively. Equation (16) gives:

(20)

with β = 1.3 × 0.65 = 0.85 or β = 1.3 × 0.83 = 1.1, depending on the radius at which the instability condition is considered (mDM is the mass of the DM within this radius). The disc becomes marginally stable when:

(21)

where mb and m* − mb are the final bulge mass and the final disc mass, respectively. We compute the final bulge fraction B/T = mb/m* by solving Eqs. (20) and (21), and find:

(22)

where . Equation (22) shows that B/T → 0 for and that:

(23)

for fd → 1. The limit for fd → 1 reduces to B/T → 0.71 for and β = 1.1.

The red solid curves are in better agreement with the symbols than the blue curves, but they systematically overestmate B/T at high masses and underestimate it at low masses. The red dashed curves show that increasing does not improve the overall agreement between Eq. (22) and the simulations because it reduces the discrepancy with the symbols at high fd, but it makes it worse at low fd.

We have also compared our results to a third model based on the conservation of energy (Cacciato et al. 2012). The specific energy of the stars in a cold disc is:

(24)

where Φ is the gravitational potential. When stars migrate from the disc to the bulge, they move to a lower energy level. The gravitational energy released by this process increases the stellar velocity dispersion σ. Cacciato et al. (2012) used the conservation of energy to study the evolution of two-component discs (gas plus stars) in a cosmological context. Here, we consider a simpler version of their model, in which the disc is isolated (there is no accretion) and purely stellar (there is no gas and thus no star formation). The internal energy of the stars, , equals the energy released by the formation of a bulge, MbΔe.

To compute Δe, one should know the initial and final radii of the stars that migrated to bulge, but one can expect within a factor of order unity (Cacciato et al. 2012) and compute B/T = Mb/M* by solving the equation:

(25)

The calculation of σ is based on the assumption that σ self-regulates so that the Toomre (1964) stability parameter is always close to critical value Qc for local stability. This assumption gives:

(26)

with Qc = 1.68 for a razor-thin sheet of stars with a Maxwellian velocity distribution (see Binney & Tremaine 2008).

Let expand Eq. (26) by using Ω = Vc/R:

(27)

is the total gravitational acceleration at radius R. The numerator measures the disc’s contribution after the formation of a bulge (the mass of the disc within R scales as ΣR2). This contribution can be decomposed into the product of two terms: the stellar contribution relative to the total stars plus DM (the fd factor, since all the stars are in the disc at t = 0) and the stellar fraction in the disc after a bulge has formed, Md/M*. The product is not exact because matter changes geometry when it migrates from the disc to the bulge, but neglecting that has a small effect. Equation (27) can thus be rewritten:

(28)

where η is a fudge factor (the disc mass within R scales as, but is not equal to, ΣR2).

By substituting Eq. (28) into Eq. (25) and using Md = M* − Mb, we find:

(29)

where all uncertainties, and notably the one on , are reabsorbed in the value of η.

It is straightforward to solve the second-degree Eq. (29) and compute B/T as a function of fd. The results are shown by the green curves on Fig. 7. The best agreement with our simulations is for η = 0.41.

The model by Cacciato et al. (2012) predicts a shallower slope for B/T versus fd than what we find in our simulations, especially at high fd and when it is applied at r = 2.2 rd. However, it is the model that works best.

5. Comparison with observations

Our simulations predict that B/T is related to Vd/Vc. Persic & Salucci (1995) and Persic et al. (1996) used 967 galaxies with I-band photometry and Hα data from Mathewson et al. (1992) to conduct one of the most careful and systematic study of the rotation curves of galaxies. The main result from Persic et al. (1996) was that Vd/Vc increases with I-band luminosity from Vd(Ropt)/Vc(Ropt) = 0.4 at MI = −18.5 to Vd(Ropt)/Vc(Ropt) = 0.96 at MI = −23.2, where Ropt = 3.2Rd (the brightest discs are close to being maximal). This finding is just another way to say that the stellar-to-total mass ratio increases with stellar mass (for example, Dekel & Silk 1986; Papastergis et al. 2012; Behroozi et al. 2013; Moster et al. 2013). We can use it to transform Eq. (19) into a prediction for how B/T grows with I-band luminosity (Fig. 8, black curve). In this section, we use three data sets to compare this prediction with observations.

thumbnail Fig. 8.

⟨B/T⟩ versus I-band magnitude MI. The black curve shows our prediction. The black circles are from Mathewson et al. (1992) when we assign to each Hubble type the ⟨B/T⟩ from Graham & Worley (2008; the uncertainty of estimating B/T determines the size of the error bars). The red squares are from Simard et al. (2011) after removing all galaxies with B/T >  0.7. The up-pointing filled blue triangles and the down-pointing empty blue triangles are from Salo et al. (2015), S4G) when we assign the bar to the bulge and the disc, respectively. The error bars on the red squares and the blue filled triangles show the standard error on the mean. We have not shown the error bars on the empty triangles so that they do not overlap with those on the filled triangles, but their width is comparable.

The catalogue in Mathewson et al. (1992) contains only 1355 galaxies and no quantitative morphologies, but it has the advantage of being the parent sample in Persic & Salucci (1995): using the morphology – luminosity relation from the same data set that we have used to pass from Vd(Ropt)/Vc(Ropt) to MI ensures that our analysis is internally consistent.

Hubble-type information can obviate the lack of quantitative morphologies. In Fig. 8, the black circles with error bars show the mean B/T from Graham & Worley (2008) and the mean I-band magnitude MI in the catalogue from Mathewson et al. (1992) for each Hubble type. Our predictions are in good agreement with the black circles within the observational uncertainties. The analysis by Hubble type has two limitations, however. First, Hubble type is not a quantitative measurement of B/T. Second, we have not considered the frequency of different Hubble types in Mathewson et al. (1992)’s catalogue (Table 4).

Table 4.

Morphological composition of the catalogue in Mathewson et al. (1992).

Simard et al. (2011) fitted the r-band surface-brightness profiles of 1.12 million galaxies from the SDSS with a Sérsic plus exponential model. Their B/T values are therefore directly comparable to ours (to the extent that r-band luminosity traces stellar mass).

To pass from r to I-band, we use Lupton’s formula5:

(30)

We retrieve Mi and Mz for galaxies in the SDSS Data Release 7 (DR7) by using the selection query:

SELECT
   objid, ra, dec, run, camcol, rerun,
   petromag_i, petromag_z,
   extinction_i, extinction_z
FROM
   PhotoPrimary
WHERE
   flags
   AND (dbo.fPhotoFlags('SATURATED')+
   dbo.fPhotoFlags('DEBLENDED_AS_PSF'))= 0
   AND (petroMag_r- extinction_r)
   BETWEEN 14.0 and 18.0
   AND Type=3

and by cross-correlating the results with the catalogue from Simard et al. (2011).

This sample selection includes elliptical galaxies, while we are interested in the B/T – MI relation for spiral galaxies. We circumvent this problem by removing all bona fide ellipticals with B/T >  0.8 and n >  2. This selection is based on the distribution of galaxies on a B/T versus n diagram (Fig. 14 of Simard et al. 2011), which shows two main populations: a spiral/S0 population with 0 <  B/T <  0.7 and 0 <  n <  7 (the mean B/T grows with n), and an elliptical population with B/T >  0.8 and n >  2. The diagram also shows a third population with B/T >  0.7 and n <  1.5 (most of these galaxies have B/T ∼ 1). Visual inspection (for example, Fig. 9) suggests that the bulges of these galaxies are misclassified discs (discs for which a Sérsic model with n <  1 worked better than an exponential fit). We therefore assign B/T = 0 to galaxies for which Simard et al. (2011) find B/T >  0.7 and n <  1.5.

thumbnail Fig. 9.

Composite image of SDSS J141007.94-002348.4 from skyserver.sdss.org. SDSS J141007.94-002348.4 is a spiral galaxy with MI = −21.8. Simard et al. (2011) assign to it B/T = 0.93 but a bulge Sérsic index of only n = 0.5.

The red squares in Fig. 8 show the B/T – MI relation for the subsample from Simard et al. (2011) selected according to the above criteria. The error bars are small because they show the error on the mean:

(31)

which is much smaller than the standard deviation σB/T (N is the number of galaxies per magnitude bin).

Our predictions are in good agreement with the red squares at MI  ≲   − 20.5. At MI >  −20.5, the B/T measured by Simard et al. (2011) are larger than those measured by Graham & Worley (2008) for galaxies with Hubble type T ≥ 5, which comprise the dominant population in this luminosity range. One explanation for this discrepancy is the challenge of fitting correctly a faint bulge component that contributes to ≲10% of the total light when the uncertainty on B/T can be as large as 20% (see Figs. 18 and 19 of Simard et al. 2011). Another explanation is that bulge/disc decompositions based on a single photometric band are not reliable (observations discussed in the next section suggest that this may be the case).

Salo et al. (2015) used the Spitzer Survey of Stellar Structure in Galaxies (S4G) to study the morphologies of 2352 with distances < 40 Mpc. The S4G’s excellent image quality allowed Salo et al. (2015) to fit three components: an exponential disc, a Sérsic bulge, and a bar modelled with a Ferrer profile. The drawback is a smaller sample and thus poorer statistics.

In Fig. 8, the filled blue triangles with error bars show B/T in the S4G when the bar is assigned to the bulge. The down-pointing open triangles show the difference when the bulge is identified with the Sérsic component only.

We note that Salo et al. (2015) published stellar masses. We converted them to baryonic masses using the relation between gas fraction and stellar mass from Boselli et al. (2014) and we used the relation between baryonic mass and MI relation for the galaxies of Persic & Salucci (1995) to obtain a final result in I-band magnitudes (the baryonic masses for the galaxies of Persic & Salucci 1995 were measured dynamically by performing a galaxy/halo decomposition, they were not inferred by fitting a stellar population synthesis model).

The higher-quality measurements of B/T from Salo et al. (2015) are in agreement with those by Simard et al. (2011) at MI <  −22 and confirm our suspicion that the latter are probably overestimated at low luminosities. Our agreement with the filled blue triangles is very good over the entire luminosity range −23 <  MI <  −20, especially given the simplifying assumptions we have made (isolated system, axisymmetric initial conditions, purely stellar disc). The fact that the agreement with the filled triangles is better than the agreement with the open ones suggests that our two-component fit effectively assigns most of the bar mass to the bulge component.

We have also checked if our predictions are consistent with the data for the Milky Way. Bland-Hawthorn & Gerhard (2016) studied the disc’s contribution fd to the gravitational acceleration at 2.2 Rd and found fd = 0.42–0.74. The main uncertainty is value of Rd. We predict B/T = 0.09–0.16 for fd = 0.53 (Rd = 2.6 kpc; Robin et al. 2003; Jurić et al. 2008) and B/T = 0.19–0.33 for fd = 0.74 (Rd = 2.15 kpc; Bissantz & Gerhard 2002; Bovy & Rix 2013). Observationally, B/T ranges between 0.19 and 0.32 depending on whether the long bar is assigned to the disc or the bulge (Bland-Hawthorn & Gerhard 2016). Hence, our predictions are consistent with the data for the Milky Way within the observational uncertainties.

6. Application to SAMs

Our key result is a model for B/T as a function of fd, the disc’s contribution to the total gravitational acceleration. To use this model predictively, we must know how fd depends on observable quantities such as luminosity or stellar mass. In Sect. 4, we followed a semi-empirical approach. The conversion from MI to fd was based on observational data. Here, Eq. (19) is incorporared into the GALICS 2.0 SAM (Cattaneo et al. 2017; Koutsouridou & Cattaneo 2019). Hence, the values of fd used to compute B/T are those predicted by our SAM.

The main difference with the semi-empirical approach is that, in our SAM, mergers cannot be switched off. Hence, we must consider two mechanisms for the formation of bulges: mergers and disc instabilities. GALICS 2.0 accounts for these two mechanisms by separating galaxies into four components: the disc, the pseudobulge, the classical bulge, and the central cusp. The cusp is composed of baryons that fell to the centre at the last major merger. It is absent in core ellipticals, which were formed in dissipationless mergers (Kormendy et al. 2009). The cusp component is irrelevant for this article. Relevant is the distinction between classical bulges formed through mergers and pseudobulges formed through disc instabilities (GALICS 2.0 is the only SAM to operate such a distinction to the best of our knowledge).

In GALICS 2.0, the distinction between major and minor mergers is based on the mass ratio ℳ1/ℳ2, where ℳ1 and ℳ2 are the total masses of the merging galaxies within their respective baryonic half-mass radii (ℳ1 > ℳ2). Cattaneo et al. (2017) assumed a sharp transition at a critical mass ratio ℳ1/ℳ2 = ϵm = 4. In major mergers (ℳ1/ℳ2 < ϵm), all stars go into a classical bulge and all gas goes into the central cusp (Toomre & Toomre 1972; Barnes 1992; Negroponte & White 1983; Barnes & Hernquist 1991; Springel et al. 2005; but also see Hopkins et al. 2009 and Hammer et al. 2018). In minor mergers (ℳ1/ℳ2 > ϵm), the merging galaxies are added component by component.

Our assumptions for minor mergers are the same as in Cattaneo et al. (2017). Based on N-body simulations by Eliche-Moral et al. (2012), we assume thay the scale-lengths of the disc and the bulge of the larger galaxy are not affected by minor-merging events. The difference is in major mergers. The model by Cattaneo et al. (2017) had the shortcoming of predicting too many pure ellipticals. The median B/T for early-type galaxies with B/T >  0.7 was close to unity in the model while it is B/T ∼ 0.85 in the observations (the red symbols with error bars in Fig. 10, data from Mendel et al. 2014). We obviate to this shortcoming by assuming that the mass fraction ftr transferrred to the bulge (or the cusp in the case of gas) is not always unity but varies gradually from ftr = 0 for ℳ1/ℳ2 = ϵm to ftr = 0 for ℳ1/ℳ2 = 1. We assume:

(32)

thumbnail Fig. 10.

B/T versus M* for all galaxies (black), spiral galaxies (B/T <  0.7, blue), and elliptical galaxies (B/T >  0.7, red). The curves are medians in bins of stellar mass. They refer to GalICS 2.0 without (left) and with (right) disc instabilities. The data points are the observations by Mendel et al. (2014). The width of the shaded areas around the curves shows the standard error on the median.

where ϵm = 3. Using ϵm = 4 would increase B/T at 1010.5 M⊙ <  M* <  1011 M⊙ but also at M* <  1011 M⊙, where there is no need for such an increase (Fig. 10).

The reader should keep in mind that the choice of a linear relation in Eq. (32) is motivated by its simplicity rather than physical considerations. The study of the dependence of ftr on ℳ1/ℳ2 would be a major project in its own right and would require an extensive campaign of numerical simulations. The GalMer database of major and minor mergers experiments (Chilingarian et al. 2010) has been the most significant effort in this direction but does not answer our question because the questions that were important for those who have analysed the GalMer database (for example, Eliche-Moral et al. 2018 for a recent study) are not those that are important for our SAM.

Besides the model for major mergers, the only other differences with the version of GALICS 2.0 described in Koutsouridou & Cattaneo (2019) are the model for disc instabilities and a detail in the calculation of disc sizes. Disc radii are important because compact discs are less stable (Sect. 3). In GALICS 2.0, rd is determined by the spin parameter λ of the DM halo (Eq. (10)), which we measure from the N-body simulation used to construct the merger trees. The gaseous disc and the stellar disc are assumed to have the same exponential scale-length. If λ fluctuates from one timestep to the next, so does rd. The effect is small for central galaxies, but it can be significant for satellites and is particularly strong in the central regions of groups and clusters, where the halo finder6 may incorrectly assign to the host system particles that in reality belong to a subhalo. Koutsouridou & Cattaneo (2019) showed that this numerical effect causes the discs of satellite galaxies to shrink at pericentric passages. Discs regain their previous sizes after moving away from their orbital pericentres, but the temporary contraction affects their stability and drives the formation of spurious bulges. In this article, we prevent this phenomenon by requiring that the sizes of discs cannot decrease.

Our model for disc instability is based on Eq. (19); fd is the contribution of the disc to the total gravitational acceleration at r = 3.2rd. The other two contributions are those of the bulge and the halo. Let Md, Mcb, and Mbar be the masses of the disc, the classical bulge inclusive of the cusp, and the bar inclusive of the pseudobulge, respectively, where Mbar = 0 initially. If the B/T computed with Eq. (19) is B/T >  Mcb/(Mcb + Mbar), then a bar will form until:

(33)

To understand how fd is computed in the case of a galaxy that already has a bar, one should consider that, in our SAM, the bar is the disc’s inner part. The formation of a bar changes the azimuthal distribution of stars in the inner disc and causes the inner disc to buckle, but does not change the total contribution of the disc-plus-bar system to the circular velocity at r = 3.2rd. We therefore compute fd as if the mass in the bar belonged to the disc (in the beginning this mass will be zero because all galaxies start as pure discs). Mbar does not change the B/T computed with Eq. (19) but affects the right-hand side of Eq. (33). If the right-hand side is larger than the left-hand side, the disc is stable and there is no need for any further increases in Mbar. A classical bulge from a previous major merger lowers fd, makes the disc more stable and can prevent the formation of a bar.

The black, blue, and red curves in Fig. 10 show B/T versus stellar mass M* in GALICS 2.0 for all galaxies, galaxies with B/T <  0.7, and galaxies with B/T >  0.7, respectively. The shaded area around each curve shows the standard error on the median. The left panel refers to a model with mergers only. The right panel shows a model that contains both mergers and disc instabilities.

We now focus on the blue curves (B/T <  0.7) and compare them to the SDSS data by Mendel et al. (2014, blue symbols with error bars). The model without disc instabilities is in reasonably good agreement with the data points for M* >  1011 M⊙ but fails completely at lower masses, where mergers make a negligible contribution to the mass growth of galaxies (Cattaneo et al. 2011; also see Cattaneo et al. 2017).

Including disc instabilities improves the SAM considerably. The agreement with the observations by Mendel et al. (2014, data points with error bars) is reasonably good at all masses but especially at M* <  1010.5 M⊙. The transition from pseudobulges formed at low masses to classical bulges formed by mergers at high masses is not only a theoretical predictions. The data themselves show that the B/T − M* relation changes slope at M* ∼ 1010 M⊙.

We remark that a model could reproduce the mean B/T galaxies with B/T <  0.7 (blue symbols) and with B/T >  0.7 (red symbols) separately and still not reproduce the mean B/T of the total galaxy population (black symbols) if the fraction of galaxies with B/T >  0.7 is not reproduced correctly. The reasonably good agreement between the black curve and the black symbols suggests that the fraction of galaxies with B/T >  0.7 is reproduced reasonably well at both low (M* ≲ 1010 M⊙) and high (M* ≳ 1011 M⊙) masses. It is only in the intermediate mass-range M* ∼ 1010.5 M⊙ that mergers seem to struggle to form enough elliptical and S0 galaxies compared to what the observations suggest they should be doing.

We note that the catalogue by Mendel et al. (2014; 660 000 galaxies) is based on the one by Simard et al. (2011) used to produce the red squares in Fig. 8. The difference is the use of additional photometric bands, which allows a more accurate fitting of the spectral energy distributions and thus more accurate mass estimates, as well as a more detailed analysis of the effects of dust absorption. Comparing the red squares in Fig. 8 with the blue circles in Fig. 10 shows that the bulge-to-total mass ratios of faint galaxies in Fig. 8 are almost certainly overestimated and highlights the danger of B/T decompositions based on one photometric band only.

7. Discussion and conclusion

Our simulations of the stability of a thin disc embedded in a static spherical halo improve the classical work by ELN because they are three-dimensional and have higher resolution, but the underlying assumptions are very similar. The only substantial difference is the density distribution of the DM. Hence, it is unsurprising and reassuring that our results and those by ELN are in substantial agreement.

However, not only do our simulations confirm that all discs with vd/vc >  0.6 are unstable unless the threshold is exceeded by a narrow margin. They also take the work by ELN one step forwards and show that vd/vc can be used to predict a galaxy’s bulge-to-total mass ratio B/T. This is the main difference between our work and that of ELN. They used fd = (vd/vc)2 to decide whether an instability would develop. We use the same parameter to compute B/T, so that our answer is more than yes or no. A pseudobulge may form, but its mass may be very small. Therefore, it is completely unrealistic to assume as Cole et al. (2000) that the entire disc collapses into a bulge whenever ELN’s stability criterion is not satisfied.

Cattaneo et al. (2017) computed the sizes and masses of pseudobulges by solving the equation:

(34)

where α was treated as an input parameter of the SAM. Discs were bulgeless when the equation admitted no solution for r. This approach has no physical justification because the ELN criterion is a global one and the stability threshold depends on the radius at which it is applied. A disc with fd(2.2rd) < 2α is likely to be globally stable and should not form any pseudobulge, statistically at least. However, it is possible that fd(r) = 2α for r <  2.2rd, in which case the SAM of Cattaneo et al. (2017) would attribute to it a pseudobulge of radius r. Cattaneo et al. (2017) compensated this overpropensity of their model to form pseudobulges by using α = 0.45 instead of α = 0.28 (it is more difficult to form pseudobulges when the instability threshold is set to a higher level) and none of their key results (stellar mass function of galaxies, disc sizes, and the Tully-Fisher relation) are sensitive to their model for disc instabilities. Still, Eq. (34) should not be regarded as a physically correct model for the formation of pseudobulges. In contrast, using Eq. (19) would represent a major step forwards with respect to the prescriptions currently used in SAMs.

Having summarised the main results of our simulations, we now discuss a number of caveats. The first is the cell size. We have run some of our simulations at lower resolution and we did not see substantial differences that would alter our conclusions.

Another potential issue is the local stability of our initial conditions. The local stability of a stellar thin disc is determined by the Toomre (1964) parameter:

(35)

where σ is the radial one-dimensional velocity dispersion, Σ is the stellar surface density and

(36)

is the epicyclic frequency. The disc is stable for Q >  1. In our razor-thin discs:

(37)

(Sect. 2.1). By substituting the σ from Eq. (37) and Ω = Vc/R with Vc from Eq. (4) into Eq. (35), we find that our discs have Q = 0.3–0.4 over most of their surface and that Q >  1 only at the centre. Hence, our initial conditions do not satisfy Toomre’s stability criterion. When Q <  1, spiral shocks heat the disc and contribute to its global stability, but a disc with Q <  1 can also develop clumps (Hockney & Hohl 1969), which can sink and merge to the centre because of dynamical friction. Relaxation effects will also lead to a redistribution of mass that could increase the central density.

To explore how this may affect our results, we have run six simulations in which we have let our initial conditions relax before we let them evolve. The relaxation procedure is as follows. We compute the initial gravitational potential and let the disc evolve in this fixed potential for one orbital time at the disc’s half-mass radius. We recompute the gravitational potential and perform a second iteration, then a third one. After three orbital times we let the disc’s potential become live and we run the simulations as usual. Table 1 lists the simulations that we have rerun in this manner and the final difference in B/T with respect to the values presented in Sect. 3.

The difference for relaxed and unrelaxed initial conditions is fairly random: Δ(B/T) = (B/T)rel − (B/T)unrel is sometimes negative and sometimes positive. Statistically, there is a systematic trend towards lower B/T when relaxed initial conditions are used, but the difference is small: −6% on average. No observer can reliably distinguish between B/T = 0.2 and B/T = 0.18, and certainly no model can claim to predict B/T with such accuracy. Hence, the effect of using unrelaxed, locally unstable initial conditions has no bearing on our conclusions within the precision that we aim to achieve (we remind the reader that the fit in Eq. (19) is accurate within to 30%).

This conclusion appears at odds with recent simulations of the dynamical evolution of cold stellar discs. Saha & Cortesi (2018) studied how such discs evolve for different stellar velocity dispersions when all other parameters are kept the same. Discs with higher σ and Q ∼ 1 formed a bar. They did not fragment into clumps. Discs with very low σ and Q <  0.5 fragmented in less than 1 Gyr. The stellar clumps migrated to the centre via dynamical friction and there coalesced into a bulge that reached B/T ≃ 0.35 at the end of the simulations. We shall argue that the simulations by Saha & Cortesi (2018) are more physical than ours, but also that their finding does not invalidate our results.

The key point is how dynamical friction works. Stellar clumps sink to the centre because they move with respect to the DM and dynamical friction transfers their angular momentum to the halo. We assume a static halo. Hence, this process cannot occur in our simulations. The disc can fragment into clumps, but the clumps do not migrate to the centre. One can see this as a limitation of our simulations. If some stellar discs really started with Q ≪ 1, their evolution into lenticular galaxies through the process described by Saha & Cortesi (2018) will not captured by our analysis. One can also see it as a justification for our assumption of a static spherical halo. We find that discs that start with Q ∼ 0.5 reach Q ∼ 0.8–1.1 in just one rotation time (measured at the half-mass radius). By assuming a static halo, we avoid that our final results are too contaminated by relaxation effects due to the Toomre instability. The latter view makes more sense if we consider that our initial conditions are artificial and that stellar discs did not start with Q ≪ 1.

To check that this explanation is correct, we have rerun four simulations with a live halo (Table 2). Orbital decay and coalescence of clumps through dynamical friction increase B/T by 12% on average. The effect is strongest in galaxies where clumps are clearly visible (those with md = 0.04 and λ = 0.1 in Figs. 2 and 3).

An effect of 10–20% (Table 2) is ultimately small. Saha & Cortesi (2018) did not report the final B/T for the simulations with Q >  0.5. Hence, we cannot make a quantitative comparison with their results. However, a glance at their images shows that the simulations with Q >  0.5 display a significant bulge component, too. Hence, there is no reason to suppose that our final B/T are inconsistent with theirs.

That is not to say that the final visual morphologies look identical in the simulations with Q <  0.5 and the simulation with Q ∼ 1 (Fig. 2 of Saha & Cortesi 2018). The simulations with Q <  0.5 form S0 galaxies with a prominent bulge but no bar or spiral structure. In the simulation with Q ∼ 1, the final galaxy is clearly barred. This is a case where a two-component bulge/disc fit and a fit with an additional bar component may give substantially different B/T and show a much stronger dependence on Q. We do not perform such an analysis, which would be more accurate in reality, because it would then be difficult to compare our results to SDSS studies based on a two-component analysis (for example, Mendel et al. 2014).

Finally, to understand the extent to which 2% gas could affect our results, we have run four simulations in which the gas fraction is exactly zero (Table 3) and we have analysed them at t = 1,1.5, 2, 2.5 and 3 Gyr. The mean variation in B/T without gas was ( − 2.7 ± 1.8)% with minimum and maximum variations with respect to the case with gas of −9.3% and +3.2%, respectively. These variations show that gas tends to concentrate in the inner regions and thus to increase B/T, but the effect is negligible for the 2% gas fraction assumed in this article, especially since the scatter in B/T for a given fd(3.2rd) is much larger (∼30%).

Based on the numerical tests above, we are confident that Eq. (19) provides the correct answer to our problem within the approximations that we have made. The adequacy of these approximations is another problem, the one that we now discuss.

Already ten years ago, Athanassoula (2008) had run more realistic simulations with a live halo and she found that supposedly stable discs (according to the ELN criterion) could still form a pseudobulge by transferring angular momentum to the DM through resonances. Conversely, supposedly unstable discs could still be stabilised by random motions in the disc or the halo or by a central DM concentration (the effects of the latter have been considered in this article). She concluded that the ELN criterion is too simplistic to describe a complex process such as the formation of pseudobulges.

Fujii et al. (2018) ran simulations with a live halo at extremely high resolution and they confirmed cases (6 out of 18) where the ELN criterion failed to predict the formation of a bar, although, in all those cases, exceeded the threshold value ϵ = 1.1 by only 30% on average. Their result is not inconsistent with ours (we, too, find a pseudobulge in a galaxy with vd/vc <  0.5 at r = 2.2rd, which should not have one according to the ELN criterion) and it does not surprise us given the additional possibilities (such as resonances) that come into play when a live halo is considered. Our simulations with a live halo (Table 2) confirms that it can make a difference, but also that the bulge-to-total mass ratios measured with a two-component fit are not too affected by it.

The comparison with Fujii et al. (2018) is interesting also in relation to two specific claims they make. Fujii et al. claim that all discs will develop a bar if one waits long enough and that what fd(2.2rd) really measures is the timescale tb for bar formation. The discs with are simply those with tb >  10 Gyr. They also claim that the shear rate:

(38)

is the physical quantity that is most relevant to determining B/T (Ω = vc/r; also see D’Onghia 2015), and demonstrate a tight correlation between B/T and Γ(2.2rd).

We cannot comment on the long-term evolution of our galaxies because we have not run any simulations for more than 3 Gyr and since our resolution is lower than that of Fujii et al. (2018) anyway (the originality of our study is rather in our systematic exploration of the parameter space). However, we have verifed that our discs evolved very little in the last gigayear before we stopped following them, in agreement with Debattista et al. (2017) and Martinez-Valpuesta et al. (2017), who find that the bar strength grows slowly after the initial rapid growth and the formation of a pseudobulge in the first 2 Gyr. Furthermore, any result for isolated discs is purely academical when used to discuss the evolution of galaxies on timescales comparable to the Hubble time7.

In relation to the second claim, our simulations confirm that B/T and Γ(2.2rd) are tightly related (Fig. 11). However, we do not consider that Γ adds substantial information with respect to fd and the reason is simple. Γ(r) is directly linked to the form of the rotation curve (Γ <  1 for a rising rotation curve; Γ >  1 for a declining one), which is the sum in quadrature of vd and vh. Substituting Eqs. (4) into (38) gives:

(39)

thumbnail Fig. 11.

B/T versus the exponent Γ of Ω(r) at r = 2.2rd. Symbols of the same shape and colour correspond to the same galaxies as in Fig. 7. The black line is a log-log linear least-square fit to the symbols with B/T >  0.01.

where:

(40)

and:

(41)

At r = 2.2rd, vd has a maximum, so that βd(2.2rd) = 0; βh(2.2rd) depends on crd, but, for reasonable values of c and rd, its value is always βh(2.2rd) ≃ 0.4 (±0.1 at most). In contrast, fd can vary by an order of magnitude from one galaxy to another. Hence, it is fd that determines Γ, as we intended to demonstrate.

Studies such as those by Athanassoula (2008) and Fujii et al. (2018) show that the formation of a bar is sensitive to the instability threshold α and that α is not universal but depends on parameters that our crude assumption of a static spherical halo cannot capture. However, predicting the stability of individual galaxies with the highest possible accuracy is not the goal of our article. Our goal is to develop a model that may be able to make statistical predictions for the most likely B/T ratios of galaxies given the halo masses, spins, concentrations, and merging histories measured in cosmological N-body simulations.

Considering the quantitative dependence of B/T on vd/vc helps to overcome the problem highlighted by the above studies. We may erroneously predict the formation of a pseudobulge in a galaxy with vd/vc >  α that should not have formed one in reality. However, if vd/vc is small, B/T will be small, and galaxies with B/T <  0.1 are, for many practical purposes, indistinguishable from galaxies with B/T = 0, particularly when one considers the difficulty of obtaining accurate quantitative morphologies from observations at high redshift. In the same way, our criterion may fail to predict the formation of a pseudobulge in a galaxy with vd/vc <  α that should develop one. However, the bulges that we miss are usually small. The clearest example is galaxy md0.3mb1 of Fujii et al. (2018). For that galaxy, the left-hand side of Eq. (1) is equal to 1.87, so the galaxy should be stable. The simulations by Fujii et al. (2018) show that it is not. However, visual inspection of the morphology after 5 Gyr does not show any massive bulge.

The assumptions under which we have derived Eq. (19) are consistent with the spirit of the semi-analytic approximation, which is to separate the evolution of the DM from that of the baryons within haloes. Cosmological hydrodynamic simulations have shown that baryonic physics can affect the structural properties of DM haloes (Pontzen & Governato 2012; Teyssier et al. 2013; Di Cintio et al. 2014; Tollet et al. 2016) and cause them to deviate from the NFW profile that SAMs assume, when they do not make the even cruder approximation of a singular isothermal sphere. Using Eq. (19) to assign morphologies to galaxies that have not experienced any mergers is no greater approximation than to model the DM with an NFW profile or to assume that major mergers transform discs into spheroids instantaneously. In fact, Robertson et al. (2006), Hopkins et al. (2009), and Hammer et al. (2009) have shown examples of discs that survived gas-rich major mergers. Yet, a lot has been learned about the Universe from this simple picture.

Finally, despite their simplifying assumptions, our simulations reproduce the magnitude – morphology relation observed in the catalogues from Mathewson et al. (1992), Simard et al. (2011), and Salo et al. (2015). SAMs based on our findings are in very good agreement with the B/T – M* relation (Mendel et al. 2014), especially at M* ≲ 1010 M⊙, where most bulges are pseudobulges formed by disc instabilities in both our SAM and observations. Our results support models in which pseudobulges grow gradually through non-axisymmetric secular disc instabilities as opposed to violent instabilities where discs collapse rapidly into bulges (à la Cole et al. 2000).

We conclude with a note on our future plans. “Pseudobulge formation through bar instabilities” might have been a more precise title for our current article, but we have chosen ‘Bulge formation through disc instabilities’ because we are planning a second article, in which explore how our results can be extended to gas-rich systems. We have seen that a 2% gas fraction increases the final B/T by ∼3% on average (Sect. 4). For high gas fractions, the instability may be not only quantitatively but also qualitatively different. Already Noguchi & Shlosman (1994) found that unstable gas-rich discs do not develop a bar. They fragment into clumps that merge into a central bulge. Bulges formed through this process may be closer to classical bulges than to peanut-shaped pseudobulges (Bournaud et al. 2007). Hence our choice of a title that is broad enough to encompass our future projects.

1

Erwin (2018) compared the frequency of bars in the S4G and the SDSS. He finds that SDSS-based studies underestimate the fraction of barred galaxies at low masses because of poor spatial resolution and the correlation between bar size and stellar mass.

2

Not all pseudobulges are peanut-shaped. Some are nuclear spirals and they are as flat as discs (Carollo et al. 1998; Kormendy & Kennicutt 2004). However, the distribution of Sérsic index, ellipticity and B/T is similar for pseudobulges in barred, oval, and normal galaxies, suggesting that the difference between classical bulges and pseudobulges is more important than the one between different types of pseudobulges (Fisher & Drory 2008). Moreover, Athanassoula (2005) has shown that disc-like bulges result from the dissipational inflow of gas into the central regions of spiral galaxies. Hence, it is reasonable for us to neglect them in an article that is about stellar discs.

3

We acknowledge that a bar and a bulge are different dynamical entities, and that it is possible to separate in a morphological decomposition (Gadotti 2009; Weinzirl et al. 2009; Salo et al. 2015; Erwin 2018, 2019). We merely state that SAMs are not an appropriate tool for such decomposition.

4

The factor of two in Eq. (14) has been introduced so that our definition of α coincides with the one in Christodoulou et al. (1995).

6

GALICS 2.0 uses a halo finder that is called HALOMAKER (Tweed et al. 2009) and is based on ADAPTAHOP (Aubert et al. 2004).

7

Even in the absence of mergers, discs still accrete gas from their environment. Slow accretion lowers the B/T ratio and could compensate the growth of pseudobulges in galaxies with tb ≳ 10 Gyr. Rapid gas accretion, whereby the gas fraction becomes significant, is beyond the scope of this article, which is on stellar discs.

Acknowledgments

This work constitutes the eight-week undergraduate thesis of T. D. Our simulations were run on the OCCIGEN supercomputer of the Centre Informatique National de l’Enseignement Supérieur and on the HPC resources of MesoPSL financed by the Région Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the Agence Nationale pour la Recherche. We also thank the anonymous referee for many remarks that have strengthened our manuscript.

References

  1. Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 597, 21 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  2. Anglés-Alcázar, D., Faucher-Giguère, C.-A., Kereš, D., et al. 2017, MNRAS, 470, 4698 [NASA ADS] [CrossRef] [Google Scholar]
  3. Athanassoula, E. 2003, MNRAS, 341, 1179 [NASA ADS] [CrossRef] [Google Scholar]
  4. Athanassoula, E. 2005, MNRAS, 358, 1477 [NASA ADS] [CrossRef] [Google Scholar]
  5. Athanassoula, E. 2008, MNRAS, 390, L69 [NASA ADS] [CrossRef] [Google Scholar]
  6. Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35 [NASA ADS] [CrossRef] [Google Scholar]
  7. Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barnes, J. E. 1992, ApJ, 393, 484 [NASA ADS] [CrossRef] [Google Scholar]
  9. Barnes, J. E., & Hernquist, L. E. 1991, ApJ, 370, L65 [NASA ADS] [CrossRef] [Google Scholar]
  10. Baugh, C. M., Cole, S., & Frenk, C. S. 1996, MNRAS, 283, 1361 [NASA ADS] [CrossRef] [Google Scholar]
  11. Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]
  12. Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143 [NASA ADS] [CrossRef] [Google Scholar]
  13. Benson, A. J. 2012, New Astronomy, 17, 175 [NASA ADS] [CrossRef] [Google Scholar]
  14. Berentzen, I., Heller, C. H., Shlosman, I., & Fricke, K. J. 1998, MNRAS, 300, 49 [NASA ADS] [CrossRef] [Google Scholar]
  15. Binney, J., & Tremaine, S. 2008, in Galactic Dynamics: Second Edition, eds. J. Binney, & S. Tremaine (Princeton, NJ USA: PrincetonUniversity Press) [Google Scholar]
  16. Bissantz, N., & Gerhard, O. 2002, MNRAS, 330, 591 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Boselli, A., Cortese, L., Boquien, M., et al. 2014, A&A, 564, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Bournaud, F., Elmegreen, B. G., & Elmegreen, D. M. 2007, ApJ, 670, 237 [Google Scholar]
  20. Bournaud, F., Elmegreen, B. G., & Martig, M. 2009, ApJ, 707, L1 [NASA ADS] [CrossRef] [Google Scholar]
  21. Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115 [NASA ADS] [CrossRef] [Google Scholar]
  22. Brook, C. B., Kawata, D., Gibson, B. K., & Freeman, K. C. 2004, ApJ, 612, 894 [NASA ADS] [CrossRef] [Google Scholar]
  23. Brook, C. B., Stinson, G., Gibson, B. K., Wadsley, J., & Quinn, T. 2012, MNRAS, 424, 1275 [NASA ADS] [CrossRef] [Google Scholar]
  24. Bryan, S. E., Kay, S. T., Duffy, A. R., et al. 2013, MNRAS, 429, 3316 [NASA ADS] [CrossRef] [Google Scholar]
  25. Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240 [NASA ADS] [CrossRef] [Google Scholar]
  26. Burkert, A., Förster Schreiber, N. M., Genzel, R., et al. 2016, ApJ, 826, 214 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cacciato, M., Dekel, A., & Genel, S. 2012, MNRAS, 421, 818 [NASA ADS] [Google Scholar]
  28. Carollo, C. M., Stiavelli, M., & Mack, J. 1998, AJ, 116, 68 [NASA ADS] [CrossRef] [Google Scholar]
  29. Cattaneo, A., Mamon, G. A., Warnick, K., & Knebe, A. 2011, A&A, 533, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Cattaneo, A., Blaizot, J., Devriendt, J. E. G., et al. 2017, MNRAS, 471, 1401 [NASA ADS] [CrossRef] [Google Scholar]
  31. Chilingarian, I. V., Di Matteo, P., Combes, F., Melchior, A. L., & Semelin, B. 2010, A&A, 518, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Christodoulou, D. M., Shlosman, I., & Tohline, J. E. 1995, ApJ, 443, 551 [NASA ADS] [CrossRef] [Google Scholar]
  33. Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, MNRAS, 319, 168 [NASA ADS] [CrossRef] [Google Scholar]
  34. Combes, F., & Sanders, R. H. 1981, A&A, 96, 164 [NASA ADS] [Google Scholar]
  35. Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A&A, 233, 82 [NASA ADS] [Google Scholar]
  36. Combes, F., García-Burillo, S., Braine, J., et al. 2013, A&A, 550, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Comerón, S., Elmegreen, B. G., Knapen, J. H., et al. 2011, ApJ, 741, 28 [NASA ADS] [CrossRef] [Google Scholar]
  38. Croton, D. J., Stevens, A. R. H., Tonini, C., et al. 2016, ApJS, 222, 22 [Google Scholar]
  39. Debattista, V. P., Ness, M., Gonzalez, O. A., et al. 2017, MNRAS, 469, 1587 [NASA ADS] [CrossRef] [Google Scholar]
  40. Dekel, A., & Silk, J. 1986, ApJ, 303, 39 [NASA ADS] [CrossRef] [Google Scholar]
  41. de Vaucouleurs, G., de Vaucouleurs, A., & Corwin, J. R. 1976, Second Reference Catalogue of Bright Galaxies (Austin: University of Texas Press), 1976 [Google Scholar]
  42. Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415 [Google Scholar]
  43. Dimauro, P., Huertas-Company, M., Daddi, E., et al. 2018, MNRAS, 478, 5410 [NASA ADS] [CrossRef] [Google Scholar]
  44. D’Onghia, E. 2015, ApJ, 808, L8 [NASA ADS] [CrossRef] [Google Scholar]
  45. Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]
  46. Efstathiou, G. 2000, MNRAS, 317, 697 [NASA ADS] [CrossRef] [Google Scholar]
  47. Efstathiou, G., Lake, G., & Negroponte, J. 1982, MNRAS, 199, 1069 [NASA ADS] [CrossRef] [Google Scholar]
  48. Eliche-Moral, M. C., González-García, A. C., Aguerri, J. A. L., et al. 2012, A&A, 547, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Eliche-Moral, M. C., Rodríguez-Pérez, C., Borlaff, A., Querejeta, M., & Tapia, T. 2018, A&A, 617, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Erwin, P. 2018, MNRAS, 474, 5372 [NASA ADS] [CrossRef] [Google Scholar]
  51. Erwin, P. 2019, MNRAS, 489, 3553 [CrossRef] [Google Scholar]
  52. Fall, S. M., & Efstathiou, G. 1980, MNRAS, 193, 189 [NASA ADS] [CrossRef] [Google Scholar]
  53. Fisher, D. B., & Drory, N. 2008, AJ, 136, 773 [NASA ADS] [CrossRef] [Google Scholar]
  54. Fisher, D. B., & Drory, N. 2011, ApJ, 733, L47 [NASA ADS] [CrossRef] [Google Scholar]
  55. Freeman, K. C. 1970, ApJ, 160, 811 [Google Scholar]
  56. Fujii, M. S., Baba, J., Saitoh, T. R., et al. 2011, ApJ, 730, 109 [NASA ADS] [CrossRef] [Google Scholar]
  57. Fujii, M. S., Bédorf, J., Baba, J., & Portegies Zwart, S. 2018, MNRAS, 477, 1451 [NASA ADS] [CrossRef] [Google Scholar]
  58. Gadotti, D. A. 2009, MNRAS, 393, 1531 [NASA ADS] [CrossRef] [Google Scholar]
  59. Gargiulo, I. D., Cora, S. A., Padilla, N. D., et al. 2015, MNRAS, 446, 3820 [NASA ADS] [CrossRef] [Google Scholar]
  60. Gonzalez-Perez, V., Lacey, C. G., Baugh, C. M., et al. 2014, MNRAS, 439, 264 [NASA ADS] [CrossRef] [Google Scholar]
  61. Graham, A. W., & Worley, C. C. 2008, MNRAS, 388, 1708 [NASA ADS] [CrossRef] [Google Scholar]
  62. Hammer, F., Flores, H., Puech, M., et al. 2009, A&A, 507, 1313 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Hammer, F., Yang, Y. B., Wang, J. L., et al. 2018, MNRAS, 475, 2754 [NASA ADS] [Google Scholar]
  64. Hatton, S., Devriendt, J. E. G., Ninin, S., et al. 2003, MNRAS, 343, 75 [NASA ADS] [CrossRef] [Google Scholar]
  65. Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hockney, R. W., & Hohl, F. 1969, AJ, 74, 1102 [CrossRef] [Google Scholar]
  67. Hohl, F. 1971, ApJ, 168, 343 [NASA ADS] [CrossRef] [Google Scholar]
  68. Hopkins, P. F., Cox, T. J., Younger, J. D., & Hernquist, L. 2009, ApJ, 691, 1168 [Google Scholar]
  69. Hubble, E. P. 1926, ApJ, 64, 321 [NASA ADS] [CrossRef] [Google Scholar]
  70. Jiang, F., Dekel, A., Kneller, O., et al. 2019, MNRAS, 488, 4801 [CrossRef] [Google Scholar]
  71. Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  72. Kalnajs, A. J. 1972, ApJ, 175, 63 [NASA ADS] [CrossRef] [Google Scholar]
  73. Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201 [Google Scholar]
  74. Kimm, T., Devriendt, J., Slyz, A., et al. 2011, ArXiv e-prints [arXiv:1106.0538] [Google Scholar]
  75. Kormendy, J. 1982, ApJ, 257, 75 [NASA ADS] [CrossRef] [Google Scholar]
  76. Kormendy, J. 1993, in Galactic Bulges, eds. H. Dejonghe, & H. J. Habing, IAU Symp., 153, 209 [NASA ADS] [CrossRef] [Google Scholar]
  77. Kormendy, J., & Kennicutt, R. C. 2004, ARA&A, 42, 603 [NASA ADS] [CrossRef] [Google Scholar]
  78. Kormendy, J., Fisher, D. B., Cornell, M. E., & Bender, R. 2009, ApJS, 182, 216 [NASA ADS] [CrossRef] [Google Scholar]
  79. Kormendy, J., Drory, N., Bender, R., & Cornell, M. E. 2010, ApJ, 723, 54 [NASA ADS] [CrossRef] [Google Scholar]
  80. Koutsouridou, I., & Cattaneo, A. 2019, MNRAS, 490, 5375 [CrossRef] [Google Scholar]
  81. Leauthaud, A., Tinker, J., Bundy, K., et al. 2012, ApJ, 744, 159 [NASA ADS] [CrossRef] [Google Scholar]
  82. Lee, J., & Yi, S. K. 2013, ApJ, 766, 38 [NASA ADS] [CrossRef] [Google Scholar]
  83. Libeskind, N. I., Yepes, G., Knebe, A., et al. 2010, MNRAS, 401, 1889 [NASA ADS] [CrossRef] [Google Scholar]
  84. Lo Faro, B., Monaco, P., Vanzella, E., et al. 2009, MNRAS, 399, 827 [NASA ADS] [CrossRef] [Google Scholar]
  85. Macciò, A. V., Stinson, G., Brook, C. B., et al. 2012, ApJ, 744, L9 [NASA ADS] [CrossRef] [Google Scholar]
  86. Martinez-Valpuesta, I., Aguerri, J. A. L., González-García, A. C., Dalla Vecchia, C., & Stringer, M. 2017, MNRAS, 464, 1502 [NASA ADS] [CrossRef] [Google Scholar]
  87. Mathewson, D. S., Ford, V. L., & Buchhorn, M. 1992, ApJS, 81, 413 [NASA ADS] [CrossRef] [Google Scholar]
  88. Meert, A., Vikram, V., & Bernardi, M. 2015, MNRAS, 446, 3943 [NASA ADS] [CrossRef] [Google Scholar]
  89. Meert, A., Vikram, V., & Bernardi, M. 2016, MNRAS, 455, 2440 [CrossRef] [Google Scholar]
  90. Mendel, J. T., Simard, L., Palmer, M., Ellison, S. L., & Patton, D. R. 2014, ApJS, 210, 3 [NASA ADS] [CrossRef] [Google Scholar]
  91. Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319 [Google Scholar]
  92. Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]
  93. Moster, B. P., Naab, T., & White, S. D. M. 2018, MNRAS, 477, 1822 [NASA ADS] [CrossRef] [Google Scholar]
  94. Muñoz-Cuartas, J. C., Macciò, A. V., Gottlöber, S., & Dutton, A. A. 2011, MNRAS, 411, 584 [NASA ADS] [CrossRef] [Google Scholar]
  95. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
  96. Negroponte, J., & White, S. D. M. 1983, MNRAS, 205, 1009 [Google Scholar]
  97. Noguchi, M., & Shlosman, I. 1994, Ap&SS, 216, 347 [CrossRef] [Google Scholar]
  98. Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467 [NASA ADS] [CrossRef] [Google Scholar]
  99. Papastergis, E., Cattaneo, A., Huang, S., Giovanelli, R., & Haynes, M. P. 2012, ApJ, 759, 138 [NASA ADS] [CrossRef] [Google Scholar]
  100. Persic, M., & Salucci, P. 1995, ApJS, 99, 501 [NASA ADS] [CrossRef] [Google Scholar]
  101. Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27 [Google Scholar]
  102. Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75 [NASA ADS] [Google Scholar]
  103. Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 [NASA ADS] [CrossRef] [Google Scholar]
  104. Porter, L. A., Somerville, R. S., Primack, J. R., & Johansson, P. H. 2014, MNRAS, 444, 942 [NASA ADS] [CrossRef] [Google Scholar]
  105. Quillen, A. C., Kuchinski, L. E., Frogel, J. A., & DePoy, D. L. 1997, ApJ, 481, 179 [NASA ADS] [CrossRef] [Google Scholar]
  106. Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74 [NASA ADS] [CrossRef] [Google Scholar]
  107. Reyes, R., Mandelbaum, R., Gunn, J. E., et al. 2012, MNRAS, 425, 2610 [NASA ADS] [CrossRef] [Google Scholar]
  108. Robertson, B., Bullock, J. S., Cox, T. J., et al. 2006, ApJ, 645, 986 [NASA ADS] [CrossRef] [Google Scholar]
  109. Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Saha, K., & Cortesi, A. 2018, ApJ, 862, L12 [CrossRef] [Google Scholar]
  111. Salo, H., Laurikainen, E., Laine, J., et al. 2015, ApJS, 219, 4 [NASA ADS] [CrossRef] [Google Scholar]
  112. Schönrich, R., & Binney, J. 2009, MNRAS, 399, 1145 [NASA ADS] [CrossRef] [Google Scholar]
  113. Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41 [Google Scholar]
  114. Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978 [NASA ADS] [CrossRef] [Google Scholar]
  115. Sheth, K., Regan, M., Hinz, J. L., et al. 2010, PASP, 122, 1397 [NASA ADS] [CrossRef] [Google Scholar]
  116. Silk, J., & Rees, M. J. 1998, A&A, 331, L1 [NASA ADS] [Google Scholar]
  117. Simard, L., Mendel, J. T., Patton, D. R., Ellison, S. L., & McConnachie, A. W. 2011, ApJS, 196, 11 [NASA ADS] [CrossRef] [Google Scholar]
  118. Springel, V., Di Matteo, T., & Hernquist, L. 2005, ApJ, 620, L79 [NASA ADS] [CrossRef] [Google Scholar]
  119. Steinmetz, M. 2012, Astron. Nachr., 333, 523 [NASA ADS] [CrossRef] [Google Scholar]
  120. Stewart, K. R., Brooks, A. M., Bullock, J. S., et al. 2013, ApJ, 769, 74 [NASA ADS] [CrossRef] [Google Scholar]
  121. Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013, MNRAS, 429, 3068 [NASA ADS] [CrossRef] [Google Scholar]
  123. Tollet, E., Macciò, A. V., Dutton, A. A., et al. 2016, MNRAS, 456, 3542 [NASA ADS] [CrossRef] [Google Scholar]
  124. Tollet, É., Cattaneo, A., Mamon, G. A., Moutard, T., & van den Bosch, F. C. 2017, MNRAS, 471, 4170 [NASA ADS] [CrossRef] [Google Scholar]
  125. Tollet, É., Cattaneo, A., Macciò, A. V., Dutton, A. A., & Kang, X. 2019, MNRAS, 485, 2511 [CrossRef] [Google Scholar]
  126. Toomre, A. 1963, ApJ, 138, 385 [NASA ADS] [CrossRef] [Google Scholar]
  127. Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
  128. Toomre, A., & Toomre, J. 1972, ApJ, 178, 623 [Google Scholar]
  129. Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. van Daalen, M. P., Schaye, J., McCarthy, I. G., Booth, C. M., & Dalla Vecchia, C. 2014, MNRAS, 440, 2997 [CrossRef] [Google Scholar]
  131. van den Bosch, F. C. 1998, ApJ, 507, 601 [NASA ADS] [CrossRef] [Google Scholar]
  132. van den Bosch, F. C. 2000, ApJ, 530, 177 [NASA ADS] [CrossRef] [Google Scholar]
  133. Villalobos, Á., & Helmi, A. 2008, MNRAS, 391, 1806 [NASA ADS] [CrossRef] [Google Scholar]
  134. Villumsen, J. V. 1985, ApJ, 290, 75 [NASA ADS] [CrossRef] [Google Scholar]
  135. Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050 [NASA ADS] [CrossRef] [Google Scholar]
  136. Weinzirl, T., Jogee, S., Khochfar, S., Burkert, A., & Kormendy, J. 2009, ApJ, 696, 411 [NASA ADS] [CrossRef] [Google Scholar]
  137. Wojtak, R., & Mamon, G. A. 2013, MNRAS, 428, 2407 [NASA ADS] [CrossRef] [Google Scholar]
  138. Yoachim, P., & Dalcanton, J. J. 2006, AJ, 131, 226 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Simulations with relaxed initial conditions: parameter values and difference in B/T with respect to the unrelaxed simulations.

In the text
Table 2.

Simulations with a live halo: parameter values and difference in B/T with respect to those with a static halo.

In the text
Table 3.

Simulations without gas: parameter values and difference in B/T with respect to those with fgas = 0.02.

In the text
Table 4.

Morphological composition of the catalogue in Mathewson et al. (1992).

In the text

All Figures

thumbnail Fig. 1.

Initial conditions and refinement regions. The black, blue, violet, green, yellow, and red curves show the isodensity contours that contain 90%, 80%, 70%, 60%, 50%, and 40% of the disc mass at t = 0, respectively. The black, blue, violet, green, yellow, and red dashed lines show the cylinders within which the cell size equals 1/8, 1/16, 1/32, 1/64, 1/128, and 1/256 of the disc exponential scale length, respectively.
In the text
thumbnail Fig. 2.

Face-on view of the simulated galaxies at the first time t when B/T has stopped growing (2 Gyr <  t <  3 Gyr). The image brightness is proportional to the logarithm of stellar surface density. This figure show the simulations with c = 5, sorted by md = Md/Mvir (lines) and halo spin λ (columns). On each image, we have shown the corresponding value of rd = Rd/Rvir. Bulge-to-total stellar mass ratios B/T have been computed by decomposing the stellar surface-density distribution into an exponential and a Sérsic (1963) profile. Galaxies that could be fitted with a single exponential profile have B/T = 0. When the fit required a bulge component, its Sérsic index n has been shown.
In the text
thumbnail Fig. 3.

Same as Fig. 2 with c = 10 instead of c = 5. To avoid overcrowding the figure, we have not shown the simulations with λ = 0.018 and λ = 0.35.
In the text
thumbnail Fig. 4.

Same as Fig. 2 with c = 15 instead of c = 5. The column λ = 0.1 is now missing because we know that simulations with λ = 0.1 will not form any (pseudo)bulge.
In the text
thumbnail Fig. 5.

Galaxy in the simulation with c = 10, md = 0.04, λ = 0.025, viewed edge-on. An X-shaped pseudobulge is clearly visible. This galaxy has been chosen as an example and is by no means atypical. The face-on image (Fig. 3) shows that this galaxy has a bar, but not a prominent one.
In the text
thumbnail Fig. 6.

Stellar surface-density for two of our simulated galaxies (black symbols). Galaxy a has been fitted with the sum of an exponential profile (blue line) and a Sérsic profile, with n = 1 in this particular case (red line). Its bulge-to-total mass ratio is B/T = 0.22. The black curve shows the sum of the two components. Galaxy b is consistent with a single exponential profile (blue line).
In the text
thumbnail Fig. 7.

Relation between bulge-to-total mass ratio B/T and initial disc fraction fd(r) for r = 2.2rd (left) and r = 3.2rd (right). Triangles pointing up, circles, and triangles pointing down correspond to simulations with c = 5, c = 10, and c = 15, respectively. The colours correspond to different disc-to-virial mass ratios: md = 0.005 (green), md = 0.01 (black), md = 0.02 (blue), and md = 0.04 (red). Symbols with B/T = 0.01 correspond to bulgeless galaxies. They have been assigned B/T = 0.01 merely to be able to show them on a logarithmic diagram. The thick solid black lines are log-log linear least-squares fits to the symbols with B/T >  0.01. They correspond to Eqs. (18) and (19). The vertical blue lines show the critical fd that corresponds to the α = 0.28 (equivalent to assuming ϵ = 1.1 in Eq. (1)) The coloured curves compare our results to previous models (Cole et al. 2000; Hatton et al. 2003; Shen et al. 2003; Cacciato et al. 2012).
In the text
thumbnail Fig. 8.

⟨B/T⟩ versus I-band magnitude MI. The black curve shows our prediction. The black circles are from Mathewson et al. (1992) when we assign to each Hubble type the ⟨B/T⟩ from Graham & Worley (2008; the uncertainty of estimating B/T determines the size of the error bars). The red squares are from Simard et al. (2011) after removing all galaxies with B/T >  0.7. The up-pointing filled blue triangles and the down-pointing empty blue triangles are from Salo et al. (2015), S4G) when we assign the bar to the bulge and the disc, respectively. The error bars on the red squares and the blue filled triangles show the standard error on the mean. We have not shown the error bars on the empty triangles so that they do not overlap with those on the filled triangles, but their width is comparable.
In the text
thumbnail Fig. 9.

Composite image of SDSS J141007.94-002348.4 from skyserver.sdss.org. SDSS J141007.94-002348.4 is a spiral galaxy with MI = −21.8. Simard et al. (2011) assign to it B/T = 0.93 but a bulge Sérsic index of only n = 0.5.
In the text
thumbnail Fig. 10.

B/T versus M* for all galaxies (black), spiral galaxies (B/T <  0.7, blue), and elliptical galaxies (B/T >  0.7, red). The curves are medians in bins of stellar mass. They refer to GalICS 2.0 without (left) and with (right) disc instabilities. The data points are the observations by Mendel et al. (2014). The width of the shaded areas around the curves shows the standard error on the median.
In the text
thumbnail Fig. 11.

B/T versus the exponent Γ of Ω(r) at r = 2.2rd. Symbols of the same shape and colour correspond to the same galaxies as in Fig. 7. The black line is a log-log linear least-square fit to the symbols with B/T >  0.01.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.