# Simulating nearby disc galaxies on the main star formation sequence. I. Bar formation and the building of the central gas reservoir

Pierrick Verwilghen, Eric Emsellem, Florent Renaud, Milena Valentini, Jiayi Sun, Sarah Jeffreson, Ralf S. Klessen, Mattia C. Sormani, Ashley T. Barnes, Klaus Dolag, et al.

► To cite this version:

Pierrick Verwilghen, Eric Emsellem, Florent Renaud, Milena Valentini, Jiayi Sun, et al. Simulating nearby disc galaxies on the main star formation sequence. I. Bar formation and the building of the central gas reservoir. Astronomy & Astrophysics - A&A, 2024, 687, 10.1051/0004-6361/202348772. insu-04644613

HAL Id: insu-04644613

https://insu.hal.science/insu-04644613v1

Submitted on 11 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Simulating nearby disc galaxies on the main star formation sequence

## I. Bar formation and the building of the central gas reservoir

Pierrick Verwilghen<sup>1,2</sup>, Eric Emsellem<sup>2,3</sup>, Florent Renaud<sup>4,5</sup>, Milena Valentini<sup>6,7,8</sup>, Jiayi Sun<sup>9</sup>, Sarah Jeffreson<sup>10</sup>, Ralf S. Klessen<sup>11,12</sup>, Mattia C. Sormani<sup>13,14</sup>, Ashley T. Barnes<sup>2</sup>, Klaus Dolag<sup>15,16</sup>, Kathryn Grasha<sup>17,18</sup>, Fu-Heng Liang<sup>2,19</sup>, Sharon Meidt<sup>20</sup>, Justus Neumann<sup>21</sup>, Miguel Querejeta<sup>22</sup>, Eva Schinnerer<sup>21</sup>, and Thomas G. Williams<sup>19</sup>

Received 28 November 2023 / Accepted 10 April 2024

## ABSTRACT

Past studies have long emphasised the key role played by galactic stellar bars in the context of disc secular evolution, via the re-distribution of gas and stars, the triggering of star formation, and the formation of prominent structures such as rings and central mass concentrations. However, the exact physical processes acting on those structures, as well as the timescales associated with the building and consumption of central gas reservoirs are still not well understood. We are building a suite of hydro-dynamical RAMSES simulations of isolated, low-redshift galaxies that mimic the properties of the PHANGS sample. The initial conditions of the models reproduce the observed stellar mass, disc scale length, or gas fraction, and this paper presents a first subset of these models. Most of our simulated galaxies develop a prominent bar structure, which itself triggers central gas fuelling and the building of an over-density with a typical scale of 100–1000 pc. We confirm that if the host galaxy features an ellipsoidal component, the formation of the bar and gas fuelling are delayed. We show that most of our simulations follow a common time evolution, when accounting for mass scaling and the bar formation time. In our simulations, the stellar mass of $10^{10} M_\odot$ seems to mark a change in the phases describing the time evolution of the bar and its impact on the interstellar medium. In massive discs ($M_* \ge 10^{10} M_\odot$), we observe the formation of a central gas reservoir with star formation mostly occurring within a restricted starburst region, leading to a gas depletion phase. Lower-mass systems ($M_* < 10^{10} M_\odot$) do not exhibit such a depletion phase, and show a more homogeneous spread of star-forming regions along the bar structure, and do not appear to host inner bar-driven discs or rings. Our results seem to be supported by observations, and we briefly discuss how this new suite of simulations can help our understanding of the secular evolution of main sequence disc galaxies.

**Key words.** hydrodynamics – galaxies: evolution – galaxies: kinematics and dynamics – galaxies: spiral – galaxies: star formation – galaxies: structure

<footnote n="1">¹ Excellence Cluster ORIGINS, Boltzmannstraße 2, 85748 Garching, Germany</footnote>

<footnote n="2">² European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany
e-mail: pierrick.verwilghen@eso.org</footnote>

<footnote n="3">³ Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="4">⁴ Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="5">⁵ University of Strasbourg Institute for Advanced Study, 5 Allée du Général Rouvillois, 67083 Strasbourg, France</footnote>

<footnote n="6">⁶ Astronomy Unit, Department of Physics, University of Trieste, Via Tiepolo 11, 34131 Trieste, Italy</footnote>

<footnote n="7">⁷ INAF – Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34131 Trieste, Italy</footnote>

<footnote n="8">⁸ ICSC – Italian Research Center on High Performance Computing, Big Data and Quantum Computing, Italy</footnote>

<footnote n="9">⁹ Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA</footnote>

<footnote n="10">¹⁰ Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA, USA</footnote>

<footnote n="11">¹¹ Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany</footnote>

<footnote n="12">¹² Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany</footnote>

<footnote n="13">¹³ Università dell'Insubria, Via Valleggio 11, 22100 Como, Italy</footnote>

<footnote n="14">¹⁴ Department of Physics, University of Surrey, Guildford GU2 7XH, UK</footnote>

<footnote n="15">¹⁵ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote n="16">¹⁶ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany</footnote>

<footnote n="17">¹⁷ Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia</footnote>

<footnote n="18">¹⁸ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

<footnote n="19">¹⁹ Sub-department of Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="20">²⁰ Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, 9000 Gent, Belgium</footnote>

<footnote n="21">²¹ Max Planck Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="22">²² Observatorio Astronómico Nacional (IGN), C/Alfonso XII 3, Madrid 28014, Spain</footnote>
# 1. Introduction

Bars and spiral arms both form naturally in galaxy disc-like systems (see Eskridge et al. 2000; Buta et al. 2015; Sellwood & Masters 2022, and references therein). They have a significant impact on the internal disc structures, induce specific torques and resonant regions, and impose, via the tumbling and varying gravitational potential, constraints on the galactic orbital skeleton. Bars have been specifically and extensively studied and modelled (Combes & Sanders 1981; Kaufmann & Contopoulos 1996; Athanassoula 1992a; Kormendy & Kennicutt 2004; Barazza et al. 2008; Gadotti 2008; Kraljic et al. 2012; Goz et al. 2015; Fragkoudi et al. 2017) and have a direct influence on the short- and long-term evolution of the system. We know that bars funnel the gas to the central regions (~100 pc) of galaxies, and have been suggested to play a role in the large-scale fuelling of the central supermassive black hole (SMBH; Fukuda et al. 1998; Shlosman et al. 1989; Ho et al. 1997).

Similarly to our own Milky Way, most nearby disc galaxies are also believed to host central SMBHs, with masses between $10^6$ and $10^9 M_\odot$, at their centres (Magorrian et al. 1998; Reines & Volonteri 2016). When gas builds into a surrounding accretion disc, those SMBHs can drive nuclear activity (Krolik 1999; Padovani et al. 2017), releasing a large amount of energy in their environment in the form of mainly radiative (quasar mode) or mechanical (radio mode) feedback (Sanders et al. 1988; Fabian 2012). Such active galactic nuclei (AGNs) may impact the dynamical evolution of their host galaxy, their chemical composition, and even modulate the central gas accretion itself (see e.g. Combes 2017, and references therein).

The exact physical mechanisms and timelines associated with the fuelling of central SMBHs are still heavily debated (see e.g. Hopkins & Quataert 2010; Cheung et al. 2015), both due to the range of spatial and timescales involved and to the complexity and non-linear nature of the physical processes involved. The duty cycle of an accretion disc more directly depends on small-scale physics and the close environment of the black hole ($10^{-6}$ to $10^{-2}$ pc), while the larger-scale kiloparsec dynamics and structures may influence its long-term evolution. Two main avenues are generally called upon for the global transport of gas in disc galaxies. An external channel or origin is when low angular momentum gas is falling from, for example, the circum-galactic medium (CGM) towards the galaxy centre (Lacey & Fall 1985; Bilitewski & Schönrich 2012). An internal channel is connected to the secular evolution and processes, such as gravitational torques and stellar-driven feedback, acting as a way to extract angular momentum from the gas (Porciani et al. 2002; Stewart et al. 2013; Übler et al. 2014; Genel et al. 2015; Pezzulli & Fraternali 2016; Schmidt et al. 2016; Valentini et al. 2017). The required amount of gas to trigger the AGN is relatively small compared with the total mass of gas available in the galaxy (Hickox et al. 2012; Padovani 2017). This means that only a very small fraction of the gas funnelled within the central few hundreds of parsecs is needed to account for the nuclear black-hole-related activity. Further considering the several orders of magnitude in spatial scale that are spanned from the galactic bar (kiloparsecs) to an accreting black hole (~$10^{-4}$–0.1 pc, see e.g. Guo et al. 2023; Izumi et al. 2023, and references therein), it is not yet clear how the bar relates to the processes that are most relevant to drive the duty cycle of AGNs.

Progress may come from the realisation that many disc galaxies host a central over-density of gas at a scale of a few hundreds of parsec from the centre (Comerón et al. 2010). While

such structures may not directly relate with the ongoing nuclear activity itself, they can play the role of a “gas reservoir” and represent an intermediate milestone that we need to understand if we wish to draw a full picture of central gas accretion and its time evolution. Those gas reservoirs are often interpreted in the context of bars as coinciding with a change in the sign of the torques (Sanders 1977; Wada 1994), that is, the location of the so-called inner Lindblad resonance (ILR). However, these models fail to reproduce the results of hydro-dynamical simulations (Sormani et al. 2015), thus calling for a new framework that can explain the full picture (see Sormani et al. 2024). In the Milky Way, a ~200 pc central molecular zone (CMZ) with a gas mass of ~$1–7 \times 10^7 M_\odot$ has been identified and extensively studied (Ferrière et al. 2007; Longmore et al. 2013; Henshaw et al. 2023). In other galaxies, such gas reservoirs (or CMZs) have also been observed and again, often interpreted in the context of the gas fuelling by bars, while details of its building and life cycle (consumption via e.g. star formation) are not well constrained.

