# How Nested Bars Enhance, Modulate, and are Destroyed by Gas Inflows

Zhi Li Shanghai Key Lab for Astrophysics, Shanghai Normal University, 100
Guilin Road, Shanghai 200234, P.R. China  Department of Astronomy, School of
Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road,
Shanghai 200240, P.R. China  Tsung-Dao Lee Institute, Shanghai Jiao Tong
University, Shanghai 200240, P.R. China  Min Du Min Du<EMAIL_ADDRESS>Department of Astronomy, Xiamen University, Xiamen, Fujian 361005, P.R. China
Victor P. Debattista Jeremiah Horrocks Institute, University of Central
Lancashire, Preston PR1 2HE, UK  Juntai Shen Juntai Shen<EMAIL_ADDRESS>Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong
University, 800 Dongchuan Road, Shanghai 200240, P.R. China  Key Laboratory
for Particle Astrophysics and Cosmology (MOE) / Shanghai Key Laboratory for
Particle Physics and Cosmology, Shanghai 200240, P.R. China  Hui Li Department
of Astronomy, Tsinghua University, Beijing 100084, P.R. China  Department of
Astronomy, Columbia University, New York, NY 10027, USA  Jie Liu Shanghai
Astronomical Observatory, Chinese Academy of Science, 80 Nandan Road,
Shanghai, 200030, P.R. China  Mark Vogelsberger Department of Physics & Kavli
Institute for Astrophysics and Space Research, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Angus Beane Center for Astrophysics $|$
Harvard & Smithsonian, Cambridge, MA, USA  Federico Marinacci Department of
Physics & Astronomy “Augusto Righi”, University of Bologna, via Gobetti 93/2,
40129 Bologna, Italy  Laura V. Sales Department of Physics and Astronomy,
University of California, Riverside, CA, 92521, USA

###### Abstract

Gas flows in the presence of two independently-rotating nested bars remain not
fully understood, which is likely to play an important role in fueling the
central black hole. We use high-resolution hydrodynamical simulations with
detailed models of subgrid physics to study this problem. Our results show
that the inner bar in double-barred galaxies can help drive gas flow from the
nuclear ring to the center. In contrast, gas inflow usually stalls at the
nuclear ring in single-barred galaxies. The inner bar causes a quasi-periodic
inflow with a frequency determined by the difference between the two bar
pattern speeds. We find that the star formation rate is higher in the model
with two bars than in that with one bar. The inner bar in our model gradually
weakens and dissolves due to gas inflow over a few billion years. Star
formation produces metal-rich/$\alpha$-poor stars which slows the weakening of
the inner bar, but does not halt its eventual decay. We also present a
qualitative comparison of the gas morphology and kinematics in our simulations
with those of observed double-barred galaxies.

galaxies: ISM — galaxies: kinematics and dynamics — galaxies: structures —
galaxies: hydrodynamics

††software:  AREPO (Springel, 2010; Weinberger et al., 2020), NumPy (Harris et
al., 2020), SciPy (Virtanen et al., 2020), Matplotlib (Hunter, 2007), Jupyter
Notebook (Kluyver et al., 2016)

##  1 Introduction

Galaxies with two nested stellar bars are termed double-barred galaxies (S2Bs,
see Laine et al., 2002; Erwin, 2004; Shlosman, 2005; Buta et al., 2015). About
20% of disk galaxies appear to be S2Bs (Laine et al., 2002; Erwin, 2011). In
such systems the two bars are observed to be randomly oriented with respect to
each other (Buta & Crocker, 1993; Friedli & Martinet, 1993), implying two
independent rotation speeds that are consistent with theoretical expectations
(e.g. Maciejewski & Sparke, 2000; Debattista & Shen, 2007; Du et al., 2015).
Inner bars in S2Bs are also dubbed as “nuclear” or “secondary” bars in some
works. These terms may be less appropriate since observations have shown that
inner bars can be as large as some short single bars (semi-major axes of
$2$-$3\;{\rm kpc}$, see Erwin, 2005; de Lorenzo-Cáceres et al., 2020). Recent
studies have suggested that inner bars are probably scaled-down replicas of
outer bars based on resolved spatial maps of the mean stellar population
properties (Méndez-Abreu et al., 2019; de Lorenzo-Cáceres et al., 2019b;
Bittner et al., 2021).

The gas inflow driven by the large-scale bar generally stalls at a starburst
nuclear ring with a typical size of a few hundred parsecs, preventing
accretion onto the supermassive black hole (SMBH, e.g. Fanali et al., 2015; Li
et al., 2017; Tress et al., 2020). Inner bars have been hypothesized to be an
important mechanism for further removing gas angular momentum, in a manner
similar to their large-scale counterparts (e.g. Shlosman et al., 1989; Friedli
& Martinet, 1993). The inflowing gas driven by the inner bar may be a possible
mechanism to contribute the growth of the SMBH in the galactic center.
However, the fuel may be choked off as the inner bar gradually dissolves for a
short time ($\sim 1\;{\rm Gyr}$) when the SMBH becomes massive enough ($\sim
0.1\%$ stellar mass of the host galaxy, Du et al., 2017; Nakatsuno & Baba,
2023), and the growth is therefore self-limiting. The remnant of the destroyed
inner bar leads to the formation of a spheroidal component that is similar to
a small classical bulge (Guo et al., 2020). Interestingly, observations have
found local S2Bs are likely to host classical instead of pseudo bulges (de
Lorenzo-Cáceres et al., 2019a), although it would be difficult to determine
the origin of these classical bulges.

One important assumption in the above scenario is that inner bars promote gas
inflows from the nuclear ring to the center. This has long been speculated to
be the case (e.g. Shlosman et al., 1989; Heller et al., 2007; Namekata et al.,
2009; Hopkins & Quataert, 2010), but studies have provided conflicting results
(e.g. Maciejewski et al., 2002; Rautiainen et al., 2002). For example, the
‘bars-within-bars’ model of Shlosman et al. (1989) and Hopkins & Quataert
(2010) supports the idea that inner bars (or small-scale bars) can reduce gas
angular momentum just like large-scale bars; however, some theoretical studies
of gas motions in S2B potentials have argued that the gas flow patterns differ
fundamentally between single bars and nested bars, and the inner bar cannot
enhance the central mass inflow rate due to the lack of shocks (Maciejewski et
al., 2002; Rautiainen et al., 2002; Shlosman & Heller, 2002).

Nevertheless, there are still large uncertainties/simplification in the
modelling of nested bars that can affect the gas behaviours. For instance the
two bars in the S2B simulations of Maciejewski et al. (2002) and Namekata et
al. (2009) are modeled as rigid (Ferrers) ellipses whose torque distributions
are different from those of real bars (Buta & Block, 2001; Li et al., 2017).
In addition, both the pattern speeds and the shapes of the bars are pulsating
in a self-consistent, dynamically decoupled double-barred system (e.g. Shen &
Debattista, 2009; Du et al., 2015; Wu et al., 2016), and this effect on gas
flows has been less explored in the literature. Moreover, star formation and
feedback in the central kiloparsec may also affect the inner bar dynamics,
including the gas inflow (Wozniak, 2015), and this requires detailed modeling
of sub-grid physics as well as sufficient numerical resolution. It remains an
important question, therefore, whether inner bars are able to drive gas
inflows and potentially activate the central SMBH, which motivates us to
perform an in-depth numerical study of gas flows in S2Bs with state-of-the-art
subgrid physics implemented. We would like to investigate: (1) whether the
inner bar can drive gas further to the center using a realistic S2B model that
qualitatively matches the observed stellar kinematics; (2) how the gas flows
respond to a pulsating inner bar; and (3) how different the gas flow pattern
is between single- and double-barred galaxies, especially in the central
region.

The paper is organized as follows: Section 2 describes our simulation setups.
Section 3 presents the evolution of our S2B models using different gas
prescriptions and isolates the effects of the inner bar with control models.
Section 4 compares the properties of our simulated S2Bs with observations. We
discuss the implications and limits of this study in Section 5. Section 6
gives a summary of our results.

##  2 Simulations

The fiducial S2B model of Du et al. (2017) nicely reproduces many observed
features in real S2B galaxies, such as $\sigma$-humps/hollows and $h_{4}$
rings(see Du et al., 2016; de Lorenzo-Cáceres et al., 2019b). The long-lived
inner bar in this model forms spontaneously from a cool inner disk via
dynamical instabilities. For our study we choose to adopt a similar setup, as
detailed below. We select the snapshot corresponding to the time $t=2.7\;{\rm
Gyr}$ of the fiducial model in Du et al. (2017) as the initial condition for
the stars in our simulations. At this time, the semi-major axis of the inner
bar is $\sim 0.75\;{\rm kpc}$, while the outer bar has a semi-major axis of
$\sim 7.5\;{\rm kpc}$. The inner bar rotates roughly three times as fast as
its outer counterpart. The two bars have reached a relatively steady state,
but the pattern speeds of the two bars are not stationary and the shape of the
bars also changes, i.e. the bars are pulsating during the evolution (see also
Debattista & Shen, 2007). The stellar disk in our model has a mass of
$6.0\times 10^{10}\;{\rm M}_{\odot}$. We use the same logarithmic rigid dark
matter halo that Du et al. (2015, 2017) used, which has a potential given by
$\Phi(r)=0.5V_{h}^{2}\ln{(r^{2}+r_{h}^{2})}$ with $V_{h}=157.4\;{\rm km}\;{\rm
s}^{-1}$ and $r_{h}=37.5\;{\rm kpc}$, where $r$ is the spherical radius with
respect to the galactic center. This rigid halo simplifies our simulation,
while the dynamics at the center where the baryonic components dominate, are
barely changed. We refer the reader to Du et al. (2015, 2017) for further
details of the model.

