© The Authors 2025

1. Introduction

Galactic bars are elongated structures of stars, gas, and dust that extend from the centres of many disc galaxies, including the Milky Way. They are common features observed in spiral galaxies across the Universe, and at least two-thirds of the local spirals have a bar (e.g. Sellwood & Wilkinson 1993; Masters et al. 2011; Gavazzi et al. 2015). Bars are also observed in high redshift galaxies, although the fraction is still uncertain (e.g. Guo et al. 2023; Le Conte et al. 2024; Guo et al. 2024). These bars form due to gravitational instabilities within the galactic disc, which lead to the redistribution of angular momentum (e.g. Athanassoula & Misiriotis 2002; Athanassoula 2003; Debattista et al. 2006; Sellwood 2012). This process could funnel interstellar gas towards the galactic centre, often resulting in an enhanced central bulge and increased star formation activity in the core regions (e.g. Spinoso et al. 2017; Donohoe-Keyes et al. 2019; George et al. 2019). Stellar bars play an active role in the secular evolution of galaxies by redistributing mass and angular momentum, thereby influencing the overall morphology and dynamics of the galaxy.

One open question is which disc galaxies form a galactic bar and when. Many studies have investigated the necessary conditions for a disc to form a bar since the early 1960s, with one example being the seminal work of Toomre (1964), which emphasises the importance of random stellar motions in maintaining the stability of galactic discs and provides a framework for understanding the dynamics of stellar systems. Toomre (1964) proposed a dimensionless parameter, the Toomre parameter Q_T, that is employed to evaluate the stability of a differentially rotating disc:

$$\begin{array}{c}\hfill Q_T=\frac{\kappa {\sigma }_{\mathrm{R}}}{3.36G\mathrm{\Sigma }},\end{array}$$(1)

where κ and σ_R denote the epicyclic frequency and velocity dispersion of the stars in the radial direction, respectively; Σ denotes the face-on stellar surface density profile; and G denotes the gravitational constant. The disc stability in response to axisymmetric density perturbations is determined by the Toomre parameter. If Q_T < 1, the disc is unstable. Ostriker & Peebles (1973) proposed that the ratio between the kinetic energy of rotation to the total gravitational energy of the galaxy has to be ≲0.14 to be stable. In addition, by studying models with exponential surface density profiles and flat rotation curves at large radii, Efstathiou et al. 1982, (ELN) investigated the formation of bars. The authors found that a critical criterion for stability in bar-like modes in disc galaxies is given by

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{ELN}}\equiv v_{\mathrm{Max}}/{(M_{\mathrm{disc}}G/r_{\mathrm{disc}})}^{1/2}\ge 1.1,\end{array}$$(2)

where v_Max is the maximum rotational velocity of the system, r_disc corresponds to the scale length of the exponential disc, and M_disc is the total disc mass. They also explored the role of a hot halo component in increasing this ratio to stabilise the galaxy against bar formation. The study highlights the importance of balancing parameters such as rotational velocity, scale length, and total disc mass to maintain stability in disc galaxies. Izquierdo-Villalba et al. (2022) have investigated if the ELN criterion is satisfied for massive galaxies with strong bars in the TNG100 and TNG50 simulations. They find that more than 70% of the barred galaxies satisfied the ELN criterion, while in galaxies that do not satisfy the criterion, it is because the bar formation is triggered by a tidal interaction rather than secular evolution. This criterion was revised by Mo et al. 1998, (MMW) with a focus on more fundamental properties, including the spin parameter (λ), the stellar disc mass fraction (m_disc = M_disc/M_h), and the angular momentum of the disc (J_disc) and halo (J_h). The authors found in their disc and halo profile that the discs are stable if λ_MMW > λ_crit, where λ_MMW is defined as

$$\begin{array}{cc}\hfill {\lambda }_{\mathrm{MMW}}& \equiv \lambda \frac{J_{\mathrm{disc}}}{J_{\mathrm{h}}}/\frac{M_{\mathrm{disc}}}{M_h}\hfill \\ \hfill & \gtrsim {\lambda }_{\mathrm{crit}}=\sqrt{2}{ϵ}_{\mathrm{ELN}}^2{m}_{\mathrm{disc}}{f}_c^{1/2}{f}_{\mathrm{R}}^{-1}{f}_{\mathrm{V}}^{-2}.\hfill \end{array}$$(3)

Here, $\lambda =\mid {\mathit{J}}_h\mid /(\sqrt{2}{M}_{200}{V}_{200}{r}_{200})$ where J_h is the angular momentum evaluated at r₂₀₀ (Bullock et al. 2001). The functions f_c, f_R, f_V weakly depend on the profile models of the disc and halo as seen in Fig. 3 by Mo et al. 1998, (MMW). Because of this, a useful but rough approximation is that λ_MMW > m_disc, which could be translated into $\lambda \frac{j_{\mathrm{disc}}}{j_{\mathrm{h}}}/\frac{M_{\mathrm{disc}}}{M_h}>1$, where j_disc and j_h are the specific angular momentum of the disc and halo.

Although these criteria are widely used to predict bar instabilities in galaxy formation models, some theoretical and observational studies have shown that bar formation could be more complex and that other ingredients may affect bar evolution. For instance, Athanassoula (2008) discusses how ELN does not take into account the (stabilising) effect of the disc velocity dispersion or the central concentration. Ghosh et al. (2023) have found that when using different contributions of the thick discs, some of their galaxies do not satisfy the ELN criterion. Romeo et al. (2023) found that for more diverse barred and unbarred disc galaxy samples, such as the galaxies of types S0–BCD from the Spitzer Photometry and Accurate Rotation Curves (SPARC) sample and the galaxies of type Im from the Local Irregulars That Trace Luminosity Extremes The HI Nearby Galaxy Survey (LITTLE THINGS), some criteria could not be applied for all the diversity of galaxies found in the local Universe. In addition, other components, such as the bulge, could aid in avoiding the formation of a bar. Saha & Elmegreen (2018) show that disc galaxies with a bulge denser than their disc can prevent bar formation even though the disc is maximal or unstable to bar instabilities. This is due to the rapid loss of angular momentum and rapid heating in the centre of initially strong bar and spiral arms. Also, high-rotating dark matter haloes could prevent the growth of bar instabilities. Saha & Naab (2013) have also showed that bar formation is favoured in haloes corotating with the disc with a dark matter spin parameter (λ_DM) in the range of 0 and 0.07. However, Long et al. (2014) have shown that the bar growth is halted in haloes with larger spins (λ_DM ≳ 0.03). Collier et al. (2018) have suggested that not only the spin parameter but also the shape of the dark matter halo is important to the evolution of bars. Yurin & Springel (2015) have considered the stability of discs by inserting already formed stellar discs in haloes using the Aquarius simulation project that consists of zoom-in dark matter-only simulations. These authors found that the 3D shape of the dark matter halo could affect the stability of the disc. Independent of their results, all of these studies show the importance of the angular momentum exchange between the dark matter and the stellar components in the formation of a bar (e.g. Sheth et al. 2012; Athanassoula et al. 2013).

The bulge and disc structural properties are related to the formation of galactic bars and are the result of the formation history of the galaxy in the first place (Mo et al. 1998; Rosito et al. 2019; Bi et al. 2022a; Joshi et al. 2024). One crucial driver that could influence the morphology of galaxies since early times is star formation feedback and active galactic nuclei (AGN) feedback. Supernova (SN) feedback refers to the impact of SN explosions on the interstellar medium (ISM). Massive stars end their lives in SN explosions, releasing vast amounts of energy, radiation, and heavy elements. This feedback process has several critical effects on the ISM: it injects energy, creating shock waves that can compress or disperse gas clouds; it enriches the ISM with heavy elements necessary for new star formation; and it regulates the star formation rate by preventing runaway star formation that could deplete the galaxy gas reservoir. AGN feedback is powered by the accretion of gas into the supermassive black holes (BHs) that are thought to be located at the centre of all (sufficiently large) galaxies, and they have been proposed to mostly affect star formation in massive galaxies (e.g. Springel et al. 2005; Bower et al. 2006; Croton et al. 2006). Both SN and AGN feedback are critical processes in all the current galaxy formation models (e.g. Springel et al. 2005; Bower et al. 2006; Croton et al. 2006). Both processes inject energy and momentum into the ISM, potentially influencing star formation and regulating the galaxy build-up and subsequently the bar formation. In particular, SN feedback has been demonstrated to be efficient in the accumulation of central low angular momentum gas, which allows for the build-up of cold discs with flat rotation curves comparable with observations (e.g. Navarro & Benz 1991).

Barred galaxies have started to be explored in a fully cosmological context recently thanks to high-resolution zoom-in simulations (e.g. Kraljic et al. 2012; Scannapieco & Athanassoula 2012; Bonoli et al. 2016). In particular, Bonoli et al. (2016) have presented a zoom-in simulation of a Milky Way-type galaxy, ErisBH, which is a sibling of the Eris simulation (Guedes et al. 2011) but with BH subgrid physics included. Bonoli et al. (2016) and Spinoso et al. (2017) find that the simulated galaxy forms a strong bar below z ∼ 1, and the authors point out that the disc in the simulation is more prone to instabilities compared to the original Eris, possibly because of early AGN feedback affecting the central part of the galaxy. Zana et al. (2019), studying a larger suite of Eris, highlights the effects of the feedback processes on the formation time and final properties of the bar. Recently, Fragkoudi et al. (2020, 2025) have studied barred galaxy formation in the Auriga suite, which consists of high-resolution, magnetohydrodynamical cosmological zoom-in simulations of galaxy formation (Grand et al. 2017). The authors find that the simulated galaxies that can reproduce many chemodynamical properties of the stellar populations seen in the Milky Way bulge have quiet merger histories. Auriga simulations have also been employed by Irodotou et al. (2022) to determine whether AGN feedback can influence the final properties of a bar but not its formation. The formation of bars in high-redshift counterpart galaxies of the local spiral galaxies formed in dense environments and with two distinct SN feedback subgrid models has been investigated by Bi et al. (2022b) using high-resolution zoom-in simulations. The authors find that the properties and evolution of galactic bars are significantly impacted by the varying SN feedback. These changes are associated with the interactions and cold accretion processes, which were influenced by the type of feedback implemented.