Some pioneering work was carried out by Athanassoula (1992b), who studied the response of gas to an idealised two-dimensional barred potential and observed the building of a central gas concentration within the central 1 kpc region. This study raised questions about the building, consumption, and overall evolution of central gas reservoirs, as well as their related physical phenomena, such as gravitational torques, shear, feedback, and the exact role of bars and spiral arms in fuelling the central regions. Those questions further triggered observational surveys and numerical works to study the fuelling rate and the size of those gas reservoirs (Combes et al. 2004; García-Burillo et al. 2005; Boone et al. 2007; Comerón et al. 2010; Sormani & Barnes 2019; Sormani et al. 2023). We still need a systematic approach to investigate numerically and characterise the different scenarios of the building and evolution of such gas reservoirs in a three-dimensional live potential (including stars, gas, and dark matter) for a set of models tailored to nearby disc galaxies.

In this paper, we focus on the design, running, and exploitation of a suite of numerical hydro-dynamical simulations of isolated disc galaxies to study the formation and evolution of bars, and the associated building of the central gas reservoirs. We motivated such a grid of simulations using the PHANGS<fnref n="1" />-ALMA galaxy sample as a guiding baseline, which provides high-resolution, high-sensitivity, and short-spacing CO spectral cubes of nearby disc galaxies and they are ideal for comparison with the models.

In Sect. 2 we introduce the grid of models. In Sect. 3 we provide details regarding the code and numerical recipes we used to perform the simulations, and in Sect. 4 we present the first results extracted from that grid of simulations. We conclude in Sect. 5 with a brief summary.

# 2. Methods

## 2.1. The PHANGS-ALMA sample

As a basis to build our simulation sample, we used the PHANGS-ALMA survey (Leroy et al. 2021) which consists of 118 low inclination nearby main sequence star-forming disc galaxies (in a distance range of ~3–30 Mpc) with varying properties, including stellar mass, gas fraction, surface density and rotation curve (Lang et al. 2020). Those input parameters can be used to extract relevant control parameters to design a grid of

<footnote n="1">¹ Physics at High Angular resolution in Nearby GalaxieS: https://sites.google.com/view/phangs/home</footnote>
initial conditions for our simulations. The most directly relevant galactic properties to build such a grid of simulations relate to the gravitational potential and associated baryonic mass distribution, hence including the stellar mass, the gas fraction and the presence of a stellar bar or not (Querejeta et al. 2021). The top panel of Fig. 1 shows the distribution of the stellar masses in the PHANGS sample, emphasising the distribution of barred and non-barred galaxies. The stellar masses range from ~$10^{9.25}$ up to ~$10^{11}$ $M_{\odot}$ and the vast majority (~75%) of that sample shows the presence of a central bar (Stuber et al. 2023).

## 2.2. The control parameters

The main objective of our grid of hydro-dynamical models is to study the impact of galaxy properties on the building and evolution of the central (few 100 pc) gas reservoirs. We thus focus primarily on the properties of galaxies that we expect would more directly impact the structure of bars and spiral arms (i.e. size, strength and pattern speed). Based on the observed properties of the PHANGS sample, we focus on five galactic parameters for our models, namely:

1. the stellar mass $M_{*}$;

2. the gas fraction $\alpha = M_{g}/(M_{g} + M_{*})$ where $M_{g} = M_{HI} + M_{H_{2}}$;

3. the typical scale length of the stellar distribution $l_{*}$;

4. the typical scale length of the gas (HI and $H_{2}$) distribution $l_{g}$;

5. and the central bulge mass fraction $\beta = M_{b}/M_{*}$.

In this work, only four of these parameters are used as control parameters since we set the value of the gas scale length to $2 l_{*}$. We decompose the stellar mass into a given disc and central ellipsoid as $M_{*} = M_{d} + M_{b}$, with $M_{d}$ and $M_{b}$ the mass of the stellar disc and the stellar ellipsoid, respectively. In the following, we will label “bulge” this central ellipsoid, emphasising the fact that such a label is often misused to describe a central excess of light above a given larger-scale disc (Gadotti et al. 2020), leading to confusion. In our semantic usage, bulge refers to a structure that bulges out of the disc and must then be both extended (with scales significantly larger than e.g. a nucleus, or the scale height of the disc) and puffed (e.g. axis ratio significantly above 0.2). Simple relations connect the gas and stellar disc masses with the $\alpha$ and $\beta$ parameters, namely:

$$M_{g} = \frac{\alpha}{1 - \alpha} \cdot M_{*}, \quad M_{d} = (1 - \beta) \cdot M_{*}. \qquad (1)$$

We select four values for the stellar mass, namely $\log_{10}M_{*}(M_{\odot}) = 9.5, 10, 10.5, 11$, to cover the PHANGS range (see Fig. 1).

We further examine the trends of all the other four input parameters as a function of the stellar mass. The right and bottom panels of Fig. 1 show the distribution of the gas fraction among the PHANGS galaxies and its values as a function of the stellar mass, respectively. The gas fraction tends to increase with decreasing stellar masses below $10^{10}$ $M_{\odot}$. We assume exponential radial profiles both for the initial stellar and gas surface density ($\Sigma_{d}, \Sigma_{g}$ accordingly) profiles and fitted all azimuthally averaged observed PHANGS profiles accordingly.

$$\Sigma_{d,*/g}(r) = \Sigma_{0,*/g} \exp\left(-\frac{r}{l_{*}/g}\right), \qquad (2)$$

where $l_{*}$ and $l_{g}$ are respectively the stellar and gas scale length, and $\Sigma_{0,g}$ is the initial bulge surface density. During the fitting of the stellar density profile, we allow the simultaneous presence of a central component with a surface brightness profile described

as an additional Sersic profile (Sérsic 1963):

$$\Sigma_{b}(r) = \Sigma_{0,b} \exp\left(-b_{n}\left[\left(\frac{r}{R_{e,b}}\right)^{1/n_{b}} - 1\right]\right), \qquad (3)$$

where $r$ is the bulge coordinate radius, $\Sigma_{0,b}$ is the initial surface density, $R_{e,b}$ is the scale length, $n_{b}$ is the Sersic index, and $b_{n}$ is a coefficient depending on $n_{b}$ calibrated to reach the half of the total luminosity at $R_{e,b}$. Since we wish that the three-dimensional ellipsoidal representation based on the one-dimensional Sersic function follows averaged properties of galactic bulges (Sérsic 1963), we forced an intrinsic ellipticity of 0.4 (axis ratio of 0.6, see Moriondo et al. 1998; Sotnikova et al. 2012), as opposed to the standard assumptions of a spheroid (i.e. a spherical bulge). As mentioned, in the following, that component will be referred to as the “central bulge”.

Figure 2 shows an example of the stellar (top panel) and gas (bottom panel) surface density fits for one specific PHANGS galaxy (NGC 5134). Fits to the observed surface brightness profiles were not always satisfactory when the data could not be well represented by a single exponential or a combined exponential plus Sersic components. We thus visually classified our fits into three categories: [1] good fits, [0] reasonable fits (exhibiting radial ranges with significant residuals), and [–1] failed fits when no good representation of the data could be found using the above-mentioned functions (see Fig. A.1 for more details). This classification is subjective and only meant as a guide for the values of reference parameters when building up the grid of initial conditions: it does not have a significant impact on the chosen grid of models we have selected. This analysis led us to fix a set of values for $l_{*}, l_{g}, R_{e,b}$ and $n_{b}$ for each galaxy in the PHANGS-ALMA sample: those are shown in Fig. 3 where they are plotted as a function of the stellar mass. We also show the values of the gas fraction $\alpha$ and those of the bulge mass fraction as derived from the obtained scale lengths.

## 2.3. The grid of models

Our grid of models relies on a selected set of values for each of the four control parameters (see previous section). As a first approximation, this set provides a fair representation of the global trends (e.g. with respect to the stellar mass) observed in the PHANGS sample. We sometimes constrained the value of a given parameter to stay constant at all masses, at the expense of missing the more relevant range of observed parameters<fnref n="2" />. This is true, for instance, for the gas fractions ($\alpha$) at the lowest stellar mass bin, where we chose to only probe 0 and 10%, thus departing from the significantly larger observed values. It is also the case for bulge mass fraction ($\beta$) for which we also kept low values (0 and 10%) while many targets exhibit values up to 40% or higher. This latter choice was motivated by the fact that we expect the central region to grow in mass during the secular evolution of the system and that most of the measured bulge indices are on the low side, hence most possibly reflecting the presence of a central (flattened) disc (not a puffed-up spheroid). In practice, those values have been used to build initial conditions for the hydro-dynamical simulations: they are illustrated by purple stars in Fig. 3 as mentioned before and are all tabulated in Table 1.

In the following, all models share a common labelling scheme, namely: each model adopts a format as

<footnote n="2">² Note that an extended grid of models, better covering those missed ranges, is planned (see Sect. 5).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Gas fraction (a) as a function of the stellar mass of the galaxies in the PHANGS sample and the corresponding histograms for the stellar mass (top panel) and gas fraction (right panel) distribution (see Sect. 2.2).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustrations of the fitted radial gas and stellar density profiles. The plots show the original data points (in black, see Sun et al. 2022, mega-tables version 3), i.e. the azimuthal-averaged surface density profiles for the stellar (top panel) and gas (bottom panel) components of the galaxy NGC 5134. The corresponding fits are superimposed and colour-coded accordingly (blue for the disc component, red for the total disc and ellipsoid (or bulge) component).</figcaption></figure>

GxxxMxxxFxxLxBxx, where “G” (stands for Galaxy) is followed by an integer (used as an internal reference), M by the log of the stellar mass multiplied by 10, F by the gas fraction, L by the typical scale length of stars, and B by the bulge mass fraction. This represents a set of 16 high-priority models, spanning typical PHANGS parameters, and tractable in terms of computing time. This paper focuses on this preliminary first subset of 16 initial conditions, with suitable values of control parameters for the comparison between all the stellar mass bins (see Table 1) that already covers a good range of properties for this sample. An extended grid of 54 models is planned but requires a significantly larger investment in computational time, and its analysis will be presented in a future paper.

## 3. Numerical simulations

### 3.1. Initial conditions

