Standing solitary waves as transitions to spiral structures in gravitationally unstable accretion disks

Astrophysical disks that are sufficiently cold and dense are linearly unstable to the formation of axisymmetric rings as a result of the disk's gravity. In practice, spiral structures are formed, which may in turn produce bound fragments. We study a nonlinear dynamical path that can explain the development of spirals in a local model of a gaseous disk on the subcritical side of the gravitational instability bifurcation. Axisymmetric equilibria can be radially periodic or localized, in the form of standing solitary waves. The solitary solutions have an energy slightly larger than a smooth disk. They are further unstable to non-axisymmetric perturbations with a wide range of azimuthal wavenumbers. The solitary waves may act as a pathway to spirals and fragmentation.


INTRODUCTION
Spirality is ubiquitous in astrophysics. The grand design spirals in galaxies and some circumstellar disks are believed to be global spiral density waves. Spirals can further regulate the formation of stars in galaxies (Roberts 1969) and collapse to form planets in circumstellar disks (see, e.g. Durisen et al. 2007;Deng et al. 2021). In this Letter we study the development of spiral structures through the gravitational instability (GI) of a massive gaseous disk orbiting in a central potential.
The physics of GI can readily be appreciated in a local patch of a thin gaseous disk (Goldreich & Lynden-Bell 1965). In this 2D local model, x and y correspond to the radial and azimuthal directions. Fluid orbiting the centre at angular velocity Ω(r)ê z is governed by the equation of motion where u is the velocity, D/Dt is the material derivative, Σ and P are the vertically integrated density and pressure, Φ is the gravitational potential due to gas selfgravity, and Φ t = −ΩSx 2 is the tidal potential. For uniform Σ and P , the steady solution u 0 = −Sxê y represents the basic orbital motion with shear rate S = hpdeng353@gmail.com gio10@cam.ac.uk −dΩ/d ln r, which equals 3Ω/2 for a Keplerian disk dominated by a central mass. Hereafter we consider only the deviation from this background shear flow, v = u − u 0 .
Axisymmetric perturbations with radial wavenumber k x oscillate at angular frequency ω given by where c is the sound speed. Long-wavelength (small k x ) perturbations are stabilized by the κ 2 = 2Ω(2Ω − S) term (squared epicyclic frequency) resulting from conservation of angular momentum, while short-wavelength perturbations are stabilized by gas pressure (Shu 2016). A band of intermediate-scale axisymmetric perturbations exists when Q = κc/πGΣ < 1 (Toomre 1964). The growth of non-axisymmetric spirals is less well understood and possibly linked to the overreflection of waves at the location that corotates with the spiral pattern. If somehow the reflected wave is deflected again toward the corotation region, a feedback loop can be established so that the wave grows exponentially (Mark 1976;Nakagawa & Sekiya 1992;Shu 2016). However, this mechanism only works efficiently for disks hovering on the brink of GI, and with enough radial structure to reflect the waves. In 3D numerical simulations of gaseous disks, spirals start to grow at Q ≈ 1.5 (Durisen et al. 2007) and ring structures are often observed as transitions to spirals (Mayer & Gawryszczak 2008;Hirose & Shi 2019;Deng et al. 2017). Inspired by these numerical simulations, we investigated nonlinear steady axisymmetric structures in the local model as transitions to spirality. We unexpectedly discovered a class of standing solitary waves in a nonlinear integro-differential equation describing the radial force balance. Self-gravity of the gas introduces nonlocality through a Hilbert transform, resembling the Benjamin-Ono equation which also admits solitons (Benjamin 1967;Ono 1975). The solitary waves are found to transition to spiral structures via a secondary instability. We solve for nonlinear axisymmetric structures in section 2, examine their energy and stability in sections 3 and 4, respectively, and draw conclusions in section 5.
For steady axisymmetric solutions (v x = 0), the PV constraint implies where σ = Σ/Σ 0 − 1 is the fractional density variation. Radial force balance reads and combines with (3) to give The dimensionless enthalpy perturbation reduces to w = ln(1 + σ) in the isothermal limit (Γ → 1). Φ is related to Σ through Poisson's equation in 3D, regarding the disk as razor-thin. This can be solved after taking a Fourier transform in x and y, with the resultΦ = − 2πG kΣ (7) for any wavenumber k = k 2 x + k 2 y = 0. In Fourier space, equation (5), when combined with (7), gives We adopt c 0 /κ as a natural unit of length. Radial force balance can be expressed in real space as where Q −1 = πGΣ 0 /κc 0 and H is the Hilbert transform (King 2009). The polytropic relation (6) makes the problem nonlinear. We look for periodic solutions with a basic wavenumber K in the form of cosine series, w n cos(nKx), (10) so that the coefficients are related by Asymptotic analysis is possible when the disk is nearly marginally stable, i.e., Q −1 and K are close to 1 with small increments δ(Q −1 ) and δK. For weakly nonlinear solutions, in which the left-hand side measures supercriticality with respect to the GI: it is positive when Q < 1 and K is sufficiently close to 1. For Γ < 5/3, which is expected, solutions exist on the subcritical side of the bifurcation with δ(Q −1 ) < 0, i.e., Q > 1. Using these approximate solutions as initial guesses, we solve equations (6) and (11) by Newton-Raphson iteration for a range of parameters Q −1 , K, and Γ. We utilize Fast Fourier Transform (FFT) algorithms with up to 5000 points to achieve good convergence when K is small. In Fig. 1, we plot the steady axisymmetric structures in various conditions. For a given (Q −1 , Γ) pair, the short-wavelength solutions (K > 1) possess a single strong peak which tends to be infinitely sharp at large K, whereas the long-wavelength (K < 1) solutions gradually approach a solitary wave. For example, in Fig. 1(a) the solitary wave is already established at K = 0.1, and we observe no change in the pattern by tripling the base wavelength, 2π/K. Hereafter the reported solitary waves for various (Q −1 , Γ), whose shapes are well established, are all calculated with a base wavenumber K = 0.08. A smaller Q −1 or a larger Γ leads to a more strongly peaked wave.
With a fixed wavenumber K = 1 and Γ < 5/3, we show the peak density perturbation of periodic ring structures as a function of Q −1 in Fig. 2. The rings become more peaked when Q −1 decreases from 1 until the solutions turn onto an upper branch at a critical Q −1 (0.712, 0.825, 0.900, 0.950 for Γ = 1.1, 1.2, 1.3, 1.4). The peak density continues to increase with Q −1 until the first trough in the periodic structure (similar to Fig. 1) reaches zero and the upper branch ceases.