With cosmological hydrodynamic simulations (Vogelsberger et al. 2014a; Schaye et al. 2015; Pillepich et al. 2018b; Nelson et al. 2018, see also the recent review of Vogelsberger et al. 2020; Crain & van de Voort 2023), it has been possible to follow the formation and evolution of the barred galaxy population in a cosmological context and with a significant statistical sample. Analysing the EAGLE simulation, Algorry et al. (2017) have found that bars slow down quickly as they evolve, expanding the inner parts of the dark matter halo. Rosas-Guevara et al. (2020, 2022) study massive barred disc galaxies at z = 0 in the TNG100 and TNG50 simulation (see also Peschken & Łokas 2019; Zhao et al. 2020; Zhou et al. 2020 for Illustris and IllustrisTNG), and they report that barred galaxies have a lower star formation rate and are more gas poor than unbarred galaxies. Following the evolution of barred galaxies back in time, Rosas-Guevara et al. (2020) find that these objects assembled most of their disc components before bar formation and at earlier times than unbarred galaxies (see also Izquierdo-Villalba et al. 2022).

Even though these simulations have shown the possibility of understanding the evolution and formation of bars in a cosmological context, many complex galaxy formation processes take place simultaneously, and it is difficult to disentangle all the relevant processes. Taking advantage of previous works, we performed zoom-in simulations of six Milky Way-like galaxies from the TNG50 simulations. These galaxies have a stable galactic bar at z = 0 and formed between z = 3 and z = 1.5 (more than 8 billion years ago), and they have quiet merger histories and isolated haloes. These zoom-in simulations were run with the AREPO code (Springel 2010) and similar initial conditions as the original TNG50 simulation (Pillepich et al. 2019; Nelson et al. 2019b), in the sense that the large-scale tidal field is very similar, but only the galaxy region was simulated with the same resolution as TNG50. Our primary goal is to study the interplay of galaxy formation physics with the formation and evolution of galactic bars. To that end, we used the same galaxy formation physics as Pillepich et al. (2018b) except that we allowed for variations in the SN and AGN feedback models in order to identify the specific impacts of these feedback channels. We also explored how these variations impact different bar instability criteria used in the field, namely, the Toomre (1964), Efstathiou et al. (1982), (ELN), and Mo et al. (1998), (MMW) criteria.

The paper is structured as follows. In Section 2, we introduce our selection of the disc galaxies and the methodology for identifying a bar and dynamics components, and we give a brief overview of the SN and AGN feedback subgrid physics of TNG50 and their variations used here. In Section 3, we study the redshift evolution of their bars and their host galaxy components. We study the criteria of bar instability and how they are affected by the different galaxy formation models in Section 4. Finally, in Sections 5 and 6, we discuss and summarise our findings.

2. Methodology

In this section, we describe the object selection criteria, the generation of the initial conditions, and the TNG galaxy formation model. We especially describe the variations of the galaxy formation model. This project accounts for 42 resimulations, with six haloes, each with seven different variations of the galaxy formation model. These simulations were performed with the moving-mesh AREPO code (Springel 2010), combining Tree-PM and Godunov/finite volume methods to discretise space. The quasi-Langragian scheme is second-order in space and time.

The cosmology parameters adopted are from the Planck Collaboration XIII (2016) Cosmology: Ω_Λ = 0.6911, Ω_m = 0.3089, Ω_b = 0.0486, σ₈ = 0.8159, h = 0.6774, and n_s = 0.9667 where Ω_Λ, Ω_m, and Ω_b are the average densities of dark energy, dark matter and baryonic matter in units of the critical density at z = 0, σ₈ is the square root of the linear variance, h is the Hubble parameter ($H_0\equiv h100{\mathrm{km}\mathrm{s}}^{-1}$) and, n_s is the scalar power-law index of the power spectrum of primordial adiabatic perturbations.

The suite of zoom-in simulations has a similar resolution to the TNG50 simulation, whose particle mass resolution is 4.5 × 10⁵ M_⊙ for dark matter particles, whereas the mean cell mass resolution is 8.5 × 10⁴ M_⊙ for gas. A comparable initial mass is passed down to stellar particles, which subsequently lose mass through stellar evolution. The spatial resolution for dark matter and stellar particles is 575 comoving pc down to z = 1, after which it remains constant at 288 pc in physical units down to z = 0. In the case of the gas component, the gravitational softening is adaptive and based on the effective cell radius, down to a minimum value of 72 pc in physical units, which is imposed at all times.

Galaxies and their haloes are identified as bound substructures using a FOF (Friends-of-Friends) and then the SUBFIND algorithm (Springel et al. 2001) and tracked over time by the SUBLINK merger tree algorithm (Rodriguez-Gomez et al. 2015). Halo masses (M₂₀₀) are defined as all matter within the radius R₂₀₀ for which the enclosed mean density is 200 times the critical density of the Universe. In each FOF halo, the central galaxy (subhalo) is the first (most massive) subhalo of each FOF group. Its satellites are the remaining galaxies within the FOF halo. The stellar mass of a galaxy is defined as all the stellar matter assigned to host subhaloes.

2.1. Object selection

The suite of our zoom-in simulations is based on barred galaxies formed in the TNG50 simulation. For this, we use the bar catalogue of Rosas-Guevara et al. (2022, herein RG22) and select central disc galaxies with a bar at z = 0. We selected those whose dark matter haloes with masses (M₂₀₀) between 10^11.5 and ${10}^{12}{\mathrm{M}}_{\odot }$ comparable to the Milky Way-like halo, yielding 31 galaxies. Furthermore, we concentrated on relatively isolated haloes. In particular, we used the tidal parameter that is defined as τ = M_200, i/M_200, tar(R_200, tar/R_i)³ (Dahari 1984), where M_200, i and M_200, tar are the halo masses of the closest massive neighbour and the target halo, respectively. R₂₀₀ is the size of the target haloes and R_i is the distance to the closest massive neighbour. We calculated τ of our possible halo candidates in order to assess the interaction strength between our candidates and their closest massive companion. This allowed us to choose MW-like haloes that are isolated in the TNG50 dark matter simulation. We calculated the tidal parameter for our halo candidates and selected those with values lower than 10⁻⁹. This condition corresponds to those haloes whose closest massive neighbour galaxy is smaller than the candidate halo (M_200, tar) and at a distance of 10R_200, tar (see Grand et al. 2017). We also verified that these galaxies have relatively quiet merger histories (no more than two major mergers in z < 2). We selected six dark matter haloes from the TNG50 simulations (Nelson et al. 2019a) whose baryonic counterpart is a disc galaxy with a stellar mass greater than ${10}^{10}{\mathrm{M}}_{\odot }$ and that have a small bulge (B/T lie between 0.12 and 0.25). These galaxies hold a stable bar at z = 0 that has formed between z = 1.5 and z = 3. The bar sizes lie between 2 and 7 kpc, and the disc sizes between 2.5 and 6 kpc. Most of the bars are strong bars (A_2, max > 0.3; see Appendix A for bar identification method). The right panel of Fig. 1 shows mock images (taken from TNG database Pillepich et al. 2019) in JWST NIRCam F200W, F115W, and F070W filters of the six barred galaxies from the TNG50 that were re-simulated in this work and Table 1 shows some properties of the bars and discs at z = 0. We note that five of the six haloes belong to the catalogue of Milky Way analogues from TNG50 (Pillepich et al. 2024). The left panel of Fig. 1 depicts the dark matter density map of the TNG50 simulation, and the magenta circles mark the locations of the selected haloes which are not located in the densest regions of the cosmic web.

Fig. 1. Left panel: Dark matter distribution of the TNG50 simulation. Magenta empty circles show the position of the selected haloes of the zoom-in simulations. Right panel: Mocked images in JWST NIRCam F200W, F115W, and F070W filters (face-on) calculated in Nelson et al. (2018), (TNG database) for the TNG50 galaxies that are re-simulated. The NIRCam blue channel highlights the young population of the galaxy, and the NIRCam red channel older populations. The original TNG50 galaxies present a massive disc with a smaller bulge and form a strong bar.

2.2. Initial conditions

The initial conditions for the zoom-in simulation suite were created by selecting all dark matter particles in the TNG50 parent simulation at z = 0 that are located within 3 × R₂₀₀ (or 5 × R₂₀₀ in cases where the former was not sufficient to get contamination free objects) of the centres of the chosen objects. The particles were then traced back to the starting redshift of the TNG50 simulation to obtain the Lagrangian region of each object. Next, a tree structure was used to produce the unperturbed particle distribution for each zoom-in simulation. Within the Lagrangian region, the tree was refined down to the same resolution as the initial particle distribution in the TNG50 simulation. Then layers of lower resolution boundary particles were added, with the resolution set based on the distance to the nearest edge of the high-resolution region, until a base resolution (corresponding to 135 particles on a side) was reached that covers the remaining part of the TNG50 volume. Perturbations were added at z = 127, the starting redshift adopted in the zoom-in simulations, using the Zeldovich approximation. Mode amplitudes and phases were drawn in Fourier space using the same random number seed and algorithm as in the N-GENIC code (Springel 2015), so that perturbations consistent with the TNG50 parent simulation were imprinted. For lower resolution boundary particles, only modes up to the Nyquist frequency that corresponds to that resolution level were included. An updated version of the employed zoom-in initial conditions code will be presented in Puchwein et al. (in prep.).

Upon starting the AREPO simulations, the particles in the initial conditions file are split into dark matter and gas components and then diagonally displaced with respect to each other by half a diagonal particle spacing, while keeping the centre of mass unchanged. A uniform magnetic seed field with a comoving field strength of 10⁻¹⁴ Gauss is assumed at the initial redshift.

2.3. The galaxy formation model and variations

We used the IllustrisTNG galaxy formation model (The Next Generation) project¹ (Nelson et al. 2018; Naiman et al. 2018; Pillepich et al. 2018b; Marinacci et al. 2018; Springel et al. 2018) and its predecessor ILLUSTRIS (Vogelsberger et al. 2013, 2014a,b; Genel et al. 2014; Nelson et al. 2015; Sijacki et al. 2015) with changes to star formation feedback, supermassive black hole growth, AGN feedback, and stellar evolution and chemical enrichment.

Gas radiative mechanisms are implemented with primordial (Katz et al. 1996) and metal line cooling, and heating by a time-dependent ultraviolet background field from stars and luminous AGN (Faucher-Giguère et al. 2009). Star formation in the dense interstellar medium is treated stochastically following an empirical Kennicutt-Schmidt relation (Springel & Hernquist 2003). Each stellar particle represents a population of stars with a common birth time following a Chabrier initial mass function. The stellar evolution is modelled in order to calculate chemical enrichment and mass expelled into the interstellar medium due to AGB stars, SNIa and SNII. Also, the evolution and production of ten elements (H, He, C, N, O, Ne, Mg, Si, Fe, & Eu) are individually tracked.

Previous studies have shown that physical processes that regulate star formation in galaxies can affect the formation of a bar (e.g. Bonoli et al. 2016; Zana et al. 2019; Irodotou et al. 2022). Taking advantage of this, combined with the parametric study done for the TNG collaboration (Pillepich et al. 2019), we examine variations in physical processes that regulate the star formation of galaxies. We can divide these variations into two categories: (1) those models that present modifications in the galactic winds generated from star formation, and (2) those models that present modifications in the BH physics of the model. The models of winds generated from star formation and AGN feedback are explained in detail in the following subsections.