Initial conditions for the RAMSES adaptive mesh refinement (AMR) code (Teyssier 2002) require information both for the particles (stars, dark matter, a central SMBH), and gas content. The grid of values gathered in Table 1 is used for the description of the distribution of each individual component. We use the Python library “pymge” based on the Multi-Gaussian Expansion (MGE; Emsellem et al. 1994a,b) method to fit the associated two-dimensional distributions (e.g. exponential discs) and deproject them into three-dimensional density distributions. The vertical axis ratios of the ellipsoid resulting from that 3D deprojection are set to $q = 0.6, 0.1, 0.05$, respectively for the bulge, stellar disc and gas disc. The dark matter is constrained by the averaged rotation curves observed via the PHANGS sample and generated using an Einasto spherical mass distribution (Einasto 1965; Ludlow & Angulo 2017) as given by:

$$ \rho_h(r) = \rho_{h,0} \exp \left( -2m \left[ \left( \frac{r}{l_h} \right)^{1/m} - 1 \right] \right), \quad (4) $$

where $\rho_h$ is the 3D density profile of the dark matter halo, $r$ is the spherical radius, $l_h$ is the scale length and $m$ is the halo index. The values we chose for the parameters of Eq. (4) are summarised in Table 2. Once those three-dimensional profiles are
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Results of the fit for different parameters of the models as a function of the stellar mass. Top panels: gas fraction (left) and stellar scale length (right). Middle panels: scale length ratio between the gas and stars (left) and scale length of the stellar bulge (right). Bottom panels: bulge index (left) and bulge mass fraction (right). Colour circles show the actual fits (see the text in Sect. 2.2 for details) to the PHANGS-ALMA sample, while the selected values for our used initial conditions (IC(model)) are shown with purple stars.</figcaption></figure>

fixed for all the components, we fix the number of particles and
derive their initial positions and velocities by solving the Jeans
equations. We assume a local anisotropy constrained by the flat-
tening of the individual components, with $\delta = 1 - \sigma_z^2/\sigma_R^2 =$
$0.6 \times (1 - q)$ where $\sigma_R, z$ the radial and vertical stellar velocity
dispersion in cylindrical coordinates (Binney et al. 2009). The
more flattened the component is, the larger its (initial) veloc-
ity anisotropy. We use a constant mass resolution of $10^4 M_\odot$ for

all stars present in the initial conditions (“old” stellar particles).
This leads to an increasing number of stellar particles, that is,
about $3 \times 10^5$, $10^6$, $3 \times 10^6$ and $10^7$, for the corresponding four
stellar masses 9.5, 10, 10.5 and 11 (in $\log_{10}[M_\odot]$), respectively.
The number of particles for the dark matter component was fixed
to $10^6$ for all models.

