Nuclear rings are the inner edge of a gap around the Lindblad Resonance

Gaseous nuclear rings are large-scale coherent structures commonly found at the centres of barred galaxies. We propose that they are an accumulation of gas at the inner edge of an extensive gap that forms around the Inner Lindblad Resonance (ILR). The gap initially opens because the bar potential excites strong trailing waves near the ILR, which remove angular momentum from the gas disc and transport the gas inwards. The gap then widens because the bar potential continuously excites trailing waves at the inner edge of the gap, which remove further angular momentum, moving the edge further inwards until it stops at a distance of several wavelengths from the ILR. The gas accumulating at the inner edge of the gap forms the nuclear ring. The speed at which the gap edge moves and its final distance from the ILR strongly depend on the sound speed, explaining the puzzling dependence of the nuclear ring radius on the sound speed in simulations.


INTRODUCTION
Gaseous nuclear rings are remarkable structures commonly found at the centres of barred galaxies.They have typical radii of 50-1000 pc (Comerón et al. 2010), total gas masses of 10 8 -10 9 M ⊙ (Sheth et al. 2005;Querejeta et al. 2021), and star formation rates spanning a wide range 0.1-10 M ⊙ yr −1 (Mazzuca et al. 2008;Ma et al. 2018).They are among the most intense star-forming regions of disc galaxies and are considered special laboratories to study star formation under extreme conditions (Moon et al. 2021;Schinnerer et al. 2023).They are sites where galactic outflows can be launched, with profound impact on the evolution of their host galaxies (Veilleux et al. 2020).They constitute cold gas reservoirs for the fuelling of central supermassive black holes.The Milky Way hosts a nuclear ring with a radius of  ≃ 120 pc that is better known as the Central Molecular Zone (Morris & Serabyn 1996;Henshaw et al. 2023).
It is well-known that nuclear rings are easy to form in simulations (e.g.Athanassoula 1992b; Kim et al. 2012;Sormani et al. 2015a, and many others).The recipe is simple: let gas flow in a non-axisymmetric rotating barred potential, and a nuclear ring will spontaneously form in the central regions.In the simplest simulations, the gas is assumed to be 2D, isothermal, non-self gravitating, and the barred potential is externally imposed, but a ring can form also if additional physics is included, for example the gas self-gravity, star formation & stellar feedback, live stellar potentials, or magnetic fields (Fux 1999;Armillotta et al. 2019;Tress et al. 2020).However, being able to watch the ring forming in simulations does not mean that we understand the underlying physical process by which it forms, which has remained elusive.
Despite the interest from several astrophysical communities, the ★ E-mail: mattiacarlo.sormani@gmail.comphysical mechanism by which nuclear rings form is not well understood.What sets the radius of the nuclear ring?What is "special" about its location?Various theories have been proposed, but we argue that they all fail to explain the formation of the rings.These previous theories are reviewed in Sect.3.2.
One of the most puzzling aspects is that the radius of the nuclear ring in isothermal simulations of gas flow in a barred potential depends very strongly on the assumed sound speed (e.g.Englmaier & Gerhard 1997;Patsis & Athanassoula 2000;Kim et al. 2012;Sormani et al. 2015a).For example, doubling the sound speed from  s = 5 km s −1 to 10 km s −1 can change the radius of the ring by a factor of two or more (see for example Figure 2 in Sormani et al. 2015a).This is surprising because the sound speed always amounts to just a small fraction of the orbital speed (typically ∼5%).The flow is always strongly supersonic.None of the currently available theories can explain the strong dependence of the ring radius on the sound speed.
In this paper, we develop a framework to understand the formation of nuclear rings.We propose that the rings are in fact the inner edge of an extensive gap that opens around the ILR due to the excitation of waves by a bar potential.These waves remove angular momentum from the gas disc, transporting the gas inwards.The nuclear ring forms due to the accumulation of gas at the inner edge of the gap.
The paper is structured as follows.In Section 2 we present some numerical experiments that illustrate the formation of nuclear rings in simulations.In Section 3 we review the constraints that we believe any plausible theory for the formation of nuclear rings should satisfy, and we review previous theories.In Section 4 we study the excitation of density waves by an external bar potential using linear theory.In Section 5 we illustrate our picture of the formation of the rings.In Section 6 we discuss various connections between this paper and previous works, in particular the works of Goldreich & Tremaine (1978, 1979) that studied the opening of the Cassini gap in Saturn's rings.We sum up in Section 7.

NUMERICAL EXPERIMENTS
We first perform some numerical experiments by letting non-selfgravitating isothermal gas flow in an external barred potential.This is useful to establish some key points and parameter dependencies that we will need later.

Numerical setup
We run a total of six 2D non-self-gravitating isothermal simulations of gas flowing in an external barred gravitational potential.The potential is described in Appendix A. Table 1 provides a summary of the simulations run.The equations of motion are: where  is the surface density, v is the gas velocity, Φ is the external gravitational potential given by Eq. (A1) and is the isothermal equation of state, where  s = constant.We will use values in the range  s = 1-20 km s −1 .We solve Eqs. ( 1) and (2) using the public grid code PLUTO (Mignone et al. 2007) on a two-dimensional static polar grid in the region  ×  = [0.1 kpc,  max ] × [0, 2].The grid is logarithmically spaced in  and uniformly spaced in  with   × 1024 cells.The resolution along the  direction is approximately Δ = 0.00529 , i.e. we have a resolution of Δ = 0.529 pc at the inner boundary of  = 100 pc.The number of cells in each direction is chosen so that the aspect ratio of the cells is approximately Δ (/Δ) ∼ 1.We use the following parameters: rk2 time-stepping, no dimensional splitting, hll Riemann solver and the default flux limiter.We solve the equations in the frame rotating at Ω p by using the rotating_frame = yes switch.Boundary conditions are outflow both on the inner boundary at  = 0.1 kpc and on the outer boundary at  =  max .
The initial density distribution is (4) Note that, since the equations of motion (1) and ( 2) are invariant under density rescaling, the density units are arbitrary.The quantity ρ therefore essentially sets the density units, and without loss of generality we set ρ = 1.The quantity   = 10 −12 ρ corresponds to the density floor imposed in the simulation to avoid crashing.We introduce the bar gradually to reduce transients (e.g.Athanassoula 1992b).We start with gas in equilibrium on circular orbits in the logarithmic axisymmetric potential Φ 0 and then linearly turn on the non-axisymmetric part of the potential Φ 1 during the first 313 Myr.

Disc with initial radius smaller than the ILR
Simulations 01-05 investigate the evolution of a uniform gas disc with an initial radius  disc = 1.2 kpc that is smaller than  ILR = 1.61 kpc (Appendix A).The only difference between these five simulations is the assumed sound speed (Table 1).Figure 1 shows the surface density as a function of time.As soon as the bar potential is turned on, trailing spiral waves are excited.These waves are clearly visible at  = 157 Myr and  = 313 Myr.The movies of the surface density as a function of time show that the waves are first excited at the outer edge of the disc, and propagate inwards.We will confirm later in Sects.4.3.4 and 4.3.5 using linear analysis that sharp edges are indeed regions where strong wave excitation takes place, and therefore play a key role in the formation of the rings.The wiggles that are visible along the spirals in some panels (for example the panel at  = 313 Myr and  s = 10 km s −1 ) are due to the wiggle instability (Wada & Koda 2004;Kim et al. 2014;Sormani et al. 2017;Mandowara et al. 2022).
Figure 2 plots a cut through the  axis of Fig. 1 at  = 157 Myr.The radial wavelength increases with increasing sound speed.This will be explained by the dispersion relation derived below (Eq.54).The amplitude of the waves decreases inward, despite the prediction of the linear analysis according to which the amplitude of density waves should increase inward due to geometric effects (see Sect 4).The reason for this behaviour is that the waves in the simulations become quickly non-linear and develop shocks.The shocks cause the waves to dissipate, decreasing their amplitude and depositing their (negative) angular momentum into the gas disc.As we will argue in Sect.5, this process is what decreases the angular momentum of the gas disc and causes it to shrink.
The final size of the ring depends very strongly on the sound speed (rightmost column in Fig. 1).This is further quantified in Fig. 3, which shows the evolution of the ring size as a function of sound speed.As can be seen in the bottom panel, increasing the sound speed by a factor of two can change the final ring size by the same factor.