ENERGY BUDGET
The energy per unit area of the axisymmetric structures, relative to the uniform disk, is where the first two terms combine the kinetic and tidal energies. Considering the isovortical constraint (3) and radial force balance (5), and after integration by parts, assuming that v y , σ, and Φ are either periodic or decaying in x, the total energy change due to the axisymmetric structure is So the mean energy change of periodic structures per unit area is We plot the energy change due to periodic ring structures and solitary waves in Fig. 3. For the K = 1 solutions in Fig. 2, the energy first increases as Q −1 decreases along the lower branch as shown in Fig. 3(a). This is expected for a subcritical bifurcation because the uniform disk is linearly stable and a local minimum of the energy. When the solution moves onto the upper branch, the energy turns over and falls with increasing Q −1 . The turnover points are stationary points of the energy (Burke & Knobloch 2006) and the upper branch solutions for Γ > 1.4 can have energy lower than the undisturbed smooth sheet. The lower-branch solutions are expected to be unstable. In Fig. 3(b), the energy change of the axisymmetric structure mostly increases with K except around the local minima near K = 1. In the long-wavelength limit (K → 0), the energy perturbation per unit area scales linearly with K, i.e., ∆E/K tends to a constant. This is a sign that the solution converges to a solitary wave with a fixed total energy change. In Fig. 3(c), the energy of solitary waves increases with Γ and Q.

NON-AXISYMMETRIC INSTABILITY
We investigate the stability of the periodic and solitary waves to non-axisymmetric perturbations (k y = 0).