Figure 4 shows a comparison between the PHANGS rotation
curves (shaded areas, encompassing individual velocity profiles
<table><caption>**Table 1.** Parameters for the initial conditions of the models of this paper.</caption><thead><tr><th>Model</th><th>log<sub>10</sub>(M<sub>★</sub>)<br>[M<sub>☉</sub>]</th><th>α<br>[%]</th><th>l<sub>★</sub><br>[kpc]</th><th>l<sub>g</sub>/l<sub>★</sub></th><th>β<br>[%]</th></tr></thead><tbody><tr><td>G001M095F10L2B00</td><td>9.5</td><td>10</td><td>2</td><td>2</td><td>0</td></tr><tr><td>G002M095F10L2B10</td><td>9.5</td><td>10</td><td>2</td><td>2</td><td>10</td></tr><tr><td>G013M095F20L2B00</td><td>9.5</td><td>20</td><td>2</td><td>2</td><td>0</td></tr><tr><td>G014M095F20L2B10</td><td>9.5</td><td>20</td><td>2</td><td>2</td><td>10</td></tr><tr><td>G037M100F10L2B00</td><td>10</td><td>10</td><td>2</td><td>2</td><td>0</td></tr><tr><td>G038M100F10L2B10</td><td>10</td><td>10</td><td>2</td><td>2</td><td>10</td></tr><tr><td>G053M100F20L2B00</td><td>10</td><td>20</td><td>2</td><td>2</td><td>0</td></tr><tr><td>G054M100F20L2B10</td><td>10</td><td>20</td><td>2</td><td>2</td><td>10</td></tr><tr><td>G105M105F10L3B00</td><td>10.5</td><td>10</td><td>3</td><td>2</td><td>0</td></tr><tr><td>G106M105F10L3B10</td><td>10.5</td><td>10</td><td>3</td><td>2</td><td>10</td></tr><tr><td>G137M105F20L3B00</td><td>10.5</td><td>20</td><td>3</td><td>2</td><td>0</td></tr><tr><td>G138M105F20L3B10</td><td>10.5</td><td>20</td><td>3</td><td>2</td><td>10</td></tr><tr><td>G161M110F10L5B00</td><td>11</td><td>10</td><td>5</td><td>2</td><td>0</td></tr><tr><td>G162M110F10L5B10</td><td>11</td><td>10</td><td>5</td><td>2</td><td>10</td></tr><tr><td>G177M110F20L5B00</td><td>11</td><td>20</td><td>5</td><td>2</td><td>0</td></tr><tr><td>G178M110F20L5B10</td><td>11</td><td>20</td><td>5</td><td>2</td><td>10</td></tr></tbody></table>

**Notes.** From left to right columns: label, total stellar mass, mass gas fraction, stellar scale length, scale length ratio (gas over stars) and bulge mass fraction.

<table><caption>**Table 2.** Values for the parameters describing the dark matter profiles (i.e. $\rho_{h,0}$, the halo density; $l_{h,0}$, the typical halo scale length; and $m$, the halo index) as a function of the stellar mass $M_{\star}$.</caption><thead><tr><th>M<sub>★</sub><br>log<sub>10</sub>([M<sub>☉</sub>])</th><th>ρ<sub>h,0</sub><br>[M<sub>☉</sub> pc<sup>-3</sup>]</th><th>l<sub>h,0</sub><br>[kpc]</th><th>m</th><th>M<sub>BH</sub><br>log<sub>10</sub>([M<sub>☉</sub>])</th></tr></thead><tbody><tr><td>9.5</td><td>3 × 10<sup>-4</sup></td><td>20</td><td>3</td><td>5.5</td></tr><tr><td>10</td><td>6 × 10<sup>-4</sup></td><td>25</td><td>3.5</td><td>6</td></tr><tr><td>10.5</td><td>7 × 10<sup>-4</sup></td><td>30</td><td>1.5</td><td>6.5</td></tr><tr><td>11</td><td>8.5 × 10<sup>-4</sup></td><td>35</td><td>2</td><td>7</td></tr></tbody></table>

**Notes.** Those values have been inspired by the work of Kun et al. (2017), but their final values have been chosen to match the PHANGS rotation curves. The central black hole masses $M_{BH}$ we use in the initial conditions are also shown in the last column.

for a given stellar mass bin, Lang et al. (2020) and our 16 models (solid curves) for the different mass ranges. Our models are, by design, in good agreement with the PHANGS galaxies, reproducing the global trend of real galaxies. We note that the limited extent in the observed rotation curves for the lower mass bins is expected considering the smaller size and shorter radial extent of the (mostly CO and HI) tracers (Leroy et al. 2019). We also emphasise that those circular velocity curves represent the initial state of our hydro-dynamical simulations, and are thus meant to evolve with time (i.e. see the black curves). It is worth mentioning that the PHANGS rotation curves are the observed gas velocity profiles. Converting such velocity profiles to actual circular velocities, as measured in our simulations, would require detailed modelling, including, for example, the effect of asymmetric drift and non-circular motions. Considering the relatively low velocity dispersion of the molecular gas, we do not expect the impact of asymmetric drift to be significant. Non-circular motions may, however, dominate in the central region. This motivated us to use those PHANGS observed rotation curves only as guidelines to build our mass models.

## 3.2. Refinement strategy

All simulations in the three lower stellar mass bins (resp. higher mass bin) are run within a 100 kpc (resp. 200 kpc) cubic box. A refinement level $l$ thus corresponds to a $100/2^l$ (resp. $100/2^l$) kpc cubic cell. We imposed a minimal refinement level of 7 (resp 8; ~780 pc) and a maximal one of 13 (resp. 14; ~12 pc). The refinement strategy we adopt for the gas is based on a Jeans polytropic approximation (see Renaud et al. 2013) allowing us to set a realistic threshold for the amount of gas contained inside a cell. The maximum mass of gas contained inside a cell at the second last level ($l$ = 12 or 13, 24 pc) before the triggering of the refinement to the last level ($l$ = 13 or 14, 12 pc) is around 14 000 $M_{\odot}$. An additional refinement criterion is constrained by the minimum number of particles per cell (stars and dark matter) that is set to 8.

## 3.3. Numerical recipes

A number of standard physical processes were implemented including gas cooling, star formation, stellar feedback and the evolution of metals such as iron and oxygen (see Agertz et al. 2013 for more detailed information). The physics of the heating follows a uniform UV background model proposed by Haardt & Madau (1996) and calibrated for redshift $z$ = 0. The physics of the cooling is implemented and divided into two different regimes. The first regime goes down two $10^4$ K and uses collisional de-excitation and atomic recombinations. The second regime, below $10^4$ K, follows the table of Sutherland & Dopita (1993). The formation of new stars is triggered above a gas density threshold via a constant efficiency, as given by the following equation:

$$ \dot{\rho}_{*} = \frac{\rho_{g}}{t_{\text{SF}}} \quad \text{for } \rho_{g} \ge m_{h} n_{*}, \qquad (5) $$

where $\dot{\rho}_{*}$ is the star-formation rate (SFR), $\rho_{g}$ is the gas density, $m_{h}$ is the mass of the hydrogen atom and $t_{\text{SF}}$ is the typical gas depletion time (~2 Gyr). The latter is given by:

$$ t_{\text{SF}} = t_{\text{ff}} / \epsilon_{\text{ff}}, \qquad (6) $$

with $t_{\text{ff}}$ the local free fall time and $\epsilon_{\text{ff}}$ the star formation efficiency per free-fall time. We set the efficiency to 2%, and $n_{*}$ to $100 \text{ cm}^{-3}$, a posteriori checking that our systems follow the trend observed in the PHANGS galaxies (see Fig. 5). Stellar feedback, which releases energy, momentum and metals in the surrounding ISM, is implemented by taking into account the contribution from supernovae (SN) SNIa, SNII and stellar winds. The energy released by the supernova depends on the local cooling radius (see Kim & Ostriker 2015). Finally, we also adopt solar metallicity as initial condition ($Z = 1 Z_{\odot}$), which does not have a significant impact on the gas processes (star formation, feedback).

# 4. Results

In this section, we first describe some general properties of the evolved systems and then illustrate the global time evolution and bar formation. We briefly discuss the bar formation times in the light of past results, and consider the time evolution in terms of specific successive phases. We finally hint at a significant difference in the phases of the building of a central gas concentration impacting the way star formation proceeds in the lowest mass bin of our sample.
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Comparison between the observed gas velocity curves of the PHANGS sample (shaded areas) and the circular velocity profiles from our models (solid curves). The vertical lines near $R = 0$ illustrate the contribution of the central SMBH. The dotted curves represent the circular velocities extracted from the first snapshot (INI Sim) of each simulation and are consistent with the analytic derivation associated with the initial conditions. The dashed curves represent the circular velocities at $\tau = 2$ (see Sect. 4.3, which corresponds to twice the bar formation time of our simulated galaxies). As mentioned in the text, the PHANGS rotation curves are meant as guidelines to build the initial mass models.

## 4.1. General properties and time evolution

### 4.1.1. The star-forming main sequence

We want our simulated galaxies to be representative of the star-forming main sequence of nearby disc galaxies, hence located on the Kennicutt-Schmidt relation (Schmidt 1959; Kennicutt 2007). Our choices for the initial properties (e.g. distribution of baryons, gas fractions, etc.) and sub-grid recipes (e.g. star formation efficiency, see Sect. 3.2) do not guarantee such a result, because star formation may partly regulate itself in such numerical experiments (Ostriker & Kim 2022). We thus computed the global SFR for all 16 simulations and compared its characteristics and trends with the instantaneous SFR of the PHANGS-ALMA sample as shown in Fig. 5. We computed the SFR by taking the time average of the SFR, only including times between $\tau = 1$ and $\tau = 2$ (see Sect. 4.3, which corresponds to the starburst phase of our simulated galaxies), thus focusing on the active star formation regime.

In Fig. 5, we confirm that our simulations are almost all following the PHANGS main sequence, except for galaxies with the lowest stellar mass bin (i.e. 9.5 in $log_{10}(M_\odot)$): those have SFR lower than the observed values in this sample. This is expected since we intentionally chose to run this first subset of simulations with only two values of the gas mass fraction, those being significant underestimating the observed gas fractions (~30-40%) for that stellar mass bin, as shown in Fig. 3. Higher gas fractions (thicker squares) naturally lead to a higher SFR, and models with a bulge (hatched squares) tend to have lower resulting SFR as expected if the self-gravity of the disc is lowered. The impact of the bulge reaches its maximum for the bulged model with the highest stellar mass, and a gas fraction of 10%. We describe in more detail the characteristics of this model in the next sections.

### 4.1.2. Building of the bar

Figure 6 illustrates the global evolution of one galaxy over time by showing the maps of gas, new and old stars. We start with an axisymmetric distribution of gas and stars, which rapidly develops low-contrast spiral arm structures. Gas follows up by cooling

<figure><img src="image_6.png" /></figure>

**Fig. 5.** Comparison of the SFR stemming from the PHANGS data sample (grey crosses) and the SFR computed from the simulated galaxies (coloured squared). The colour of the squares corresponds to the different stellar masses. The red dashed line represents the star formation main sequence at $z = 0$ from Leroy et al. (2019).

and forming new stars (top left panel). After a few hundred Myr we start to observe the formation of a bar and spiral arms in the old and new stars (top right panel). Once the bar has formed, the gas concentration in the central region starts to increase, leading to the build-up of an early gas reservoir (bottom left panel). This could be induced by the emergence of an inner Lindblad resonance (ILR: Lin et al. 2008; Sormani et al. 2024). Finally, we observe the building and growth of a central gas reservoir (bottom right panel).

This scenario is similar for most of the 16 simulated galaxies and the evolution over time of the 16 simulations is shown in Appendix B where the different times are chosen to emphasise the evolution of the building of the gas reservoir within the central 1 kpc region (see Sect. 4.1.3).

One of the main features developed by almost all our simulations is a stellar bar (except for models G162M110F10L5B10 and G178M110F20L5B10). We present
G053M100F20L2B00

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution over time (in millions of years) of the gas and stars of the model labelled G053M100F20L2B00 (one of the sixteen simulations presented in this paper). The two big panels illustrate four main stages, chronologically ordered from top left (0–400 Myr), to top right (500–900 Myr), bottom left (1000–1800 Myr) and bottom right (2000–2800 Myr). Top left: first spiral structures, cooling of the gas, onset of star formation, initial bar structure emerging. Top right: bar strengthening and active local star formation. Bottom left: building of a central concentration of gas (and new stars). Bottom right: growth of the gas reservoir. In each big panel, from top to bottom, we present maps for the gas mass density (GAS), the mass of young stars (≤50 Myr, NS), of old stars (OS) and finally the cumulative mass of new stars (formed since the beginning of the simulation, CNS). Each panel shows a 10 × 10 kpc² region and the colour scaling is adapted for each panel.</figcaption></figure>

model G053M100F20L2B00 in Fig. 6 to illustrate the growth of the bar witnessed in our simulations (see third row). To more quantitatively characterise the evolution of the bar over time, we estimated the amplitude of the bar via a polar (R, θ) Fourier decomposition of the face-on surface stellar mass density, and specifically used the traditional A₂ Fourier coefficient (Efstathiou et al. 1982; Athanassoula & Misiriotis 2002; Athanassoula et al. 2013) given by

$$A_2 = \frac{\sqrt{a_2^2 + b_2^2}}{a_0}, a_n = \sum_{i=1}^{N} m_i \cos(n\theta_i), b_n = \sum_{i=1}^{N} m_i \sin(n\theta_i), \quad (7)$$

where $m_i$ is the mass of the ith star and $\theta_i$ its corresponding position angle. The full radial $A_2(R)$ profile is then used, for example, to track its maximum over time, giving us a quantitative assessment of the evolution of the bar strength, as shown in Fig. 7 for the model labelled G053M100F20L2B00. In that figure, we see that the bar grows rapidly during the first 500–600 Myr, roughly following an exponential behaviour (Binney 2020; Bland-Hawthorn et al. 2023). Then, it reaches a first maximum after ~750 Myr and a second after ~1000 Myr. This characteristic shape of $A_2$ was already reproduced by Athanassoula (2002). At later times, $A_2$ exhibits some oscillations until the end of the simulation. These modulations coincide with the relative cycling
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Evolution of the bar strength (measured through the $A_2$ Fourier coefficient) for the model G053.</figcaption></figure>

(and alignment) between the bar and the outer spiral arms which have different pattern speeds (Hilmi et al. 2020).

We emphasise that the resulting $A_2$ values should depend on the tracer to compute it, for instance, using particle masses, accounting only for old or newly formed stars, or using luminosity-weighted quantities as is naturally the case for observations. It may also be influenced by the presence of dust and by additional processing stages (e.g. deprojection of observed photometry). Still, the $A_2$ values found in our simulations are in the range [0.2–0.5], consistent with the ones observed by Díaz-García et al. (2016) for star-forming late-type galaxies, and also consistent with the majority of barred galaxies in the PHANGS sample. Stuber et al. (2023) report significantly larger values ($A_2$ > 0.4), but those are representative of the high tail of the bar strength distribution. The fact that we do not reproduce this tail of high $A_2$ values may naturally arise from the early bar phase that we are probing with our simulations (i.e. up to about 3 Gyr). Some of the strong bars observed in the PHANGS sample could be at a later evolutionary stage for which the bar has experienced a second (secular) growth (see Athanassoula 2013).

### 4.1.3. The central 1 kpc

In the PHANGS-ALMA galaxy sample, we observe inner molecular rings with typical sizes ranging from ~100 pc to ~1 kpc (Leroy et al. 2021). Our simulations suggest their size and mass evolve with time. We thus monitored the evolution of the central gas mass concentration and star formation within the central 1 kpc, and their relation with the bar using our 16 simulations. In the rest of this section, we refer to the (central 1 kpc) “gas reservoir” to describe the mass of gas inside the central 1 kpc region.

In Fig. 8, the left (resp. right) column presents time evolution profiles for models with the lower (9.5 and 10 in $log_{10}(M_\odot)$) (resp. higher, 10.5 and 11 in $log_{10}(M_\odot)$) stellar masses. The middle panels show the evolution of the total gas mass contained inside a cylindrical radius of 1 kpc and a thickness of 1 kpc. The initial gas mass within 1 kpc correlates with the initial stellar mass and gas fraction as expected. After a generic plateau lasting about 300 Myr, the gas concentration increases steadily. This increase of the gas mass inside the central 1 kpc region is coincident with an increase of the SFR: those roughly correspond to $t_{bar}$, the time at which the $A_2$ coefficient reaches the value of

~0.2, whatever the stellar mass, gas fraction or the presence of a bulge. Since bulges tend to increase $t_{bar}$ (see Sects. 4.1.1 and 4.2), models with bulges also show a corresponding delay in the start of the increase of the gas mass.

G162M110F10L5B10 seems to be a peculiar case among the two models that do not form a bar: $A_2$ never reaches a value of 0.2. Gas is consequently not funnelled towards the centre as for other models, and there are no star formation bursts with SFR staying around $10^{-3} M_\odot$ yr$^{-1}$. In that simulation, the bulge has a stabilising effect, enough in this specific case to prevent any bar to form for the duration of our simulation. This is a well-known result, and has been already studied and recently emphasised in various papers (see e.g. Sellwood & Evans 2001; Saha & Naab 2013; Kataria & Das 2018; Fujii et al. 2019). In that context, the bar in G162M110F10L5B10 possibly has a formation time $t_{bar}$ significantly larger than a Hubble time (see also Bland-Hawthorn et al. 2023).

### 4.2. Typical bar formation time

Using simulation G053M100F20L2B00 as an example, we witness the onset of bar formation as early as $t = 300$ Myr (Fig. 6). However, the bar seems to more robustly emerge between $t = 400$ and $t = 500$ Myr and this corresponds with the time when $A_2$ reaches a value of ~0.2. This is also true for the other simulations in our sample. In the following, we therefore use $t_{bar}$, the time corresponding to $A_2 = 0.2$, as the reference time for the bar formation (see Sect. 4.1.3). Such a reference measurement is commonly used in the literature (Athanassoula & Misiriotis 2002; Fujii et al. 2018).

All simulations, except two (models G162M110F10L5B10 and G178M110F20L5B10) form bars. The values for $t_{bar}$ are shown for the 14 simulations forming a bar in Fig. 9. To further probe the stability of our models against bar formation, we have applied the criterion established by Efstathiou et al. (1982) for all 16 models. We indeed find that only the two above-mentioned models do overshoot the threshold for stability. While this is an interesting result, suggesting that we should indeed not expect the formation of a bar in models G162M110F10L5B10 and G178M110F20L5B10, we caution the reader regarding the use of that stability criterion as a robust quantitative diagnostic. As pointed out by Athanassoula (2008), this criterion does not take into account the interaction between the disc and the halo, nor the velocity dispersion, which plays a major role in disc stability. They also suggested that this parameter is not designed for the stability analysis of multi-component discs. This has been also shown in a recent work carried out by Romeo et al. (2023), whose conclusion is that this parameter fails at separating barred and non-barred galaxies in 55% of the cases.

We further emphasise the fact that $t_{bar}$ should be considered a relative indicator that exhibits significant systematics. We have conducted a series of tests using the simulation G037M100F10L2B00 as a reference, varying the initial number of particles and the maximum level of refinement, and even shutting off the gas cooling and star formation recipes. As long as we keep the gas disc live within the simulation, we do not witness a significant change in the evolution of the $A_2$ profile and thus in $t_{bar}$. When the number of disc particles is changed (from 0.5 to 1, 2, 5 or 10 million), the time increase (i.e. increasing rate) of $A_2$ looks similar. However, the starting time of that initial increase, tracing the growth of $m = 2$ modes typically varies from one simulation to the next by ±75 Myr. That timing offset does not seem to depend monotonically on the number of particles, as it is, for instance, slightly lower for 2 and 10 mil-
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Evolution of the maximum of the $A_2$ coefficient (top panels), the mass of gas within 1 kpc (middle panels), and the SFR (bottom panels) over time for the 16 simulations (evolved until max 3 Gyr). The colour gradient represents the four stellar masses (9.5, 10, 10.5 and 11 [log($M_\odot$)]) from the lightest to the most massive. The left column represents the less massive galaxies (9.5 and 10, in green and cyan, respectively) and the right column shows the most massive ones (10.5 and 11, in blue and dark blue, respectively). The solid and dashed curves illustrate the models without and with a bulge, respectively and the thickness of the lines accounts for the gas fraction (10 and 20% for the thinnest and thickest, respectively). The typical bar formation time $t_{bar}$ is also shown (red diamonds) within a time interval of 200 Myr (red part of the curves).</figcaption></figure>

lion disc particles, and thus higher for 1 and 5 million disc parti-
cles. We interpret it as partly due to the contribution of the spiral
modes that may interfere with the bar mode, as well as varia-

tions in the initial relaxation phase. While it is beyond the scope
of this study to further probe such a dependence, it does show
that a significant systematic uncertainty should be included when
discussing $t_{\text{bar}}$ measured in such a way. In the rest of the paper, we added a (conservative) systematic uncertainty of ±100 Myr on our $t_{\text{bar}}$ measurements and checked that it did not impact our results.

The most obvious trend associated with $t_{\text{bar}}$ is connected with the impact of the initial ellipsoid (or bulge). When an ellipsoid or bulge (hatched squares in the figure) is added to the model, a delay of a few 100's of Myr in the bar formation is induced, leading to a delay in the building and fuelling of the central 1 kpc (see also Sect. 4.1.3). Note that the induced time delay is significant, and much larger than the above-mentioned systematic uncertainty in $t_{\text{bar}}$. For models without bulges (empty squares), a larger gas fraction (thick squares) tends to decrease $t_{\text{bar}}$ confirming early results by Athanassoula et al. (2013) who found a similar outcome for systems with up to 50% of gas. For more gas-rich systems (e.g. gas fractions equal or larger than 75%) Athanassoula et al. (2013) witnessed an inverse trend considering the initial growth of the bar. They also find that such trends significantly depend on the halo triaxiality (see their Fig. 7), something we do not test with our present simulations. We note that the delay in the early bar growth due to higher gas fractions in our simulations does not simply hold for our bulged models, suggesting, together with the above-mentioned results from Athanassoula et al. (2013) that a more relevant parameter may drive the initial bar growth.

Fujii et al. (2018) recently suggested to use of $f_{\text{disc}}$, the ratio between the disc mass and the total galaxy mass within $2.2 R_{\text{disc}}$ as a relevant driving parameter for the bar growth timescale. This parameter $f_{\text{disc}}$ is admittedly one way to quantify the importance of the stellar disc within the global potential (baryonic + dark matter) that may be connected with the onset of bar instabilities. Fujii et al. (2018) used a series of pure N-body simulations of systems with stellar masses of a few $10^{10} M_{\odot}$, with varying bulge to disc ratios, using 500 million particles (8 million for the baryonic disc) and 10 pc softening length to study timescales associated with bar and spiral growth (see also Valencia-Enríquez et al. 2017). They suggested a relation between $f_{\text{disc}}$ and $t_{\text{bar}}$. Bland-Hawthorn et al. (2023) reviewed the Fujii et al. (2018) results by running simulations both using dry runs (stars and dark matter particles only) and a few wet runs (including gas but excluding cooling and star formation) using the AMR code RAMSES, with spatial sampling of 6 and 12 pc: they mostly confirmed the existence of a relation between $f_{\text{disc}}$ and $t_{\text{bar}}$, the “Fujii relation”, albeit with a large scatter for $f_{\text{disc}}$ larger than 0.5 (their Fig. 4). Bland-Hawthorn et al. (2023) further suggested that the addition of gas tends to decrease the value of $t_{\text{bar}}$, as also witnessed in our bulgeless simulations (but not in our bulged models; and see Athanassoula et al. 2013).

We take the opportunity of our sample of simulations to test further the validity range of the “Fujii relation”. Our simulations have similar resolution and particle numbers, and use a very similar code than Bland-Hawthorn et al. (2023) but include cooling and star formation. As compared to the original Fujii et al. (2018) simulations that are pure N-body via a tree-code (Barnes & Hut 1986), our simulations have 8 times fewer particles in the disc (and a mass resolution for the dark matter halo about 100 times larger) and again include gas and star formation. Figure 10 specifically illustrates the trend of $t_{\text{bar}}$ as a function of $f_{\text{disc}}$ for our set of simulations, emphasising the fact that we are so far only probing values of $f_{\text{disc}}$ larger than 0.5 by design.

The left (resp., right) panel shows the Fujii et al. (2018) relation for $f_{\text{disc}}$ (dashed lines) computed with the star and gas (resp. only star) component. The inclusion of the gas as part of the disc

<figure><img src="image_10.png" /><figcaption>Fig. 9. Typical bar formation time $t_{\text{bar}}$ as a function of the stellar mass for all 16 simulations. The colour of the squares corresponds to the different stellar masses. The size of each square represents 200 Myr (i.e. ±100 Myr), hence illustrating the uncertainty in $t_{\text{bar}}$ (see Sect. 4.2).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Typical bar formation time $t_{\text{bar}}$ as a function of $f_{\text{disc}}$ as defined in Bland-Hawthorn et al. (2023). The dashed lines show the relation given by Fujii et al. (2019). $f_{\text{disc}}$ is derived using the sole stellar disc (right panel) or the total baryonic content (i.e. stars and gas; left panel). In the right panel, $f_{\text{disc}}$ (stars only) decreases when the gas fraction increases as expected (the relative contribution of the stellar disc gets smaller), while in the left panel, $f_{\text{disc}}$ (total) increases with the gas fraction (as the total baryonic contribution increases). As in Fig. 9, the size of each square represents 200 Myr (i.e. ±100 Myr), hence illustrating the uncertainty in $t_{\text{bar}}$ (see Sect. 4.2).</figcaption></figure>

potential (left panel) increases $f_{\text{disc}}$ only slightly, as expected: for galaxies without a bulge, increasing the gas fraction also systematically decreases $t_{\text{bar}}$: that can be understood in terms of the destabilising effect of the additional disc of gas. For the more massive simulations with bulges, the trend seems the opposite, illustrating that the addition of an ellipsoid dominates the budget for $t_{\text{bar}}$. Finally, the simulation at the lowest stellar mass bin and with a bulge shows a smaller $t_{\text{bar}}$ when the gas fraction is increased.

At fixed $f_{\text{disc}}$, our sample of simulations seems to overall lie a factor of two below the Fujii et al. (2019) relation, with differences of about 500 Myr or more, well beyond the expected uncertainty associated with the measured $t_{\text{bar}}$ values. As emphasised, one difference with the results from Fujii et al. (2019) is that our simulations contain gas. While, as mentioned above, the bar formation time may depend on the gas fraction, this discrepancy does not change significantly if we include the gas component in the value of $f_{\text{disc}}$. Another significant difference between our simulations and those by Fujii et al. (2019) is that the latter used pure N-body simulations with a large number of particles (with a softening length of about 10 pc). Still, we are in a regime where we should have enough mass resolution, at least for the
3 larger mass bins (leading to a minimum of 2 million particles), and we do not detect a trend in stellar mass in the context of the $t_{\text{bar}}$ discrepancy with Fujii et al. (2019). Even considering the added systematic error of ±100 Myr (see the beginning of the present section), we do not retrieve the Fujii et al. (2018) relation. We also need to emphasise the fact that simulations performed by Bland-Hawthorn et al. (2023) use a similar setup as ours (and the same AMR code, RAMSES), and do not see such a discrepancy with Fujii et al. (2019, their Fig. 4).

More fundamentally, other differences could be at the root of this discrepancy, including more concentrated discs and dark matter distributions in Bland-Hawthorn et al. (2023), different initial velocity distributions (anisotropies), or the lack of cooling, star formation and feedback in the latter experiments. As emphasised by Athanassoula et al. (2013), other parameters, including the shape of the halo and the initial dynamical set up for the disc, can strongly influence the growth of unstable modes, including the bar and spirals.

Our results confirm and extend previous studies showing that $t_{\text{bar}}$ seems to depend on many details in the initial conditions, including the steepness of the inner potential (e.g., via the index of the Sersic bulge), the overall contribution of the disc components (w.r.t. other baryonic and non-baryonic components; see also Athanassoula 2013). This may naively suggest that $t_{\text{bar}}$ is not fully described by a single parameter such as $f_{\text{disc}}$. However, we cannot yet conclude, as we would need to probe a more extended set of simulations, including, for example, models with varying properties (e.g. anisotropies) at fixed $f_{\text{disc}}$, and models with value of $f_{\text{disc}}$ below 0.5 to detect the exponential increase of $t_{\text{bar}}$ (and where values of $t_{\text{bar}}$ are significantly larger than about 2 Gyr).

### 4.3. Evolutionary phases of isolated barred discs

Since the characteristic bar formation time seems to play a major role in the formation of the central gas reservoir, we have normalised the time by $t_{\text{bar}}$ and have introduced the dimensionless parameter $\tau = t/t_{\text{bar}}$. We have also normalised the gas mass by the initial amount of gas contained within 1 kpc as shown in Fig. 11. In this figure, we can see that the bars and SFRs share a similar time evolution with a slight trend of increasing amplitude as a function of the stellar mass, for all simulations with initial stellar masses equal or above $10^{10} M_{\odot}$ (see Sect. 4.5 for a discussion on simulations in the lower stellar mass bin). These similarities between barred systems confirm the relevance of the choice of $\tau$ allowing us to naturally account for the delay caused by the presence of a bulge.

Except for the models G162M110F10L5B10 and G178M110F20L5B10, the normalised gas mass always peaks with $\tau$ between 1 and 1.5. This peak also roughly coincides with the maximum of the bar strengths and the SFRs. Those are associated with a central starburst leading to an increase in the number of new stars and a decrease in the gas mass inside the central 1 kpc. For simulations with $M_{\star} \ge 10^{10} M_{\odot}$, we can roughly decompose the time evolution of the gas mass inside 1 kpc in three main phases referenced by the bar formation timescale $t_{\text{bar}}$:

1. From $\tau = 0$ to ~1, the bar forms and the gas mass inside the central 1 kpc region is nearly constant or increase weakly. The star formation rate slowly increases with time (by about 1 order of magnitude in 500 Myr).

2. From $\tau = 1$ to ~1.5, the bar is strong (by definition, $A_2 > 0.2$) and significantly influences the redistribution of gas. Gas is being funneled towards the centre, leading to a more rapid

increase of the gas mass, that triggers a starburst inside that region. The star formation evolution then reaches a plateau.

3. From $\tau = 1.5$ to ~2, the bar keeps its high amplitude, and a steeper decrease in the amount of gas in the central 1 kpc region is associated with a slow decline in star formation rate. This gas depletion phase marks the emergence of a small gaseous and stellar central mass overdensity, a discy structure that is fuelled by the bar. Follow-up star formation tends to predominantly occur within this inner disc region building up further beyond $\tau = 2$.

The fact that those phases are similar for all barred models having an initial stellar mass bigger or equal to $10^{10} M_{\odot}$ suggests that the fuelling and the consumption of gas inside the central 1 kpc region is driven by the same physical phenomena and evolution, and only depends weakly on the initial gas fraction and the presence or absence of a central ellipsoid (within the range of parameters probed by our simulations). While there are clear local differences between simulations, it is still an interesting and relevant result as it sets the stage of the evolution of barred systems which could serve as an “isolated case” reference for further studies. Note that we are not here discussing the evolution beyond 3 Gyr, which could, for instance, lead to a secondary growth of the stellar bar, and to the emergence of large bars (see e.g. Schinnerer et al. 2023; Bland-Hawthorn et al. 2023).

### 4.4. Gas structure in the central 1 kpc

In Sect. 4.1.3, we briefly discussed the property and global evolution in the central kpc. In Fig. 12, we now present the gas distributions and reservoir by providing a gas density map for all 16 simulations at $\tau = 2$, in order to illustrate the emergence of distinct inner gas structures. In that figure, the four lines correspond to the four different stellar mass bins (colour-coded labels), the models having a lower (resp. higher) gas fraction ($\alpha$) being in the two left (resp. right) columns. The models without bulges ($\beta = 0$) are given by the first and third columns, whereas the models with bulges (here $\beta = 10$) are given by the second and the fourth columns.

We observe a distinct central concentration of gas for almost all barred models having a stellar mass larger than $10^{9.5} M_{\odot}$. For these stellar mass ranges, only model G162M110F10L5B10 does not show any bar structure and central gas concentration. We see that models with a bulge have a more extended gas reservoir. While the bulge delays the bar formation, once the bar is formed, the bulge does not prevent the formation of a gas reservoir and actually allows the emergence of a prominent inner structure. We discuss the detailed properties and growth of those structures in a subsequent paper, but we can already suggest here that the extent of such an inner disc (or rings, see Fig. 12) closely follows the change of the inner mass concentration.

### 4.5. The conditional onset of inner stellar discs

As illustrated in Fig. 8, we do not witness strong differences between simulations above and below the $10^{10} M_{\odot}$ values in terms of the amplitude of the bar ($A_2$) or the global SFR. The main difference between those two sub-samples shows up in the evolution of the gas in the central 1 kpc. The decrease of the central gas mass leading to a long-term gas depletion mentioned as Phase 3 in the previous section is observed for all galaxies with an initial mass of $10^{10} M_{\odot}$ and above. For the less massive galaxies (i.e. models G001M095F10L2B00, G002M095F10L2B10, G013M095F20L2B00 and G014M095F20L2B10, having an
initial stellar mass of $10^{9.5} M_\odot$), the gas mass stays roughly constant or even increases during that phase. This means that the channelling of gas towards the central region (and subsequent star formation) is still occurring but seems to proceed in a different manner.

To understand this better, we turn to the spatial distribution of star formation (and new stars). Figure B.1 shows (second row in each panel) that for the lowest stellar mass models, star formation occurs mainly along and inside the bar until the end of the simulation without showing the emergence of a central gas reservoir. For the higher stellar mass models (Figs. B.2–B.4), star formation occurs also along and inside the bar until $\tau = 1.5–2$, when we start to see the emergence of a central gas overdensity. Subsequent star formation is mainly distributed inside this central gas reservoir. This result also connects with the fact that distinct inner streaming lanes, as well as disc or ring structures form and grow for the more massive end systems, while the lower mass systems exhibit more heterogeneous gas distribution within the bar region. This is quite an intriguing result as it means that star formation is either prevented or triggered in significantly different regions for the lowest mass simulation bin.

This trend has been reported by Fraser-McKelvie et al. (2020) using emission-line mapping via integral-field spectroscopic observations for a large sample of nearby galaxies (but see also Díaz-García et al. 2020). The stellar mass value of $10^{10} M_\odot$ at which this occurs also coincides with the one in our simulation set. It is worth mentioning that most of the galaxies presented in the above survey have a much higher gas fraction compared with the PHANGS sample. Their typical gas fraction<fnref n="3" /> in the $10^{9.5–10} M_\odot$ stellar mass bin is above ~30% and goes up to ~60–70%. Since we tuned our initial conditions to an observed sample of star-formation main sequence galaxies, this may mean either that the change is built in our set of morphological parameters (e.g. scale lengths) or that it is related to the varying relative contributions of physical processes at play, for example, the influence of feedback versus the strength of the gravitational potential, associated with a change of stellar mass.

The physical origin of this change in the regime between the lowest and highest stellar masses is not yet fully understood. Fraser-McKelvie et al. (2020) extrapolated from earlier simulations (Emsellem et al. 2015) that shear may be an important ingredient in setting up such differences. While this may play a role, we note that two of our simulations, one at $10^{9.5} M_\odot$ and one at $10^{10} M_\odot$, share almost exactly the same radial mass gradient (only the mass scaling is different), and those two do present the above-mentioned change in the evolved morphology. We thus suggest that the main driver for such a difference lies with the relative contribution of the stellar-driven feedback within those scaled gravitational potential (see e.g. Collins & Read 2022, and references therein). This will be specifically discussed in a subsequent paper.

# 5. Summary and conclusion

In this work, we have performed a set of 16 three-dimensional high-resolution hydro-dynamical simulations (see Table 1) using the RAMSES AMR code to study and characterise the building and evolution of central gas reservoir in nearby main sequence disc galaxies. We have designed this grid based on the PHANGS-ALMA sample. We made use of four control parameters (i.e. stellar mass, gas fraction, scale length for the star dis-

tribution, and bulge mass fraction), and for each model started from axisymmetric initial conditions including gas, stars and dark matter.

We have further quantified the characteristic formation time for the bar $t_{\text{bar}}$ in our simulations using a reference value of 0.2 for $A_2$, and compared them with the relation suggested by Fujii et al. (2019). We have found that our simulations are located significantly below the relation (have smaller $t_{\text{bar}}$ at fixed $f_{\text{disc}}$). We could not conclude robustly on the origin of such a discrepancy as it may both reflect the intrinsic scatter of such a relation and the imposed variety in the initial conditions (i.e. halo concentration, anisotropy). A larger set of simulations including gas and star formation covering lower values of $f_{\text{disc}}$ (hence larger values of $t_{\text{bar}}$) are needed to confirm this trend.

We have found that models G162M110F10L5B10 and G178M110F20L5B10 are expected to be significantly more stable against bar formation (according to the criterion established by Efstathiou et al. 1982; but see Sect. 4.2). Note that the criterion we use to decide if a bar is formed or not (i.e. $A_2 = 0.2$) is met for model G178M110F20L5B10 despite the lack of an apparent proper bar structure: we interpret that case as $A_2$ capturing the evolution of the strong spiral arms.

We have studied the impact of three control parameters on the evolution of the central 1 kpc region. The mass inside the central gas reservoir naturally increases with the initial stellar and gas mass. The presence of a bulge delays the formation of the bar (i.e. values of $t_{\text{bar}}$ are larger) and thus the formation of the gas reservoir, but does not prevent its formation.

The global evolution of the 12 models having a stellar mass $\ge 10^{10} M_\odot$ can be roughly described using a dimensionless bar formation time parameter $\tau = t/t_{\text{bar}}$, including three subsequent phases:

1. A formation phase: from $\tau = 0$ to ~1, the bar forms and the gas mass inside the central 1 kpc region is nearly constant or increases weakly. The star formation rate slowly increases with time.

2. A fueling and growth phase: from $\tau = 1$ to ~1.5, the bar is strong enough (i.e. $A_2 > 0.2$) and starts to transport gas towards the centre, leading to a steeper increase of the gas mass and a starburst inside that region. The star formation reaches a plateau.

3. A depletion phase: from $\tau = 1.5$ to ~2, the bar stays strong, and a steep decrease in the amount of gas in the central 1 kpc region is associated with a slow decline in star formation rate. This phase witnesses the emergence of a central stellar mass seed growing into a more extended inner stellar and gas structure (discs and rings). The sizes of the inner discy structure seem to vary with e.g. the initial stellar mass and the presence or absence of an ellipsoid.

Simulated galaxies with initial stellar masses below $10^{10} M_\odot$, falling in the lowest mass bin, exhibit differences with respect to the more massive galaxies in the sample as for their inner gas structures, the spatial distribution of star forming regions for $\tau \ge 2$ and the gas depletion timescales in the central kpc. More specifically, the two first above-mentioned phases are also witnessed for the lower stellar mass models ($M_\star = 10^{9.5} M_\odot$), but we do not observe Phase 3, that is, a steep decrease of the gas mass inside the 1 kpc central region at $\tau > 1.5$. The gas mass in the central kpc thus stays roughly constant and is associated with a relatively constant SFR over time (until at least $\tau = 4$). We have also shown evidence for two distinct star formation distributions for $\tau \ge 2$ below and above the $10^{10} M_\odot$ initial stellar mass of our models. For models with $M_\star < 10^{10} M_\odot$ (hence only for the lower mass bin), we observe that the star

<footnote n="3">³ Value added Catalogs (SDSS): https://data.sdss.org/sas/dr17/-env/MANGA_HI</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Evolution of the maximum of the $A_2$ coefficient (top panels), the mass of gas within the central 1 kpc normalised by the initial mass of gas within the same radius (middle panels), and the SFR (bottom panels). The evolution is shown through the dimensionless parameter $\tau$, which is the ratio between the time of the simulations and the corresponding time when $A_2$ reaches the value of 0.2. The dashed vertical lines and shaded coloured areas show peculiar values of $\tau$ we use to describe the phases of the fuelling (i.e., $\tau \in 0-1$ (blue area); $\tau \in 1-1.5$ (purple area); $\tau \in 1.5-2$ (orange area)). The colour code and the meaning of the different lines are the same as in Fig. 8.</figcaption></figure>

formation mainly occurs along the bar, whereas star formation is mostly triggered inside the inner gas concentration for models with $M_{\star} \ge 10^{10} M_{\odot}$. This trend echoes the reporting made by Fraser-McKelvie et al. (2020) using $H\alpha$ two-dimensional map-

ping of a sample of nearby galaxies. We suggest that this relates to the relative contribution of stellar-driven feedback within a given gravitational potential: this will be probed and discussed in detail in a subsequent paper (Verwilghen, in prep.). The subset
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Density map of gas for the 16 simulations at the corresponding time for which $\tau = 2$. The first three rows show a box of 10 kpc side while the last row shows a box of 20 kpc side. The model numbers are colour-encoded according to the corresponding stellar mass (from the lightest to the most massive from the top to the bottom, i.e. 9.5, 10, 10.5, and 11 in $\log_{10}(M_{\odot})$). The two left columns illustrate the model with a gas fraction of 10% while the two right columns illustrate the models with a gas fraction of 20%. The odd (even) reference number accounts for models without (with) a bulge.</figcaption></figure>

of simulations presented here only covers a restricted range of observational properties and are constrained by a fixed and limited set of initial structural parameters. We can thus already presume that stellar mass may not be the only or even the prime driver of the differences we observe in our simulations (see e.g. Díaz-García et al. 2020).

In this paper, we provide a first pilot study at the structures and time evolution of a set of simulations probing the star formation main sequence of nearby disc galaxies. Subsequent papers in this series will focus on examining in detail the advent and growth of the above-mentioned inner gas structures (and their associated stellar content), as well as study in more detail the origin of the change of regime around $10^{10} M_{\odot}$ observed in Fraser-McKelvie et al. (2020) and reproduced with this first subset of models. We have already planned for an extended sam-

ple of 54 models better covering the range of observed properties of the PHANGS-ALMA sample in this stellar mass range: a more full-fledged account of this “complete” set of hydrodynamical simulations will be presented when available. Such studies can then serve as a benchmark for simulations embedded in a more comprehensive environment, including, for example, gas accretion, interactions or evolution in a cosmological context.

*Acknowledgements.* We thank the anonymous referee for a report that helped clarify this manuscript. P.V. acknowledges support from the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (German Research Foundation) under Germany’s Excellence Strategy - EXC-2094-390783311. The simulations in this paper have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP). We are grateful for the support by Alexey Krukau and Margarita Petkova through
C2PAP. F.R. acknowledges support provided by the University of Strasbourg Institute for Advanced Study (USIAS), within the French national programme Investment for the Future (Excellence Initiative) IdEx-Unistra. M.V. is supported by the Fondazione ICSC National Recovery and Resilience Plan (PNRR), Project ID CN-00000013 "Italian Research Center on High-Performance Computing, Big Data and Quantum Computing" funded by MUR - Next Generation EU. M.V. also acknowledges partial financial support from the INFN Indark Grant. R.S.K. acknowledges financial support from the European Research Council via the ERC Synergy Grant "ECOGAL" (project ID 855130), from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) "STRUCTURES", and from the German Ministry for Economic Affairs and Climate Action in project "MAINN" (funding ID 50002206). R.S.K. thanks for computing resources provided by the Ministry of Science, Research and the Arts (MKW) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and also for data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. K.D. acknowledges financial support from "BiD4BEST" - European Innovative Training Network (ITN) funded by the Marie Skłodowska-Curie Actions (860744) in Horizon 2020 and by the COMPLEX project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. J.N. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694343). K.G. is supported by the Australian Research Council through the Discovery Early Career Researcher Award (DECRA) Fellowship (project number DE220100766) funded by the Australian Government. K.G. is supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. M.C.S. acknowledges financial support from the European Research Council under the ERC Starting Grant "GalFlow" (grant 101116226) and from the Royal Society (URF) R\221118).