2.3.1. The galactic wind model and its variations

For our suite of zoom-in simulations, we use the TNG model for galactic-scale, star formation-driven kinetic feedback (Pillepich et al. 2018b). In this model, wind particles absorb the thermal energy of the inherited gas cell and are injected isotropically. The initial orientation of the wind particles is random and are launched with an initial speed that scales with the local, one-dimensional dark matter velocity dispersion σ_DM (as in Eq. (14) of Oppenheimer & Davé 2006, 2008; Vogelsberger et al. 2013), calculated with a weighted kernel over the N = 64 nearest DM particles. In addition, the wind velocity has a dependency on redshift and a minimum value ${\mathit{v}}_{\mathit{w},\min }$ given by

$$\begin{array}{c}\hfill v_w=\max [k_w{\sigma }_{\mathrm{DM}}{(H_0/H(z))}^{1/3},v_{w,\min }],\end{array}$$(4)

where k_w is a dimensionless parameter. We note that the wind velocity depends on redshift, this ensures that the wind velocity and virial halo mass increase in scale with redshift in the same way (see Pillepich et al. 2018b). This choice of redshift-independent wind velocities at constant halo mass is motivated by semi-analytical model findings, where a comparable approach was required to recreate observed stellar mass functions and rest frame B- and K-band brightness functions throughout redshift (Henriques et al. 2013). Once the wind injection velocity is known, the wind mass loading factor is determined by the specific energy available, e_w, which is connected to the energy released by SNII per produced stellar mass in the model. The TNG model includes two categories that affect the available wind energy: (i) some fraction of this energy is thermal, as determined by a parameter τ_w; and (ii) wind energy is dependent on the metallicity of the star-forming gas cells such that galactic winds are weaker in higher metallicity environments. The mass loading factor upon injection is calculated as follows:

$$\begin{array}{c}\hfill {\eta }_w={\dot{M}}_w/{\dot{M}}_{\mathrm{SFR}}=\frac{2}{v_w^2}{e}_w(1-{\tau }_w),\end{array}$$(5)

where ${\dot{M}}_{\mathit{w}}$ corresponds to the rate of gas mass to be converted into wind particles and ${\dot{M}}_{\mathrm{SFR}}$ the instantaneous, local, star formation rate. Then the wind energy available from the star-forming gas cells depends on metallicity Z as follows:

$$\begin{array}{c}\hfill e_w=\overline{e_w}[f_{w,Z}+\frac{1-f_{w,Z}}{1+Z/{(Z_{w,\mathrm{ref}})}^{{\gamma }_w,Z}}]N_{\mathrm{SNII}}{E}_{\mathrm{SNII},51}{10}^{51}\mathrm{etg}{M}_{\odot }^{-1},\end{array}$$(6)

where $\overline{e_{\mathit{w}}}$ is a dimensionless free parameter of the model, E_SNII, 51 denotes the available energy per core-collapse SNe in units of 10⁵¹ erg. N_SNII is the number of SNII per formed stellar mass (in solar mass units) and depends on the shape of the IMF and the assumed minimum mass of core-collapse SNe. f_w, Z is a parameter that reduces the energy at injection for gas cells with metallicities much larger than a given value $Z_{\mathit{w},\mathrm{ref}}$.

We focus on the energy available from star-forming cells of the galactic winds over time in our suite of zoom-in simulations, which limits the intensity of the galactic wind. The energy is regulated by the dimensionless free parameter, $\overline{e_{\mathit{w}}}$, which is between $\overline{e_{\mathit{w}}}=0$ and $\overline{e_{\mathit{w}}}=7.2$. These values correspond to the No Wind (NW) model and the Strong Wind (SW) model, respectively. The value $\overline{e_{\mathit{w}}}=3.6$ is employed in the TNG50 simulation, and we refer to it here as the TNG50-like model (TNG50-like). In total, we have five modifications, including weaker winds (WW) and medium winds (MW). These variations are summarised in Table 2.

We remark that some of these variations have been explored to see their impact on the complete galaxy population in Pillepich et al. (2018b). We can see in their Fig. B1 how the global star formation history of a galaxy population can be affected by the strength of the winds. It is also worth mentioning that we have tried other variations at low resolution, such as a higher fraction of the thermal energy in the wind (τ_w) or higher wind velocity (κ_w), and we found that the effects on the star formation of a galaxy were similar as varying $\overline{e_{\mathit{w}}}$. However, at high redshift, the variations in the strength of the winds $\overline{e_{\mathit{w}}}$ affect the star formation density of the galaxy population (see Fig. B1 of Pillepich et al. 2018b).

2.3.2. The BH model and its variations

The other model used in our suite of zoom-in simulations is the subgrid physics of supermassive black holes (BHs) that was presented in Weinberger et al. (2017). The BHs are formed in massive haloes ($M>7.38\times {10}^{10}{\mathrm{M}}_{\odot }$) with an initial black hole mass of $1.18\times {10}^6{\mathrm{M}}_{\odot }$, and can grow via two growth channels: BH mergers and gas accretion. Gas accretion is Eddington limited and allowed to accrete at the Bondi-Hoyle-Lyttleton accretion rate. To define high accretion and low accretion states, we checked whether the Bondi-Hoyle-Lyttleton accretion rate exceeds a specific fraction of the Eddington limit using

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Bondi}}/{\dot{M}}_{\mathrm{Edd}}\ge f_{\mathrm{Edd},\mathrm{th}},\end{array}$$(7)

where the Bondi-Hoyle-Lyttleton accretion rate and the Eddington limit are defined, respectively, as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Bondi}}=\frac{4\pi G^2{M}_{\mathrm{BH}}^2\rho }{c_s^3}\text{and}{\dot{M}}_{\mathrm{Edd}}=\frac{4\pi G{M}_{\mathrm{BH}}{m}_{\mathrm{p}}}{{ϵ}_{\mathrm{r}}{\sigma }_{\mathrm{T}}c},\end{array}$$(8)

where c is the speed of light in vacuum, m_p the mass of the proton, and σ_T the Thompson cross-section. ϵ_r represents the radiative accretion efficiency. M_BH is the mass of the black hole, while ρ is the gas density of the surrounding gas; c_s is the effective speed of sound near the black hole and corresponds to (c_s, therm² + B²/4πρ)^1/2. The f_Edd, th scales with the BH mass to promote the transition from high to low accretion rates for the most massive BHs at late times such that

$$\begin{array}{c}\hfill f_{\mathrm{Edd}.\mathrm{th}}=\min [f_{\mathrm{Edd},0}{\left(\frac{M_{\mathrm{BH}}}{{10}^8{\mathrm{M}}_{\odot }}\right)}^{\beta },0.1],\end{array}$$(9)

where β and f_Edd, 0 are both free parameters. It is worth noting that the values of f_Edd, 0, β, and the threshold mass 10⁸M_⊙ do not have independent values. These values were selected so that the most massive BHs can attain the low accretion rate state, while the less massive BHs struggle to do so.

In the model, there are two modes of AGN feedback, and the transition between the modes is given by Eq. (9). The two modes are (i) thermal quasar mode that heats the surrounding gas of the BH at high accretion rates (e.g. Springel et al. 2005; Di Matteo et al. 2005). The feedback energy in the quasar mode is released continuously as thermal energy. This thermal energy is injected into the surrounding gas and defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_{\mathrm{high}}={ϵ}_{f,\mathrm{high}}{ϵ}_r\dot{M}{c}^2\mathrm{\Delta }t,\end{array}$$(10)

where ϵ_r is the radiative efficiency, taking typical values between 0.04 − 0.4, for optically thin, geometrically thin, and radiatively efficient accretion discs depending on the SMBH spin. In the original TNG runs, it has been set to 0.2, and we do not change it to prevent additional degeneracy. The parameter ${ϵ}_{f,\mathrm{high}}$ is the fraction of energy in quasar mode that couples with the surrounding gas, set to be 0.1 in the original TNG simulation. (ii) The second mode is the kinetic wind mode that produces winds, typically when the SMBH is massive enough and is subject to low accretion rates.

For our suite of zoom-in simulations, we consider two variations in the quasar mode feedback: (1) the model NBH, which is a simulation performed without BH physics, and (2) BHlowEff, where the energy fraction that couples to the surrounding gas in the AGN feedback in quasar mode, ${ϵ}_{f,\mathrm{high}}$, is set to a lower value or 0.05. Both models are compared to the TNG50-like model, which takes the same parameters as the parent simulation TNG50. The values adopted for the modified parameters are found in Table 2.

2.4. Kinematic decomposition and identification of bars

To identify the different morphological components of each galaxy, we employ the kinematic decomposition algorithm MORDOR (Zana et al. 2022) to determine more specific galaxy components. The decomposition is based on the circularity (ϵ) and binding energy (E) phase space, where a minimum in E is identified. The methodology identified five components: classical bulge, pseudobulge, thin disc, thick disc and stellar halo. To calculate the disc size, we use the face-on stellar surface density profiles of the total disc (thin and thick components), fitting an exponential profile.

A bar is identified using the Fourier decomposition of the face-on stellar surface profile. We use the maximum value reached in the ratio between the second and zero terms of the Fourier expansion (A₂(R)) as a proxy for the bar strength. The length of the bar is defined as the maximum radius where the phase (Φ(R)) is constant. If the bar strength exceeds 0.2, the galaxy is classified as barred. Appendix A.1 provides detailed information and is based on Rosas-Guevara et al. (2022). The thickness of the disc (using both thin and thick components) was calculated by fitting the edge-on surface density of the disc component with a square arc secant function in terms of the z-direction. In the case of the bulge size, we employ the face-on stellar surface density profiles, which are simultaneously fitted to the sum of a Sersic and an exponential profile for unbarred galaxies and the sum of two Sersic and an exponential profile where a bar is identified using the Fourier decomposition of the face-on stellar surface profile. Examples and more details of kinematic decomposition, identification of bars and the method of three-component decomposition of surface face-on density profiles can be found in Appendix A. Table 3 shows the median values of the bulge and disc structural properties for the different models. It also includes the median values of the bar strength and length for the models. We do not find any particular trend for the bar strength when we vary models since the relation of the bar length with stellar mass is flat, indicating no direct dependence on stellar mass. From hereafter, we refer to a disc as the sum of the thin and thick disc components and a bulge as a classical bulge.

3. Impact on bar and galaxy properties

In this section, we investigate the impact of altering the strength of the winds and BH physics on the properties of the galaxies and their structural parameters as a function of redshift.

3.1. Properties of bars and galaxies at z = 0