Disc with initial radius larger than the ILR
Simulations 04 and 04_Large only differ in the size of the initial gas disc,  disc = 1.2 kpc vs.  disc = 5 kpc.Thus, simulation 04 includes only the flow inside the ILR ( ILR = 1.61 kpc), while simulation 04_Large comprises the large-scale flow in the entire bar region.
Figure 4 shows that the final ring size is approximately the same in both simulations, and therefore that ring size does not depend on the large-scale flow outside the ILR.This implies that the mechanism determining the radius of the ring must be "local" (see point 5 in Sect.3.1).
Figure 5 illustrates the evolution of the axisymmetrised surface density as a function of radius in the simulation 04_Large.In particular, we can see that an extensive gap of low surface density is opened around the ILR.The nuclear ring is the inner edge of this gap, where the material that once was in the gap has accumulated.

Conditions that a plausible theory must satisfy
We introduce the conditions that we believe any plausible theory for the formation of nuclear rings must satisfy.We take the approach that numerical experiments, such as those in Sect.2, guide us on how the ring properties should depend on the underlying parameters.We summarise the insights obtained from simulations into the following five conditions: 1.The radius of the ring must depend on the circular rotation curve.Athanassoula (1992b) and Li et al. (2015) have shown that the radius of the ring in simulations changes if we change the circular velocity curve of the underlying gravitational potential, i.e. if we change the axisymmetric part of the gravitational potential (see in particular Fig. 4 in Li et al. 2015), while keeping everything else fixed.
2. The radius of the ring must depend on the non-axisymmetric part of the underlying potential.Sormani et al. (2015b) has shown that the radius of the ring can change significantly if we change the quadrupole of the potential while keeping the monopole (and therefore the rotation curve) fixed.Hence, a theory aiming to explain why rings form at a certain location must take into account a dependence on the non-axisymmetric part of the potential.
3. The radius of the ring must depend on the bar pattern speed.Many authors (e.g.Athanassoula 1992b;Li et al. 2015;Sormani et al. 2015b, among others) have shown that the radius of the ring depends on the rotation speed of the bar.
4. The radius of the ring must depend on the equation of state of the gas.Many authors (e.g.Englmaier & Gerhard 1997;Patsis & Athanassoula 2000;Kim et al. 2012;Sormani et al. 2015a, among many others) have shown that the size of the ring strongly depends on the sound speed.This is confirmed by the numerical experiments we conducted in Section 2.2 (see in particular Figs. 1 and 3).Thus,  the radius of the ring does not depend solely on the gravitational potential but must involve the equation of state of the gas.
5. The radius of the ring must be determined "locally", i.e. the final ring size should not depend on the larger-scale flow at  >  ILR .This is demonstrated by the numerical experiments in Fig. 4, which shows two simulations that differ only for the extent of the simulated gas disc.The 04_Large simulation (left) covers the entire "bar region", out to  disc = 5 kpc.It includes the usual bar-driven accretion flow from the disc to the ring.The 04 simulation (right) is the same shown in the second row of Fig. 1, and it only simulates a gas disc of  disc = 1.2 kpc, which is all contained within the ILR at  ILR = 1.61 kpc (Fig. A1).The final size of the ring is essentially the same in the two simulations (it is slightly larger in the 04_Large simulation because fresh gas is continuously brought from outside the ILR, which takes longer to lose angular momentum).This shows that removal of angular momentum also happens in the vicinity of the ring.Material from outside the ILR that crosses the ILR continually loses angular momentum up to when it settles on the ring.

Previous theories for the formation of nuclear rings
Here we briefly summarise previous theories for the formation of the ring and for determining its location.We argue that none of them provides a satisfactory explanation for the formation of the rings by showing that each of them fails to satisfy at least one of the conditions outlined in Section 3.1.
• The resonant theory (Combes 1988(Combes , 1996;;Buta & Combes 1996).This is perhaps the most widely accepted theory, especially in the extragalactic community.It states that the ring forms at the Lindblad resonance under the continuous action of gravity torques from the bar potential.
× Refutation: This theory satisfies conditions 1, 3, 5, and if the notion of ILR is generalised to include strongly barred potentials (van Albada & Sanders 1982;Athanassoula 1992a), it may satisfy condition 2. However, since the position of the resonance does not depend on the equation of state of the gas, it does not satisfy condition 4.Moreover, the numerical experiments shown in Fig. 1 show that the radius of the ring forms at a radius  that is much smaller than  ILR = 1.61 kpc (on this point, see also Regan & Teuben 2003).
• The minimum shear theory (Lesch et al. 1990;Krumholz & Kruĳssen 2015).This theory states that the ring forms at the radius at which the shear, as calculated from the axisymmetric rotation curve, is minimum.This conclusion stems from an analogy with accretion disc theory, in which transport is more efficient where shear is higher, so gas is expected to pile up and form a ring at the point of minimum shear.
× Refutation: This theory satisfies conditions 1 and 5, but does not satisfy conditions 2, 3, 4. Sormani & Li (2020) demonstrate in detail that simulations are inconsistent with this theory.
• The reverse shear theory (Sormani et al. 2018b).This theory states that the ring forms in a region where a family of nonaxisymmetric closed periodic orbits called  2 orbits displays "reverse shear" that prevents viscous spreading, making it possible to confine a stable ring.
× Refutation: This theory satisfies conditions 1, 2, 3, and 5, but it does not satisfy condition 4 since it predicts that the radius of the ring depends exclusively on the gravitational potential.
In conclusion, none of the currently available theories satisfies all the criteria introduced in Sect.3.1, highlighting the need for a new theory.

LINEAR DISC DYNAMICS
The simulations in Sect. 2 suggest that density waves are important for the removal of angular momentum of the disc and the opening of the gap.To gain insight into this process, in this section we study the excitation of waves by the bar potential in the linear regime, and estimate the amount of angular momentum that they remove from the gas disc as a function of the unperturbed density profile  0 and of the sound speed  s .As we shall see, the excitation of waves happens primarily at two locations: i) near the ILR.This regime has been studied in detail by Goldreich & Tremaine (1979), and we will not repeat their calculations here; ii) at sharp edges (i.e., strong gradients) in the unperturbed density distribution  0 .The insight gained in this section will be used in Sect. 5 to develop a new picture for the formation of nuclear rings.
Consider a 2D axisymmetric differentially rotating fluid disc in equilibrium in an external gravitational potential.Our goal is to study the propagation of small perturbations (waves) and their excitation by a "small" external potential.These "density waves" are conceptually similar to sound waves in air, but the rotation makes the dynamics richer and mathematically more complex.In Appendix B we present a 1D toy problem that can be solved fully analytically and provides a mathematically simpler analog to the more complicated problem studied in this section.
We start from the disc's unperturbed steady state and linearize the equations of motions around it.We ignore the self-gravity of the gas.The equations of motion are the continuity and Euler equations: where  is the surface density, v =   ê +   ê is the velocity,  is the pressure, and Φ(x, ) is the external gravitational potential.We assume a polytropic equation of state where  ≥ 1 and  is a constant.To simplify the calculations it is convenient to introduce the enthalpy ℎ defined by: substituting ( 7) into (8) and integrating we find: Using (8), the equations of motion ( 5) and ( 6) can be expanded in polar coordinates (, ) as:

Unperturbed state
We assume that the density, velocity and gravitational potential of the unperturbed steady-state are: Substituting these into (10)-( 12) and assuming steady-state and axisymmetry (  =   = 0), we see that the continuity equation ( 10) and the azimuthal Euler equation ( 12) are already satisfied, while the radial Euler equation ( 11) gives: In the following, ℎ 0 , Φ 0 and Ω are prescribed functions of  that satisfy Eq. ( 17).Note that given Φ 0 (), there formally exists an equilibrium solution ℎ 0 () for any arbitrary rotation profile Ω().
However not all possible profiles are physical.To avoid instability, the unperturbed state must satisfy the Rayleigh stability criterion, which states that a necessary and sufficient condition for the local axisymmetric stability of an inviscid differentially rotating fluid disc is that the specific angular momentum monotonically increases with In this paper we will assume mainly two types of density profiles.
The first is a constant density profile The second is a truncated disc profile, i.e. a density that is roughly constant at  ≪  edge , has a relatively sharp transition at an edge  edge during which it drops at a much lower value, and is then roughly constant again at  >  edge .Note that the edge cannot be made too thin, otherwise it would violate the Rayleigh criterion (18).When later in the paper it will be necessary to assume a specific truncated profile for numerical calculations, we will use the following: where where  edge is the position of the edge and Δ controls its width.
The quantity ρ is a constant that, as noted in Sect.2.1, essentially defines the units used for density.Physically meaningful results do not depend on the particular value of this quantity since the equations of motion ( 5) and ( 6) are invariant under density rescaling.Without loss of generality, we set ρ = 1.

Linearised equations
To study the propagation of small waves on top of the unperturbed state described in the previous section, we expand all quantities as: Substituting equations ( 22)-( 25) into ( 10)-( 12) and linearising by keeping only first-order terms in the quantities with subscript 1, we obtain: where we have defined the convective derivative of the unperturbed state and the Oort parameter Without loss of generality, we can write all the "small" subscript-1 quantities as: where ρ1 , ṽ1 etc. are complex, and the "physical" quantity is the real part.The general solution of equations ( 26)-( 28) can always be decomposed in such modes because the equations are linear and the superposition principle applies.Each mode evolves independently from the others in the linear approximation.In this paper, we will only be concerned with  = 2 as this is the only non-zero term in the expansion of the external potential described in Appendix A.
Hereafter, we drop the ˜symbol to avoid cluttering.With these substitutions, we have   = − and   = .Substituting Eq. ( 31) into (9) we have: where we have introduced the sound speed of the unperturbed medium: Eq. ( 36) is valid for  ≥ 1 (including equality).We also define This is the angular frequency with which each mode appears to rotate, as can be understood by noting that In this paper, we will always take Ω p to be the same as the pattern speed of the bar described in Appendix A, since only modes at this frequency can be excited by the external potential in the linear approximation.Substituting ( 31)-( 35) into ( 26)-( 28) we obtain: Isolating  1 and  1 from ( 41) and ( 42) we find: where we have defined The points where  = 0 define the Lindblad resonances,2 while the point where Ω = Ω p defines the Corotation resonance.Now we can substitute ( 43) and ( 44) into (40) and use (36) to eliminate  1 to obtain an equation in the variable ℎ 1 : where Eq. ( 47) coincides with Eq. ( 13) of Goldreich & Tremaine (1979).
In order to eliminate the first order derivative from Eq. ( 47), it is convenient to define a new variable  1 such that Substituting Eq. ( 52) into Eq.( 47), one finds  53) for a uniform disc (full black line) and a truncated disc (dashed black line) in the case  s = 10 km s −1 .
In the WKB approximation, this represents the wavenumber versus radius of free density waves that rotate with the same pattern speed of the bar (see Sect. 4.3.2).The cyan lines compare with the wavenumber given by the Lin-Shu dispersion relation (59).Middle: the forcing term  () in Eq. ( 53).
Bottom: the uniform ( 0 = 1) and truncated disc (Eq.20 with  edge = 0.6 kpc and Δ = 0.03 kpc) density profiles assumed in this figure.The red vertical dashed line marks the ILR.The green vertical dashed line marks  ★ , which is defined as the radius where  () = 0 (see Sect. 4.3.1). where Eq. ( 53) is the fundamental equation that governs linear modes in the disc.It is similar to Eq. (B10) in the toy problem in Appendix B, but is more complicated because  is not constant.To follow the calculations in the following section more easily, it is useful to note that Eq. ( 53) is equivalent to that of a forced harmonic oscillator,   +  2 () = (), where  replaces ,  is the mass,  () is a timedependent spring constant, and () is a time-dependent external force.

Analysis of Eq. (53)
Equation ( 53) describes the dynamics of the most general linear perturbation in the presence of an external potential.To calculate the amplitude of waves excited by the bar potential we need to solve this equation with appropriate boundary conditions.Unfortunately, no general analytic solution is available, so we need to resort to various approximations that are valid in different radial ranges.Figure 7 provides an overview of the various regimes that we analyse.53) apply."WKB is valid" denotes where the general solution of Eq. ( 53) is well approximated as the sum of the WKB solution ( 56) and the equilibrium solution   given by (60)."WKB not valid" denotes the region near the edge where Eq. ( 67) is more appropriate.The shaded "matching regions" denote where both solutions are simultaneously valid and we can apply the method of matched asymptotic expansions.The region within approximately one wavelength  from the ILR is where the analysis of Goldreich & Tremaine (1979) is appropriate."Stage 1" and "Stage 2" denote the regions correponding to the two stages in our picture of the formation of the rings described in Sect. 5.
This section is structured as follows.In Sect.4.3.1 we identify special points where the treatment of Eq. ( 53) require special care because the coefficient  either vanishes or diverges.In Sect.4.3.2we derive the WKB solution of the homogeneous equation associated with (53), and show that it is generally very accurate away from the special points and away from sharp edges (see Fig. 7).In Sect.4.3.3we derive a particular solution of the non-homogeneous (53) that is approximately valid when  s is sufficiently low and away from special points and sharp edges.In Sect.4.3.4we obtain exact numerical solutions of Eq. ( 53) in a few selected cases, to illustrate that truncated discs with sharp edges excite much stronger waves than uniform discs.In Sect.4.3.5 we present an approximated analytical solution of Eq. ( 53) that is valid near sharp edges and estimate the flux of angular momentum at sharp edges.

Special points
There are two types of points where Eq. ( 53) requires special attention: • Turning points.These are the points  ★ where  ( ★ ) = 0.At these points, the character of the solutions changes from oscillatory to exponential.
Figure 6 shows the coefficients of Eq. ( 53) for a uniform and a truncated disc profile in the case  s = 10 km s −1 .In the region of interest for this paper there is typically one turning point  ★ and one singular point at  ILR , with  ★ <  ILR .As we shall see below,  ★ is where the medium becomes absorbing and leading waves incident from  <  ★ are reflected into trailing waves that subsequently travel inwards.The position of  ★ depends on both the sound speed  s and the shape of the unperturbed density profile  0 ().In the limit  s → 0 we have  ★ →  ILR .However, for a finite value of the sound speed, the two points are distinct.