## References

Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25

Athanassoula, E. 1992a, MNRAS, 259, 328

Athanassoula, E. 1992b, MNRAS, 259, 345

Athanassoula, E. 2002, in Disks of Galaxies: Kinematics, Dynamics and Peturbations, eds. E. Athanassoula, A. Bosma, & R. Mujica, ASP Conf. Ser., 275, 141

Athanassoula, E. 2008, MNRAS, 390, L69

Athanassoula, E. 2013, in Secular Evolution of Galaxies, eds. J. Falcón-Burroso, & J. H. Knapen, 305

Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35

Athanassoula, E., Machado, R. E. G., & Rodionov, S. A. 2013, MNRAS, 429, 1949

Barazza, F. D., Jogee, S., & Marinova, I. 2008, ApJ, 675, 1194

Barnes, J., & Hut, P. 1986, Nature, 324, 446

Bilitewski, T., & Schönrich, R. 2012, MNRAS, 426, 2266

Binney, J. 2020, MNRAS, 496, 767

Binney, J., Tremaine, S., & Freeman, K. 2009, Phys. Today, 62, 56

Bland-Hawthorn, J., Tepper-Garcia, T., Agertz, O., & Freeman, K. 2023, ApJ, 947, 80

Boone, F., Baker, A. J., Schinnerer, E., et al. 2007, A&A, 471, 113

