Standing Solitary Waves as Transitions to Spiral Structures in Gravitationally Unstable Accretion Disks

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 center at an angular velocity ${\rm{\Omega }}(r)\,\hat{{{\boldsymbol{e}}}_{z}}$ is governed by the equation of motion:

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{u}}}{{Dt}}+2{\rm{\Omega }}\,\hat{{{\boldsymbol{e}}}_{z}}\times {\boldsymbol{u}}=-\displaystyle \frac{{\boldsymbol{\nabla }}P}{{\rm{\Sigma }}}-{\boldsymbol{\nabla }}({\rm{\Phi }}+{{\rm{\Phi }}}_{{\rm{t}}}),\end{eqnarray} \tag{ 1 }$$

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 self-gravity, and Φ_t = −ΩSx² is the tidal potential. For uniform Σ and P, the steady solution ${{\boldsymbol{u}}}_{0}=-{Sx}\,{\hat{{\boldsymbol{e}}}}_{y}$ represents the basic orbital motion with shear rate $S=-d{\rm{\Omega }}/d\mathrm{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 ₀.

Axisymmetric perturbations with radial wavenumber k_x oscillate at an angular frequency ω given by

$$\begin{eqnarray}&&{\omega }^{2}={\kappa }^{2}-2\pi G{\rm{\Sigma }}| {k}_{x}| +{c}^{2}{k}_{x}^{2},\end{eqnarray} \tag{ 2 }$$

where c is the sound speed. Long-wavelength (small k_x ) perturbations are stabilized by the κ² = 2Ω(2Ω − S) term (squared epicyclic frequency) resulting from the 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 nonaxisymmetric 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 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.

Two material invariants of the ideal fluid model are the specific entropy and the potential vorticity (PV) (2Ω − S + ∂_x v_y − ∂_y v_x )/Σ. Throughout this paper, we consider accessible solutions with the same uniform entropy, PV, and mean density as a uniform sheet (denoted by a subscript zero). For a perfect gas of (2D) adiabatic index Γ > 1, isentropy implies the polytropic relations P = AΣ^Γ and dP/Σ = dW, with specific enthalpy W = AΓΣ^Γ−1/(Γ − 1) = c²/(Γ − 1). Hydrostatic 3D disks have an effective Γ between (3γ − 1)/(γ + 1) and 3 − 2/γ (Gammie 2001), i.e., 1.33 < Γ < 1.57 for a cold disk of molecular hydrogen with γ = 1.4. Lower effective values of Γ may be relevant when radiative processes are taken into account, with Γ = 1 corresponding to the isothermal limit of instantaneous thermal relaxation.

For steady axisymmetric solutions (v_x = 0), the PV constraint implies

$$\begin{eqnarray}&&{\partial }_{x}{v}_{y}=\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\sigma ,\end{eqnarray} \tag{ 3 }$$

where σ = Σ/Σ₀ − 1 is the fractional density variation. The radial force balance reads

$$\begin{eqnarray}&&-2{\rm{\Omega }}{v}_{y}=-{\partial }_{x}W-{\partial }_{x}{\rm{\Phi }}\end{eqnarray} \tag{ 4 }$$

and combines with Equation (3) to give

$$\begin{eqnarray}&&{\partial }_{x}^{2}(W+{\rm{\Phi }})={\kappa }^{2}\sigma .\end{eqnarray} \tag{ 5 }$$

The dimensionless enthalpy perturbation

$$\begin{eqnarray}&&\displaystyle \frac{W-{W}_{0}}{{c}_{0}^{2}}=w=\displaystyle \frac{{(1+\sigma )}^{{\rm{\Gamma }}-1}-1}{{\rm{\Gamma }}-1}\end{eqnarray} \tag{ 6 }$$

reduces to $w=\mathrm{ln}(1+\sigma )$ 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 that

$$\begin{eqnarray}&&\tilde{{\rm{\Phi }}}=-\displaystyle \frac{2\pi G}{k}\tilde{{\rm{\Sigma }}}\end{eqnarray} \tag{ 7 }$$

for any wavenumber $k=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}\ne 0$. In Fourier space, Equation (5), when combined with Equation (7), gives

$$\begin{eqnarray}&&(2\pi G{{\rm{\Sigma }}}_{0}| {k}_{x}| -{\kappa }^{2})\tilde{\sigma }={c}_{0}^{2}{k}_{x}^{2}\tilde{w}.\end{eqnarray} \tag{ 8 }$$

We adopt c₀/κ as a natural unit of length. Radial force balance can be expressed in real space as

$$\begin{eqnarray}&&{w}_{{xx}}+2{Q}^{-1}{ \mathcal H }{\sigma }_{x}=\sigma ,\end{eqnarray} \tag{ 9 }$$

where Q⁻¹ = π GΣ₀/κ c₀ and ${ \mathcal H }$ is the Hilbert transform (King 2009). The polytropic relation (Equation (6)) makes the problem nonlinear. We look for periodic solutions with a basic wavenumber K in the form of a cosine series,