WKB solution of the homogeneous equation
Consider the homogeneous equation associated with Eq. ( 53), i.e. the equation obtained setting  = 0.This equation describes the propagation of "free" density waves on top of the unperturbed disc in the absence of the external bar potential.In this case, Eq. ( 53) is of the same form of Eq. (C1) and it can be solved in the WKB approximation.The general solution is given by Eq. (C9), which adapted to the notation used here reads: where  1 and  2 are arbitrary complex constants and  0 is an arbitrary radius.
The two terms on the right-hand side of Eq. ( 56) represent two waves travelling in opposite directions, analogously to the two sound waves that are possible in a uniform medium at a given frequency (see Appendix B).The quantity  is the wavenumber, which varies with radius.When  2 > 0, the solution (56) has oscillatory character and waves can travel, while when  2 < 0 it has exponential character and the medium is absorbing.Thus, as can be seen from Fig. 6, travelling waves can exist only at  <  ★ .Equation ( 54) implicitly contains , and therefore for fixed  this expression can be seen as a dispersion relation  =  ().
The direction of propagation of the waves can be understood from the group velocity.In Appendix D we calculate the group velocity of the WKB waves and we find that trailing waves ( 1 ≠ 0 and  2 = 0) propagate inwards, while leading waves ( 1 = 0 and  2 ≠ 0) propagate outwards.
The angular momentum flux associated with the WKB waves ( 56) is calculated in Appendix E and is given by Eq. (E10): Since  1 and  2 are constant for a given WKB wave, this equation shows that the flux of angular momentum is constant as a function of .It can be shown that the angular momentum flux corresponds to the adiabatic invariant associated with the general WKB solution (C9).Equation ( 57) also shows that the trailing wave has   > 0, while the leading wave has   < 0. Thus, trailing (leading) wave packets remove (increase) the amount of angular momentum in the region where they travel.What is the range of validity of the WKB approximation?The WKB approximation is expected to work well when the following parameter is small (see Eq. C3): Figure 8 shows that the WKB approximation works exceptionally well at  <  ★ , but breaks down near  =  ★ .The WKB approximation will also fail near sharp edges, because d/d becomes large (e.g.Fig. 6).
The WKB approximation used here is not completely equivalent to the more well-known Lin-Shu approximation.The Lin-Shu dispersion relation in the absence of self-gravity ( = 0) is given by (Eq.6.55 of Binney & Tremaine 2008): where  is given by Eq. ( 45).The top panel in Fig. 6 53) and the WKB approximation.The full black line shows the solution obtained numerically integrating Eq. ( 53) from  = 0.1 kpc with initial conditions  1 = 1, d 1 /d = 0 (full black line).The cyan dashed line shows the WKB approximation (56).We assumed  s = 10 km s −1 and constant unperturbed density  0 () = 1.Bottom: the "small" parameter of the WKB approximation (Eq.58).At  ★ it diverges.The WKB approximation works well only at  <  ★ .
given by Eq. ( 54).The two are similar at  <  ★ , but differ considerably around  ★ and  ILR .In particular, in the Lin-Shu approximation the turning point (which is the point where waves are absorbed) coincides with the ILR, while it is at a smaller radius ( ★ ) according to Eq. ( 54).This is because the Lin-Shu dispersion relation assumes very small sound speed, while the dispersion relation (54) takes into account the effect of finite sound speed.Indeed, in the limit of vanishing sound speed we recover the Lin-Shu dispersion relation from our dispersion relation (54).This can be shown by noting that in this limit  () ≃ −/ 2 s (Eq.49), while  2 ≪  and d/d ≪  (Eq.54).

Approximate non-oscillatory solution of the non-homogeneous equation
An approximate particular solution of Eq. ( 53) is: In the analogy with the harmonic oscillator, this solution corresponds to following the "instantaneous" equilibrium position of the oscillator as the external force slowly varies.It is expected to be valid when the "force" () varies slowly enough compared to the frequency of the harmonic oscillator .More formally, one can substitute  1 =   in Eq. ( 53), and impose that the first term on the left-hand side is small, i.e. d 2  1 /d 2 ≪  2  1 .This gives the following condition: This condition is verified in particular at low sound speed, since  → ∞ as  s → 0 at fixed  (see Eq. 54).Eq. ( 60) is equivalent to Eq. ( 15) of Goldreich & Tremaine (1979).It is a non-wave solution which is the analogue of the black dashed solution for the toy problem in Fig. B1 in Appendix B.

Excitation of density waves in uniform and truncated discs
We numerically solve Eq. ( 53) in a few selected cases.The goal is to calculate the amplitude of density waves excited by the bar potential in a region [ 1 ,  2 ] where  1 <  2 <  ★ , to illustrate how the amplitude depends on  s and on the unperturbed density profile  0 .A sharp edge might be present inside [ 1 ,  2 ].The appropriate boundary conditions are "radiation" boundary conditions (see Fig. 9).Causality requires that waves propagate away from the region [ 1 ,  2 ], because a solution in which waves come towards it would require a source of waves outside this region.Therefore, the correct solution to our problem contains waves propagating inwards at  =  1 , and outwards at  =  2 .These are schematically shown as the two waves  1 and  2 in Fig. 9.The goal is to calculate the amplitude of  1 and  2 .
Although the numerical solutions of Eq. ( 53) described in these section are exact, to impose the boundary conditions we need to use the results of the WKB analysis as we need to identify the direction of propagation of the waves.We proceed as follows.Let  1 be an exact solution of Eq. ( 53) that satisfies the radiation boundary condition.We assume that at  1 and  2 the conditions ( 58) and ( 61) are valid, so in a neighbourhood of these points we can decompose the solution as the sum of the WKB solution (56) and of the "equilibrium" solution (60).Therefore in a neighbourhood of  1 we can write and in a neighbourhood of  2 where  1 and  2 are trailing and leading WKB waves respectively (Eq.56): To find  1 and  2 , we use the shooting method.We start from  =  1 with an initial guess for  1 (which is a complex number, so the guess involves two real numbers) and initial conditions given by Eq (62), and integrate until  2 .At  2 we decompose  1 −   into its WKB components.This decomposition is unique and can be found by equating  1 −   and d( 1 −   )/d with Eq. ( 56) and its derivative and solving the resulting algebraic system of two equations in the two unknowns that give the amplitude of the two waves.We then vary the initial guess for  1 until the solution at  2 only contains an outgoing wave.The amplitude of the latter gives  2 .
Figure 10 shows the result of this procedure applied to [ 1 ,  2 ] = [0.1,1.0] for a uniform (left) and truncated disc profile (right), and for two different values of the sound speed.The truncated disc profile is chosen so that the width of the edge is comparable to the wavelength ( = 2/) at the edge.The edge cannot be much smaller than this without violating the Rayleigh criterion (Sect.4.1).The top panel shows the solution  1 , the middle panel the corresponding density profiles, and the bottom panel the flux of angular momentum associated with the waves  1 and  2 obtained performing the WKB decomposition as a function of .The amplitude of the excited waves is given by the oscillations around the equilibrium solution   (dashed line).
The figure illustrates the following points: • Waves excited in uniform discs are weak (right panels).The amplitude is essentially zero at  s = 1 km s −1 , while it is small but visible at  s = 10 km s −1 .This reflects the fact that the approximate solution   is very accurate at low sound speed, but is less accurate when the sound speed is slightly larger.Physically, the reason why stronger waves are excited for larger  s is that the wavelength  of WKB increases with  s (Fig. 11), so that waves couple more effectively to the forcing term  which varies on large scales.As we shall see in Sect.5, waves excited in uniform discs are too weak to remove the angular momentum necessary to open the gap.
• Waves excited at the edge of a truncated disc are strong (left panels).The total density becomes negative near the edge ( 1 +  0 < 0), indicating that the linear approximation breaks down.The waves become highly non-linear and in reality they will develop shocks very quickly near the edge, as indeed seen in the simulations of Sect. 2. The amplitude of the waves is similar at  s = 1 km s −1 and  s = 10 km s −1 , but the flux of angular momentum is much larger for  s = 10 km s −1 .This will be explained by Eq. ( 72) below.
There is a simple explanation for why strong waves are excited at sharp edges.When the background density  0 () varies rapidly, such as at the edge of a disc, the forcing term () on the right-hand side of Eq. ( 53) will have a localised bump on the same scale (see dashed line in the middle panel of Fig. 6).This localised bump acts like an impulsive force.In the analogy with the harmonic oscillator, this force will impart a finite amount of "momentum" equal to the integral of the force.The amplitude of the resulting oscillations gives the amplitude of the waves excited at the edge.