Figure 1: Initial gas surface density profiles used in our models. The solid
curve shows the ‘HOLE’ profile, with the two arrows denoting the inner and
outer boundary at $R_{1}=3.75\;{\rm kpc}$ and $R_{2}=10.0\;{\rm kpc}$; the
dashed line shows the ‘EXP’ profile with the central density
$\Sigma_{0}=25\;{\rm M}_{\odot}\;{\rm pc^{-2}}$ and the scale radius
$R_{d}=5\;{\rm kpc}$.

We insert a lightweight gas disk into this snapshot to study gas flows under
the live S2B potential. We start the simulation time from when we insert the
gas disk. There are two different gas surface density profile at $t=0$ used in
our models as shown in Figure 1.: one is set to an annulus with inner and
outer radii at $R_{1}=3.75\;{\rm kpc}$ and $R_{2}=10.0\;{\rm kpc}$,
respectively, where $R$ is the cylindrical radius with respect to the galactic
center (labeled as ‘HOLE’ profile). $R_{1}$ is chosen to be close to the
middle of the outer bar. This annulus configuration aims to protect the inner
bar from the rapid formation of a gaseous central massive concentration (CMC)
that may dissolve the inner bar quickly (e.g. Wozniak, 2015; Du et al., 2017;
Guo et al., 2020; Nakatsuno & Baba, 2023). The initial inner boundary $R_{1}$
also help us focusing on the inflowed gas during the evolution of the model.
The surface density of the annulus gas disk between $R_{1}$ and $R_{2}$ is
uniform111A uniform gas disk is commonly used in studying bar-driven gas flows
because it results in no pressure gradient in the radial direction (e.g.
Athanassoula, 1992; Englmaier & Gerhard, 1997; Kim et al., 2012; Sormani et
al., 2015). with $\Sigma_{\rm gas}=5\;{\rm M}_{\odot}/{\rm pc}^{2}$ and is
tapered off by a cubic spline function with a characteristic width of
$0.2\;{\rm kpc}$ outside $R_{1}$ and $R_{2}$ (see also a similar initial gas
disk setup in Beane et al., 2023). The total gas mass is therefore $1.5\times
10^{9}\;{\rm M}_{\odot}$; the other is set to have an exponential surface
density profile (labeled as ‘EXP’ profile). The scale radius of the
exponential gas disk is $R_{d}=5\;{\rm kpc}$ and the central gas surface
density is $\Sigma_{0}=25\;{\rm M}_{\odot}\;{\rm pc^{-2}}$. The total gas mass
is therefore $3.9\times 10^{9}\;{\rm M}_{\odot}$. Although this gas disk
changes both the gas surface density profile and the gas mass at the same
time, it primarily tests the effects of a massive CMC on bars in a relatively
extreme case. The vertical density profile is chosen to be in hydrostatic
equilibrium following Springel et al. (2005) and Wang et al. (2010) for both
gas disk types. We start the gas disk on circular orbits obtained from the
axisymmetrized galactic potential. The initial temperature of the gas disk is
set to $1.5\times 10^{4}{\;{\rm K}}$. All gas cells at $t=0$ have solar
metallicity and alpha-elements as described in Asplund et al. (2009).

Table 1: Model Setup

Model | $\gamma$ | Gas Self-Gravity | SMUGGLE | Gas Profile | Gas Fraction | Galaxy Type
---|---|---|---|---|---|---
S2BISO | $1.0$ | Off | Off | HOLE | $2.5\%$ | S2B
S2BSMU | $5/3$ | On | On | HOLE | $2.5\%$ | S2B
S1BSMU | $5/3$ | On | On | HOLE | $2.5\%$ | SB
S0BSMU | $5/3$ | On | On | HOLE | $2.5\%$ | SA
S2BEXP | $5/3$ | On | On | EXP | $6.5\%$ | S2B

Note. — List of models in this work. Column 1: model names; Column 2:
adiabatic index $\gamma$, $\gamma=1$ corresponds to the isothermal EoS while
$\gamma=5/3$ is for ideal monoatomic gas; Column 3: switch for gas self-
gravity; Column 4: switch for SMUGGLE model; Column 5: initial surface density
profile of the gas disk, HOLE means an initial uniform gas disk with an inner
boundary of $R_{1}=3.75\;{\rm kpc}$ and outer boundary of $R_{2}=10.0\;{\rm
kpc}$, while EXP means an initial exponential gas disk without the inner
boundary (see also Figure 1); Column 6: gas mass fraction relative to the
stellar disk. Column 7: galaxy morphology of the model.

Gas cooling and heating, star formation, and stellar feedback are included via
the ISM and stellar feedback model model SMUGGLE (Stars and MUltiphase Gas in
GaLaxieEs model, Marinacci et al., 2019). This model reproduces well the
observed feedback-regulated star formation and has been widely used for
studying galaxy formation and evolution (e.g. Burger et al., 2022; Li et al.,
2022; Sivasankaran et al., 2022; Tacchella et al., 2022). We adopt a local
star formation efficiency (SFE) $\epsilon_{\rm sf}=0.01$ and a star formation
density threshold $\rho_{\rm th}=100\;{\rm cm}^{-3}$ in SMUGGLE. We use the
same stellar yields and evolution model as the original Illustris simulation
(Vogelsberger et al., 2013).

Figure 2: Evolution of model S2BISO. Top: stellar surface density. Bottom: gas
surface density. White lines represent the contours of the stellar surface
density placed at $0.16,0.54,1.85$, and $6.31\times 10^{3}\;{\rm
M}_{\odot}\;{\rm pc^{-2}}$. The two stellar bars are clearly outlined by the
contours. The first two columns show the early evolution when the two bars are
misaligned, while they are parallel in the third column and perpendicular in
the fourth column. The evolution time is labeled at the top. The disk rotates
counterclockwise.  Figure 3: Similar to Figure 2 but for S2BSMU.

The simulations are performed using the moving-mesh code AREPO (Springel,
2010; Weinberger et al., 2020). The simulation box size is $60\;{\rm kpc}$ so
the boundary effects on the gas disk are negligible. There are 3.8 million
stellar particles with a fixed softening length of $25\;{\rm pc}$ and 1.3
million initial gas cells with adaptive softening enabled. The minimum
adaptive softening is set to $0.25\;{\rm pc}$ and the target mass of gas cells
is $\approx 1100\;{\rm M}_{\odot}$. The galaxy properties modeled with SMUGGLE
nicely converge at this resolution (Marinacci et al., 2019). The minimum gas
cell radius in our simulations reaches $\sim 0.14\;{\rm pc}$.

We present 5 models with configurations summarized in Table 1. The fiducial
S2B model is referred to as S2BSMU, with SMUGGLE enabled. S2BISO has the same
initial conditions as S2BSMU but uses the isothermal equation of state (EoS)
and excludes gas self-gravity to avoid gravothermal catastrophe. The effective
sound speed in S2BISO is $c_{s}\approx 10\;{\rm km}\;{\rm s}^{-1}$ (obtained
from initial gas temperature). This sound speed reflects the average velocity
dispersion between molecular clouds rather than the microscopic temperature of
the diffuse gas. We further run two additional control models (S1BSMU and
S0BSMU) to highlight the gravitational effects of bars: We generate S1BSMU by
azimuthally scrambling222Scrambling means the in-plane position $(x,y)$ and
velocity $(v_{x},v_{y})$ of each stellar particle is rotated by a random angle
$\phi$ with respect to the galactic center. $\phi$ follows a uniform
distribution on the interval $[0,2\pi)$ (see also Brown et al., 2013; Li et
al., 2014). the stellar particles at $R\leq 0.75\;{\rm kpc}$ of S2BSMU in the
initial condition. S1BSMU is thus a single-barred (SB) galaxy whose bar
structure is almost the same as the outer bar of S2BSMU. Similarly, we
generate S0BSMU by scrambling stellar particles at $R\leq 7.5\;{\rm kpc}$,
thus producing an unbarred galaxy. These 4 models have the exactly same
azimuthally-averaged mass profile and rotation curve at $t=0$. The gas mass
fraction relative to the stellar disk in these 4 models is $2.5\%$; these 4
models are run for $\sim 1.9\;{\rm Gyr}$. S2BEXP has nearly identical initial
conditions as S2BSMU, with the exception of an initial exponential gas disk
rather than a uniform gas disk. Consequently, the gas mass fraction is
increased to $6.5\%$; this model is run for $\sim 3.8\;{\rm Gyr}$.