Buta, R. J., Sheth, K., Athanassoula, E., et al. 2015, ApJS, 217, 32

Cheung, E., Trump, J. R., Athanassoula, E., et al. 2015, MNRAS, 447, 506

Collins, M. L. M., & Read, J. I. 2022, Nat. Astron., 6, 647

Combes, F. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, P. Di Matteo, F. Herpin, et al., 223

Combes, F., & Sanders, R. H. 1981, A&A, 96, 164

Combes, F., García-Burillo, S., Boone, F., et al. 2004, A&A, 414, 857

Comerón, S., Knapen, J. H., Beckman, J. E., et al. 2010, MNRAS, 402, 2462

Díaz-García, S., Salo, H., & Laurikainen, E. 2016, A&A, 596, A84

Díaz-García, S., Moyano, F. D., Comerón, S., et al. 2020, A&A, 644, A38

Efstathiou, G., Lake, G., & Negroponte, J. 1982, MNRAS, 199, 1069

Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87

Emsellem, E., Monnet, G., & Bacon, R. 1994a, A&A, 285, 723

Emsellem, E., Monnet, G., Bacon, R., & Nieto, J. L. 1994b, A&A, 285, 739

Emsellem, E., Renaud, F., Bournaud, F., et al. 2015, MNRAS, 446, 2468