The evolution of linear perturbations (primed variables) is governed by
Dv y Dt where D/Dt = ∂ t + (−Sx + v y )∂ y is the material derivative on the axisymmetric equilibrium. The perturbation can be expressed as a ladder of shearing waves (Goldreich & Lynden-Bell 1965), e.g., with k x,n = nK + k y St. We only consider isovortical perturbations with zero potential vorticity change, i.e., The density perturbation can thus be directly related to the velocity perturbations, leaving two coupled linear ordinary differential equations (ODEs) for each n: Here, the Fourier coefficients of v y and ∂ x v y for the axisymmetric equilibrium can be easily obtained from the cosine series. The Fourier coefficients of Σ can be related to the velocity perturbation coefficients by the isovortical condition. Φ n is readily obtained from Σ n by the Poisson equation, and W n can be calculated via a convolution between Σ n and (c 2 /Σ) n . The ODEs have t-dependent coefficients but the system is periodic in t, because the k x,n ladder shifts one place in n after the recurrence time T = K/k y S. Hence Floquet's method can be used to determine the growth rates of various perturbations with |n| not exceeding a truncation order N (Vanon & Ogilvie 2016). To that end, a 2(2N + 1) × 2(2N + 1) monodromy matrix was produced by applying 2(2N + 1) sets of initial conditions where all variables but one are set to zero, with a different nonzero variable in each case, to equation (21) and integrating it for one recurrence time. To ensure the v ±N entries are properly evaluated, we extended the Fourier series, employing 2(N +N 0 )+1 wavenumbers for each variable in equation (21). Here N 0 is the number of nonzero (> 10 −14 ) modes in the equilibrium axisymmetric structures so that the convolutions in equation (21) are properly handled for all modes (−N → N ) in the monodromy matrix. We used up to N = 800 wavenumbers for some strongly peaked solitary waves with N 0 up to 700. FFT and inverse FFT were utilized for the convolutions.
The growth rate s is determined by the eigenvalues λ of the monodromy matrix (i.e., Floquet multipliers): The corresponding eigenvectors give the velocity perturbations in the Fourier space. We note that there is typically only one growing mode for a given k y . We focus on the stability of the solitary waves since they are localized and not affected by the radial boundary condition. In Fig. 4(c), we plot the growth rates of non-axisymmetric modes with different k y in isothermal disks. The sharper solitary waves at smaller Q −1 are susceptible to more vigorous instabilities over a wider k y range. We note that a larger value of Γ narrows the window of unstable modes in k y , but the fastest growth rate over all k y is only slightly affected.
The v x structure for the fastest growing mode with Q −1 = 0.95 is shown in Fig. 4(b). We verified our Floquet analysis by carrying out direct hydrodynamic simulations with the ATHENA code (Stone et al. 2008) 1 We initialized the solitary wave in a square sheet of size 2π/K and added random velocity noise of 0.1% of the sound speed. The sheet is resolved by 8192 cells, i.e., 104 cells per c 0 /κ, to minimize numerical dissipation, and we note that the instability grows slower in low-resolution simulations. In Fig. 4(a), the simulation reproduces the velocity structure in the analytical calculation, and the measured fastest growing mode in the simulation has the predicted wavelength and growth rate. Fig. 4(a) shows contamination by other modes close to k y ≈ 2.5 because they have comparable growth rates.
The solitary waves are expected to transition quickly to spiral structures. In our direct simulations of isothermal disks, the growing non-axisymmetric perturbation leads to runaway collapse of the solitary wave on scales of a few tenths of c 0 /κ (Deng et al. 2021). This is expected in isothermal environments and resolving the collapsing clumps becomes increasingly challenging. Unfortunately, solving the energy equation in the Γ = 1 case with our current hydrodynamical code (Mullen et al. 2021) leads to inaccuracy, preventing us from directly simulating the expected transition to spirals. Further analysis and simulations on the effects of vertical stratification and the equation of state are desirable.

CONCLUSIONS.
We investigated the structure of nonlinear axisymmetric equilibria on the subcritical side of the GI bifurcation. Periodic ring structures of a wide range of length scales exist. A class of solitary waves with slightly higher energy than the smooth disk is of particular interest, which may be induced by external perturbations. The solitary waves are unstable to non-axisymmetric perturbations, providing a new way of understanding the generation of spirals in subcritical self-gravitating accretion disks 6. ACKNOWLEDGMENTS