Figure 4: Comparison of the inner gas surface density between S2BISO (upper)
and S2BSMU (lower). The four columns represent 4 different epochs when the two
bars have different angles with respect to each other (indicated by the red
and blue sticks at the bottom-right corners of the upper row). The snapshots
are taken at $t=1.271,1.287,1.310,1.334\;{\rm Gyr}$ for S2BISO and
$t=1.256,1.272,1.295,1.317\;{\rm Gyr}$ for S2BSMU. The dashed and dotted-
dashed circles denote two characteristic cylindrical radii of $R_{\rm
outer}=3.75\;{\rm kpc}$ and $R_{\rm inner}=0.75\;{\rm kpc}$. The gas disk
rotates counterclockwise.  Figure 5: Bar properties and gas inflows in S2BISO
(left) and S2BSMU (right) models. First row: outer and inner bar amplitudes.
The Fourier components ${\rm A_{2}}$ and ${\rm A_{0}}$ are measured in the
radial range of $2.4-7.5\;{\rm kpc}$ for the outer bar and $0-0.75\;{\rm kpc}$
for the inner bar; Second row: enclosed mass of gas (S2BISO) and gas plus
newly formed stars (S2BSMU) inside $R_{\rm inner}$ and $R_{\rm outer}$; Third
row: mass inflow rate obtained from the time derivative of the lines in the
second row. Vertical grey solid and dashed lines indicate the moments when the
two bars are aligned and perpendicular, respectively. Last row: enclosed star
formation rate in S2BSMU. This quantity is calculated based on the total mass
of new stars inside the given radius that are formed in the past $10\;{\rm
Myr}$.

##  3 Results

###  3.1 Comparing S2B Model with Different Gas Assumptions

In this subsection, we show how stellar feedback affects the gas morphology as
well as the mass inflow rate in the same S2B potential by comparing the fully
star-forming model S2BSMU with the isothermal model S2BISO.

####  3.1.1 Gas Morphology in the Inner and Outer Region

We first present gas flow patterns in S2BISO with the isothermal EoS. This EoS
assumes the specific internal energy of the ISM is the result of a balance
between heating and cooling, and can be regarded as a first-order
approximation to the observed cold gas in galactic disks. Although it is a
simplification of the true multiphase nature of the ISM, previous studies have
shown that the isothermal EoS reproduces the observed gas dynamics in bars
reasonably well (e.g. Weiner et al., 2001; Lin et al., 2013; Fragkoudi et al.,
2017). We therefore start with this simple model to highlight the
gravitational effects of bars. It also helps us to better compare our model
with previous S2B simulations that use the isothermal EoS (e.g. Maciejewski et
al., 2002; Shlosman & Heller, 2002; Namekata et al., 2009).

Figure 2 plots the surface density of the stellar (top) and the gas (bottom)
disk in S2BISO at 4 different times. The models are rotated such that the
outer bar is always aligned with the $x$-axis for easier comparison. The two
stellar bars are clearly outlined by the white contours of stellar surface
density. Note the inner bar has different orientations with respect to the
outer bar due to its different (higher) pattern speed. The first two columns
show the earlier evolution of the model. The initial circular and uniform gas
disk is gradually distorted by the bars, forming a pair of strong shocks at
the leading side of the outer bar together with a two-arm spiral around
$10\;{\rm kpc}$. Gas flows inward to the inner boundary soon after the
simulation starts (first column), then accumulates at a high-density nuclear
ring with a size of $\sim 1.5\;{\rm kpc}$ (second column). The inner bar lies
within the gaseous nuclear ring, as seen also in observations (Buta et al.,
2015). Overall, the large-scale flow pattern at $R\gtrsim 4\;{\rm kpc}$
resembles the results of previous single-barred simulations (e.g. Kim et al.,
2012; Li et al., 2015; Sormani et al., 2015). The morphology and kinematics of
gas in this region is mainly dominated by the outer bar potential, and is not
sensitive to the inner bar orientation, as shown in the third (parallel bars)
and fourth (perpendicular bars) columns.

Figure 6: Stellar and gas surface density comparison using S2BSMU, S1BSMU, and
S0BSMU at $t=1.8\;{\rm Gyr}$, from left to right as labelled. The top row is
the stellar surface density and the bottom row is the gas surface density. The
outer bar is aligned with the $x$-axis for S2BSMU and S1BSMU. The dashed
circles in the bottom row denote the two characteristic cylindrical radii
$R_{\rm outer}=3.75\;{\rm kpc}$ and $R_{\rm inner}=0.75\;{\rm kpc}$.

We next present S2BSMU, which resolves the multiphase gas with explicit star
formation and stellar feedback included. When comparing the surface densities
in Figure 3 with those in Figure 2, the stellar component is almost the same.
The smooth gas features (e.g. the shocks, spirals, and the nuclear ring) in
S2BISO become more flocculent and clumpy in S2BSMU due to local stellar
feedback, but their overall shapes remain mostly unchanged. This is expected
since the gas mass fraction is only $2.5\%$ in our models. The global star
formation rate (SFR) is around $0.1-0.2{\;{\rm M}_{\odot}\;{\rm yr}^{-1}}$
during the entire simulation, so the gas evolution on the large scale is still
governed by the non-axisymmetric stellar potential.

The major differences between S2BISO and S2BSMU emerge in the central
$R\lesssim 2\;{\rm kpc}$, which is shown in Figure 4. In this figure, we
present gas surface densities in 4 snapshots with different bar orientations
indicated by the red (outer bar) and blue (inner bar) sticks at the lower
right of the first row. The shape of the nuclear ring in S2BISO clearly
depends on the relative orientations between the two bars: it tends to be more
elliptical when the bars turn from parallel to perpendicular (second and third
columns), and vice versa. Furthermore, the ring slightly leads the inner bar
during the rotation. The cyclically changing shape of the nuclear ring is
consistent with earlier S2B simulations and orbit analysis (Maciejewski &
Sparke, 2000; Laine et al., 2002). We also see straight shocks (which would be
traced by dust lanes) at the edges of the inner bar (more obvious in the
second and third columns) that can drive gas to the center. This flow pattern
is different from those found by Maciejewski et al. (2002) and Shlosman &
Heller (2002), and we discuss the possible reasons in §5.1. On the other hand,
the ring in S2BSMU roughly maintains the same shape regardless of the bar
relative orientations. This can be understood in terms of a simple estimation
of the value of the gas kinetic energy and the energy injected by stellar
feedback: the gas in the ring has a typical rotation velocity of $\sim
200\;{\rm km}\;{\rm s}^{-1}$ and the gas mass of the ring is $\sim 2.0\times
10^{8}\;{\rm M}_{\odot}$, which gives a kinetic energy of $8.0\times
10^{55}\;{\rm erg}$; the energy released by the supernova (SN) feedback in the
ring during one ring rotation ($\sim 60\;{\rm Myr}$) is $\approx 6-9\times
10^{55}\;{\rm erg}$. The two comparable energies suggest feedback becomes as
important as gravity in about one dynamical timescale, and the effect of the
inner bar on the ring shape is therefore less prominent. We also find that a
small gas bar is formed and co-rotates with the inner bar in S2BISO, but this
structure seems to be erased by stellar feedback in S2BSMU.

Figure 7: A zoom-in view of Figure 6 to the central $2.5\;{\rm kpc}$. The
first and second columns are the stellar and gas surface densities. The third
and fourth column plot the metallicity (represented by [Fe/H]) and $\alpha$
abundance (represented by [Mg/Fe]) of the newly formed stars.

####  3.1.2 Bar Amplitudes and Mass Inflow Rate

A more quantitative analysis of the evolution of the two S2B models is
presented in Figure 5. The two columns display the properties of the
isothermal (S2BISO) and the multiphase (S2BSMU) models. We first compare the
bar amplitude measured by the normalized Fourier component ${\rm A_{2}/A_{0}}$
in the top row. The bars in S2BISO are relatively stable, although the bar
amplitudes pulsate as the inner bar rotates with respect to the outer one,
which is also seen in previous simulations (Debattista & Shen, 2007; Wozniak,
2015; Du et al., 2015; Wu et al., 2016). On the other hand, the bar amplitudes
exhibit a weakening trend in S2BSMU. This is because gas self-gravity is
included in this model, and the accumulating gas forms a CMC that can weaken
the bars (Athanassoula, 2003; Shen & Sellwood, 2004; Du et al., 2017). In the
second row, we plot the enclosed gas mass at two characteristic radii ($R_{\rm
outer}=R_{1}=3.75\;{\rm kpc}$ and $R_{\rm inner}=0.75\;{\rm kpc}$), which can
be seen by the white dashed and dotted-dashed circles in the first column of
Figure 4. $R_{\rm outer}$ is equal to $R_{1}$, which corresponds to the inner
boundary of the initial gas disk (see §2), and most of the gas inflow driven
by the outer bar accumulates inside this radius; $R_{\rm inner}$ is
approximately the size of the inner bar. We see that the enclosed mass quickly
increases with time: it reaches $\sim 3.0\times 10^{8}\;{\rm M}_{\odot}$
inside $R_{\rm outer}$ within $1\;{\rm Gyr}$ in both models. This value is
about $20\%$ of our initial gas disk mass and is about $0.5\%$ of the stellar
disk mass. Note that the enclosed mass includes both gas and newly formed
stars in S2BSMU. It is also worth noting that when the inner bar in S2BSMU
starts to weaken at $t\gtrsim 1.3\;{\rm Gyr}$, the mass enclosed within
$R_{\rm inner}$ reaches $\sim 0.1\%$ of the disk mass. This is consistent with
the inner bar dissolution criterion predicted in Du et al. (2017). The outer
bar in this model is also weakened over time.