Eskridge, P. B., Frogel, J. A., Pogge, R. W., et al. 2000, AJ, 119, 536

Fabian, A. C. 2012, ARA&A, 50, 455

Ferrière, K., Gillard, W., & Jean, P. 2007, A&A, 467, 611

Fragkoudi, F., Di Matteo, P., Haywood, M., et al. 2017, A&A, 606, A47

Framer-McKelvie, A., Aragón-Salamanca, A., Merrifield, M., et al. 2020, MNRAS, 495, 4158

Fujii, M. S., Bédorf, J., Baba, J., & Portegies Zwart, S. 2018, MNRAS, 477, 1451

Fujii, M. S., Bédorf, J., Baba, J., & Portegies Zwart, S. 2019, MNRAS, 482, 1983

Fukuda, H., Wada, K., & Habe, A. 1998, MNRAS, 295, 463

Gadotti, D. A. 2008, in Formation and Evolution of Galaxy Bulges, eds. M. Bureau, E. Athanassoula, & B. Barbuy, 245, 117

Gadotti, D. A., Bittner, A., Falcón-Burroso, J., et al. 2020, A&A, 643, A14

García-Burillo, S., Combes, F., Schinnerer, E., Boone, F., & Hunt, L. K. 2005, A&A, 441, 1011

