Bulge formation through disc instability -- I

We use simulations to study the growth of a pseudobulge in an isolated thin exponential stellar disc embedded in a static spherical halo. We observe a transition from later to earlier morphological types and an increase in bar prominence for higher disc-to-halo mass ratios, for lower disc-to-halo size ratios, and for lower halo concentrations. We compute bulge-to-total stellar mass ratios $B/T$ by fitting a two-component S\'ersic-exponential surface-density distribution. The final $B/T$ is strongly related to the disc's fractional contribution $f_{\rm d}$ to the total gravitational acceleration at the optical radius. The formula $B/T=0.5\,f_{\rm }^{1.8}$ fits the simulations to an accuracy of $30\%$, is consistent with observational measurements of B/T and f_d as a function of luminosity, and reproduces the observed relation between $B/T$ and stellar mass when incorporated into the GalICS~2.0 semi-analytic model of galaxy formation.


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. 2015Meert et al. , 2016Dimauro 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 (2018Erwin ( , 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 (S 4 G; 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 resolution 1 .
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 > ∼ 10 11 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 < 10 11 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 = 10 9 M to ∼ 100% at M = 10 10.7 M .
Most of the bulges in spirals with M < 10 11 M are pseudobulges with different kinematics than elliptical galaxies and do not follow the fundamental plane (Kormendy 1982(Kormendy , 1993Kormendy & 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).
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 M d to the halo mass M vir . 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 scaleheight 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 andGargiulo 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 M d decreases until the disc becomes stable again. In discs where V c (2.2R d ) 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, Frenk, & White (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 Section 2, we describe our initial conditions, which are entirely specified by three parameters (the ratio r d = R d /R vir of the disc scale-length R d to the virial radius R vir , the ratio m d = M d /M vir of the disc mass M d to the total mass M vir 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 Section 3, we present our findings for the dependence of B/T on r d , m d , and c. In Sections 4 and 5, we compare our results with previous models and the observed morphology-luminosity relation, respectively. In Section 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.

Initial conditions
We make our problem dimensionless by expressing all lengths in units of the virial radius R vir , all masses in units of the virial mass M vir and all speeds in units of the virial velocity V vir = GM vir /R vir . 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 lowercase letters refer to dimensional and dimensionless quantities, respectively. As M vir is the total mass within R vir , m d = M d /M vir and 1 − m d 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: 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: 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 Section 4), we run all our simulations for h/r d = 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/r d is in the range 0.088 -0.18. Fig. 1 shows the isodensity contours that contain 40%, 50%, 60 %, 70%, 80%, and 90% of the disc mass for our initial configuration. Eq. (3) is used to generate random coordinates for a million stellar particles within the optical radius r opt = 3.2r d 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 10 4 K and is not allowed to form stars. The total mass of the stellar particles is 0.98 × 0.83 m d in dimensionless virial units.
The circular velocity v c (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: where v d (r) is computed as in Freeman (1970) and Stars have velocity dispersion: determined from Eq. (3) through the requirement that our initial condition should be in equilibrium (albeit unstable). Hence, their velocities: 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 v rot (r) = v φ (r) equals the circular velocity v c (r) only for σ = 0. For σ > 0, v rot and v c are linked by the condition: which gives: since ∆v φ = 0 and ∆v 2 φ = σ 2 . Our rotation speeds and thus our particle velocities are computed using Eq. (9) everywhere except in a small central region where we set v rot = 0 because v c < σ. This central region has r r d by construction, since the normalisation of σ(r) is determined by the disc scale-height (Eq. 6) and we have assumed that h r d . However, the fact that σ is computed considering the vertical equilibrium only and that 3σ 2 > v 2 c at r < ∼ 0.07 r d 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 (Section 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 Section 4.

Parameter space
Having fixed the disc scale height h = 0.044 r d and the gas fraction f gas = 0.02, our initial conditions are entirely determined by three parameters: r d , m d , and c.
The dimensionless radius of a thin exponential disc of dimensionless mass m d embedded in an NFW halo with concentration c, ) is completely determined by the disc's spin parameter which is identical to the halo's spin parameter if we follow Mo, Mao, & White (1998) and we assume that the disc Eq. (10) implies that we can use the spin distribution of DM haloes from cosmological N-body simulations to derive plausible values for r d . This is true even if conservation of angular momentum is a crude approximation (Kimm et al. 2011;Stewart et al. 2013;Jiang et al. 2018) because the model by Mo, Mao, & White (1998) reproduces the correct mass -relation for disc galaxies (for example, Cattaneo et al. 2017).
For a flat rotation curve with v c = 1, Eq. (10) gives r d = λ/2. We do not make this approximation. We substitute Eq. (4) into Eq. (10) and solve Eq. (10) numerically for each combination of λ, m d , and c that we wish to consider. However, we find that the approximation r d λ/2 is accurate at the 20% level. Muñoz-Cuartas et al. (2011) found that λ follows a log-normal distribution withλ = 0.044 and σ ln λ = 0.57. Burkert et al. (2016) found a similar distribution with λ = 0.052 and σ ln λ = 0.46. So did Cattaneo et al. (2017) withλ = 0.049 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 r d 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 r d 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. 2013Behroozi et al. , 2019Moster et al. 2013Moster et al. , 2018Tollet 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 r d that we probe in reality.
In contrast to λ, which appears to be a true random quantity, m d strongly correlates with M vir . Considerable observational evidence shows that the stellar-to-halo mass ratio increases with M vir 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 M vir = 10 12 M , the typical stellar mass of the central galaxy is M = 4 × 10 10 M (m d = 0.04 if we assume that the galaxy started as a pure disc). In a halo with M vir = 10 11 M , the typical stellar mass is two orders of . 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 M vir > 10 11 M .
Six values for λ, four values for m d 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 (Section 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, 2, and 3, respectively. The simulations run for this study are 34 + 14 = 48 in total. B/T=0.22 Fig. 5. Galaxy in the simulation with c = 10, m d = 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.

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 128 3 cells and side length 32 r d . Hence, the resolution on the scale of the entire computational volume is l = 0.25 r d . The resolution is increased within six nested cylinders by a factor of two each time (Fig. 1). The radii r/r d 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/r d = 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 r d ) large enough to ensure that the disc remains well resolved when it buckles up, since the scalelength (and therefore the scale-height) of a pseudobulge is usually several times smaller than r d (Gadotti 2009).

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 ). Figs. 2 to 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.
Figs. 2, 3, 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 edgeon. In most bulges, an X-shaped peanut is clearly visible (Fig. 5).
The galaxy with m d = 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).
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 Section 2.1 after Eq. (9) coupled to the lower resolution in the central region. The radius r ∼ 0.1 r d 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 m d = 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 m d correspond to earlier-type morphologies in the Hubble sequence and to higher B/T (Figs. 2 to 4). We can also increase B/T by reducing λ at constant c and m d . At constant λ and m d , B/T usually decreases with c. Hence, if a simulation does not form a pseudobulge, simulations with lower m d , 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 . 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 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: where m d (2.2r d ) = 0.65m d is the disc mass within 2.2r d (Freeman 1970). The first equality in Eq.(12) implies that Eq.
(1) can be rewritten as: By introducing the new parameter α = 0.31/ , Eq. (13) can be written in the alternative form: 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 is the disc's fractional contribution to the gravitational acceleration at r = 2.2r d and is related to disc's mass fraction within 2.2r d by the equation: .
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 = 3r d rather than at r = 2.2r d for the reason that the rotation curve v c (r) peaks at larger radii when the halo's contribution is considered. We have verified that Eq. (16) keeps holding at r = 3.2r d (the factor 1.3 does not change to the first decimal digit).
Eqs. (12) and (14) explain why B/T increases with m d and why it decreases with λ and c (Figs. 2, 3, and 4). The higher m d , 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 semiqualitative considerations and to explore the relation between B/T and f d in quantitative detail. Fig. 7 shows B/T versus f d at r = 2.2r d (where the rotation curve of an exponential disc peaks) and at the optical radius r opt = 3.2r d , although we have investigated the relation at other radii down to r = 0.5r d . The vertical blue lines correspond to: for α 0.28. The qualitative picture is the same at r = 2.2r d and r = 3.2r d . At f d f crit d , 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 f d f crit d , all galaxies develop a pseudobulge and display a tight correlation between B/T and f d . At intermediate f d , 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 f d that separates the two populations remains fundamentally valid.  Fig. 7. Relation between bulge-to-total mass ratio B/T and initial disc fraction f d (r) for r = 2.2r d (left) and r = 3.2r d (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: m d = 0.005 (green), m d = 0.01 (black), m d = 0.02 (blue), and m d = 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 f d 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).
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 f d , we find: and 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 f d is lower at 3.2r d than it is at 2.2r d . This makes sense because f d (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 f d 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 f d (r) (that is, the rootmean-square residuals) decreases slowly but systematically from 0.12 dex at r = 0.5r d to 0.11 dex at r = 3.2r d . Hence, the ELN criterion is fairly insensitive to the radius at which it is applied, but a criterion at r = 3.2r d appears to be preferable.
If the condition f crit d 0.31 is applied at r = 3.2r d , then all discs with f d (3.2r d ) > f crit d 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 f d (3.2r d ) < f crit d 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.

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 r d and 3.2 r d are 0.65 m * and 0.83 m * , respectively. Eq. (16) gives: 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 (m DM is the mass of the DM within this radius). The disc becomes marginally stable when: where m b and m * − m b are the final bulge mass and the final disc mass, respectively. We compute the final bulge fraction B/T = m b /m * by solving Eqs. (20) and (21), and find: where f crit d ≤ for f d → 1. The limit for f d → 1 reduces to B/T → 0.71 for f crit d = 0.31 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 f crit d does not improve the overall agreement between Eq. (22) and the simulations because it reduces the discrepancy with the symbols at high f d , but it makes it worse at low f d .
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: 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 twocomponent 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, 3 2 M * σ 2 , equals the energy released by the formation of a bulge, M b ∆e.
To compute ∆e, one should know the initial and final radii of the stars that migrated to bulge, but one can expect ∆e ∼ V 2 c within a factor of order unity (Cacciato et al. 2012) and compute B/T = M b /M * by solving the equation: The calculation of σ is based on the assumption that σ self-regulates so that the Toomre (1964) stability parameter is always close to critical value Q c for local stability. This assumption gives: with Q c = 1.68 for a razor-thin sheet of stars with a Maxwellian velocity distribution (see Binney & Tremaine 2008). Let expand Eq. (26) by using Ω = V c /R: V 2 c (R)/R 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 ΣR 2 ). This contribution can be decomposed into the product of two terms: the stellar contribution relative to the total stars plus DM (the f d factor, since all the stars are in the disc at t = 0) and the stellar fraction in the disc after a bulge has formed, M d /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. Eq. (27) can thus be rewritten: where η is a fudge factor (the disc mass within R scales as, but is not equal to, ΣR 2 ). By substituting Eq. (28) into Eq. (25) and using M d = M * − M b , we find: where all uncertainties, and notably the one on ∆e/V 2 c , are reabsorbed in the value of η.
It is straightforward to solve the second-degree Eq. (29) and compute B/T as a function of f d . 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 f d than what we find in our simulations, especially at high f d and when it is applied at r = 2.2 r d . However, it is the model that works best.

Comparison with observations
Our simulations predict that B/T is related to V d /V c . 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) where R opt = 3.2R d (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.
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 V d (R opt )/V c (R opt ) to M I 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 M I 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). 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 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. 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 formula 5 : We retrieve M i and M z 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 -M I 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. The red squares in Fig. 8 show the B/T -M I 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: 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 M I < ∼ −20.5. At M I > −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 (S 4 G) to study the morphologies of 2352 with distances < 40 Mpc. The S 4 G'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 S 4 G 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 M I 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 M I < −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 < M I < −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 f d to the gravitational acceleration at 2.2 R d and found f d = 0.42-0.74. The main uncertainty is value of R d . We predict B/T = 0.09-0.16 for f d = 0.53 (R d = 2.6 kpc; Robin et al. 2003;Jurić et al. 2008) and B/T = 0.19-0.33 for f d = 0.74 (R d = 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.

Application to SAMs
Our key result is a model for B/T as a function of f d , the disc's contribution to the total gravitational acceleration. To use this model predictively, we must know how f d depends on observable quantities such as luminosity or stellar mass. In Section 4, we followed a semi-empirical approach. The conversion from M I to f d was based on observational data. Here,Eq. (19) is incorporared into the GalICS 2.0 SAM Koutsouridou & Cattaneo 2019). Hence, the values of f d 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 ). 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 M 1 /M 2 , where M 1 and M 2 are the total masses of the merging galaxies within their respective baryonic half-mass radii (M 1 > M 2 ). Cattaneo et al. (2017) assumed a sharp transition at a critical mass ratio M 1 /M 2 = m = 4. In major mergers (M 1 /M 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 andHammer et al. 2018). In minor mergers (M 1 /M 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 f tr transferrred to the bulge (or the cusp in the case of gas) is not always unity but varies gradually from f tr = 0 for M 1 /M 2 = m to f tr = 0 for M 1 /M 2 = 1. We assume: where m = 3. Using m = 4 would increase B/T at 10 10.5 M < M * < 10 11 M but also at M * < 10 11 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

B/T
Mergers and disc instabillities 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.
of f tr on M 1 /M 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 (Section 3). In GalICS 2.0, r d 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 r d . 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 finder 6 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); f d is the contribution of the disc to the total gravitational acceleration at r = 3.2r d . The other two contributions are those of the bulge and the halo. Let M d , M cb , and M bar be the masses of the disc, the classical bulge inclusive of the cusp, and the bar inclusive of the pseudobulge, respec-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). tively, where M bar = 0 initially. If the B/T computed with Eq. (19) is B/T > M cb /(M cb + M bar ), then a bar will form until: To understand how f d 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.2r d . We therefore compute f d 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). M bar 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 M bar . A classical bulge from a previous major merger lowers f d , 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 * > 10 11 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 * < 10 10.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 * ∼ 10 10 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 * < ∼ 10 10 M ) and high (M * > ∼ 10 11 M ) masses. It is only in the intermediate mass-range M * ∼ 10 10.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;660000 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.

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 v d /v c > 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 v d /v c 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 f d = (v d /v c ) 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: 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 f d (2.2r d ) < 2α is likely to be globally stable and should not form any pseudobulge, statistically at least. However, it is possible that f d (r) = 2α for r < 2.2r d , 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: where σ is the radial one-dimensional velocity dispersion, Σ is the stellar surface density and is the epicyclic frequency. The disc is stable for Q > 1. In our razor-thin discs: (Section 2.1). By substituting the σ from Eq. (37) and Ω = V c /R with V c 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 Section 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 m d = 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 f d (3.2r d ) 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, V c (2.2R d )/ GM d /R d 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 v d /v c < 0.5 at r = 2.2r d , 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 f d (2.2r d ) really measures is the timescale t b for bar formation. The discs with f d (2.2r d ) < f crit d 0.31 are simply those with t b > 10 Gyr. They also claim that the shear rate: is the physical quantity that is most relevant to determining B/T (Ω = v c /r; also see D'Onghia 2015), and demonstrate a tight correlation between B/T and Γ(2.2r d ).
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 10 0 7 × 10 −1 8 × 10 −1 9 × 10 −1 Γ(2.2r d ) 10 −2 10 −1 B/T Fig. 11. B/T versus the exponent Γ of Ω(r) at r = 2.2r d . 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.
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 time 7 . In relation to the second claim, our simulations confirm that B/T and Γ(2.2r d ) are tightly related (Fig. 11). However, we do not consider that Γ adds substantial information with respect to f d 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 v d and v h . Substituting Eq. (4) into Eq. (38) gives: where: and: At r = 2.2r d , v d has a maximum, so that β d (2.2r d ) = 0; β h (2.2r d ) depends on cr d , but, for reasonable values of c and r d , its value is always β h (2.2r d ) 0.4 (±0.1 at most).
In contrast, f d can vary by an order of magnitude from one galaxy to another. Hence, it is f d that determines Γ, as we intended to demonstrate.
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 t b > ∼ 10 Gyr. Rapid gas accretion, whereby the gas fraction becomes significant, is beyond the scope of this article, which is on stellar discs. 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 v d /v c helps to overcome the problem highlighted by the above studies. We may erroneously predict the formation of a pseudobulge in a galaxy with v d /v c > α that should not have formed one in reality. However, if v d /v c 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 v d /v c < α 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), andHammer 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), andSalo 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 * < ∼ 10 10 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