The median mass of the halo in the original TNG50 simulation is $M_{\mathrm{crit}200}={10}^{10.91}{\mathrm{M}}_{\odot }$ (see Table 1) and the median of the zoom-in haloes for the different models ranges between $M_{200}={10}^{10.92}{\mathrm{M}}_{\odot }$ (BHlowEff model) and $M_{200}={10}^{10.98}{\mathrm{M}}_{\odot }$ (NW model) at z = 0 which is at most a difference of 0.06 dex in halo mass. This is not the case for the rest of the properties when higher dispersion is found in varying SNe feedback. Table 3 shows some of the median properties of the galaxies and their bars at z = 0. The median stellar mass ranges between $M_{\ast }={10}^{10.43}{\mathrm{M}}_{\odot }$ (SW model) and $M_{\ast }={10}^{11.06}{\mathrm{M}}_{\odot }$ (NW model) at z = 0. When this is compared to the stellar masses of the haloes in the original TNG50 simulations ($M_{\ast }={10}^{10.61}{\mathrm{M}}_{\odot }$) or TNG50-like model ($M_{\ast }={10}^{10.72}{\mathrm{M}}_{\odot }$), the difference varies between 0.4 and 0.2 dex in stellar mass. The most massive galaxies are those without SN feedback (NW model), whereas the least massive galaxies are those in the SW model. This is expected since SN processes are known to be physical processes that eject gas or prevent gas from forming stars, and could influence the angular momentum of the gas in the galaxy. It is interesting to note that the morphology roughly remains the same as in the original TNG50 simulations at z = 0 in the sense that galaxies have a high disc dominant component, although the stellar mass and sizes of the disc and the bulge are different in all the models (see Table 3). The galaxies with the lowest D/T = 0.48 ((D/T)_fit = 0.51, using Fitting surface profiles) are those in the NW model, whereas the highest D/T = 0.74 ((D/T)_fit = 0.88) are those in the MW model. In the case of the bulges, the lowest B/T = 0.09 ((B/T)_fit = 0.10) is for the galaxies in the SW model, whereas the highest B/T = 0.21 is for the galaxies in the WW model (and in the NW model using (B/T)_fit = 0.55). The compactness also changes. The most compact and massive bulges and thinner and massive discs are those in galaxies in the NW model (R_bulge = 0.44 kpc, h_disc = 1.33 kpc, h_z = 0.55 kpc). We note, however, that the smallest bulge radii correspond to galaxies in the weak wind model. The least compact bulge (R_bulge ≈ 1 kpc) is found in the galaxies in SW and MW models. Also, the galaxies in SW and MW models exhibit the most extended and thicker discs found in galaxies with h_disc = 4.47 kpc, h_z = 2.00 kpc and h_disc = 5.19 kpc, h_z = 1.71 kpc, respectively. Regarding the galaxy properties in the models varying BH physics, the galaxy properties are pretty similar to those in the TNG50-like model. The parent TNG50 galaxies exhibit a bar, while in the zoom-in simulations, 83% of them present a bar at z = 0 (5 of 6 galaxies) in the TNG50-like model. This galaxy, in particular, does not form a bar in any of the models, including the weak wind model, where the bar is present in the other galaxies. The median bar sizes lie between 2 kpc and 3 kpc, and the disc sizes all of which are strong bars (A_2, max ≥ 0.4) as in the parent TNG50 simulations.

The top panel of Fig. 2 illustrates two galaxies at z = 0 that are subject to the various variations in the wind model from NW (no wind model) to SW (strong wind model). The face-on (top panels) and edge-on (bottom panels) stellar density maps are represented in each column, with the energy injection per SN event intensifying from left to right. The disc-like morphology of the galaxy is generally independent of the model, as mentioned previously. The figure shows that the bar does not form in either the NW model or the SW model. The galaxy appears to be more concentrated in the no-wind model, while the galaxy appears to be less concentrated in the SW model. The two objects exhibit this behaviour. The bar forms in the TNG50-like model and the WW model. For the MW model, only one of the galaxies forms a bar. It is essential to acknowledge, however, that the stellar bar properties, including its length and strength, are different for different models.

Fig. 2. Face-on stellar density maps of two distinct galaxies at z = 0. The top figure shows the wind model variations, and each column corresponds to a distinct wind model. The left column shows a galaxy without SN feedback, and the right columns show galaxies with stronger SN feedback. The bottom figure shows different quasar black hole models but with the same SNe feedback model as indicated by Table 2. The black star symbolises the location of the BH. In both the no wind and strong wind models, the disc galaxy does not exhibit the formation of a bar. Quasar BH physics appear to have a limited influence on the formation of a bar.

The two objects are also depicted in the bottom panel of Fig. 2 in two distinct variations in the BH physics: a model without BH physics (NBH) and a BH with less efficient quasar AGN feedback (BHlowEff) when the BH is accreting at higher rates compared to the TNG50-like model but with the same stellar feedback assumed as in the TNG50-like model. We observe that galaxy morphology is not significantly altered by the variations in the BH model, as opposed to the variations in the wind models. However, the final properties of the bar and the evolution of the BH itself will differ.

3.2. Galactic bar and galaxy evolution

We focus now on the evolution of the bar and galaxy properties. Fig. 3 illustrates the median redshift evolution of the bar properties and morphology of the galaxies for each of the four galactic winds models, the three BH models, and the galaxy counterparts in the TNG50 simulation. We observed that only the galaxies in the WW and TNG50-like models develop a stable bar in the galactic wind models. The galaxies in the other variations do not form a stable bar. We also observe that the weak wind model develops the strongest bars, which are also shorter in comparison to the galaxies in the TN50-like model. It is important to note, however, that galaxies in the no wind model develop non-axisymmetric instabilities, but they do not grow over time. This suggests that the disc may be unstable (see Sect. 4); however, refrain from forming a bar. We observed that the galactic bar formed earlier (see cyan vertical lines for the weak wind model, z_bf = 2.5) compared to the TNG50-like models (black vertical lines, with z_bf = 2). However, the scatter in the time of bar formation is large in all the models (see vertical colour lines in Fig. 3), including TNG50 galaxies. This could be attributed to the chaotic-like behaviour caused by feedback processes in the TNG model, affecting the scatter in some galaxy properties (e.g. Genel et al. 2019; Joshi et al. 2024).

Fig. 3. Median evolution of the bar and galaxy morphology for simulations with different galactic wind (top) and BH feedback (bottom) models. Each column corresponds to a different model, as indicated in the figures. From top to bottom rows, the bar strength, bar extent, the stellar mass, and the disc-to-total and bulge-to-total mass fraction are shown. The medians and the distribution between 20th and 80th percentiles correspond to solid lines and shaded regions, respectively. Thinner solid lines correspond to the evolution of individual galaxies. Vertical solid lines and shaded region represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. The stronger the wind, the lower the probability of the galaxy developing a bar, except for the no-wind galaxy. The build-up of a well-defined massive cold disc is delayed for stronger galactic wind models. No significant change is seen in the black hole variation models.

The morphology of the galaxies in the NW and WW models has been well-defined since early times. For instance, the galaxies in the no wind model have already assembled a substantial disc by z_disc = 3 and are also the most massive galaxies in comparison to the other models. When we move to later times, the massive disc component is delayed when the wind intensity generated by SN feedback is greater, whereas the bulge forms earlier and is less massive in comparison to the disc. It should be noted that the bar formation seems to be linked to the earlier build-up of the massive disc (see Zana et al. 2019; Rosas-Guevara et al. 2020; Izquierdo-Villalba et al. 2022; Khoperskov et al. 2024). However, galaxies do not form a stable bar in the no wind model. Conversely, galaxies in the strong wind models developed a substantial disc that was later formed, but it never developed a stable bar. In the BH variations models, the bar is present in all simulations, but the bar formation time varies dramatically. Additionally, we can observe some variations in the length and strength of the bars for galaxies across the various models. The BHlowEff variation, in particular, displays the bars with the highest length and strength.

Fig. 4 illustrates the evolution of the disc and bulge masses and their sizes in the galactic wind and the quasar BH physics models. The galaxies in all models possess a massive disc which is larger than ${10}^{10}{\mathrm{M}}_{\odot }$ and becomes dominant, but their build-up is distinct: the disc build-up is slower as the winds from SN feedback become stronger. In order to quantify this, we calculate the formation time (redshift) of the massive disc by determining the time at which it reaches 50% of its z = 0 mass. This is depicted as vertical dotted lines in the figures, with the 20th and 80th percentiles for all objects. These measured times (redshifts), $t_{50,\mathrm{disc}}(z_{50,\mathrm{disc}})$, show that the disc formation occurs later for galaxies in the strong wind model ($z_{50,\mathrm{disc}}=0.8$) and earlier for galaxies in the no wind model ($z_{50,\mathrm{disc}}=2.25$). We note that the scatter of $z_{50,\mathrm{disc}}$ in the variations of the BH model is higher than that in the TN50-like model, oscillating between z = 0.8 and z = 1.3 in the case of galaxies in the no black holes model. However, the growth of the disc is not significantly affected.

Fig. 4. Mass and size evolution of discs and bulges. Galactic winds model variations (top) and quasar BH physics model variations (bottom). The last column shows the original TNG50. Solid lines and shaded regions represent median values and the 20th to 80th percentile distribution. Thinner lines correspond to the evolution of each galaxy. Vertical solid lines denote bar formation median redshift, and shaded regions its scatter. Arrows with small vertical lines correspond to the median redshift formation of the disc/bulge ($z_{50,\mathrm{disc}}$/zb) and its scatter. Stronger SN feedback delays disc buildup and makes discs more extended.

The second row of Fig. 4 illustrates the disc length. This generally increases with increasing time for all the models; however, the final disc length increases as the intensity of the galactic winds increases. The most extended and least dense discs are those in the strong wind model, while the most compact and dense are in the no wind model. This is important because these properties play a role in the bar formation criterion, as we show in the next section.

The bulge growth for variations in the wind and BH models is illustrated in the third row of Fig. 4. The figures demonstrate that the bulge evolution is more complicated, particularly in light of the wind model variations. Additionally, it is important to observe that the fluctuations between galaxies in each model are greater. In the case of NW and WW models and TNG50 galaxies, the bulge forms later in time than the disc, suggesting that the stellar bar may facilitate the expansion of the bulges.

The SN feedback intensity is high in models with moderate (MW) and (SW) strong winds, resulting in a more complex bulge assembly, subject to significant fluctuations. Initially, the bulge experiences an increase in mass as time progresses, reaching up to 50% of its current mass even earlier. Subsequently, the bulge mass decreases over time. We track the stellar particles that are part of the bulge at the moment of their mass peak. We find that a fraction of the stellar particles that were not part of the bulge later in time present higher velocities and at larger distances, associated with the halo component. This might be due to a fast ejection of matter.