Genel, S., Fall, S. M., Hernquist, L., et al. 2015, ApJ, 804, L40

Goz, D., Monaco, P., Murante, G., & Curir, A. 2015, MNRAS, 447, 1774

Guo, M., Stone, J. M., Kim, C.-G., & Quataert, E. 2023, ApJ, 946, 26

Haardt, F., & Madau, P. 1996, ApJ, 461, 20

Henshaw, J. D., Barnes, A. T., Battersby, C., et al. 2023, in Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, ASP Conf. Ser., 534, 83

Hickox, R. C., Wardlow, J. L., Smail, I., et al. 2012, MNRAS, 421, 284

Hilmi, T., Minchev, I., Buck, T., et al. 2020, MNRAS, 497, 933

Ho, L. C., Filippenko, A. V., & Sargent, W. L. W. 1997, ApJ, 487, 591

Hopkins, P. F., & Quataert, E. 2010, MNRAS, 407, 1529

Izumi, T., Wada, K., Imanishi, M., et al. 2023, Science, 382, 554

Kataria, S. K., & Das, M. 2018, MNRAS, 475, 1653

Kaufmann, D. E., & Contopoulos, G. 1996, A&A, 309, 381

Kennicutt, R. C. 2007, in Triggered Star Formation in a Turbulent ISM, eds. B. G. Elmegreen, & J. Palous, 237, 311

Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 802, 99

Kormendy, J., & Kennicutt, R. C. 2004, ARA&A, 42, 603

Kraljic, K., Bournaud, F., & Martig, M. 2012, ApJ, 757, 60

Krolík, J. H. 1999, Active Galactic Nuclei: From the Central Black Hole to the Galactic Environment (Princeton: Princeton University Press)

Kun, E., Keresztes, Z., Simkó, A., Szűcs, G., & Gergely, L. Á. 2017, A&A, 608, A42

Lacey, C. G., & Fall, S. M. 1985, ApJ, 290, 154

Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122

Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24

Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, ApJS, 257, 43

Lin, L.-H., Yuan, C., & Buta, R. 2008, ApJ, 684, 1048

Longmore, S. N., Bally, J., Testi, L., et al. 2013, MNRAS, 429, 987

Ludlow, A. D., & Angulo, R. E. 2017, MNRAS, 465, L84

Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285

Moriondo, G., Giovanardi, C., & Hunt, L. K. 1998, A&AS, 130, 81

Ostriker, E. C., & Kim, C.-G. 2022, ApJ, 936, 137

Padovani, P. 2017, Front. Astron. Space Sci., 4, 35

Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&ARV, 25, 2

Pezzulli, G., & Fraternali, F. 2016, Astron. Nachr., 337, 913

Porciani, C., Dekel, A., & Hoffman, Y. 2002, MNRAS, 332, 325

Queirejeta, M., Schinnerer, E., Meidt, S., et al. 2021, A&A, 656, A133

Reines, A. E., & Volonteri, M. 2016, Am. Astron. Soc. Meet. Abstr., 227, 119.01

Renaud, F., Bournaud, F., Emsellem, E., et al. 2013, MNRAS, 436, 1836

Romeo, A. B., Agertz, O., & Renaud, F. 2023, MNRAS, 518, 1002

Saha, K., & Naab, T. 2013, MNRAS, 434, 1287

Sanders, R. H. 1977, ApJ, 217, 916

Sanders, D. B., Soifer, B. T., Elias, J. H., Neugebauer, G., & Matthews, K. 1988, ApJ, 328, L35

Schinnerer, E., Emsellem, E., Hensaw, J. D., et al. 2023, ApJ, 944, L15

Schmidt, M. 1959, ApJS, 129, 243

Schmidt, T. M., Bigiel, F., Klessen, R. S., & de Blok, W. J. G. 2016, MNRAS, 457, 2642

Sellwood, J. A., & Evans, N. W. 2001, ApJ, 546, 176

Sellwood, J. A., & Masters, K. L. 2022, ARA&A, 60, 36

Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41

Shlosman, I., Frank, J., & Begelman, M. C. 1989, Nature, 338, 45

Sormani, M. C., & Barnes, A. T. 2019, MNRAS, 484, 1213

Sormani, M. C., Binney, J., & Magorrian, J. 2015, MNRAS, 451, 3437

Sormani, M. C., Barnes, A. T., Sun, J., et al. 2023, MNRAS, 523, 2918

Sormani, M. C., Sobacchi, E., & Sanders, J. L. 2024, MNRAS, 528, 5742

Sotnikova, N. Y., Reshetnikov, V. P., & Mosenkov, A. V. 2012, Astron. Astrophys. Trans., 27, 325

Stewart, K. R., Brooks, A. M., Bullock, J. S., et al. 2013, ApJ, 769, 74

Stuber, S. K., Schinnerer, E., Williams, T. G., et al. 2023, A&A, 676, A113

Sun, J., Leroy, A. K., Rosolowsky, E., et al. 2022, AJ, 164, 43

Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253

Teyssier, R. 2002, A&A, 385, 337

Übler, H., Naab, T., Oser, L., et al. 2014, MNRAS, 443, 2092

Valencia-Enríquez, D., Puerari, I., & Chaves-Velasquez, L. 2017, Rev. Mex. Astron. Astrofis., 53, 257

Valentini, M., Murante, G., Borgani, S., et al. 2017, MNRAS, 470, 3167

Wada, K. 1994, PASJ, 46, 165
**Appendix A: Estimation of the fit values**

As introduced in Sect. 2 and illustrated in Fig. 2, we used a subjective rating process to sort the different results coming from the exponential and Sersic fits of the 1D stellar and gas density profiles. Fig. A.1 briefly illustrates this rating for two profiles that were attributed the values 0 and -1, and thus considered as ‘satisfying’ and ‘bad’, respectively. The left panel (NGC 7476) represents the case for which we considered the fit to be satisfying (value 0) and the right panel (NGC 3239) represents the case we considered the fit to be bad (value -1). The ‘bad’ case shows one example of a density profile for which the adopted functional form (exponential plus Sersic) is clearly not appropriate. While this classification is subjective, we were mostly interested in the global trends, and it did not strongly influence the choice of the control parameter values themselves (see Fig. 3).

<figure><img src="image_14.png" /><figcaption>**Fig. A.1.** Example of fits for which we attributed the value 0 (satisfying but with significant residuals; top panel) and -1 (bad; bottom panel).</figcaption></figure>

**Appendix B: Evolution of the models as a function of the parameter τ**

As mentioned in Sect. 4.3, we have identified three distinct phases in the building and evolution of the gas reservoir by using the dimensionless parameter $\tau = t/t_{bar}$, where $t_{bar}$ is the typical bar formation time. The following figures illustrate the time evolution of the system for the 16 simulations providing snapshots at four different values of $\tau$ (i.e. 1, 1.5, 2, and End, for the end of the simulation). For each simulation, we gather the four models belonging to the same initial stellar mass into one figure. Each panel corresponds to one labelled simulation: each include a zoom (box of 10 kpc on a side) and 4 row with the surface densities of (from top to bottom) the gas, of the (new) stars formed within 50 Myr of the corresponding snapshot time, of all old (initial) stars and of all (new) stars (formed since the beginning of the simulation). Fig. B.1, B.2, B.3, B.4 refer to the models with a stellar mass of 9.5, 10, 10.5 and 11 $\log_{10}(M_{\odot})$, respectively.

In the top rows of the most massive initial stellar mass bins (Fig. B.2, B.3, B.4) we witness the formation of a clearly marked inner gas structure (except for model G162M110F10L5B10, and for model G178M110F20L5B10), while simulations presented in Fig. B.1 do not exhibit such gas structures. We also observe a change in the distribution of new stars that tend to emerge along the bar structure for the lowest stellar mass models, and are more localised (inner region) for the higher stellar mass models from $\tau \ge 2$. There is an associated difference in the structure of the bar with the more massive galaxies presenting a more distinct and rounder stellar concentration (within the bars). This central concentration could be caused by the emergence of an ILR as discussed in Sect. 4.1.2 and 4.5.
<figure><img src="image_15.png" /><figcaption>**Fig. B.1.** Time evolution as a function of τ of the gas surface density (top row, GAS), formed stars' surface density (second row, NS), old stars' surface density (third row, OS), and cumulative formed stars' surface density (bottom row, CNS) of the four lowest stellar mass models (10⁹.⁵ M☉). Each panel shows a box with a side length of 10 kpc.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. B.2. Same as Fig. B.1 for the four $10^{10} M_{\odot}$ stellar mass models.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. B.3. Same as Fig. B.1 for the four $10^{10.5} M_{\odot}$ stellar mass models.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. B.4.** Same as Fig. B.1 for the four $10^{11} M_{\odot}$ stellar mass models. The box size has been extended to 20 kpc for this most massive stellar bin.</figcaption></figure>