About $80\%$ of the gas in the nuclear ring is accumulated inside $R_{\rm
inner}$ in S2BISO, but the fraction drops to $\sim 30\%$ in S2BSMU. There are
probably two reasons for this: (1) Rapid star formation consumes gas in the
nuclear ring at $R\sim 2\;{\rm kpc}$ of S2BSMU, which makes less gas available
for the inner bar to drive inwards to $R_{\rm inner}$. The gravitational
dynamical timescale (i.e. free-fall timescale) for a gas cell to form stars in
the ring with a typical gas density of $0.5\;{\rm M}_{\odot}\;{\rm pc^{-3}}$
is around $10\;{\rm Myr}$, while the gas rotation period in the ring is $\sim
60\;{\rm Myr}$, and it usually takes a few rotation periods for the inner bar
to remove gas angular momentum. Since the dynamical timescale is longer than
the local star formation timescale, a large fraction of gas will form stars in
the ring instead of being driven inwards by the inner bar. By contrast, there
is no mechanism to consume gas in S2BISO, and the inner bar can therefore
drive most gas in the ring down to the center. (2) The ring in S2BISO is
smaller in size, and periodically changes its shape (Figure 4). This helps gas
in the ring to be more easily affected by the inner bar torque in S2BISO
compared to S2BSMU.

Figure 8: Comparison of the enclosed mass (top row), mass inflow rate (second
row), molecular gas fraction (third row), and SFR (bottom row) at two
characteristic radii ($R_{\rm outer}=3.75\;{\rm kpc}$ and $R_{\rm
inner}=0.75\;{\rm kpc}$) of models S2BSMU, S1BSMU, and S0BSMU. The molecular
gas fraction is defined as the ratio of molecular gas mass and total gas mass.
Note the range of the $y$-axis is not the same in the left and right columns.

The third row of Figure 5 shows the mass inflow rate at the two radii,
obtained by taking the time derivative of the lines in the second row. We
focus on $0.6-1.6\;{\rm Gyr}$ when the inner bars are relatively stable in
both models. The mass inflow rates oscillate with time and seem to be quasi-
periodic. We label the moments when the inner bar and outer bar are parallel
(solid vertical lines) and perpendicular (dashed vertical lines) to each
other. The inflow rate tends to peak when the bars are perpendicular,
consistent with more clearly shock features in the top middle two columns of
Figure 4. The oscillation period ($\sim 78\;{\rm Myr}$) is therefore
determined by the pattern speed difference (or the beat frequency) of the two
bars. This bar-modulated periodic gas inflow is also reported in some of the
previous S2B simulations (Shlosman & Heller, 2002; Namekata et al., 2009), but
here we show for the first time that the inflow is strongest when the two bars
are perpendicular, possibly due to more clear shocks inside the inner bar at
this angle. Interestingly, this trend remains in S2BSMU, although the spatial
distribution of gas seems to be independent with respect to the bar
orientations due to local stellar feedback (see the bottom row in Figure 4).
It is worth noting that there appears to be no regular pattern in the enclosed
SFR curve of S2BSMU within $R_{\rm inner}$, despite the fact that the model
has a nearly periodic gas supply in this region (but also see the results in
Moon et al., 2022).

Figure 5 demonstrates that the inner bar does promote gas inflow to the
center. However, the mass inflow rate at $R_{\rm inner}$ in S2BSMU is
considerably lower than that in S2BISO. It is therefore natural to ask whether
the amplitude of the central mass inflow rate is still dominated by the inner
bar potential when stellar feedback is included.

Figure 9: Solid lines: total SFR as a function of time in S2BSMU, S1BSMU, and
S0BSMU. Dashed lines: SFR at $R\geq R_{\rm outer}=3.75\;{\rm kpc}$ as a
function of time in the models.

###  3.2 Comparing S2B Model with SB and SA Models

In this subsection we compare three models with different bar numbers but all
with the SMUGGLE model (S2BSMU, S1BSMU, S0BSMU) to isolate the effects of the
inner bar on gas dynamics from (stochastic) stellar feedback.

####  3.2.1 Mass Inflow Rate and Star Formation Rate

Figure 10: KS relation in models S2BSMU, S1BSMU, and S0BSMU. Solid black
curves are the contours of the scatter points. Only gas and newly formed stars
in the disk region ($\left\lvert z\right\rvert\leq 2\;{\rm kpc}$) are
considered in this plot. For each model, we include 11 snapshots separated by
$\sim 100\;{\rm Myr}$ to increase the sample size. The black dashed line is
the global fit from Querejeta et al. (2021) with a slope of $0.97$. Two black
dotted lines representing the depletion time of $0.1\;{\rm Gyr}$ and $10\;{\rm
Gyr}$ are included for comparison. The subplot in the middle panel illustrates
the placement of the hexagonal bins in our models. The bins are colored
according to their galactocentric radius, and the underlying plot is the same
gas surface density plot of S2BSMU shown in Figure 7.

Figure 6 shows the stellar and gas surface densities of S2BSMU, S1BSMU, and
S0BSMU at $t=1.8\;{\rm Gyr}$. These three models have the same azimuthally-
averaged mass distribution at $t=0$. Gas flow patterns in S2BSMU and S1BSMU
are very similar outside $R\sim 1\;{\rm kpc}$; they both form a pair of spiral
arms and a bright nuclear ring produced by the outer bar. S0BSMU only shows
flocculent spirals and little gas can flow into $R_{\rm outer}$ as no bar is
present333The faint gas blob at the center of S0BSMU is probably due to
numerical diffusion.. The gas distribution of this model looks similar to NGC
2775 (e.g. Leroy et al., 2021).

Inside $R\sim 1\;{\rm kpc}$, the gas behaves differently in S2BSMU and S1BSMU
as shown in Figure 7. The left two columns are zoom-in views of the stellar
and gas surface densities in Figure 6. Inflow signatures can be spotted in
S2BSMU as multiple gas streams formed inside the nuclear ring down to the
galactic center. On the other hand, the nuclear ring has a relatively well-
defined inner boundary in S1BSMU, implying the inflow has stalled at this
radius. This is not surprising since the nuclear ring in single-barred
galaxies serves as a barrier that prevents gas from reaching the nucleus, even
if stellar feedback is included (Fanali et al., 2015; Li et al., 2017; Tress
et al., 2020). The right two columns illustrate the spatial distribution of
the newly formed stars color-coded by their chemistries ([Fe/H] and [Mg/Fe]).
Although we do not observe an inner gas bar in S2BSMU, newly formed stars
inside the ring follow the shape and rotation of the stellar inner bar,
similar to the cascade bar formation scenario proposed in previous studies
(e.g. Shlosman et al., 1989, 1990; Hopkins & Quataert, 2010). This bar-shaped,
newly formed stellar structure is more metal rich and less $\alpha$-enhanced
compared to the surroundings, consistent with the observations of the inner
bar in NGC 5850 (de Lorenzo-Cáceres et al., 2019b). The newly formed stars in
S1BSMU show a similar pattern, but not as strongly as in S2BSMU.

In Figure 8 we plot the enclosed mass, mass inflow rate, and SFR for these
three models at $R_{\rm outer}$ (left column) and $R_{\rm inner}$ (right
column). The accumulated masses at $R_{\rm outer}$ are almost identical for
S2BSMU and S1BSMU, but the two models diverge inside $R_{\rm inner}$. The
inner bar in S2BSMU drives $\sim 3$ times more mass inside $R_{\rm inner}$
compared with the single-barred model. The mass inflow rate of S2BSMU at
$R_{\rm inner}$ is quasi-periodic and is related to the inner bar dynamics, as
shown in Figure 5. By contrast, the mass inflow rate of S1BSMU at $R_{\rm
inner}$ is mostly dominated by stellar feedback, and the periodicity
disappears. We obtained a time-averaged mass inflow rate of $0.144{\;{\rm
M}_{\odot}\;{\rm yr}^{-1}}$ at $R_{\rm outer}$ and $0.013{\;{\rm
M}_{\odot}\;{\rm yr}^{-1}}$ at $R_{\rm inner}$ for S1BSMU, within the time
interval of $t=0.6-1.6\;{\rm Gyr}$; the ratio of these two inflow rates is
$\sim 10.8$, similar to the results in Tress et al. (2020). S2BSMU increases
the central inflow rate to $0.044{\;{\rm M}_{\odot}\;{\rm yr}^{-1}}$, reducing
the ratio to $\sim 3.4$. The control model S0BSMU has weak inflows (probably
caused by diffusion and numerical effects), which are negligible compared with
the other two models. We therefore conclude that the inner bar in S2Bs indeed
helps gas flowing inwards.

Figure 11: Evolution of ${\rm A_{2}/A_{0}}$ (top) and pattern speed (bottom)
in the inner bar region ($R\leq 0.75\;{\rm kpc}$) of two models with different
gas fraction. At $t\gtrsim 1.7\;{\rm Gyr}$ in S2BEXP, the pattern speed
calculated using the ${\rm A_{2}}$ phase angle time derivatives may be less
accurate due to the small ${\rm A_{2}}$ amplitude.  Figure 12: Kinematic maps
of S2BEXP, viewed face-on. From left to right are 4 different epochs. The
second ($\sigma$) and fourth ($h_{4}$) Gauss–Hermite moments of the vertical
velocity $v_{z}$ distribution are shown in the top and bottom rows,
respectively. White lines represent the contours of the stellar surface
density. The outer bar is aligned with the $x$-axis.