The bulge size is presented in the bottom panels of the figure. As the bulge mass, their evolution is more complicated, and it is more affected by the wind model variations. The bulge size increases with time for NW, WW, and TNG50-like models; the MW and SW models are more complicated and show similar fluctuations, as seen in the bulge mass attributed to the intensity of SN feedback. The most compact and massive bulges are those in the no wind model. The most extended and least massive bulges are found in the strong wind models.

In the strongly barred galaxies, discs and bulges are formed at a more rapid pace than those in unbarred galaxies. Our results suggest that gas is swiftly consumed to form stars in galaxies in the models with NW, WW, and TNG50-like models, resulting in more dense discs and bulges. In contrast, galaxies in models with medium and strong feedback discharge gas at earlier times as a result of the less intense gravitational potential of the halo and the intense SN feedback. This would have an impact on the properties of the gas and the formation of stars from early times. This is in agreement with the result from Bi et al. (2022a). The authors, using zoom-in simulations with a mass resolution compatible with the TNG50 simulation, find that SN feedback impacts the cold accretion flow (see their Figure 1) in high-redshift Milky Way galaxies.

These galaxies properties have been used to assess the stability of a disc against gravitational instabilities, such as bars, which we study in the next section.

4. Bar formation

We study whether the galaxies with different galactic winds and quasar BH physics model variations satisfy bar instability criteria, which are contingent upon the structural parameters. In particular, we focus on three criteria: the Toomre 1964 (Eq. (1)), Efstathiou et al. 1982, (ELN) (Eq. (2)), and Mo et al. 1998, (MMW) (Eq. (3)) criteria.

The evolution of the parameters ϵ_ELN, λ_MMW and Q_T is illustrated in Fig. 5 for the galactic wind models (top figure) and BH physics models (bottom figure) models. Except for the no wind model, the parameters ϵ_ELN, λ_MMW roughly determine whether the disc will be stable to bar formation. In particular, for WW and TNG50-like models, both criteria have been satisfied since the time of bar formation, since their values lie in the magenta-shaded region below the horizontal dotted line in the figure. In the case of the Toomre criterion, the Toomre parameters present a higher scatter. However, the parameter Q_T is roughly able to determine if the galaxies are stable to bar formation. The exceptions are galaxies in the no wind and weak wind models.

Fig. 5. Evolution of relevant parameters for bar formation. From top to bottom, the ELN criterion (Eq. (2)), MMW criterion (Eq. (3)), and Toomre parameter (Eq. (1)) for different galactic wind models (top figure) and BH physics models (bottom figure) and the original TNG50 are shown. Horizontal dashed lines represent the threshold that defines if a disc is unstable, and magenta shaded region represents the values where the disc is unstable. Thinner lines correspond to the evolution of each halo. Vertical solid lines and shaded region represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. For TNG50-like, MW and SW models, the three instability criteria determine if a bar forms or not. For NW model, ELN and MMW criteria fail to predict a stable disc, whereas Toomre criterion fails to predict bar formation for WW model.

4.1. ELN criterion

In this section, we focus on the evolution of the disc mass and size evolution and the maximum circular velocity (V_max) of the dark matter halo that concerns the ELN criterion Efstathiou et al. 1982, (ELN). In the last section, in Fig. 4, we have studied the mass and size growth of the disc for the different galactic wind and BH model variations. We observe that the evolution of the disc scale length and the disc mass have been affected by the strength of the galactic winds (kinetic energy released into the ISM) since early times. The most extreme case is the variation where there is no galactic wind or, opposite, the galaxy experiences a strong galactic wind. In both NW and SW models, the galaxy forms a massive disc ($M_{\ast }>{10}^{10}{\mathrm{M}}_{\odot }$) and a relatively small bulge (B/T < 0.4, see Table 3), and neither of both models forms a stable bar. The other models show an intermediate effect between these two examples, and only the galaxies in the TNG50-like and weak wind models form a bar. In contrast, the evolution of the maximum circular velocity (V_max) of the dark matter halo almost remains constant across all scenarios, as seen in Fig. 6. There is no distinction when the dark matter haloes are the sole factor considered. Therefore, the size of the disc is a determinant of the ELN criterion and, with less significance, the stellar mass disc. The ELN criterion fails to predict the disc stability of the galaxies in the no wind model, while the ELN criterion is met by the galaxy in the WW and TNG50-like models, at least during the time of bar formation.

Fig. 6. Median evolution of the maximum circular velocity of the dark matter halo. Each column represents the variations of galactic winds (upper Fig.) and BH physics models (bottom Fig.) in order from left to right, and vertical solid lines represent the bar formation redshift. The original galaxy in TNG50 is represented by the final right column. Thinner lines correspond to the evolution of each halo. Vertical solid lines and shaded region represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. No significant variation is found in the evolution of vmax among the different models.

The BH physics variations show some differences in disc size evolution, but not in disc mass. However, the values and evolution of the ϵ_ELN parameter remain consistent, as seen in Fig. 5. The ELN criterion emphasises that a compact, dense disc is a prerequisite for forming a bar; however, our findings show this condition alone is not sufficient.

4.2. MMW criterion

In order to conduct a more thorough investigation of the MMW criterion (Mo et al. 1998), the middle panels of Fig. 5 show the values of the parameter λ_MMW, and the magenta shaded region corresponds to the values of λ_crit as seen in Eq. (3). For details in the calculation of λ_crit, see Eqs. (23), (32), (34) in Mo et al. (1998). Therefore, if (λ_MMW/λ_crit) > 1, the disc is approximately stable. For the galaxies in the models that form a bar, the MMW criterion correctly determines that the disc is prone to form a bar since the parameter (λ_MMW/λ_crit)≲1² lies in the magenta-shaded region, especially at the time of bar formation. An exception are the galaxies in the models without SN feedback, where the MMW criterion fails to predict that the disc is stable since a bar does not form. On the other hand, the MMW criterion is satisfied for galaxies in the models with medium and strong winds since they do not form a bar and have (λ_MMW/λ_crit) > 1.

To understand why the criterion fails in some cases and others not, we follow the evolution of the galaxy properties that the parameter λ_MMW depends on: the spin parameter (λ), the ratio of the angular momentum of the disc and the halo (J_disc/J_h), and the ratio of the disc mass and the dm halo mass (M_disc/M_h). The evolution of these galaxy properties is presented in Fig. 7. The evolution of the spin parameter remains almost constant regardless of the galactic wind model or the varied BH physics, with values between 0.02 and 0.04, with the exception of the no wind model. This model presents spin parameters smaller than 0.03. On average, the spin parameter in all the models presents low values in agreement with those studied in N-body simulations in favour of the formation of a bar (e.g. Long et al. 2014). We observe the spin parameter to be almost unaffected by SN feedback and BH physics model variations. However, J_disc/J_h is subject to change as the intensity of the galactic winds is modified. This ratio increases over time, transitioning from values less than 0.02 at z = 2 to values just smaller than 0.05 at z = 0, except for the galaxies in the no wind model. These galaxies present high ratios and reach the highest value (J_disc/J_h ∼ 0.06 at z = 0). The galaxies with the lowest values are those in the medium and strong wind models, nearly reaching 0.04 at z = 0. This is expected since SN feedback allows for gradual gas cooling and accretion from a hot and diffuse gas halo that, in some cases, galaxies might form stable discs (e.g. Sales et al. 2012). Similarly, M_disc/M_h is affected by SN feedback with the highest values presented for the no wind model (≳0.05) and the lowest values for the strong wind models (≲0.02) at z = 0. Apart from the no-wind model, this ratio increases with time. Strong feedback can expel baryonic matter from galaxies since early times, cycling it through a hot and diffuse gas before it can be re-accreted and converted into stars (e.g. Sales et al. 2012). We note that for the WW, TNG50-like models, where galaxies form a bar, present higher values of high angular momentum and higher disc content with respect to the halo. This result agrees with the ELN criterion that a massive disc is a necessary condition for a disc to form a bar, but not sufficient. A good example is the galaxies in the NW model. In this case, the parameter (λ_MMW/λ_crit) determines that the discs should form a bar that is not there. As we show in Sect. 5.1, both the ELN and MMW criteria do not take into account the effect of the presence of a large and compact bulge.

Fig. 7. Median evolution of spin parameter, the ratio of the angular momentum of the disc over the halo, and the disc mass fraction (solid thicker lines). From left to right, different galactic wind models (top Fig.) and BH physics models (bottom Fig.) are shown. The rightmost column corresponds to the original galaxy in TNG50. Thinner solid lines correspond to individual haloes. Vertical solid lines and shaded regions represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. The angular momentum and mass of the disc with respect to the halo are smaller with stronger SN feedback. As the ELN criterion, the MMW criterion demonstrates the essential condition of the early buildup of a massive disc to favour bar formation. However, it is insufficient as it fails to predict stable discs in the NW model.

4.3. Toomre criterion and swing parameter

As discussed in the introduction, the Toomre (1964, Q_T, see Eq. (1)) is one of the criteria for bar instability. This parameter is responsible for disc stability in axisymmetric density perturbations. If Q_T < 1, the disc is unstable. In conjunction with the Toomre parameter, the swing parameter (X) evaluates the amplification of the density perturbations through the swing amplification effect. This parameter is defined by X = rκ²/(4πG)Σ, where κ is the epicyclic frequency, G is the gravitational constant, and Σ is the stellar surface density. Both parameters Q_T and X depend on the kinematic properties of the disc as a function of the radius. The Toomre parameter should be approximately 1, but it should be at most 2 to achieve substantial amplification (e.g. Binney & Tremaine 2008). To the extent that Q_T is excessively high, the disc is too stable to generate large-scale perturbations. In addition, a value in the range 1 ≤ X ≤ 3 suggests a regime in which swing amplification can substantially amplify non-axisymmetric perturbations, such as spirals or bars.

The bottom panels of Fig. 5 display the average Q_T and X for the various models within a stellar half-mass radius and a z-slide of 2 kpc, which is comparable to the thickness of a disc (refer to Table 3). We can observe that the values of Q_T and X are located in the bar instability region in the TNG50-like and BH variation models, particularly during the bar formation phase. It is interesting to note that the values of Q_T are just above 2 at later times. This is in agreement with the values presented in recent observations, which show that some local barred discs can reach values of Q_T of ≈3 (e.g. Romeo et al. 2023). The swing parameter, X, is located in the instability regions of galaxies for all the models, especially in the medium and strong wind models, while Q_T is located above the instability region with values larger than 3 in both the medium and strong wind models (see Cols. 4 and 5 of top figure of Fig. 5). The high values of Q_T determine that the disc is highly stable against bar instabilities. Interestingly, the Toomre parameter is substantially higher in the no/weak wind model for galaxies (Q_T > 4) than in the other models. We focus on Fig. 8 to better understand this behaviour. We note that the galaxies in the no/weak wind models have higher κ values. They also present higher stellar surface densities, which could be used to offset the highest values in κ. Consequently, we obtained similar values in X at z ≲ 2. However, the radial velocity dispersion, which is significantly influenced by the SN feedback, is the primary cause of the discrepancy in Q_T between the no/weak wind models and the other models. In addition, this could be the reason for the opposite predictions between this instability criterion, which predicts stable bars for galaxies in NW and WW models, and the ELN and MMW criteria, which predict unstable discs. Even more galaxies in the NW model do not form a stable bar, whereas galaxies in WW galaxies do. The origin of the discrepancy may be because the ELN and MMW criteria do not consider σ_R (Eqs. (2) & (3)) while Q_T does depend on it (Eq. (1)), indicating that the Toomre criterion is sensitive to the precise value of σ_R. We discuss this in the next section.