Analytical estimate of the waves excited at the edge of a truncated disc
In this section, we derive an analytical estimate for the amplitude of waves excited at a sharp edges.
Consider an edge at  edge of width  out −  in = Δ, where  in and  out >  in are the two extremities of the region over which the edge extends (see Fig. 7).The shape of the edge can be arbitrary.The edge is assumed to be thin but not too thin, otherwise the unperturbed density profile would violate the Rayleigh criterion (18) and become unstable.In practice, considering the Lin-Shu approximation (Eq.59) this means that the edge should not be thinner than approximately one wavelength,  = 2/.
Away from the edge and from turning and singular points, , i.e. at  <  in and  out <  <  ★ , the general solution of Eq. ( 53) can be approximated as the sum of the WKB solution (56) and of the particular solution (60).We assume that waves are excited only near the edge, i.e. at  in <  <  out .We impose radiation boundary conditions, so that at  <  in we have only the trailing wave and at  >  out only the leading wave.Therefore outside the edge we Waves excited by the bar potential in the region [ 1 ,  2 ] = [0.1,1.0] kpc in uniform and truncated discs in the linear approximation.Waves excited in uniform discs are small, while waves excited at the sharp edge of a truncated disc are much larger.For each of the four cases shown, the three panels from top to bottom display the following.Top: the solution of Eq. ( 53) with radiation boundary conditions at  = 0.1 kpc and  = 1.0 kpc (full black line), and the approximate "instantaneous equilibrium" solution   given by Eq. (60) (dashed black line).The oscillations of  1 around   give the amplitude of the excited waves.Middle: the thin full black line shows the unperturbed density profile  0 , which can be either a uniform disc ( 0 = 1) or a truncated disc given by Eq. ( 20) with  edge = 0.6 kpc and Δ = 0.03 kpc (for  s = 10 km s −1 ) or Δ = 0.003 kpc (for  s = 1 km s −1 ).The dashed and full thick black lines show the total density (unperturbed + perturbation) that corresponds to the solutions  1 and   shown in the top panel.Bottom: the flux of angular momentum associated to the two WKB waves into which  1 −   can be decomposed.See Section 4.3.4 for more details. write: The constants  in and  out will be determined by solving Eq. ( 53) near the edge and matching the two solutions.
Near the edge, the forcing term  varies rapidly, violating condition (61), and the equilibrium solution (60) fails (dashed line in the middle panel of Fig. 6).To solve Eq. ( 53) near the edge, we proceed as follows.We assume that  is approximately constant across the edge, i.e. in the range  in <  <  out .This is justified by our assumption that the edge is relatively sharp (see also the black dashed line in the top panel of Fig. 6).Under this assumption, Eq. ( 53) can be solved using the method of variation of parameters.The general solution is: where  ≃  ( in ) ≃  ( out ).The constants  1 and  2 are determined by the condition that the solution contains only waves travelling inwards at radii  <  in , and waves travelling outwards at radii  out < .These calculations are reported in Appendix F1.We find: where  out = ( out ) and  in = ( in ).
Both Eq. ( 66) and Eq. ( 67) are valid solutions of Eq. ( 53) in a neighbourhood of  in and in a neighbourhood of  out (shaded regions in Fig. 7).Matching these two solutions gives (see Appendix F1): The coefficients  in and  out give the amplitude of density waves excited at the edge.It can be shown that the absolute values | in | and | out | are independent of  0 , as they should since the angular momentum flux at the edge should be independent of this arbitrary radius.
The rapid variation of  near the edge is what generates density waves.The coupling between the forcing term  and the density waves is expected to be maximum when the scale-length over which  varies is comparable to the wavelength of the waves  = 2/ (similarly to the toy problem in Appendix B), i.e. when Δ ≃ .
We can use Eqs.( 57) and ( 70) to calculate the angular momentum flux carried by the waves excited at the edge.To obtain a closed formula it is necessary to make some further assumptions on the edge.If the edge of the disc is marginally stable to the Rayleigh criterion ( 18), one has | out − in | ∼  s /Ω ∼ 1/.This is the sharpest edge that can be constructed without making the unperturbed density distribution unstable.Then the exponential exp(− ) in the integral of Eq. ( 70) is nearly constant.As shown in Appendix F2, in this case Eq. ( 70) reduces to: Note that this is essentially the impulse approximation, i.e. we have assumed that the force  gives an instantaneous "kick" at  =  edge .Using Eq. (E10), the flux of angular momentum of waves excited at a sharp edge is then Eq. ( 72) is correct when the distance of the edge from the inner Lindblad resonance is larger than approximately one wavelength, i.e.

THE FORMATION OF NUCLEAR RINGS
We are now ready to put everything together and describe our picture of the formation of nuclear rings.For simplicity, we illustrate our scenario by starting from a uniform density distribution that extends from  = 0 to  ≫  ILR .This is essentially the same situation as in simulation 04_Large shown (Fig. 5 and Sect.2.3).The formation of the ring can be schematically divided into two stages, which depend on the distance of the edge of the gas disc from the ILR.The regions corresponding to the two stages are marked in Fig. 7.The simulations 01-05 shown in Fig. 1 only include the second stage.