In addition, we find S2BSMU has the highest SFR in both regions as shown in
the last row of Figure 8. This is expected as a consequence of more
accumulated gas in the center. Observations have shown that bars can enhance
central star formation (e.g. Lin et al., 2017, 2020), and our results imply
multiple bars may further amplify this effect. In fact, the total SFR in
S2BSMU is larger than the other two models during most of the time in our
simulations (solid lines in Figure 9). This indicates that bars may not
significantly reduce the SFR within a short period in our gas-poor models (but
see also Khoperskov et al., 2018). The SFR in S2BSMU and S1BSMU at $R>R_{\rm
outer}$ is slightly lower than that in S0BSMU (dashed lines in Figure 9). This
is probably due to a smaller gas surface density in the outer disk region, as
a significant fraction of gas has been driven inside $R_{\rm outer}$.

####  3.2.2 Kennicutt–Schmidt Relation

We investigate also the molecular Kennicutt–Schmidt (KS) relation (Schmidt,
1959; Kennicutt, 1998) that connects the molecular gas surface density
($\Sigma_{\rm mol}$) to the star formation rate surface density ($\Sigma_{\rm
SFR}$) in our models. Motivated by Querejeta et al. (2021) who used PHANGS-
ALMA data (Leroy et al., 2021) to explore the KS relation in different
galactic environments, we dissect our simulations into $\sim 1.3\;{\rm
kpc}$-wide hexagonal bins444See the subplot in the middle panel of Figure 10,
which is similar to Fig. 3 in Querejeta et al. (2021).. The gas and SFR
surface densities are then estimated inside each bin and are plotted as
colored dots in Figure 10. We divided the bins into three groups based on
their radius: inside the inner bar (orange); outside the inner bar but roughly
covering the nuclear ring (blue); and the disk region (grey). In general, the
KS relation in our model is similar to the global fit in Querejeta et al.
(2021) (black dashed line), which implies a weak environmental (or radial)
dependence of the star formation efficiency. Centers (orange points) in barred
models (S2BSMU and S1BSMU) have shorter depletion times than the rest of the
disk, also consistent with the findings in Querejeta et al. (2021). In
addition, the depletion time in the center of S2BSMU seems to be even smaller
than that in S1BSMU. This is possibly due to a higher molecular gas fraction
in the central regions of S2BSMU, as shown in the third row of Figure 8. On
the other hand, the low gas surface density as well as SFR in the center of
the unbarred model S0BSMU is simply due to our centrally-holed initial gas
disk setup, and thus cannot reflect the real central environments in unbarred
galaxies. We therefore conclude that SMUGGLE model can reproduce the latest
observations reasonably well. In SMUGGLE model we adopt a constant star
formation efficiency per free-fall time ($\epsilon_{\rm sf}=0.01$) and compute
the SFR based on the molecular gas only, the KS relation is then captured by
the self-regulated star formation and feedback processes. However, there are
hints that $\epsilon_{\rm sf}$ may vary systematically with $\Sigma_{\rm mol}$
and $\Sigma_{\rm SFR}$ (Sun et al., 2023). We leave a detailed comparison of
the star formation laws between simulations and observations for future
studies.

Figure 13: Comparison between an observed S2B galaxy NGC 1433 and our S2BSMU
model. We adopt a distance of $18.63\;{\rm Mpc}$, a systemic velocity of
$1057.4\;{\rm km}\;{\rm s}^{-1}$, a disk position angle of $199.7\degree$, and
an inclination angle of $28.6\degree$ for this galaxy (Buta et al., 2001; Lang
et al., 2020). Top row: from left to right are S4G $3.6\;\mu m$ image (unit:
${\rm L_{\odot}/pc^{2}}$), PHANGS-ALMA CO(2-1) integrated intensity map
(moment 0, unit: ${\rm K\;{\rm km}\;{\rm s}^{-1}}$), intensity weighted mean
velocity map (moment 1, unit: ${\rm km\;s^{-1}}$), and residual velocity map
obtained by subtracting circular motions (unit: ${\rm km\;s^{-1}}$). Note the
spatial scales are different for S4G image and ALMA maps. Second row: from
left to right are projected stellar surface density map (unit: ${\rm
M_{\odot}/pc^{2}}$), projected molecular gas surface density map (unit: ${\rm
M_{\odot}/pc^{2}}$), molecular gas line-of-sight velocity map (unit: ${\rm
km\;s^{-1}}$), and molecular gas residual velocity map (unit: ${\rm
km\;s^{-1}}$) for the model at $t=818.2\;{\rm Myr}$. Bottom row: same as the
second row but at $t=892.7\;{\rm Myr}$. The green plus sign marks the center
of the image $(x=0,y=0)$.

###  3.3 Long Term Evolution

Figure 5 has shown that the inner bar in S2BSMU gradually decays due to
accumulated central mass. We therefore use S2BEXP to test how long the inner
bar will survive if there is an even higher gas mass at the center from the
beginning. The gas mass inside $R_{\rm inner}$ is $\sim 0.1\%$ stellar mass at
$t=0$ for this model, which meets the inner bar dissolution criterion in Du et
al. (2017). The evolution of the Fourier component $A_{2}/A_{0}$ and inner bar
pattern speed is shown in Figure 11. We find that the independent pattern
speed inside $R_{\rm inner}$ disappears at $t\sim 2\;{\rm Gyr}$. At this
moment, the contours of the stellar surface density in the central $1\;{\rm
kpc}$ are nearly round (see the third column of Figure 12). We conclude the
inner bar fully dissolves at $t\sim 2\;{\rm Gyr}$ in S2BEXP. The second and
third columns of Figure 12 are two snapshots before and after the dissolution
of the inner bar. We note the outer bar also weakens considerably in this
model.

We plot the evolution of 2D stellar kinematic maps of S2BEXP in Figure 12. The
top and bottom row shows the second ($\sigma$) and fourth ($h_{4}$)
Gauss–Hermite moments (Gerhard, 1993; van der Marel & Franx, 1993) of the
vertical velocity ($v_{z}$) distribution. We include both the old stars and
the newly formed stars in this figure. The $\sigma$-humps (de Lorenzo-Cáceres
et al., 2008; Du et al., 2016) along the minor axis of the inner bar is
clearly seen in the first column, but this feature disappears when the inner
bar has fully been destroyed. The central region then remains roughly
axisymmetric till the end of the simulation (see also Fig. 8 in Guo et al.,
2020). The $h_{4}$ maps reveal a ring-like structure and a central blob that
peak better after $t\sim 1.7\;{\rm Gyr}$, which are possibly associated with
the newly formed stars in the simulation. Even though the star formation is
high in the inner bar region, the newly formed stars mainly follow chaotic
orbits scattered by the CMC, and can hardly support the inner bar structure as
suggested by Nakatsuno & Baba (2023). Strong stellar feedback of young stars
may also partially contribute to the dissolution of the inner bar, as it
prevents star forming clouds settling on the inner bar supporting orbits (see
also §3.1.1 and Figure 13). These kinematic signatures could be potentially
useful for detecting dissolved inner bars in future integral field unit (IFU)
observations.

##  4 Comparison with Observations

Inflow signatures related to the inner bar have been reported in observations
(e.g. Fathi et al., 2006; Schinnerer et al., 2006, 2007; González-Alfonso et
al., 2021). Based on the recent PHANGS-ALMA data, we find 6 S2B candidates
(NGC 1068, NGC 1087, NGC 1317, NGC 1433, NGC 4304, and NGC 4321) according to
the classifications in Erwin (2004) and Buta et al. (2015). NGC 1433 probably
has the most similar bar shapes (i.e. stellar surface density contours)
compared to our S2BSMU model. Although our model S2BSMU is not designed to
match any particular galaxy, we still compare it with NGC 1433 in Figure 13.
The top row shows the S4G $3.6\;\mu m$ image, CO(2-1) integrated intensity
map, CO intensity-weighted mean line-of-slight (LOS) velocity map, and the
residual velocity map obtained by subtracting the circular motion (i.e. the
spider diagram) from the $V_{\rm los}$ map of NGC 1433. The first two CO maps
are the same as in the PHANGS-ALMA survey paper (moment 0 and moment 1 maps in
Leroy et al., 2021). The circular velocity curve used to derive the residual
velocity map is from Lang et al. (2020) (see their Eq. 10 and Table 4, with
$V_{0}=204.5\;{\rm km}\;{\rm s}^{-1}$ and $R_{t}=0.59\;{\rm kpc}$). The middle
and bottom rows show the projected stellar surface density map, projected
molecular gas surface density map, molecular gas LOS velocity map, and
residual velocity map in S2BSMU at $t=818\;{\rm Myr}$ and $892\;{\rm Myr}$.
These plots are obtained by projecting the stellar and molecular gas disk in
the simulation with an inclination of $28.6^{\circ}$ and a disk position angle
of $199.7^{\circ}$ (Buta et al., 2001; Lang et al., 2020). The circular
velocity curve is derived by DiskFit (Sellwood & Spekkens 2015) that applies
the tilted ring method555The tilted ring method may fail to extract the true
circular velocity curve when bar-induced non-circular motions dominate the LOS
velocity field. A more detailed analysis will be presented in Liu et al. (in
prep.). to the molecular gas $V_{\rm los}$ field ($V_{0}=270.0\;{\rm km}\;{\rm
s}^{-1}$ and $R_{t}=0.30\;{\rm kpc}$). We assume purely circular motion and do
not consider radial or bisymmetric flows during the fitting, similar to Lang
et al. (2020).