Fig. 8. Evolution of radial velocity dispersion, stellar surface density, epicyclic, and angular frequency in a stellar half-mass radius aperture. From left to right, the different galactic wind models (top figure) and BH physics models (bottom figure) are shown. The last column on the right corresponds to the original galaxy in TNG50. Thinner solid lines correspond to individual haloes. Vertical solid lines and shaded regions represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. The three galaxy properties are highly affected by the SN feedback process, especially the radial velocity dispersion. The highest velocity dispersion is presented in the NW and WW models, causing the Toomre criterion to determine stable discs in contrast to the predictions of ELN and MMW criteria.

5. Discussion

In this section, we discuss our findings with respect to previous works and future prospects.

5.1. Effect of the central bulge on bar formation

The ELN and MWW criteria predict that the disc is unstable to bar-like perturbations for galaxies in the no wind model; however, a stable bar does not form. These galaxies exhibit the highest stellar mass content in the disc since early times. For example, the findings of Bonoli et al. (2016) and Spinoso et al. (2017), which analyse the twin zoom-in simulations of the Milky Way, ErisBH and Eris, suggest that ErisBH forms a bar as a result of AGN feedback, which diminishes the bulge size in comparison to Eris. This enables the activation of bar instabilities. In addition, Saha & Elmegreen (2018) demonstrate that cold stellar discs surrounded by an initial relative dense bulge (ρ_bulge/ρ_disc > 1 for compact bulges) do not form a bar even though the disc was susceptible to bar instabilities (lower values of the Toomre parameter). Similarly, Kataria & Das (2018) found that massive and more concentrated bulges can inhibit bar formation and growth, particularly when the bulge-to-disc mass fraction is higher than 0.3 for concentrated bulges and even higher (0.6) for less concentrated bulges. This was achieved using N-body simulations.

Fig. 9 illustrates the bulge-to-total-mass fraction (B/D) as a function of bulge-to-disc sizes ratio and colour-coded by redshift. Galaxies in the no wind model show an increase of B/D with values higher than 0.3 and bulge-to-disc size ratios smaller than 0.5 after z = 1. The upper limit suggested by Kataria & Das (2018) for compact bulges to prevent bar formation is B/D = 0.3. Their simulations indicate that bars can only form in disc galaxies when the radial force due to the bulge is less than 0.35 in comparison to the total force at the disc scale length. This is an outcome of the bulge’s increasing mass, which increases the velocity dispersion of the disc and the radial force. At high redshift (> 1.5), the galaxies in the no-wind model, generally exhibit B/D < 0.3; however, the bulge-to-disc size ratios take values higher than 0.5 (less compact bulges). Consequently, the disc is unstable in principle, but the non-axisymmetric instabilities are unable to grow (see Fig. 3) since the bulge mass continues to increase constantly (see Fig. 4).

Fig. 9. Bulge-to-disc mass fractions as a function of bulge-to-disc size ratios. Left to right, galactic wind models (top Fig) and BH physics model variations (bottom Fig) are shown. The colour code represents the redshift. Horizontal grey dashed lines represent the limit given by Kataria & Das (2018) for dense bulges that do not favour the formation of a bar, even if the disc could be unstable. Galaxies in the NW model present a large and compact bulge since early times, preventing the disc from forming a bar.

At the time of bar formation (z_bf ∼ 2) and until z ≈ 0.5, the galaxies in the WW model have bulge-to-disc size ratios that are less than 0.5 and their B/D values are less than 0.3. This suggests the disc is still capable of forming a bar even though the galaxies in the WW model have concentrated bulges. At low redshift (z < 0.5), the bulge continues to grow, but the bar length does not increase to the same extent as in the TNG50-like and BH variations models. This may suggest that the bulge could be influencing bar evolution. In contrast, the medium (MW) and strong (SW) wind models exhibit galaxies that show similar behaviour to those in the TNG50-like model, with B/D < 0.3 at low redshift. The galaxies present higher B/D > 0.3, and the bulge-to-disc size ratio takes values higher than 0.5 at higher redshift. Despite the fact that the bulge may have a minor effect, the discs in the MW and SW models do not form a bar, and according to the three instability criteria, these discs are stable against bar formation in general.

5.2. Impact of feedback processes on bar formation

As seen in Section 3, feedback processes regulate star formation and affect the galaxy build-up and its components. In particular, stronger SN feedback delays the assembly of the disc (Fig. 4). As a result, discs exhibit less mass (Fig. 7), lower radial velocity dispersion (Fig. 8) and larger sizes (Fig. 4), as the SN feedback strength increases. These properties are tightly related to determining whether a disc is stable or not (Eqs. (1), (2) and (3)). In particular, galaxies in the strong wind model experience more gradual gas cooling and accretion from a hot and diffuse gas halo, than galaxies in the other models (e.g. Sales et al. 2012). As a result, galaxies form stable discs. In the case of the galaxies in the No-wind model, insufficient feedback may lead to runaway star formation and rapid consumption of gas, leading to the rapid formation of a massive compact bulge (e.g. Navarro & Benz 1991). This has repercussions in the conditions that promote bar formation.

We also note that the galaxies with weak or no-SN feedback are not realistic, and we have only studied them to gain an understanding of the effects of feedback processes on bar formation, but we already know that these models predict galaxies that overestimate the baryon content of galaxies (see Pillepich et al. 2018a).

Varying quasar BH physics models does not significantly affect the time of disc or the bulge formation (Fig. 4). As a result, the mass (Fig. 6), size (Fig. 4) and radial dispersion velocity (Fig. 8) are pretty similar. However, the scatter in the time of the formation of the bulge and bar is higher as the efficiency of AGN feedback in quasar mode is lower (Fig. 4). This could be because this mode of AGN feedback affects the central gas of the galaxy. Indeed, Irodotou et al. (2022), which studies a galaxy of the Auriga simulations, has found that bars form with/without quasar mode feedback. This mode of AGN feedback can affect the final bar properties of the bar.

Our findings suggest that the formation of bars may serve as a constraint on the subgrid feedback models of galaxy formation, thereby contributing significantly to our understanding of galaxy evolution. We caution, however, that our galaxy sample, as well as the variations of the subgrid physics, are limited. Therefore, bar formation could be affected by other factors beyond our study, such as galactic surroundings (e.g. Peschken & Łokas 2019) which, foster bar formation under given conditions.

5.3. Lessons learned from disc instability criteria

As mentioned in Section 4, most galaxies among the different models studied here satisfy the three criteria of bar instability. An exception is the galaxies in the NW and WW models. Galaxies in the NW model do not form a bar, while galaxies in the WW model do. Nevertheless, the Toomre criterion predicts stable discs, whereas the ELN and WWM criteria forecast unstable discs. It is possible that a missing condition is required for bar formation, as the ELN and WWM criteria have opposite predictions to the Toomre criteria. Indeed, the galaxies in NW and WW models have formed the most massive compact bulges of the models studied here, and these components may influence bar formation (Fig. 9, Section 5.1). The effect of the bulge is not taken into account in the ELN and WWM criteria. This has been studied by Athanassoula (2008) who find that the ELN criterion does not take into account the effect of the central concentration to stabilise the disc. Kataria & Das (2018) using N-body simulations propose a limit on bar formation in discs with bulges by assessing the ratio of radial force due to bulge and disc components at disc scale length. This is because the velocity dispersion of disc stars becomes large enough to prevent bar-type instability. Similarly, the Toomre parameter also fails to predict a bar in the presence of a massive bulge in galaxies in the weak wind model. The Toomre parameter does not adequately account for the varying effects of different bulge masses (compact versus less dense) and their interaction with the halo, which can lead to scenarios where bars form or do not form despite being predicted as unstable or stable discs based on a parameter alone. When a bulge is strong enough, also a disc may not form a bar due to non-linear responses that are not captured by the Toomre criterion. For instance, Jang & Kim (2023), using isolated simulation, propose a new condition that includes the central mass concentration (CMC) together with the Toomre parameter to predict if the disc is stable or not.

In addition to this, more ingredients have not been investigated in this study due to the nature of our sample, which comprises galaxies that exhibit bar formation in TNG50 at high redshift (z > 1.5) and subsequently have low halo spin parameters. For instance, Izquierdo-Villalba et al. (2022) find that a small percentage of barred galaxies in TNG50 do not satisfy the ELN criterion. The potential reason for this is that they had a close encounter with a galaxy during the bar formation process (see also Rosas-Guevara et al. 2024). Conversely, other studies have found that close encounters do not promote bar formation, but rather delay bar formation (e.g. Zana et al. 2018). Algorry et al. (2017) also find that the ELN criterion is not sufficient to predict stable bars in the EAGLE simulation and proposed an additional parameter which evaluates the dynamical contribution of the baryons at a half mass radius scale with respect to the global dynamics of the dark matter halo (maximum circular velocity of the halo).

In conclusion, there is no consensus on the importance and the limits of secondary conditions needed to promote/prevent bar formation, maybe reflecting the complex processes that are involved in bar formation.

5.4. Effect of stochasticity on feedback processes in bar formation

Something to be aware of when trying to understand our results is to take into account the stochasticity provided by subgrid physics in simulations. In particular, Genel et al. (2019) found that simulations with feedback mechanisms could exhibit a stronger butterfly effect than those without feedback. This suggests that the dynamics at galactic scales could significantly influence how small differences in initial conditions can evolve into larger discrepancies in galaxy properties (see also Borrow et al. 2023 for other galaxy formation models). It is essential to acknowledge that the morphology of the models under investigation here is not significantly affected; however, the radial velocity dispersion and the extent of the disc could be. These qualities are employed to determine the stability of a disc, as seen in Figs. 4 and 8.