$$\begin{eqnarray}&&\sigma =\sum _{n=1}^{\infty }{\sigma }_{n}\cos ({nKx}),\quad w=\sum _{n=0}^{\infty }{w}_{n}\cos ({nKx}),\end{eqnarray} \tag{ 10 }$$

so that the coefficients are related by

$$\begin{eqnarray}&&{(2{Q}^{-1}{nK}-1){\sigma }_{n}=({nK})}^{2}{w}_{n}.\end{eqnarray} \tag{ 11 }$$

Asymptotic analysis is possible when the disk is nearly marginally stable, i.e., Q⁻¹ and K are close to 1 with small increments δ(Q⁻¹) and δ K. For weakly nonlinear solutions,

$$\begin{eqnarray}&&2\delta {({Q}^{-1})-(\delta K)}^{2}\approx -\displaystyle \frac{(2-{\rm{\Gamma }})(5-3{\rm{\Gamma }})}{8}{\sigma }_{1}^{2},\end{eqnarray} \tag{ 12 }$$

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⁻¹) < 0, i.e., Q > 1. Using these approximate solutions as initial guesses, we solve Equations (6) and (11) by the Newton–Raphson iteration for a range of parameters Q⁻¹, K, and Γ. We utilize Fast Fourier Transform (FFT) algorithms with up to 5000 points to achieve good convergence when K is small.

In Figure 1, we plot the steady axisymmetric structures in various conditions. For a given (Q⁻¹, Γ) pair, the short-wavelength solutions (K > 1) possess a single strong peak that tends to be infinitely sharp at large K, whereas the long-wavelength (K < 1) solutions gradually approach a solitary wave. For example, in Figure 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⁻¹, Γ), whose shapes are well established, are all calculated with a base wavenumber K = 0.08. A smaller Q⁻¹ or a larger Γ leads to a more strongly peaked wave.

Zoom In Zoom Out Reset image size

With a fixed wavenumber K = 1 and Γ < 5/3, we show the peak density perturbation of periodic ring structures as a function of Q⁻¹ in Figure 2. The rings become more peaked when Q⁻¹ decreases from 1 until the solutions turn onto an upper branch at a critical Q⁻¹ (0.712, 0.825, 0.900, and 0.950 for Γ = 1.1, 1.2, 1.3, and 1.4). The peak density continues to increase with Q⁻¹ until the first trough in the periodic structure (similar to Figure 1) reaches zero and the upper branch ceases.

Zoom In Zoom Out Reset image size

The energy per unit area of the axisymmetric structures, relative to the uniform disk, is

$$\begin{eqnarray}&&{\rm{\Sigma }}\left(\displaystyle \frac{1}{2}{v}_{x}^{2}+\displaystyle \frac{2{{\rm{\Omega }}}^{2}}{{\kappa }^{2}}{v}_{y}^{2}\right)+\displaystyle \frac{1}{2}({\rm{\Sigma }}-{{\rm{\Sigma }}}_{0}){\rm{\Phi }}+\displaystyle \frac{P-{P}_{0}}{{\rm{\Gamma }}-1},\end{eqnarray} \tag{ 13 }$$

where the first two terms combine the kinetic and tidal energies. Considering the isovortical constraint (Equation (3)) and radial force balance (Equation (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

$$\begin{eqnarray}&&\iint \displaystyle \frac{{{\rm{\Sigma }}}_{0}{c}_{0}^{2}}{{\rm{\Gamma }}}\left[w+\left(1-\displaystyle \frac{{\rm{\Gamma }}}{2}\right)\sigma w\right]{dx}\,{dy}.\end{eqnarray} \tag{ 14 }$$

So the mean energy change of periodic structures per unit area is

$$\begin{eqnarray}&&{\rm{\Delta }}E=\displaystyle \frac{{{\rm{\Sigma }}}_{0}{c}_{0}^{2}}{{\rm{\Gamma }}}\left[{w}_{0}+\left(\displaystyle \frac{1}{2}-\displaystyle \frac{{\rm{\Gamma }}}{4}\right)\sum _{n=1}^{\infty }{\sigma }_{n}{w}_{n}\right].\end{eqnarray} \tag{ 15 }$$

We plot the energy change due to periodic ring structures and solitary waves in Figure 3. For the K = 1 solutions in Figure 2, the energy first increases as Q⁻¹ decreases along the lower branch as shown in Figure 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⁻¹. 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.

Zoom In Zoom Out Reset image size

In Figure 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 Figure 3(c), the energy of solitary waves increases with Γ and Q.

We investigate the stability of the periodic and solitary waves to nonaxisymmetric perturbations (k_y ≠ 0). The evolution of linear perturbations (primed variables) is governed by

$$\begin{eqnarray}&&\displaystyle \frac{{{Dv}}_{x}^{{\prime} }}{{Dt}}-2{\rm{\Omega }}{v}_{y}^{{\prime} }=-{\partial }_{x}({W}^{{\prime} }+{{\rm{\Phi }}}^{{\prime} }),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{Dv}}_{y}^{{\prime} }}{{Dt}}+(2{\rm{\Omega }}-S+{\partial }_{x}{v}_{y}){v}_{x}^{{\prime} }=-{\partial }_{y}({W}^{{\prime} }+{{\rm{\Phi }}}^{{\prime} }),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D{{\rm{\Sigma }}}^{{\prime} }}{{Dt}}+{v}_{x}^{{\prime} }{\partial }_{x}{\rm{\Sigma }}=-{\rm{\Sigma }}({\partial }_{x}{v}_{x}^{{\prime} }+{\partial }_{y}{v}_{y}^{{\prime} }),\end{eqnarray} \tag{ 18 }$$

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.,