In Figure 13, the two snapshots of S2BSMU have similar bar orientations as in
the S4G image (first column). In the second column, multiple gas streams
extending from the nuclear ring to the central part are seen both in the
observation and in our model. This indicates a strong gas inflow that is also
evident in the other five observed S2Bs. The presence of clumpy gas features
around the nuclear ring suggests that these regions are undergoing intense
star formation. The $V_{\rm los}$ maps in the third column also match
qualitatively. Note here we show the velocity information in the whole region
of the simulation, which presents a larger view of gas flows. We argue that
reproducing exactly the observed CO $V_{\rm los}$ map of S2Bs would be quite
challenging, since the effects of (stochastic) stellar feedback on gas
morphology and kinematics are as important as the gravitational effect of the
inner bar (§3.1.1). We also demonstrated this in the second and third rows of
Figure 13: the stellar surface densities are almost identical at $t=818\;{\rm
Myr}$ and $892\;{\rm Myr}$, but the gas morphology and $V_{\rm los}$ maps have
a relatively large difference. In the fourth column we show the residual
velocity maps of the observation and simulation. The residual velocity map in
the simulation look more chaotic, possibly also due to strong stellar
feedback. We note the structures in the residual velocity map somewhat relies
on whether people can measure the inclination and position angle of disks
accurately (e.g. Kolcu et al., 2023), thus a large uncertainty may exist to
make a reasonable comparison.

Another way to compare our model with real S2Bs is to calculate the mass
inflow rate at different radii. We predict that the ratio of the mass inflow
rate from the outer bar to the nuclear ring and from the nuclear ring to the
center is around 3.4 for S2Bs, while this number increases to $\gtrsim 10$ for
single-barred galaxies (§3.2). Wu et al. (2021) obtained a mass inflow rate of
$12{\;{\rm M}_{\odot}\;{\rm yr}^{-1}}$ along the dust lanes of the outer bar
in the S2B galaxy NGC 3504 based on ALMA CO data. It would be interesting to
measure the radial motions of gas in this galaxy (as well as other S2Bs) at
different radii, and compare the results with single-barred galaxies (e.g.
similar to Haan et al., 2009). However, this requires a detailed analysis of
the observational data, which is the subject of future works.

##  5 Discussion

###  5.1 Comparison with previous S2B models

Our models suggest that inner bars are efficient at driving gas inflows, which
does not support the previous S2B simulations in Maciejewski et al. (2002) and
Rautiainen et al. (2002). We first note that the results in Maciejewski et al.
(2002) may suffer from a now-recognized bug in the code CHMOG (see Kim et al.,
2012). Indeed Namekata et al. (2009) adopted a similar analytical S2B
potential but a different nuclear bulge component compared with Maciejewski et
al. (2002) and found significant inflows caused by the inner bar. The authors
speculated that the shallow central potential in Maciejewski et al. (2002) may
be the reason for the lack of inflow. While the model of Rautiainen et al.
(2002) consisted of a live S2B model, its resolution may have been too low to
resolve the central gas behavior since it employed only 40000 sticky gas
particles within a $4.5\;{\rm kpc}$ radius gas disk. Thus, it is possible that
the absence of strong gas inflows in Maciejewski et al. (2002) and Rautiainen
et al. (2002) may mostly be attributed to numerical issues.

The inflow due to the inner bar is quasi-periodic in our S2B models. Shlosman
& Heller (2002) and Namekata et al. (2009) found a similar cyclic inflow in
their S2B simulations, although the pattern speeds and the shapes of the bars
in their models are constant along time. This may imply a pulsating inner bar
is not necessary to produce a periodic inflow. The periodicity is likely due
to the fact that the inner bar sweeps gas more efficiently (or shocks are more
easily formed) when the two bars tend to be perpendicular with each other,
i.e. a pattern speed differences for the two bars would be sufficient. We also
find the inner bar in our S2BEXP model does not revive after dissolution,
which appears to differ from Wozniak (2015) who also included star formation
and stellar feedback. This may imply that star formation is not the only key
parameter for generating a recurrent inner bar, or that the gas fraction is
still too low in our models. It could also be due to the different
implementations of the sub-grid physics in the two studies.

###  5.2 For how long does the inner bar promote gas inflows?

According to S2BEXP model, the inner bar is fully dissolved within $\sim
2\;{\rm Gyr}$ with a gas disk of $\sim 5\%$ stellar mass. The gas inflow then
stalls at the nuclear ring, similar to the case in S1BSMU. This timescale is
longer than that obtained in Du et al. (2017) ($\lesssim 1\;{\rm Gyr}$),
probably because the central potential of the CMC is more ‘soft’ than a nearly
Keplerian one used in Du et al. (2017).

One important aspect that is not considered in our study is the feedback due
to the active galactic nucleus (AGN). The huge energy released by the central
engine may help to clear the accumulated gas at the center. The inner bar is
then expected to be more robust against gas inflows than what we have found
here. It is therefore reasonable to suggest that the inner bar can promote
central gas inflows on a timescale of a few Gyrs, although the inflows may be
episodic. On the other hand, Irodotou et al. (2022) suggested there is an
anti-correlation between AGN feedback and bar strength in cosmological
simulations, i.e. bars tend to be stronger and shorter if AGN feedback is not
included. Future studies of the relation between AGN feedback and bar
stability would be very useful to better understand this issue.

###  5.3 Do we expect a clear correlation between S2Bs and AGNs?

In our S2B model, gas is funneled close to the Bondi radius (a few tens of
parsecs) if a $1.0\times 10^{8}\;{\rm M}_{\odot}$ SMBH (Schwarzschild radius:
$R_{S}=9.57\times 10^{-6}\;{\rm pc}$, Eddington accretion rate: $\dot{M}_{\rm
Edd}=2.2{\;{\rm M}_{\odot}\;{\rm yr}^{-1}}$) had been present at the galatic
center. The accretion flows around a SMBH are generally classified into the
radiative mode, also known as quasar mode, which operates when the SMBH has
high accretion rates; and kinetic mode, also known as radio mode, which is
typically important when the SMBH has low accretion rates. The boundary
between the two modes is usually $\sim 2\%\dot{M}_{\rm Edd}=0.043{\;{\rm
M}_{\odot}\;{\rm yr}^{-1}}$ (summarized in Yuan et al., 2020). Our fiducial
S2B model with a lightweight gaseous disk has the same order of gas inflow
rate, and this level of fueling can be easily achieved in local disk galaxies
that host quasars (Zhao et al., 2019). If the SMBHs in S2Bs are in quasar
mode, the timescale for the gas to flow from the accretion disk
($R=1000R_{S}\approx 0.01\;{\rm pc}$) to feed the SMBH is $0.44\;{\rm Gyr}$
based on Eq. 6 of Yuan et al. (2018). With an additional free-fall timescale
(a few tens of Myr) for gas to flow into the accretion disk from the Bondi
radius, the overall timescale of gas inflows from the Bondi radius to feed the
central SMBH can be estimated at a few hundred Myr, which is longer than the
rotation period of the inner bar in our model. Though gas inflows are more
prominent when the inner bar is perpendicular to its outer counterpart, it is
not likely to result in a clear correlation between the relative orientation
of the two bars and any AGN activity.

Do we expect to find a higher frequency of AGNs in S2B galaxies? Our
simulations suggest that SMBHs may grow rapidly when an inner bar exists.
About $10^{7}\;{\rm M}_{\odot}$ of gas accumulates at the central $100\;{\rm
pc}$ of S2BSMU within $0.5\;{\rm Gyr}$ even with a rather lightweight initial
gaseous disk. As suggested in Du et al. (2017) and Guo et al. (2020), and also
from the results of this paper, the inner bar will be destroyed once enough
gas accumulates in the galaxy center, no matter whether the gas has fallen
into the SMBH. The inner bar may be in the process of dissolution, considering
a potential delay of AGN activity by around a few hundred Myr estimated above.
We also note that any possible link between S2Bs and AGNs may disappear if the
first inner bar that powers the AGN decays and a new one forms at later times,
similar to the scenario proposed for large-scale bars (e.g. Sellwood & Moore,
1999; Li et al., 2017). It is therefore understandable that no clear relation
between inner bars and AGNs is reported in observations (Laine et al., 2002;
Erwin & Sparke, 2002; Erwin, 2011). On the other hand, recent studies seem to
suggest that AGNs are preferentially found in barred galaxies (Alonso et al.,
2018; Silva-Lima et al., 2022), although the S2B fraction in the samples
remains unclear.

###  5.4 Are inner bars long-lived or recurrent?