The TNG50 simulation predicted the presence of six strongly barred galaxies; however, only five of these galaxies present a strong bar in the TNG50-like model at z = 0. For this particular galaxy, we note that the original TNG50 galaxy has a strong bar with the weakest strength among the six galaxies studied here ($A_{2,\max }0.38$, see Table 1) at z = 0. In the TNG50-like model, the galaxy forms a bar that appears to weaken over time, as it evolves. We suspect that the cause of this weakening is the fact that the galaxy in the TNG50-like model has higher gas fractions in the centre than the original TNG50 galaxy. This could be because the black hole is less massive, and thus a smaller amount of energy is released in radio mode. This is an example of how the stochasticity of subgrid physics could affect bar formation, but in general we find minor variations in the characteristics of strong bars in our study. On the other hand, when the feedback mechanisms differ, the merging histories are very similar to the original TNG50. The stochastic uncertainty that is incorporated into simulations limits the physical conclusions that can be derived. The broad causal tendencies we have identified suggest that this constraint is of lesser significance.

Finally, we emphasise that our study is limited to five galaxies with a strong bar, which does not allow us to make general conclusions regarding the effects of the stochasticity of feedback processes, so a more comprehensive study in this matter using larger cosmological simulations with a larger galaxy statistical sample is needed for future projects, especially for the weak bars.

5.5. Outlook

The zoom-in hydrodynamic simulations enabled us to investigate the influence of feedback processes with the advantage of previous works using the TNG50 simulation. We re-simulated six Milky Way-like galaxies from the TNG50 simulations. These galaxies have a stable galactic bar at z = 0 which was formed between z = 3 and z = 1.5 (more than 8 gigayears ago) and with quiet merger histories and isolated haloes. Even so, it is crucial to recognise the limitations of this investigation. The physics above scales of approximately 288 pc can be resolved by these zoom-ins. Although this is sufficient for a comprehensive examination of the formation and evolution of bars within the most massive galaxies in the local universe, simulations with a higher resolution are required to correctly track the initial phases of bar build-up at high redshifts. Recent observations from JWST have shown that bars exist at high redshifts (e.g. Guo et al. 2023; Costantin et al. 2023; Guo et al. 2024).

The conditions that galaxies are subject to at such high redshifts are notably distinct from lower redshifts. It is anticipated that the frequency of encounters and mergers is higher at such early times (e.g. Conselice et al. 2008; Bi et al. 2022a,b). A cosmological context would be an interesting setting in which to conduct this experiment with an extensive statistical galaxy sample and in different environments at higher redshifts.

It is thought that bars drive flows into the galaxy centre, which in turn affects central star formation or the accretion onto the supermassive black hole (e.g. Fanali et al. 2015; Spinoso et al. 2017). The fate of the gas in the nuclear zone, affected by the bar and other physical processes like SN and AGN feedback, should be studied. In particular, it would be interesting to see gaseous bars formed by stellar bar torques (e.g. Englmaier & Shlosman 2004).

The effects of the star formation processes on the bar properties, especially the bar pattern speed, which is the least constrained, need to be studied. The bar pattern speed is most affected by the interaction between the bar, disc, and the central dark matter distribution (e.g. Cuomo et al. 2020). The main effect is a slowdown of the bar as a function of time. This seems to be in contradiction with observations, which have found that bar pattern speeds are rapid. Using the zoom-in simulations, we can see the effect of the star formation processes at high redshift that affect the initial bar pattern speed when the bar forms (e.g. Semczuk et al. 2024).

6. Summary

In this work, we have investigated the impact of feedback processes on the formation and evolution of bars in galaxies similar to the Milky Way. We studied a suite of zoom-in cosmological hydrodynamical simulations of Milky Way-sized haloes selected from the TNG50 simulations in which the original galaxies form a strong bar between z = 1.5 and z = 3 In particular, we explored variations of the SN feedback and BH physics models within the TNG galaxy formation model. Specifically, our investigation focused on five different intensity levels within the galactic wind models (energy injected into the galactic wind per SN event) and two models of BH physics affecting the AGN feedback in the quasar mode regime. The results of our investigation are as follows:

• Our analysis revealed that the feedback from SN impacts the formation of a bar, while the effect of the model variations of the quasar mode BH physics studied here is negligible (Figs. 1 and 3). However, BH physics can affect the ultimate characteristics of a bar, such as strength and length. In order to gain a deeper understanding, we examined the characteristics of the various components of the galaxy.

• We observed that the morphology at low redshift remained consistent in the sense that in all the models, galaxies developed a massive disc ($>{10}^{10}{\mathrm{M}}_{\odot }$), which is dominant in comparison with the bulge mass (Fig. 3).

• Our findings indicate that SN feedback has a significant impact not only on the stellar mass content but also on the structural properties of the disc and bulge, including the disc size, thickness, and bulge size (see Table 3 and Fig. 4). SN feedback also affects the formation time of the disc and bulge and then bar formation, which aligns with the finding by Zana et al. (2019). In extreme cases where the SN feedback is strong, bar formation is stopped. Our analysis revealed that the bulge and disc in galaxies lacking SN feedback are noticeably more compact and concentrated than those with stronger feedback at z = 0. In such galaxies, bar formation also stops.

• Based on these findings, it is evident that the feedback from SNe plays a crucial role in shaping the various components found within galaxies. When the feedback strength is increased, it causes a decrease in the rate at which the disc grows. As a result, this affects the (thin) disc and bulge characteristics. The models that vary the BH physics in quasar mode do not influence the disc or the bulge formation, as also seen in Irodotou et al. (2022) using the Auriga simulations.

• We thoroughly reviewed the bar instability criteria in our different models, particularly during the process of bar formation. Specifically, we focused on three criteria: the Toomre 1964, Efstathiou et al. (1982, ELN), and Mo et al. (1998, MMW) criteria. Our analysis revealed that the sizes of the discs (Fig. 4), the stellar mass of the disc over halo mass (Fig. 7), and the radial velocity dispersion (Fig. 8) are strongly influenced by SN feedback and used to determine if a disc is prone to instability and bar formation (Fig. 5). Most of the models satisfy the instability criteria at the moment of bar formation, as in the barred galaxies in TNG50 (Izquierdo-Villalba et al. 2022). The exception is the no wind model, which presents a higher radial velocity dispersion, leading to high values of the Toomre parameter, whereas the ELT and MMW criteria predict an unstable disc.

• In galaxies in the weak winds and TNG50-like models, discs exhibit instability and give rise to bar formation, while galaxies in the strong and medium wind models do not form a bar. On the other hand, galaxies in the no wind model do not form a bar even though the disc satisfies some instability criteria. This could be explained by the presence of a massive compact bulge (see Fig. 9). This is in agreement with Kataria & Das (2018), who found that galaxies with B/D > 0.3 for concentrated bulges (B/D > 0.6 for less concentrated bulges) prevent bar formation.

The comparison of galaxy formation models in our study offers a deeper understanding of the relationship between the emergent structures of galaxies and the input physics assumptions during their formation and evolution. This understanding could inspire new approaches to the prescription of galaxy formation models, thus offering hope for further advancements in the field. Our research can also be employed to evaluate the predictive capabilities of analytic methods employed to identify the presence of a bar. This would enhance our comprehension of bar formation and enhance the modelling of barred galaxies in more empirical models.

Acknowledgments

The authors thank Volker Springel for providing them with the AREPO code and the TNG model and his suggestions regarding the project. The authors would like to thank the anonymous referee for their comments which improved the quality of this paper. The authors thankfully acknowledge the computer resources at MareNostrum and the technical support provided by Barcelona Supercomputing Centre (RES-AECT-2023-2-0015). YRG acknowledges the support of the“Juan de la Cierva Incorporacion” Fellowship (IJC2019-041131-I). S.B. & YRG acknowledge support from the Spanish Ministerio de Ciencia e Innovación through project PID2021-124243NB-C21. SC acknowledges the support of the ‘Juan de la Cierva Incorporacíon’ fellowship (IJC2020-045705-I). S.B., SC & YRG acknowledge support from the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement number 101086388 (LACEGAL). Technical and human support provided by the DIPC Supercomputing Centre is gratefully acknowledged.

References

Appendix A: Kinematic decomposition and bar identification

In order to identify more precise galaxy components, we adopt the kinematic decomposition algorithm MORDOR³ presented in Zana et al. (2022) to determine more specific galaxy components. The decomposition is predicated on the circularity (ϵ) and binding energy (E) phase space, where a minimum in E is identified for each galaxy, E_cut. All four components listed below were identified for the simulations employed in the paper:

• Classical Bulge: is defined as the set of stellar particles that exhibit the highest degree of binding, with a value that is contingent upon the E_cut of each galaxy. These particles also exhibit counterrotation, which is characterised by a negative value of ϵ (< 0). Next, Monte Carlo sampling is employed to select a distribution that is equivalent to the component of the distribution that exhibits positive circularity and allocate it to the Bulge.

• Thin/Cold disc: is defined as the set of stellar particles that have the highest degree of binding and are not allocated to the Bulge, but have positive values of ϵ (ϵ > 0.7).

• Pseudobulge: is defined as the remaining stellar particles that exhibit a high degree of binding without being assigned to the Bulge or the thin disc.

• Thick/Warm disc: is defined as the set of stellar particles that have positive ϵ values (ϵ > 0.7) and are not assigned to the bulge/pseudobulge, despite exhibiting a reduced degree of binding.

The stellar density maps components of a galaxy are illustrated in Fig. A.1 for the TNG50-like model, in which the galaxy forms a bar, and for the strong wind model, in which the galaxy does not form a bar. MORDOR is employed to generate these maps. The distinctions between all of the components are readily apparent. For this example, the thick disc seems to trace the bar, while the pseudobulge appears to be indistinguishable from the bulge. This also seems to be the case at the time of bar formation.

Fig. A.1. Stellar density maps of the galaxy components identified by MORDOR. From left to right, the columns show the total galaxy, thin disc, thick disc, pseudobulge, bulge, and halo. Top panels correspond to a galaxy in the TNG50-like model. The bottom panels correspond to the same galaxy in the strong wind model.

A.1. Identification of bars

The face-on stellar surface density is Fourier decomposed to identify non-axisymmetric structures (e.g. Athanassoula & Misiriotis 2002; Zana et al. 2018; Rosas-Guevara et al. 2020). We concentrate on A₂(R), which is defined by the ratio between the second and zero terms of the Fourier expansion and its phase Φ(R) (refer to equations 1 and 2 in RG22). A bar structure is characterised by A₂(R) and Φ(R), with the bar strength defined as A_2, max, the peak of A₂(R). Inside the bar, the phase should be constant. To establish a constant phase, we calculate the standard deviation (σ) of Φ(R) every time a new cylinder shell is added and set σ ≤ 0.1. The bar extent (r_bar) is the greatest radius where σ ≤ 0.1, and the A₂ profile first dips below 0.15 or the minimum value of A₂(R). Large values of A₂(R) may be attributable to transient events like mergers and interactions. Thus, we assumed the bar represents a long-lasting characteristic if the following criteria were met: the maximum of A₂, A_2, max is greater than 0.2; r_bar > r_min, where r_min = 1.38 × ϵ_*, z is an imposed minimum radius and ϵ_*, z corresponds to the proper softening length for stellar particles; and the estimated age of the bar is larger than the time between the analysed output and two previous simulation outputs. Although this filter may omit newer bars at lower redshifts, it aligns with the bar selection criteria employed by RG20; RG22 when the bar structure was found in previous outputs.;