First stage (|𝑅
In the first stage, a trailing spiral wave is excited near the ILR by the external bar potential.This is the regime analysed by Goldreich & Tremaine (1979).The wave travels inwards but for realistic strengths of the bar potential it very quickly becomes non-linear and develops into a shock.The wave then dissipates, depositing its (negative) angular momentum into the gas disc (i.e., removing angular momentum from the gas disc). 3This reduces the angular momentum of the disc, causing the gas to move inward.A gap opens around the ILR.
The width of the gap opened in the first stage is the range of validity of the calculations of Goldreich & Tremaine (1979), which is approximately one wavelength, i.e. | edge −  ILR | ∼ , where  = 2/.Using the Lin-Shu approximation (Eq.59) and approximating  ≃ ( −  ILR )(d/d) (recall that  = 0 at the resonance) we have  ∼  s /| | 1/2 ∼  s /|( edge −  ILR )(d/d)| 1/2 and therefore The size of the gap is therefore much smaller than the radius of the resonance, and increases for increasing sound speed.
The velocity at which the edge of the gap moves can be estimated by dividing the flux of angular momentum of the waves,  A , by the amount of angular momentum per unit radius in the unperturbed disc, 2 0  3 Ω: The angular momentum flux  A during the first stage can be estimated using Eq. ( 46) of Goldreich & Tremaine (1979), bearing in mind that these calculations are valid in the linear approximation and should not be expected to be too accurate for the highly non-linear waves excited by a strong bar potential considered here.Taking into account that  = 2, we find Inserting the numbers of our gravitational potential (Appendix A) into Eq.( 74) we obtain d edge /d ≃ 6 km s −1 .The duration of the first stage can be estimated by dividing the size of the gap by the velocity of the edge.Using | edge −  ILR | ∼ ( s / 0 ) 2/3  ILR ,  0 = 220 km s −1 ,  ILR = 1.6 kpc and the value of d edge /d found above we obtain: This time is relatively short compared to the total timescales involved (see Fig. 3).In reality, the evolution is likely to be even faster than Eq. ( 75) suggests because of the non-linearity of the process.The evolution of the gap after the first stage and its final size are determined by the second stage.

Second stage (|𝑅
At the beginning of the second stage there is a gap around the ILR, and the distance between the inner edge of the gap and the ILR is approximately one wavelength.Since the width of the edge at this point can be at most one wavelength (because the edge tail cannot extend beyond the ILR), the edge is "sharp" by definition and strong waves will be excited at its location according to the analysis in Sect.4.3.Similarly to the waves excited near the ILR in the first stage, the waves excited near the edge will become quickly non linear and dissipate, removing the angular momentum from the gas disc and causing the edge to move inwards.Gas will accumulate at the edge, forming a ring.This process can be seen in action in Figs. 1 and 2. The first figure shows trailing waves excited by the bar potential (see for example  = 157 Myr).The pitch angle of these waves is in good agreement with that predicted by the WKB analysis of Sect.4.3.2,indicating that these are indeed trailing waves of the same type studied in the linear analysis.Fig. 2 shows that the amplitude of the waves decreases inwards, contrary to what would be predicted in the linear approximation (e.g.Fig. 10).This is because when the waves become strongly non-linear and develop shocks, they quickly dissipate and decrease their amplitude.This dissipation is what allows the wave to deposit their angular momentum into the gas disc.
The speed at which the edge moves during the second stage can be estimated using Eq. ( 73), where we use Eq. ( 72) to estimate the flux of angular momentum   of waves excited at sharp edges.Using the Lin-Shu approximation (Eq.59) to write  ≃ | | 1/2 / s , and  = 2, we find: where the factor of 2 takes into account that the outward-travelling leading wave excited at the edge will be reflected at  =  ★ into an inward-travelling trailing wave.Note that Eq. ( 74) and ( 76) only differ for the factor in the first round parentheses on the right-handsides, and this factor coincides in the two equations at a distance of approximately one wavelength from the ILR.This is the point where we transition from the analysis of Goldreich & Tremaine (1979) to the analysis in Sect.4.3.5, and from the first to the second stage.
Fig. 3 compares the size of the ring as a function of time predicted by Eq. ( 76) to that measured in the numerical experiments of Sect.2.2.We find that the equation captures the general trends in the figure, including the fact that the edge moves faster for larger sound speed, but it tends to underestimate the speed at which it moves, especially at large sound speed.That the analytic prediction is not quantitatively accurate is not surprising considering that Eq. ( 76) is derived in the linear approximation, but the waves excited at the edge are strongly non-linear (Figs. 2 and 10).The flux of angular momentum generated in the case of a uniform disc is too small to move the edge significantly over the course of several Gyr (Fig. 10).
When does the edge stop moving?The process above continues until waves can be effectively excited at the edge, which happens when both of the following conditions are satisfied: (i) the edge is "sharp", i.e. the edge width is smaller than a few times the wavelength of density waves  = 2/; (ii) the gravitational potential Φ 1 is sufficiently strong.The distance between the edge and the ILR poses an upper limit to the width of the edge since the edge cannot cross the ILR, Δ < | ILR −  edge |.4 Therefore, when the edge is not sufficiently far from the ILR, it must be sharp.In particular, we can expect the edge to keep moving until it is located a few wavelengths away from the ILR.Since  increases linearly with  s (Eq.59 and Fig. 11), we expect the edge to move farther at larger sound speed, which explains why the ring radius depends on the sound speed.
Predicting exactly where the edge will stop, and therefore the final radius of the ring, is a difficult task.The process is highly non-linear, and the unperturbed density profile changes in a way that cannot be calculated in the linear approximation.Empirically, we find from the numerical experiments in Sect. 2 that for our assumed gravitational potential the ring stops when one can fit approximately 7 wavelengths  between  edge and  ILR .For weaker barred potential, the edge might stop sooner if Φ 1 is too small to generate sufficient flux of angular momentum at the edge.
Finally, we note that our theory satisfies all the 5 conditions that we laid out in Sect.3.1.Conditions 1-3 are satisfied because clearly the radius of the ring depends on the rotation curve, on the nonaxisymmetric part of the gravitational potential, and on the pattern speed of the bar which sets the location of the ILR.All these dependencies are also evident in Eq. ( 76).Condition 4 is satisfied because the radius of the ring depends on the sound speed of the gas in two ways: first because the speed at which the edge moves away from the ILR increases as a function of  s (see Eq. 76), and second because the final ring size is determined by the condition that the edge should be a few wavelengths away from the ILR, which results in smaller rings at larger sound speed since as discussed above the wavelength increases with the sound speed.Condition 5 is satisfied since the excitation of density waves at the edge is a local process.

Comparison with the works of Goldreich & Tremaine
During the late 70's and early 80's, Peter Goldreich and Scott Tremaine published a series of papers in which they studied the dynamics of planetary rings.The calculations and physical processes studied in these works have much in common with those presented in the present paper.Here we highlight the main similarities and differences between them and the present paper.Goldreich & Tremaine (1978) (hereafter GT78) developed a picture for the formation of the Cassini division in Saturn's ring that has several similarities with our picture for the formation of nuclear rings described in Sect. 5.In both cases: (i) a gap opens near the 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Lindblad resonance due to waves excited at the resonance; (ii) subsequent excitation of waves at the edge of the gap continues to widen the gap.The main differences are: (i) The bar potential considered here is a perturbation many orders of magnitude stronger than the one from Saturn's satellite Mimas; (ii) the sound speed is negligibly small in Saturn's problem, while the effects of finite sound speed are important in our problem (Sect.5); (iii) self-gravity is negligible for our case, but it is not negligible in Saturn's problem.In particular, gravity is the main means of transport of angular momentum in Saturn's problem, while advective transport through pressure is the main mechanism for transport in our problem.Goldreich & Tremaine (1979) (hereafter GT79) studied the excitation of density waves in a differentially rotating gas disc by a rigidly rotating external potential.Their calculations have similarities with those presented in Sect.4, but there are three key differences: (i) GT79 assume that  s → 0, while we take into account the effects of finite sound speed; (ii) GT79 assume that  0 varies slowly; (iii) GT79 included the self-gravity of the gas disc, which we have neglected.The combination of (i) and (ii) is why GT79 find that waves can be excited only at the resonances, while for example in Sect.4.3.4we find waves excited away from the resonance.The physical explanation is the following.The external potential can couple effectively to density waves only when the wavelength of WKB waves ( = 2/) is comparable to the typical scalelength over which the forcing term  in Eq. ( 53) varies, i.e. when  ∼ /(d/d).The quantity /(d/d) is determined by the external potential Φ 1 and by the unperturbed density distribution  0 , and is therefore typically very large unless there are sharp edges in  0 .In the limit of vanishing sound speed, /(d/d) ≫  everywhere except near the turning point  ★ at which  → ∞ (see Fig. 11).In the limit  s → 0 the turning point merges with the ILR (Sect.4.3.1)and  is small everywhere else.Thus, in this limit waves can be excited only at the resonance.For finite sound speed instead  can become large and comparable to /(d/d) away from the resonance and waves can be excited (see Fig. 11).As we have seen in Sect.5, the effects of finite sound are important in the formation of nuclear rings.

Relation to 𝑥 2 orbits
Several works have suggested a connection between nuclear ring and  2 orbits (e.g.Regan & Teuben 2003;Li et al. 2015;Sormani et al. 2018b).The  2 orbits are a family of non-circular closed orbits that can exist in the central regions of a bar potential, and are elongated in the direction perpendicular to the major axis of the bar (e.g.Contopoulos & Grosbol 1989; Athanassoula 1992a).Figure 12 illustrates the relation between these orbits and the present paper.The streamlines of the "equilibrium" solution Eq. ( 60) are very similar to closed  2 orbits in the same bar potential.Therefore, the WKB waves excited by the bar potential that we studied in Sect.4.3 travel on top of an  2 gas disc.Our picture for the formation of the rings is therefore consistent with the idea that gas in nuclear rings flows on  2 orbits.

Relation with the resonant theory
Our theory is somewhat the "opposite" of the resonant theory, which states that the ring forms at the ILR (Combes 1988(Combes , 1996;;Buta & Combes 1996).In our theory the gas is pushed away from the ILR rather than accumulating at it.Our theory is more consistent with the fact that the rings are typically inside the ILR in simulations (e.g.Englmaier & Gerhard 1997;Patsis & Athanassoula 2000;Kim et al. 2012;Sormani et al. 2015a;Li et al. 2015) and with observations that show that for example in the Milky Way the radius of the nuclear ring is ≃ 100-200pc while the ILR is at  > 500 pc (Henshaw et al. 2023).A key difference between our theory and all previous theories, including the resonant theory, is that we can explain the puzzling dependence of nuclear ring size on the sound speed seen in simulations.

Brief considerations on magnetic fields and turbulent pressure
The mechanism for the formation of rings described in Sect. 5 relies on waves propagated through pressure.We would therefore expect that adding magnetic fields, which create magnetic pressure in the gas, could have a similar effect as increasing the sound speed, and would therefore lead to smaller rings.Turbulent pressure seems to have a smaller effect than "real" microscopic pressure on the size of nuclear rings.Salas et al. (2020) performed some numerical experiments with external turbulence driving.Their Figs. 1 and 2 show that turbulence driving changes the size of the nuclear ring by a smaller amount than an increase in the sound speed when the injected turbulent energy is comparable to the corresponding increase in thermal energy.We attribute this to the fact that, due to the presence of inelastic collisions, in a gas with turbulent pressure sound waves do not propagate as efficiently as in a gas that has the same amount of microscopic pressure.

On the assumption of an isothermal equation of state
Throughout this work, we have assumed an isothermal equation of state.This follows a tradition of works adopting the isothermal prescription to study the dynamics of the ISM on galactic scales and in nuclear rings (e.g.Roberts 1969;Cowie 1980;Athanassoula 1992b;Englmaier & Gerhard 1997;Fux 1999;Maciejewski 2004;Kim et al. 2012;Sormani et al. 2015a;Fragkoudi et al. 2017;Li et al. 2022).However, the real ISM is multi-phase, turbulent and highly inhomogeneous.It is therefore natural to ask whether the isothermal prescription can capture the basic mechanism for the formation of nuclear rings.
Numerical simulations that include a multi-phase medium via thermal instabilities (Sormani et al. 2018a) as well as gas self-gravity, star formation and supernova feedback (Armillotta et al. 2019;Tress et al. 2020) show that the large-scale morphology of nuclear rings is only moderately affected by the presence of this additional physics.The main differences are observed at small scales (smaller or comparable to the width of the rings), where gas condenses into molecular clouds and collapses to make star formation.The large-scale properties of the ring, such as its radius and width, can be often mimicked by using an 'effective' isothermal sound speed.For example, we found that the morphology, width and radius of the nuclear ring in the simulations of Sormani et al. (2018a), which include a non-equilibrium chemical network that produces a two-phase medium via the thermal instability, are very similar to those obtained by replacing the non-equilibrium network and the associated cooling function with an isothermal equation of state with a low sound speed of  s ≃ 1 km s −1 .This low value is because the gas in the ring is very cold in these simulations, as they did not include any sources of heating or turbulence such as stellar feedback.When star formation and stellar feedback are added to the simulations (e.g.Armillotta et al. 2019;Tress et al. 2020), they heat up the gas and generate turbulent pressure, and the morphology of the rings can still be crudely mimicked by raising the isothermal sound speed (although as noted in Sect.6.4 by less than the turbulent velocity dispersion, which would be the naive way of doing it).Similarly, the effects of magnetic fields can be crudely mimicked by increasing the sound speed by summing in quadrature the typical Alfvén speed.
In conclusion, the isothermal prescription should be viewed as an 'effective' equation of state that takes into account in a phenomenological way the additional physics via a single parameter that can be easily controlled.Ultimately, the key property that needs to be captured in this approach is the ability of the medium to propagate waves through pressure.Thus, the sound speed does not correspond to the actual kinetic temperature of the gas, but to an 'effective' temperature that takes into account in a crude way averaging over different phases, turbulent motions on unresolved scales, and other effects such as magnetic pressure.Reassuringly, numerical simulations suggest that the basic mechanism for the formation of the ring is well captured using this approach.

CONCLUSION
We have used both hydrodynamical simulations and analytical calculations of linear disc dynamics to construct a new theory for the formation of nuclear rings in barred galaxies.According to this theory, nuclear rings are an accumulation of gas at the inner edge of a gap that forms around the Inner Lindblad Resonance (ILR) of a bar potential.The gap initially opens because the bar potential excites strong trailing waves around the ILR, which remove angular momentum from the gas disc and push the gas inwards.The gap then continues to widen because the bar potential excites trailing waves at the inner edge of the gap, until the edge stops at a distance of several wavelengths from the ILR.The gas accumulates at the inner edge of the gap, forming a ring.The speed at which the gap edge moves and its final distance from the ILR, which determine the radius of the nuclear ring, depend on the gas sound speed through the dispersion relation.
Our theory has much in common with the picture for the formation of the Cassini gap in Saturn's ring proposed by Goldreich & Tremaine (1978).The most important differences are that (i) the effects of finite sound speed are important in our problem, while the sound speed can be assumed to be vanishingly small in the planetary problem; (ii) we have neglected the effects of self-gravity, which are typically less important in the nuclear ring problem, but cannot be neglected in the planetary rings problem.The equations of motion of this system are the same as Eqs.( 1), ( 2) and (3) where the gradient is replaced by d/d since the problem is one dimensional.We linearise these equations around the background state by writing (, ) =  0 +  1 (, ) and (, ) =  1 (, ) and keeping only the first-order terms in the quantities with subscript 1.We obtain: Without loss of generality, we can write all variables as Where F are complex quantities.We use complex notation for mathematical convenience, but it is understood that the physical quantities are given by the real part.Substituting all perturbation variables in the form (B3) into (B1) and (B2), and omitting the symbol ˜hereafter for simplicity of notation, we obtain Isolating  1 from (B5) and introducing the variable Substituting (B6) into (B4) we obtain an ODE for  1 : where Eq. ( B7) is the equation of a forced harmonic oscillator.Now consider an oscillating Gaussian potential of the form: where Φ 1 is the strength of the potential,  0 is the width of the Gaussian perturbation,  0 is the oscillation frequency.Since in the linear approximation there is no coupling between modes at different frequencies, only modes with frequency  =  0 will be excited by this potential.Hence we assume  =  0 hereafter.Introducing the dimensionless coordinate  = / 0 and using (B9), Eq. (B7) becomes: where and we have introduced the following dimensionless parameters: The parameter  is the inverse of the sound speed, normalised with the typical scale-length and frequency of the problem.The parameter  is the normalised strength of the external potential.

B2 Analytical solution
The general solution of Eq. ( B10) is where  1 and  2 are arbitrary constants and Here, erf () = 2 √  ∫  0  − 2 d is the error function, which is defined for complex argument  (to evaluate the integral, you can choose any integration path in the complex plane that leads to ).Note however that the functions  and  are real because the erf function has the properties erf () = erf () and erf (−) = erf (), where the bar denotes the complex conjugate.
Figure B1 plots the function  (, ) for various values of .In the limit  → ±∞ we have that erf( + /2) → ±1 for any fixed , so compared to the others.Then Eq. (C6) becomes: This equation can be integrated and the solution is Plugging (C8) into (C4) we find that the general solution of Eq. (C1) in the WKB approximation is: where  1 and  2 are arbitrary constants.We can also write the following approximate expression for the derivative by neglecting the small terms  and : Note that the total energy of a simple harmonic oscillator is not in general conserved when () is not constant.However, if we calculate  using the approximate solution (C9) and (C10), we find that the following quantity is constant: ) This means that the amplitude of the oscillation becomes a function of .If we increase  slowly then we slowly decrease it to its original value, at the end of the process the amplitude will be the same as it was at the start.It is easy to see that this is violated if () does not change slowly (think for example of abruptly changing  when the system passes through  = 0: in this case the energy does not change, but  does).The quantity  is an example of an adiabatic invariant (see for example Arnold 1978 andLandau &Lifshitz 1976 for more on adiabatic invariants).

APPENDIX D: GROUP VELOCITY
In this Appendix, we calculate the group velocity of the WKB waves.Consider a WKB solution (56) with  1 ≠ 0 and  2 = 0.This is of the form: where Equation ( D1) is of the same form of Equation (1) of Toomre (1969) or (1.26) of Whitham (1974).The analysis in these references shows that the group velocity, i.e. the velocity at which a wave packet travels, can be defined by isolating  from the dispersion relation ( 54) and then taking the derivative with respect to : Figure D1 shows the group velocity for the case  s = 10 km s −1 .The group velocity of solutions with  1 ≠ 0 and  2 = 0 is negative, meaning that these waves travel inward, while waves with  1 ≠ 0 and  2 = 0 travel outwards.The group velocity of the two types of waves has the same magnitude but different sign.The group velocity loses meaning and becomes imaginary at  >  ★ , when the medium becomes absorbing.
The cyan line compares our group velocity with that obtained from the Lin-Shu dispersion relation, which is given by (see Equation 20of Goldreich & Tremaine 1979): where  Lin−Shu is given by Eq. ( 59).The two group velocities are similar away from  ★ .

APPENDIX E: ANGULAR MOMENTUM TRANSPORT
An equation for the angular momentum transport in a fluid disc can be obtained from Eq. ( 2).Multiplying the azimuthal component of this equation by , using standard cylindrical coordinates (, , ) and rearranging gives: where The quantity   is the angular momentum per unit volume, while F  is the flux of angular momentum, which is the sum of contributions due to bulk motions of the gas and pressure forces.The term Φ/ is a source term representing the changes in angular momentum due to torques from the external potential.When Φ/ = 0, the total angular momentum of the system is conserved.Indeed, the only agent that can change the total angular momentum in our problem is the external bar potential.Integrating Eq. (E1) over the volume  of a cylinder of radius  0 and using the divergence theorem, 5 we obtain the following equation for the rate of change of the total angular momentum contained within the cylinder: where is the total  angular momentum contained inside the cylinder, and are the fluxes of angular momentum in and out of the cylinder.Eq. (E4) states that the change in the total angular momentum of the contained within the cylinder is the sum of two contributions:   , the angular momentum flux due to advection, and  Φ , the gravitational torques from the external bar potential.
The quantity   appears in Goldreich & Tremaine (1979) as their Eq.( 26). A > 0 means that material inside the cylinder is losing angular momentum.Notice that even in a steady-state, in which single fluid elements neither gain nor lose angular momentum on average, it is nevertheless possible that  A ≠ 0. This can happen if fluid elements carry more angular momentum on their outward journey (as they are exiting the cylinder) than on their return.This type of transport has been named lorry transport by Lynden-Bell & Kalnajs (1972), who explained how fluid elements can "transport angular momentum just as a system of lorries can transport coal without accumulating a growing store on the lorries themselves".This is similar to a plane sound wave transporting linear momentum in a steady-state situation.
5 The divergence theorem states that for any vector-valued function F(x): where  out = ( out ).Since the waves are travelling outwards at  =  out , the term proportional to    should vanish.This condition gives Eq. ( 68).
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .Figure 2 .
Figure 1.Surface density of simulations 01-05 (see Table1), illustrating the formation of nuclear rings for various values of the sound speed  s .The only difference between these simulations is the assumed  s .Time increases from left to right.Sound speed decreases from top to bottom (top row is the 05 simulation, bottom row is the 01 simulation).The black dashed circle indicates the ILR.The white dotted circle indicates the instantaneous ring radius according to the definition used in Fig.3.All panels are rotated so that the major axis of the bar potential (i.e. the  = 0 line in Eq.A1) coincides with the  axis.The sense of rotation is clockwise.Trailing density waves excited by the bar potential are visible (see in particular the second and third column from the left).The radius of the ring at the end of the simulation (rightmost column) strongly depends on the sound speed.Regions with densities  < 10 −2 are shown white.

Figure 3 .
Figure 3. Top: The full lines show the radius of the ring as a function of time in the simulations 01-05.The radius of the ring is calculated as  ring =√︁     , where   and   are the locations of the density maxima along the  and  axis respectively.The dashed lines show the prediction according to Eq. (76) obtained in the linear approximation (see Sect. 5.2).Bottom: The radius of the ring mediated over simulation time  = 1152-1252 Myr as a function of the sound speed.The radius strongly depends on the sound speed.The red dashed line indicates the inner boundary of the computational grid.

Figure 4 .
Figure 4. Surface density of the 04_Large (left and zoom-in panel) and 04 (right) simulation at the end of the simulation ( = 1252 Myr).The only difference between the two simulations is that in simulation 04_Large the initial gas disc extends to  max = 5.0 kpc, while in simulation 04 only to  disc = 1.2 kpc.Simulation 04 is the same as shown in the second row of Fig. 1.The dashed circle indicates the ILR.All panels are rotated so that the major axis of the bar potential (i.e. the  = 0 line in Eq.A1) coincides with the  axis.The sense of rotation is clockwise.Comparison between the two simulations shows that the ring always reaches the same final size, regardless of the larger-scale flow outside the ILR, demonstrating that the physical process determining the radius of the ring must be "local".

Figure 5 .
Figure 5. Axisymmetrised surface density ⟨⟩  =∫ (,  )d/(2  ) as a function of cylindrical radius  for the simulation 04_Large at three different times.An extensive gap opens around the ILR.The material that once was in the gap is transported inwards and accumulates at the inner edge of the gap, forming a nuclear ring.

Figure 6 .
Figure 6.Top: the coefficient  () in Eq. (53) for a uniform disc (full black line) and a truncated disc (dashed black line) in the case  s = 10 km s −1 .In the WKB approximation, this represents the wavenumber versus radius of free density waves that rotate with the same pattern speed of the bar (see Sect. 4.3.2).The cyan lines compare with the wavenumber given by the Lin-Shu dispersion relation (59).Middle: the forcing term  () in Eq. (53).Bottom: the uniform ( 0 = 1) and truncated disc (Eq.20 with  edge = 0.6 kpc and Δ = 0.03 kpc) density profiles assumed in this figure.The red vertical dashed line marks the ILR.The green vertical dashed line marks  ★ , which is defined as the radius where  () = 0 (see Sect. 4.3.1).

Figure 7 .
Figure 7. Schematic diagram of where the various approximate solutions of Eq. (53) apply."WKB is valid" denotes where the general solution of Eq. (53) is well approximated as the sum of the WKB solution (56) and the equilibrium solution   given by (60)."WKB not valid" denotes the region near the edge where Eq. (67) is more appropriate.The shaded "matching regions" denote where both solutions are simultaneously valid and we can apply the method of matched asymptotic expansions.The region within approximately one wavelength  from the ILR is where the analysis ofGoldreich & Tremaine (1979) is appropriate."Stage 1" and "Stage 2" denote the regions correponding to the two stages in our picture of the formation of the rings described in Sect. 5.

Figure 9 .
Figure 9. Schematic diagram of excitation of waves in the region [ 1 ,  2 ].  1 and  2 are the waves excited by the barred potential in this region.Arrows indicate the direction of propagation.See Sect.4.3.4 for more details.

Figure 12 .
Figure 12.Black full lines: Streamlines of the non-wave solution (60) in the case  s → 0 and constant  0 .Cyan dashed lines: ballistic closed  2 orbits in the barred potential described in Appendix A. The two are similar, showing that the density waves excited by the bar potential discussed in Sect. 4 are essentially perturbations travelling on top of an  2 disc.

Figure A1 .
Figure A1.Top: the circular velocity of our potential,  c = (dΦ 0 /d) 1/2 .Middle: the curves Ω and Ω ± /2, where Ω =  c / and  is the epyciclic frequency (Eq.46).The intersection of these curves with the horizontal line at Ω p gives the position of the resonances, indicated by vertical dashed lines.Bottom: the quadrupole (Eq.A3).

Figure D1 .
Figure D1.Group velocity of WKB waves for  = 10 km s −1 and constant unperturbed density  0 () = 1 (see Eq. D4).The cyan line shows the group velocity according to the Lin-Shu dispersion relation (D5).The ILR and  ★ are marked by vertical dashed lines.

Table 1 .
Summary of the simulations run in this paper.The parameters are defined in Sect. 2.
−   2 ∫  out  in ()   d + +  −   in 2 2    in −     out 2 2  −  out −We approximate Eq. (70) as follows.First, we neglect the terms proportional to  in and  out because  varies rapidly at radii  in <  <  out .We obtain()d =  1 +  2 +  3 +  4 ,Since | out / in − 1| ≪ 1 and far from the ILR the integrand is bounded, we have  1 ≃  4 ≃ 0.We calculate  2 and  3 below.The idea is to integrate by parts in order to isolate the integral of a bounded function.We have 2 = −2 in () −  d + + in () −  d .(F5) Second, the exponential exp(− ) is nearly constant as we have assumed | out −  in | ∼  ∼ 1/, so we can write | in | ≃ | out | ≃