The high frequency of observed S2Bs (Erwin, 2011) indicates that inner bars
are either long-lived or easily reformed after dissolution. de Lorenzo-Cáceres
et al. (2019a) performed a detailed photometric analysis of 17 nearby S2Bs and
found that most bulges in their sample have intermediate Sérsic index
($n\approx 2$) and lay at the top sequence of the Kormendy relation (Kormendy,
1977), i.e. they are similar to classical bulges. On the other hand, the
simulation of Guo et al. (2020) show that the remnant of the dissolved inner
bar share such properties with those of observed bulges in S2Bs. This makes
one wonder if the inner bar in S2B can reform after dissolution instead of
being a long-lived structure. If inner bars have dissolved and reformed over
time, we should see these classical bulges as relics in observed S2Bs, which
seems to be consistent with de Lorenzo-Cáceres et al. (2019a). However, since
there are multiple ways to form classical bulges, this does not necessarily
mean that inner bars have to be recurrent. It is possible that two kinds of
inner bars exist in the universe as de Lorenzo-Cáceres et al. (2020) found
there are two distinct groups of inner bars in terms of their in-plane length
and ellipticity. Whether the properties of long-lived inner bars are
systematically different from those of recurrent ones needs further
investigation.

##  6 Summary

In this study, we have used high-resolution numerical simulations to explore
gas flow patterns in a galaxy with two independently rotating bars (S2Bs).
Different gas physics setups are tested, and the gravitational effect of the
inner bar is highlighted with control models. Our main findings are summarized
as follows:

(1) We find the inner bar can drive gas flow inwards. The inflow is periodic,
and is stronger when the two bars are perpendicular with respect to each
other, possibly due to more clear shocks in the inner bar at this orientation
(Figure 5). The time-averaged central inflow rate in our S2B model is about 3
times higher than that in the single-barred one (Figure 8).

(2) Gas forms multiple streams down to the center from the nuclear ring in our
S2B model, while this is not observed in the single-barred simulation (Figure
7). Outside the nuclear ring, the flow patterns are very similar between
single- and double-barred galaxies. A gas inner bar is formed in the S2B model
without stellar feedback, but it is erased when feedback is included (Figure
6). The gas flow pattern in our S2B model with SMUGGLE enabled is
qualitatively consistent with a few observed S2Bs (Figure 13).

(3) The S2B model has the highest star formation rate compared with the
single- and non-barred models (Figure 9) at a given stellar mass. This may be
explained by a higher star formation efficiency (or a shorter depletion time),
together with a higher accumulated gas mass in the center of the S2B model. We
also find that the depletion time in the central region is shorter than that
in the disk, similar to recent observations (Figure 10).

We thank the anonymous referee for suggestions that help to improve the
presentation of the paper. We thank Feng Yuan and Mou-Yuan Sun for insightful
discussions. The research presented here is supported by the National Natural
Science Foundation of China under grant No. 12103032, and is partially
supported by the National Key R&D Program of China under grant No.
2018YFA0404501; by the National Natural Science Foundation of China under
grant Nos. 12025302, 11773052, 11761131016; by the “111” Project of the
Ministry of Education of China under Grant No. B20019; and by the Chinese
Space Station Telescope project. M.D. acknowledges the support of the
Fundamental Research Funds for the Central Universities. J.S. acknowledges the
support of a Newton Advanced Fellowship awarded by the Royal Society and the
Newton Fund. A.B. was supported by the Future Investigators in NASA Earth and
Space Science and Technology (FINESST) award number 80NSSC20K1536 during the
completion of this work. This work made use of the Gravity Supercomputer at
the Department of Astronomy, Shanghai Jiao Tong University, the facilities of
the Center for High Performance Computing at Shanghai Astronomical
Observatory, and the High-performance Computing Platform of Peking University.