To identify a bar structure, the formation time, t_bar, must also be determined. We followed the growth of A_2, max up to the moment when A_2, max ≤ 0.2 and r_bar(t_bar)≤r_min for several snapshots. Additionally, the difference between A_2, max at a snapshot and the two prior snapshots must not exceed 0.45 during this time frame. This confirms that the bars we find are stable.

A.2. Method of three-component decomposition of surface face-on density profile

In this section, we outline the process of decomposing the surface brightness profiles of the galaxies to determine the length scale of the disc, the effective radius of the bulge, and the properties of the bar when it is identified through Fourier decomposition. The surface density profiles are calculated in a face-on perspective within concentric annuli that are 0.16 dex in width and centred on the stellar component minimal potential.

We employ the particle swarm optimisation (PSO) code PSOBacco ⁴ (Aricò et al. 2021) to determine the minimumχ² of a fit of the sum of an exponential profile and two Sersic profiles (one of which corresponds to the bar). This value is expressed in terms of the stellar mass as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Sigma }(R)& ={\mathrm{\Sigma }}_{\mathrm{d},0}\exp (-R/h_{\mathrm{disc},\ast })+{\mathrm{\Sigma }}_{\mathrm{b},\mathrm{e}}\exp (-b_n[{(r/r_{\mathrm{b},\mathrm{eff}})}^{1/n}-1])\hfill \\ \hfill & +{\mathrm{\Sigma }}_{\mathrm{bar},\mathrm{e}}\exp (-b_n[{(r/r_{\mathrm{bar},\mathrm{eff}})}^{1/n_{\mathrm{bar}}}-1])\hfill \end{array},\end{array}$$(A.1)

where Σ_d, 0 is the central surface density of the disc component, h_disc, * is the disc scale length, $r_{b,\mathrm{eff}}$ (r_bar) the effective radius that encloses half of the stellar mass of the one predicted from the Sersic profile, Σ_b, e (Σ_bar, e) corresponds to the surface density at $r_{b,\mathrm{eff}}$, and n is the Sersic index (of the bulge/bar). The value of b_n is contingent upon the complete gamma function, as Γ(2n) = 2γ(2n, b_n).

We apply the method of fitting three components in the region at radii smaller than $r_{\mathrm{fit}}=\max ({\log }_{10}(\mathrm{\Sigma }(\mathrm{r})/{\mathrm{\Sigma }}_{\max })\ge -2.6,10\mathrm{kpc})$ when a bar is detected by Fourier decomposition. If not, we fit two components. This condition guarantees to suit both high-redshift discs, which are more compact (smaller disc scale lengths for a given stellar mass), and their analogues at low redshifts, whose discs are less compact (larger disc scale lengths for a given stellar mass).

The surface density profiles and the fits for a galaxy in the TNG50-like model, where a bar is formed, and the same galaxy in the strong wind model, where a bar is not formed, are illustrated in Fig. A.2. A black line denotes the total fitted profile. The disc, bar, and bulge components are represented by blue, pink, and red lines, respectively. As a point of reference, the scale-length of a disc, h_disc, *, the effective radius of the bulge, r_b, eff, and the length of the bar in the case of barred galaxies are represented by coloured vertical lines. The vertical lines that are dashed represent the radius at which the fit is performed.

Fig. A.2. Face-on stellar surface density profiles in terms of R/R*, 50 where R*, 50 is the stellar half-mass radius, for a galaxy with a bar in the TNG50-like model (top) and unbarred galaxy in the strong wind model (bottom) at z = 0, 1. The profiles (black lines) were obtained by fitting simultaneously a Sersic profile for a bulge (red lines) and a Sersic for a bar (pink lines), and an exponential profile (blue lines). Purple lines represent the surface density profile from the predicted Sersic profile and exponential profile. Horizontal dashed lines correspond to the radius at which the fit is done. Dotted horizontal lines correspond to the bar length (pink), the effective radius of the bulge (red), and the scale length of the disc(blue).

All Tables

All Figures

Fig. 1. Left panel: Dark matter distribution of the TNG50 simulation. Magenta empty circles show the position of the selected haloes of the zoom-in simulations. Right panel: Mocked images in JWST NIRCam F200W, F115W, and F070W filters (face-on) calculated in Nelson et al. (2018), (TNG database) for the TNG50 galaxies that are re-simulated. The NIRCam blue channel highlights the young population of the galaxy, and the NIRCam red channel older populations. The original TNG50 galaxies present a massive disc with a smaller bulge and form a strong bar.

In the text

Fig. 2. Face-on stellar density maps of two distinct galaxies at z = 0. The top figure shows the wind model variations, and each column corresponds to a distinct wind model. The left column shows a galaxy without SN feedback, and the right columns show galaxies with stronger SN feedback. The bottom figure shows different quasar black hole models but with the same SNe feedback model as indicated by Table 2. The black star symbolises the location of the BH. In both the no wind and strong wind models, the disc galaxy does not exhibit the formation of a bar. Quasar BH physics appear to have a limited influence on the formation of a bar.

In the text

Fig. 3. Median evolution of the bar and galaxy morphology for simulations with different galactic wind (top) and BH feedback (bottom) models. Each column corresponds to a different model, as indicated in the figures. From top to bottom rows, the bar strength, bar extent, the stellar mass, and the disc-to-total and bulge-to-total mass fraction are shown. The medians and the distribution between 20th and 80th percentiles correspond to solid lines and shaded regions, respectively. Thinner solid lines correspond to the evolution of individual galaxies. Vertical solid lines and shaded region represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. The stronger the wind, the lower the probability of the galaxy developing a bar, except for the no-wind galaxy. The build-up of a well-defined massive cold disc is delayed for stronger galactic wind models. No significant change is seen in the black hole variation models.

In the text

Fig. 4. Mass and size evolution of discs and bulges. Galactic winds model variations (top) and quasar BH physics model variations (bottom). The last column shows the original TNG50. Solid lines and shaded regions represent median values and the 20th to 80th percentile distribution. Thinner lines correspond to the evolution of each galaxy. Vertical solid lines denote bar formation median redshift, and shaded regions its scatter. Arrows with small vertical lines correspond to the median redshift formation of the disc/bulge ($z_{50,\mathrm{disc}}$/zb) and its scatter. Stronger SN feedback delays disc buildup and makes discs more extended.

In the text

Fig. 5. Evolution of relevant parameters for bar formation. From top to bottom, the ELN criterion (Eq. (2)), MMW criterion (Eq. (3)), and Toomre parameter (Eq. (1)) for different galactic wind models (top figure) and BH physics models (bottom figure) and the original TNG50 are shown. Horizontal dashed lines represent the threshold that defines if a disc is unstable, and magenta shaded region represents the values where the disc is unstable. Thinner lines correspond to the evolution of each halo. Vertical solid lines and shaded region represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. For TNG50-like, MW and SW models, the three instability criteria determine if a bar forms or not. For NW model, ELN and MMW criteria fail to predict a stable disc, whereas Toomre criterion fails to predict bar formation for WW model.

In the text

Fig. 6. Median evolution of the maximum circular velocity of the dark matter halo. Each column represents the variations of galactic winds (upper Fig.) and BH physics models (bottom Fig.) in order from left to right, and vertical solid lines represent the bar formation redshift. The original galaxy in TNG50 is represented by the final right column. Thinner lines correspond to the evolution of each halo. Vertical solid lines and shaded region represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. No significant variation is found in the evolution of vmax among the different models.

In the text

Fig. 7. Median evolution of spin parameter, the ratio of the angular momentum of the disc over the halo, and the disc mass fraction (solid thicker lines). From left to right, different galactic wind models (top Fig.) and BH physics models (bottom Fig.) are shown. The rightmost column corresponds to the original galaxy in TNG50. Thinner solid lines correspond to individual haloes. Vertical solid lines and shaded regions represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. The angular momentum and mass of the disc with respect to the halo are smaller with stronger SN feedback. As the ELN criterion, the MMW criterion demonstrates the essential condition of the early buildup of a massive disc to favour bar formation. However, it is insufficient as it fails to predict stable discs in the NW model.

In the text

Fig. 8. Evolution of radial velocity dispersion, stellar surface density, epicyclic, and angular frequency in a stellar half-mass radius aperture. From left to right, the different galactic wind models (top figure) and BH physics models (bottom figure) are shown. The last column on the right corresponds to the original galaxy in TNG50. Thinner solid lines correspond to individual haloes. Vertical solid lines and shaded regions represent the median redshift and the distribution between 20th and 80th percentiles of the bar formation time. The three galaxy properties are highly affected by the SN feedback process, especially the radial velocity dispersion. The highest velocity dispersion is presented in the NW and WW models, causing the Toomre criterion to determine stable discs in contrast to the predictions of ELN and MMW criteria.

In the text

Fig. 9. Bulge-to-disc mass fractions as a function of bulge-to-disc size ratios. Left to right, galactic wind models (top Fig) and BH physics model variations (bottom Fig) are shown. The colour code represents the redshift. Horizontal grey dashed lines represent the limit given by Kataria & Das (2018) for dense bulges that do not favour the formation of a bar, even if the disc could be unstable. Galaxies in the NW model present a large and compact bulge since early times, preventing the disc from forming a bar.

In the text

	Fig. A.1. Stellar density maps of the galaxy components identified by MORDOR. From left to right, the columns show the total galaxy, thin disc, thick disc, pseudobulge, bulge, and halo. Top panels correspond to a galaxy in the TNG50-like model. The bottom panels correspond to the same galaxy in the strong wind model.
In the text

Fig. A.2. Face-on stellar surface density profiles in terms of R/R*, 50 where R*, 50 is the stellar half-mass radius, for a galaxy with a bar in the TNG50-like model (top) and unbarred galaxy in the strong wind model (bottom) at z = 0, 1. The profiles (black lines) were obtained by fitting simultaneously a Sersic profile for a bulge (red lines) and a Sersic for a bar (pink lines), and an exponential profile (blue lines). Purple lines represent the surface density profile from the predicted Sersic profile and exponential profile. Horizontal dashed lines correspond to the radius at which the fit is done. Dotted horizontal lines correspond to the bar length (pink), the effective radius of the bulge (red), and the scale length of the disc(blue).

In the text