$$\begin{eqnarray}&&{{\boldsymbol{v}}}^{{\prime} }=\sum _{n=-\infty }^{\infty }{{\boldsymbol{v}}}_{n}^{{\prime} }(t)\exp ({{ik}}_{x,n}x+{{ik}}_{y}y),\end{eqnarray} \tag{ 19 }$$

with k_x,n = nK + k_y St. We only consider isovortical perturbations with zero potential vorticity change, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }_{x}{v}_{y}^{{\prime} }-{\partial }_{y}{v}_{x}^{{\prime} }}{2{\rm{\Omega }}-S}=\displaystyle \frac{{\rm{\Sigma }}^{\prime} }{{{\rm{\Sigma }}}_{0}}.\end{eqnarray} \tag{ 20 }$$

The density perturbation can thus be directly related to the velocity perturbations, leaving two coupled linear ordinary differential equations (ODEs) for each n:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dv}}_{x,n}^{{\prime} }}{{dt}}+{{ik}}_{y}\displaystyle \sum _{m=-\infty }^{\infty }{v}_{y,m}{v}_{x,n-m}^{{\prime} }\\ \quad -\,2{\rm{\Omega }}{v}_{y,n}^{{\prime} }=-{{ik}}_{x,n}({W}_{n}^{{\prime} }+{{\rm{\Phi }}}_{n}^{{\prime} }),\end{array}\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dv}}_{y,n}^{{\prime} }}{{dt}}+{{ik}}_{y}\displaystyle \sum _{m=-\infty }^{\infty }{v}_{y,m}{v}_{y,n-m}^{{\prime} }+(2{\rm{\Omega }}-S){v}_{x,n}^{{\prime} }\\ \ +\,\displaystyle \sum _{m=-\infty }^{\infty }{\left({\partial }_{x}{v}_{y}\right)}_{m}{v}_{x,n-m}^{{\prime} }=-{{ik}}_{y}({W}_{n}^{{\prime} }+{{\rm{\Phi }}}_{n}^{{\prime} }).\end{array}\end{eqnarray} \tag{ 21b }$$

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 ${\rm{\Sigma }}^{\prime}$ can be related to the velocity perturbation coefficients by the isovortical condition. ${{\rm{\Phi }}}_{n}^{{\prime} }$ is readily obtained from ${{\rm{\Sigma }}}_{n}^{{\prime} }$ by the Poisson equation, and ${W}_{n}^{{\prime} }$ can be calculated via a convolution between ${{\rm{\Sigma }}}_{n}^{{\prime} }$ and ${({c}^{2}/{\rm{\Sigma }})}_{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 (21a) and integrating it for one recurrence time. To ensure the ${{\boldsymbol{v}}}_{\pm N}^{{\prime} }$ entries are properly evaluated, we extended the Fourier series, employing 2(N + N₀) + 1 wavenumbers for each variable in Equation (21a). Here, N₀ is the number of nonzero (>10⁻¹⁴) modes in the equilibrium axisymmetric structures so that the convolutions in Equation (21a) 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₀ 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):

$$\begin{eqnarray}&&s=\displaystyle \frac{1}{T}\mathrm{maxRe}(\mathrm{ln}\lambda ).\end{eqnarray} \tag{ 22 }$$

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 because they are localized and not affected by the radial boundary condition. In Figure 4(c), we plot the growth rates of nonaxisymmetric modes with different k_y in isothermal disks. The sharper solitary waves at smaller Q⁻¹ 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.

Zoom In Zoom Out Reset image size

The ${v}_{x}^{{\prime} }$ structure for the fastest-growing mode with Q⁻¹ = 0.95 is shown in Figure 4(b). We verified our Floquet analysis by carrying out direct hydrodynamic simulations with the ATHENA code (Stone et al. 2008). ³ We initialized the solitary wave in a square sheet of size 2π/K and added a random velocity noise of 0.1% of the sound speed. The sheet is resolved by 8192 cells, i.e., 104 cells per ${c}_{0}/\kappa$, to minimize numerical dissipation, and we note that the instability grows slower in low-resolution simulations. In Figure 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. Figure 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 nonaxisymmetric perturbation leads to the runaway collapse of the solitary wave on scales of a few tenths of c₀/κ (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.

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 nonaxisymmetric perturbations, providing a new way of understanding the generation of spirals in subcritical self-gravitating accretion disks.

This research was funded by the Isaac Newton Trust (University of Cambridge) through research grant 21.07(d). H.D. also acknowledges support from the Swiss National Science Foundation via a postdoctoral fellowship. We thank the anonymous referee for suggestions that improved the clarity of the paper.

Software: SciPy (Virtanen et al. 2020), ATHENA (Stone et al. 2008).