## References

  * Alonso et al. (2018) Alonso, S., Coldwell, G., Duplancic, F., Mesa, V., & Lambas, D. G. 2018, A&A, 618, A149, doi: 10.1051/0004-6361/201832796
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Athanassoula (1992) Athanassoula, E. 1992, MNRAS, 259, 345
  * Athanassoula (2003) —. 2003, MNRAS, 341, 1179, doi: 10.1046/j.1365-8711.2003.06473.x
  * Beane et al. (2023) Beane, A., Hernquist, L., D’Onghia, E., et al. 2023, ApJ, 953, 173, doi: 10.3847/1538-4357/ace2b9
  * Bittner et al. (2021) Bittner, A., de Lorenzo-Cáceres, A., Gadotti, D. A., et al. 2021, A&A, 646, A42, doi: 10.1051/0004-6361/202039505
  * Brown et al. (2013) Brown, J. S., Valluri, M., Shen, J., & Debattista, V. P. 2013, ApJ, 778, 151, doi: 10.1088/0004-637X/778/2/151
  * Burger et al. (2022) Burger, J. D., Zavala, J., Sales, L. V., et al. 2022, MNRAS, 513, 3458, doi: 10.1093/mnras/stac994
  * Buta & Block (2001) Buta, R., & Block, D. L. 2001, ApJ, 550, 243, doi: 10.1086/319736
  * Buta & Crocker (1993) Buta, R., & Crocker, D. A. 1993, AJ, 105, 1344, doi: 10.1086/116514
  * Buta et al. (2001) Buta, R., Ryder, S. D., Madsen, G. J., et al. 2001, AJ, 121, 225, doi: 10.1086/318024
  * Buta et al. (2015) Buta, R. J., Sheth, K., Athanassoula, E., et al. 2015, ApJS, 217, 32, doi: 10.1088/0067-0049/217/2/32
  * de Lorenzo-Cáceres et al. (2008) de Lorenzo-Cáceres, A., Falcón-Barroso, J., Vazdekis, A., & Martínez-Valpuesta, I. 2008, ApJ, 684, L83, doi: 10.1086/592145
  * de Lorenzo-Cáceres et al. (2019a) de Lorenzo-Cáceres, A., Méndez-Abreu, J., Thorne, B., & Costantin, L. 2019a, MNRAS, 484, 665, doi: 10.1093/mnras/sty3520
  * de Lorenzo-Cáceres et al. (2020) —. 2020, MNRAS, 494, 1826, doi: 10.1093/mnras/staa780
  * de Lorenzo-Cáceres et al. (2019b) de Lorenzo-Cáceres, A., Sánchez-Blázquez, P., Méndez-Abreu, J., et al. 2019b, MNRAS, 484, 5296, doi: 10.1093/mnras/stz221
  * Debattista & Shen (2007) Debattista, V. P., & Shen, J. 2007, ApJ, 654, L127, doi: 10.1086/511264
  * Du et al. (2016) Du, M., Debattista, V. P., Shen, J., & Cappellari, M. 2016, ApJ, 828, 14, doi: 10.3847/0004-637X/828/1/14
  * Du et al. (2017) Du, M., Debattista, V. P., Shen, J., Ho, L. C., & Erwin, P. 2017, ApJ, 844, L15, doi: 10.3847/2041-8213/aa7ecb
  * Du et al. (2015) Du, M., Shen, J., & Debattista, V. P. 2015, ApJ, 804, 139, doi: 10.1088/0004-637X/804/2/139
  * Englmaier & Gerhard (1997) Englmaier, P., & Gerhard, O. 1997, MNRAS, 287, 57
  * Erwin (2004) Erwin, P. 2004, A&A, 415, 941, doi: 10.1051/0004-6361:20034408
  * Erwin (2005) —. 2005, MNRAS, 364, 283, doi: 10.1111/j.1365-2966.2005.09560.x
  * Erwin (2011) —. 2011, Memorie della Societa Astronomica Italiana Supplementi, 18, 145
  * Erwin & Sparke (2002) Erwin, P., & Sparke, L. S. 2002, AJ, 124, 65, doi: 10.1086/340803
  * Fanali et al. (2015) Fanali, R., Dotti, M., Fiacconi, D., & Haardt, F. 2015, MNRAS, 454, 3641, doi: 10.1093/mnras/stv2247
  * Fathi et al. (2006) Fathi, K., Storchi-Bergmann, T., Riffel, R. A., et al. 2006, ApJ, 641, L25, doi: 10.1086/503832
  * Fragkoudi et al. (2017) Fragkoudi, F., Athanassoula, E., & Bosma, A. 2017, MNRAS, 466, 474, doi: 10.1093/mnras/stw3023
  * Friedli & Martinet (1993) Friedli, D., & Martinet, L. 1993, A&A, 277, 27
  * Gerhard (1993) Gerhard, O. E. 1993, MNRAS, 265, 213, doi: 10.1093/mnras/265.1.213
  * González-Alfonso et al. (2021) González-Alfonso, E., Pereira-Santaella, M., Fischer, J., et al. 2021, A&A, 645, A49, doi: 10.1051/0004-6361/202039047
  * Guo et al. (2020) Guo, M., Du, M., Ho, L. C., Debattista, V. P., & Zhao, D. 2020, ApJ, 888, 65, doi: 10.3847/1538-4357/ab584a
  * Haan et al. (2009) Haan, S., Schinnerer, E., Emsellem, E., et al. 2009, ApJ, 692, 1623, doi: 10.1088/0004-637X/692/2/1623
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Heller et al. (2007) Heller, C. H., Shlosman, I., & Athanassoula, E. 2007, ApJ, 657, L65, doi: 10.1086/513104
  * Hopkins & Quataert (2010) Hopkins, P. F., & Quataert, E. 2010, MNRAS, 407, 1529, doi: 10.1111/j.1365-2966.2010.17064.x
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Irodotou et al. (2022) Irodotou, D., Fragkoudi, F., Pakmor, R., et al. 2022, MNRAS, 513, 3768, doi: 10.1093/mnras/stac1143
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Khoperskov et al. (2018) Khoperskov, S., Haywood, M., Di Matteo, P., Lehnert, M. D., & Combes, F. 2018, A&A, 609, A60, doi: 10.1051/0004-6361/201731211
  * Kim et al. (2012) Kim, W.-T., Seo, W.-Y., Stone, J. M., Yoon, D., & Teuben, P. J. 2012, ApJ, 747, 60, doi: 10.1088/0004-637X/747/1/60
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Kolcu et al. (2023) Kolcu, T., Maciejewski, W., Gadotti, D. A., et al. 2023, MNRAS, 524, 207, doi: 10.1093/mnras/stad1896
  * Kormendy (1977) Kormendy, J. 1977, ApJ, 218, 333, doi: 10.1086/155687
  * Laine et al. (2002) Laine, S., Shlosman, I., Knapen, J. H., & Peletier, R. F. 2002, ApJ, 567, 97, doi: 10.1086/323964
  * Lang et al. (2020) Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122, doi: 10.3847/1538-4357/ab9953
  * Leroy et al. (2021) Leroy, A. K., Hughes, A., Liu, D., et al. 2021, ApJS, 255, 19, doi: 10.3847/1538-4365/abec80
  * Li et al. (2022) Li, H., Vogelsberger, M., Bryan, G. L., et al. 2022, MNRAS, 514, 265, doi: 10.1093/mnras/stac1136
  * Li et al. (2017) Li, Z., Sellwood, J. A., & Shen, J. 2017, ApJ, 850, 67, doi: 10.3847/1538-4357/aa9377
  * Li et al. (2015) Li, Z., Shen, J., & Kim, W.-T. 2015, ApJ, 806, 150, doi: 10.1088/0004-637X/806/2/150
  * Li et al. (2014) Li, Z.-Y., Shen, J., Rich, R. M., Kunder, A., & Mao, S. 2014, ApJ, 785, L17, doi: 10.1088/2041-8205/785/1/L17
  * Lin et al. (2017) Lin, L., Li, C., He, Y., Xiao, T., & Wang, E. 2017, ApJ, 838, 105, doi: 10.3847/1538-4357/aa657a
  * Lin et al. (2020) Lin, L., Li, C., Du, C., et al. 2020, MNRAS, 499, 1406, doi: 10.1093/mnras/staa2913
  * Lin et al. (2013) Lin, L.-H., Wang, H.-H., Hsieh, P.-Y., et al. 2013, ApJ, 771, 8, doi: 10.1088/0004-637X/771/1/8
  * Maciejewski & Sparke (2000) Maciejewski, W., & Sparke, L. S. 2000, MNRAS, 313, 745, doi: 10.1046/j.1365-8711.2000.03270.x
  * Maciejewski et al. (2002) Maciejewski, W., Teuben, P. J., Sparke, L. S., & Stone, J. M. 2002, MNRAS, 329, 502, doi: 10.1046/j.1365-8711.2002.04957.x
  * Marinacci et al. (2019) Marinacci, F., Sales, L. V., Vogelsberger, M., Torrey, P., & Springel, V. 2019, MNRAS, 489, 4233, doi: 10.1093/mnras/stz2391
  * Méndez-Abreu et al. (2019) Méndez-Abreu, J., de Lorenzo-Cáceres, A., Gadotti, D. A., et al. 2019, MNRAS, 482, L118, doi: 10.1093/mnrasl/sly196
  * Moon et al. (2022) Moon, S., Kim, W.-T., Kim, C.-G., & Ostriker, E. C. 2022, ApJ, 925, 99, doi: 10.3847/1538-4357/ac3a7b
  * Nakatsuno & Baba (2023) Nakatsuno, N., & Baba, J. 2023, arXiv e-prints, arXiv:2308.15029, doi: 10.48550/arXiv.2308.15029
  * Namekata et al. (2009) Namekata, D., Habe, A., Matsui, H., & Saitoh, T. R. 2009, ApJ, 691, 1525, doi: 10.1088/0004-637X/691/2/1525
  * Querejeta et al. (2021) Querejeta, M., Schinnerer, E., Meidt, S., et al. 2021, A&A, 656, A133, doi: 10.1051/0004-6361/202140695
  * Rautiainen et al. (2002) Rautiainen, P., Salo, H., & Laurikainen, E. 2002, MNRAS, 337, 1233, doi: 10.1046/j.1365-8711.2002.05955.x
  * Schinnerer et al. (2007) Schinnerer, E., Böker, T., Emsellem, E., & Downes, D. 2007, A&A, 462, L27, doi: 10.1051/0004-6361:20066711
  * Schinnerer et al. (2006) Schinnerer, E., Böker, T., Emsellem, E., & Lisenfeld, U. 2006, ApJ, 649, 181, doi: 10.1086/506265
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243, doi: 10.1086/146614
  * Sellwood & Moore (1999) Sellwood, J. A., & Moore, E. M. 1999, ApJ, 510, 125, doi: 10.1086/306557
  * Sellwood & Spekkens (2015) Sellwood, J. A., & Spekkens, K. 2015, arXiv e-prints, arXiv:1509.07120, doi: 10.48550/arXiv.1509.07120
  * Shen & Debattista (2009) Shen, J., & Debattista, V. P. 2009, ApJ, 690, 758, doi: 10.1088/0004-637X/690/1/758
  * Shen & Sellwood (2004) Shen, J., & Sellwood, J. A. 2004, ApJ, 604, 614, doi: 10.1086/382124
  * Shlosman (2005) Shlosman, I. 2005, in American Institute of Physics Conference Series, Vol. 783, The Evolution of Starbursts, ed. S. Hüttmeister, E. Manthey, D. Bomans, & K. Weis, 223–240, doi: 10.1063/1.2034990
  * Shlosman et al. (1990) Shlosman, I., Begelman, M. C., & Frank, J. 1990, Nature, 345, 679, doi: 10.1038/345679a0
  * Shlosman et al. (1989) Shlosman, I., Frank, J., & Begelman, M. C. 1989, Nature, 338, 45, doi: 10.1038/338045a0
  * Shlosman & Heller (2002) Shlosman, I., & Heller, C. H. 2002, ApJ, 565, 921, doi: 10.1086/324403
  * Silva-Lima et al. (2022) Silva-Lima, L. A., Martins, L. P., Coelho, P. R. T., & Gadotti, D. A. 2022, A&A, 661, A105, doi: 10.1051/0004-6361/202142432
  * Sivasankaran et al. (2022) Sivasankaran, A., Blecha, L., Torrey, P., et al. 2022, MNRAS, 517, 4752, doi: 10.1093/mnras/stac2759
  * Sormani et al. (2015) Sormani, M. C., Binney, J., & Magorrian, J. 2015, MNRAS, 449, 2421, doi: 10.1093/mnras/stv441
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Springel et al. (2005) Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776, doi: 10.1111/j.1365-2966.2005.09238.x
  * Sun et al. (2023) Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2023, ApJ, 945, L19, doi: 10.3847/2041-8213/acbd9c
  * Tacchella et al. (2022) Tacchella, S., Smith, A., Kannan, R., et al. 2022, MNRAS, 513, 2904, doi: 10.1093/mnras/stac818
  * Tress et al. (2020) Tress, R. G., Sormani, M. C., Glover, S. C. O., et al. 2020, MNRAS, 499, 4455, doi: 10.1093/mnras/staa3120
  * van der Marel & Franx (1993) van der Marel, R. P., & Franx, M. 1993, ApJ, 407, 525, doi: 10.1086/172534
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vogelsberger et al. (2013) Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031, doi: 10.1093/mnras/stt1789
  * Wang et al. (2010) Wang, H.-H., Klessen, R. S., Dullemond, C. P., van den Bosch, F. C., & Fuchs, B. 2010, MNRAS, 407, 705, doi: 10.1111/j.1365-2966.2010.16942.x
  * Weinberger et al. (2020) Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS, 248, 32, doi: 10.3847/1538-4365/ab908c
  * Weiner et al. (2001) Weiner, B. J., Sellwood, J. A., & Williams, T. B. 2001, ApJ, 546, 931, doi: 10.1086/318289
  * Wozniak (2015) Wozniak, H. 2015, A&A, 575, A7, doi: 10.1051/0004-6361/201425005
  * Wu et al. (2016) Wu, Y.-T., Pfenniger, D., & Taam, R. E. 2016, ApJ, 830, 111, doi: 10.3847/0004-637X/830/2/111
  * Wu et al. (2021) Wu, Y.-T., Trejo, A., Espada, D., & Miyamoto, Y. 2021, MNRAS, 504, 3111, doi: 10.1093/mnras/stab1087
  * Yuan et al. (2020) Yuan, F., Ostriker, J. P., Yoon, D., et al. 2020, in Perseus in Sicily: From Black Hole to Cluster Outskirts, ed. K. Asada, E. de Gouveia Dal Pino, M. Giroletti, H. Nagai, & R. Nemmen, Vol. 342, 101–107, doi: 10.1017/S174392131800371X
  * Yuan et al. (2018) Yuan, F., Yoon, D., Li, Y.-P., et al. 2018, ApJ, 857, 121, doi: 10.3847/1538-4357/aab8f8
  * Zhao et al. (2019) Zhao, D., Ho, L. C., Zhao, Y., Shangguan, J., & Kim, M. 2019, ApJ, 877, 52, doi: 10.3847/1538-4357/ab1921