Schrödinger evolution of self-gravitating discs

Abstract

An understanding of the long-term evolution of self-gravitating discs ranks among the classic outstanding problems of astrophysics. In this work, we show that the secular inclination dynamics of a geometrically thin quasi-Keplerian disc, with a surface density profile that scales as the inverse square-root of the orbital radius, are described by the time-dependent Schrödinger equation. Within the context of this formalism, nodal bending waves correspond to the eigenmodes of a quasi-particle's wavefunction, confined in an infinite square well with boundaries given by the radial extent of the disc. We further show that external secular perturbations upon self-gravitating discs exhibit a mathematical similarity to quantum scattering theory. Employing this framework, we derive an analytic criterion for the gravitational rigidity of a nearly-Keplerian disc under external perturbations. Applications of the theory to circumstellar discs and Galactic nuclei are discussed.

methods: analytical, protoplanetary discs, Galaxy: nucleus

1 INTRODUCTION

Astrophysical discs are among the most ubiquitous objects in the known universe. Generically understood to be a consequence of energy dissipation within gravitationally bound rotating systems, these structures span a staggering assortment of scales, ranging from galaxies to protoplanetary nebulae and circumplanetary rings (Goldreich & Tremaine 1982). Accordingly, characterizing the long-term evolution of self-gravitating discs constitutes one of the key challenges of dynamical astronomy.

Owing to their inherent diversity, astrophysical discs arise in nature with variable compositions, and can occupy rather distinct physical regimes. For example, active galactic nuclei and protoplanetary discs are predominantly composed of hydrogen and helium gas, and are representative of fluid discs. As such, their evolution is governed by gravitational as well as (magneto-)hydrodynamic effects. Conversely, planetesimal/debris discs, as well as discs of stars that orbit supermassive black holes in the centres of galaxies (i.e. so-called particle discs), are subject to essentially pure gravitational dynamics¹ (Latter et al. 2017).

Observed instances of astrophysical discs often encircle central objects that are much more massive than the discs themselves. The resulting dominance of the central body's gravitational potential leads to quasi-Keplerian particle motion that resembles planetary orbits on short (orbital) time-scales, but can exhibit non-trivial behaviour over much longer periods of time, due to self-gravity (Tremaine 2001; Touma & Tremaine 2014). Describing the long-term exchange of angular momentum within such systems is the primary aim of this paper.

The specific goals of our analysis are essentially two-fold. Our first aim is to formulate a tangible theory for the secular inclination evolution of dynamically cold self-gravitating discs. Our second goal is to derive a dimensionless number, somewhat akin to Toomre's Q, that characterizes the propensity of discs to warp under external perturbations. In other words, we seek to obtain an analytical criterion for the gravitational rigidity of astrophysical discs.

Quantifying the self-gravitational evolution of discs and their global response to external perturbations is essential to interpreting modern observations, as well as to the development of planet formation theory. In particular, the dynamical origins of the warped stellar disc at the centre of our Galaxy have been attributed to a complex interplay between self-gravitational effects and torques from surrounding stellar clusters (Kocsis & Tremaine 2011). Likewise, large-scale warps and spiral morphology of young circumstellar discs² (see e.g. Golimowski et al. 2006; Backman et al. 2009 and the references therein) are routinely ascribed to interactions between the discs themselves and perturbing stellar or planetary companions (Nesvold, Naoz & Vican 2016; Nesvold et al. 2017). Even the seemingly unrelated question of the provenance of spin-orbit misalignments in exo-planetary systems may be deeply rooted in the process of gravitational torquing of protoplanetary nebulae by bound or passing stars (Bate et al. 2010; Spalding & Batygin 2014).

Importantly, all of these phenomena arise on the so-called secular time-scale – one that greatly exceeds the characteristic orbital period, but is significantly shorter than the physical lifetime of the system (e.g. ∼10⁵ yr). By operating in between the aforementioned temporal extrema, the dynamical mechanisms at play are poorly represented by either an explicit (N-body) description of orbital motion, or a diffusive model of angular momentum transfer within discs (Pringle 1981). As a result, secular behaviour of continuous self-gravitating systems remains imperfectly understood, and warrants the development of a simple theoretical model.

At present, there exist three primary means of analysing the dynamics of self-gravitating discs. The most straight-forward route is to employ direct N-body simulations. In spite of remarkable advances in computation that have transpired within the last decade (Gaburov et al. 2009; Moore & Quillen 2011), this method remains too computationally expensive and specialized for most problems of interest. A second, more compact technique utilizes the collisionless Boltzmann equation, which, instead of granting N individual trajectories as the solution, yields the evolution of the system's distribution function in phase-space (Binney & Tremaine 1987; Sridhar & Touma 2017). A third approach relies on Gauss's averaging method of celestial mechanics, to replace the individual bodies with a series of massive wires that interact gravitationally among one another on time-scales much longer than the orbital period (Touma et al. 2009; Batygin 2012).

In this work, we focus exclusively on the third, mean-field model of self-gravitating discs. Specifically, working within the context of the Lagrange–Laplace secular perturbation theory, we demonstrate that the N-ring description of self-gravitating discs can be reduced to an evolution governed by Schrödinger's equation. Correspondingly, this simplified framework allows for the computation of gravitational rigidity of astrophysical discs, and yields new insight into their long-term dynamical evolution.

The paper is structured as follows. In Section 2, we obtain Schrödinger's equation as a continuum limit of Hamilton's equations. In Section 3, we apply this calculation to the inclination dynamics of nearly circular self-gravitating discs, and derive the corresponding eigenfunctions and eigenfrequencies. Importantly, we note that while the eigenfunctions of the disc can be approximately obtained by only considering nearest-neighbour interactions (Section 3.2.3), the computation of the system's eigenfrequencies requires accounting for collective effects within the disc (Section 3.2.5). In Section 4, we extend this formalism to account for external perturbations, and deduce an analytic criterion for gravitational rigidity of perturbed systems. We summarize and discuss our results in Section 5.

2 SCHRÖDINGER's EQUATION FROM HAMILTON's EQUATIONS

Prior to considering self-gravitating discs explicitly, let us deleniate a general framework that the calculations will follow. Recall that quantum and classical mechanics are routinely thought to be governed by Schrödinger's and Hamilton's equations, respectively. Being macroscopic in nature, the latter is generally considered to be a limiting case of the former, where the action quantum ℏ → 0. This notion is strongly related to Ehrenfest's theorem, and is often referred to as the correspondence principle (see e.g. Sakurai 1985). In this section, we demonstrate the converse to also be true in a specific case. That is, for a particular Hamiltonian, Schrödinger's equation naturally arises as a limiting case of Hamilton's equations.

Consider an infinite sequence of coupled objects lying on the x-axis, whose individual dynamical state is described by a pair of action-angle variables ... (Φ_j − 1, ϕ_j − 1), (Φ_j, ϕ_j), (Φ_j + 1, ϕ_j + 1) ... . Let the spacing between the objects be equidistant, and denote it δx (Fig. 1). Finally, let the purely classical Hamiltonian describing the evolution of object j be

\begin{eqnarray} \mathcal {H}&=& 2\,c_1\,\Phi _j+2\,c_2\,\sqrt{\Phi _j\,\Phi _{j+1}} \cos (\phi _j-\phi _{j+1}) \nonumber \\ &&+\, 2\,c_2\,\sqrt{\Phi _j\,\Phi _{j-1}} \cos (\phi _j-\phi _{j-1}), \end{eqnarray}

(1)

where c₁ and c₂ are quantities that can depend on x and t. Note that in the above expression, only interactions between the nearest neighbours are considered.

Figure 1.

An infinite chain of harmonically coupled oscillators. The positions of the objects are fixed on the x-axis, such that the distance between neighbouring sites is δx. The dynamical state of each oscillator is described by a pair of action-angle coordinates, (Φ, ϕ).

Next, we define canonical Cartesian analogues (ξ, ζ) to the action-angle variables (Φ, ϕ):

\begin{eqnarray} \xi =\sqrt{2\,\Phi }\cos (\phi )\quad \zeta =\sqrt{2\,\Phi }\sin (\phi ), \end{eqnarray}

(2)

and organize them into a single complex coordinate:

\begin{eqnarray} \Psi = \frac{\xi +\imath \,\zeta }{\sqrt{2}} = \sqrt{\Phi } \exp (\imath \,\phi ), \end{eqnarray}

(3)

where |$\imath =\sqrt{-1}$|⁠. The Hamiltonian then reads:

\begin{eqnarray} \mathcal {H}&=& 2\,c_1\,\Psi _j \,\Psi ^*_j+2\,c_2\,\big (\Psi _j \, \Psi ^{*}_{j+1} + \Psi ^{*}_{j} \Psi _{j+1}+ \Psi _j \, \Psi ^{*}_{j-1} \nonumber \\ &&+\, \Psi ^{*}_{j} \Psi _{j-1}\big ), \end{eqnarray}

(4)

where Ψ* is a complex conjugate to Ψ.

In terms of complex variables, Hamilton's equations take on a rather succinct form (Strocchi 1966):

\begin{eqnarray} \frac{\text{d}\Psi _j}{\text{d}t} = \imath \,\frac{\mathrm{\partial} \mathcal {H}}{\mathrm{\partial} \Psi _j^{*}} = \imath \, \big (2\,c_1\,\Psi _j +c_2 \left(\Psi _{j+1} + \Psi _{j-1} \right) \big ). \end{eqnarray}

(5)

If we envision Ψ_j − 1, Ψ_j, Ψ_j + 1, ... to be a discrete representation of a continuous, complex field Ψ, we may write down the central difference approximation for its second derivative as follows:

\begin{eqnarray} \frac{\mathrm{\partial} ^2\,\Psi _{j}}{\mathrm{\partial} \,x^2} \simeq \frac{\Psi _{j+1} - 2\,\Psi _j + \Psi _{j-1}}{\left(\delta x\right)^2}. \end{eqnarray}

(6)

We note that a similar approximation is routinely made in the elastic continuum representation of harmonic crystals (Animalu 1977), and corresponds to the limit where the wave numbers of interest greatly exceed the spacing between atomic sites (i.e. 2π/k ≫ δx).

Plugging in equation (6) into equation (5) and multiplying both sides by −ℏ/ı, we obtain:

\begin{eqnarray} \imath \,\hbar \, \frac{\mathrm{\partial} \Psi _j}{\mathrm{\partial} t} = -\hbar \,\bigg (2(c_1+c_2)\Psi _{j} + c_2 \left(\delta x \right)^2 \frac{\mathrm{\partial} ^2\,\Psi _j}{\mathrm{\partial} \,x^2} \bigg ). \end{eqnarray}

(7)

Note that here, we have taken advantage of the fact that the dynamical objects are fixed in their x-coordinate to set dΨ/dt = ∂Ψ/∂t. Adopting the following expressions for the interaction constants:

\begin{eqnarray} c_1 = -\frac{V(x,t)}{2}- \frac{\hbar }{2\,\mu \left(\delta x \right)^2}\quad c_2 = \frac{\hbar }{2\,\mu \left(\delta x \right)^2}, \end{eqnarray}

(8)

and taking the limit as δx → 0 (which renders the relationship in equation 6 exact), we obtain the time-dependent Schrödinger's equation:

\begin{eqnarray} \imath \,\hbar \,\frac{\mathrm{\partial} \, \Psi (x,t)}{\mathrm{\partial} t} = \left[ -\frac{\hbar ^2}{2\,\mu } \frac{\mathrm{\partial} ^2}{\mathrm{\partial} \,x^2} + V(x,t) \right] \Psi (x,t). \end{eqnarray}

(9)

This derivation is trivially generalizable to three spatial dimensions. Moreover, a simple modification of the same procedure naturally leads to the non-linear variant of Schrödinger's equation. That is, addition of a non-linear action term to equation (1), such that

\begin{eqnarray} \mathcal {H}^{\prime } = \mathcal {H}+ \frac{\kappa }{2}\,\Phi _j^2,=\mathcal {H}+\frac{\kappa }{2}\,(\Psi _j\,\Psi _j^*)^2 \end{eqnarray}

(10)

yields

\begin{eqnarray} \imath \,\hbar \frac{\mathrm{\partial} \, \Psi (x,t)}{\mathrm{\partial} t} &=& \left[ -\frac{\hbar ^2}{2\,\mu } \frac{\mathrm{\partial} ^2 }{\mathrm{\partial} \,x^2} + V(x,t) \right] \Psi (x,t) \nonumber \\ &&+\, \kappa \,|\Psi (x,t)|^2\,\Psi (x,t). \end{eqnarray}

(11)

It is worth noting that a Hamiltonian of the form (10) is often referred to as the ‘second fundamental model of resonance,’ and possesses a phase-space topology similar to that of a harmonic oscillator or a pendulum, depending on the assumed values of the constants (Henrard & Lemaitre 1983).

This work is by no means the first to derive Schrödinger's equation from purely classical principles. For example, Nelson (1966) showed that Schrödinger's equation can be obtained as a description of a stochastic system, where particles are subjected to Brownian motion with a diffusion coefficient ℏ/2m, in addition to Newton's laws of motion. Extension of this formalism to an effective field theory has been considered by Guerra (1981), and has led to the so-called stochastic interpretation of quantum mechanics.

In contrast with these works, here we did not draw on stochastic fluctuations to augment classical mechanics. At the same time, we emphasize that even though the above derivation is self-consistent, it is not general, and only applies to a particular Hamiltonian that describes a continuous chain of forced harmonic oscillators with specific interaction coefficients. Incidentally, the governing equations of Largrange–Laplace secular theory of celestial mechanics can be cast into this form.

3 SECULAR EVOLUTION OF SELF-GRAVITATING DISCS

Having outlined the general course of action in the previous section, we now proceed to show how the Schrödinger equation can be used to understand the long-term evolution of self-gravitating discs. However, prior to computing the dynamical evolution itself, it is necessary to define a number of basic properties of our model.

3.1 Disc profile

Consider a disc of material, comprised of a large number of point masses in orbit around a single central body of mass M. Envision that the gravitational potential of the central body dominates, such that the trajectories of the individual bodies follow Kepler's third law:

\begin{eqnarray} n = \sqrt{\frac{\mathcal {G}M}{a^3}}, \end{eqnarray}

(12)

where n is the mean motion, and a is the semimajor axis. Further, assume that the orbital eccentricities as well as mutual inclinations of neighbouring orbits are not exceedingly large, meaning that the inherent velocity dispersion of a population of objects occupying the same semimajor axis range is modest compared to the Keplerian velocity. Correspondingly, we define the aspect ratio of the disc, β, as an intrinsic small parameter of the problem, and take it to be constant:

\begin{eqnarray} \beta \equiv \frac{h}{a} = \rm {const.} \ll 1, \end{eqnarray}

(13)

where h is the characteristic scale height. Crucially, the value of β sets the gravitational softening length of the disc³ (Adams et al. 1989; Touma 2002).

Let us now replace this collection of secondary bodies with a series of N massive ‘streams,’ thus turning the aggregate of orbiting particles into a disc comprised of nested elliptical wires, with normalized thickness β (Fig. 2). Importantly, we assume that the wires do not cross and that the epicyclic motion of constituent particles (i.e. distribution of eccentricities at a given semimajor axis) is fully encapsulated by softening parameter (Binney & Tremaine 1987). The spacing of the wires is taken to be geometric, such that the semimajor axis ratio of neighbouring orbits, α, is constant, and roughly corresponds to a single softening length:

\begin{eqnarray} \alpha \equiv \frac{a_{j-1}}{a_j} = \frac{a_{j}}{a_{j+1}} = \frac{1}{1+\beta } \approx 1. \end{eqnarray}

(14)

In light of this definition, it is useful to scale the semimajor axis by the disc's inner truncation radius, a_in, and introduce a dimensionless logarithmic radial coordinate

\begin{eqnarray} \rho \equiv \log \left( \frac{a}{a_{\rm {in}}} \right)\quad \mathcal {L} \equiv \log \left( \frac{a_{\rm {out}}}{a_{\rm {in}}} \right), \end{eqnarray}

(15)

where a_out denotes the outer boundary of the disc. In terms of this quantity, the wires are spaced equidistantly, and the disc boundaries are given by |$\rho \in [0,\mathcal {L}]$|⁠. To this end, we note that any realistic system will have |$\mathcal {L}$| on the order of a few, rarely exceeding 10. Coincidentally, for a typical protoplanetary nebula (e.g. a_in ∼ 0.05 au and a_out ∼ 50 au), |$\mathcal {L}\sim 2\,\pi$|⁠.

Figure 2.

Geometric setup of our model. A quasi-Keplerian disc of N ≫ 1 particles is modelled as a sequence of geometrically spaced massive wires. The gravitational potential of each wire is softened by the aspect ratio of the disc, accounting for the inherent velocity dispersion of the constituent particles. A surface density profile that scales inversely with the square root of the semimajor axis is assumed.

Qualitatively speaking, the process of ‘smearing out’ individual bodies into massive rings, whose line densities are inversely proportional to orbital velocity, is equivalent to canonically averaging over the rapidly varying orbital angles (i.e. mean longitudes). This procedure yields orbital semimajor axes, a, that are frozen in time, because the action conjugate to mean longitude is proportional to |$\sqrt{a}$| (Morbidelli 2002). In turn, this conservation implies that the two-body (Keplerian) energy of each object is preserved, and the wires only exchange angular momentum. Physically, this process leads to slow precession of the apsidal and nodal lines of the orbits as well as variations in eccentricities and inclinations.

A conventional parametrization of density variations in astrophysical discs assumes that the surface density, Σ, scales as some negative power of the orbital radius (Armitage 2010). Here, we follow this standard prescription and adopt an inverse square root surface density profile of the disc:

\begin{eqnarray} \Sigma = \Sigma _0 \, \bigg (\frac{a_0}{a}\bigg )^{1/2}, \end{eqnarray}

(16)

where Σ₀ is the surface density at a reference semimajor axis, a₀. We note that while this profile is somewhat shallower than, say, a classical Mestel (1963) disc which takes Σ ∝ 1/a, it is routinely employed in numerical studies of protoplanetary discs (e.g. Bitsch & Kley 2011; Lambrechts & Lega 2017). Strictly speaking, this approximation is not necessary. However, here we choose to adopt expression (16) for illustrative purposes, as it will simplify some of the calculations down the line.

With the above definitions in place, we obtain the following expression for the mass of an individual wire (of index j):

\begin{eqnarray} m_j = \oint \int _{a(1-\beta /2)}^{a(1+\beta /2)} \Sigma \,a\,\text{d}a\,\text{d}\varphi \approx 2\,\pi \,\beta \,\Sigma _0\sqrt{a_0\,a_j^3}+\mathcal {O}(\beta ^3), \end{eqnarray}

(17)

where φ denotes the azimuthal coordinate. Correspondingly, assuming that a_in ≪ a_out, the total disc mass is given by

\begin{eqnarray} m_{\rm {\text{disc}}} = \oint \int _{a_{\rm {in}}}^{a_{\rm {out}}} \Sigma \,a\,\text{d}a\,\text{d}\varphi \approx \frac{4\,\pi }{3}\,\Sigma _0\,\sqrt{a_0\,a_{\rm {out}}^3}. \end{eqnarray}

(18)

The approximations employed above are only sensible if the dynamical evolution that ensues on the orbital time-scale does not deviate significantly from purely Keplerian motion. Quantitatively, this statement corresponds to the requirement of the gravitational stability of the system (Safronov 1960; Toomre 1964):

\begin{eqnarray} Q = \frac{h\,n^2}{\pi \,\mathcal {G}\,\Sigma } \gtrsim 1. \end{eqnarray}

(19)

In practical terms, this restriction translates to an upper limit on |$m_{\rm {\text{disc}}}$|⁠:

\begin{eqnarray} m_{\rm {\text{disc}}} \lesssim \beta \, \frac{M}{2}, \end{eqnarray}

(20)

and yields a connection between the inherent velocity dispersion of constituent matter and surface density of the disc. With the preliminary specifications of the model now in place, we now continue on to compute the dynamical evolution.

3.2 Inclination dynamics

Lagrange–Laplace secular theory constitutes one of the earliest, and best-known results of perturbation theory in celestial mechanics. Within the framework of this model, the phase-averaged gravitational potential of the interacting bodies (i.e. the negative disturbing function) is expanded as a Fourier series in the orbital angles and as a power-series in eccentricities and inclinations (Ellis & Murray 2000). Truncating the expansion at second order in both quantities yields secular vibrations of inclinations that are decoupled from the oscillations of eccentricities.⁴ Taking advantage of this disconnect between the degrees of freedom, here we treat the inclination dynamics, simply assuming that the eccentricities are small.

3.2.1 Governing equation

Consider the inclination dynamics of a wire (labeled by index j ≠ 1, N) embedded within a radially extended disc. As a first approximation, let us restrict the range of interactions of wire j to its nearest neighbours, j + 1 and j − 1. This approximation has the obvious shortcoming of reduced coupling within the disc, and we will revisit (and fix) this limitation later in the manuscript. The relevant (scaled) disturbing function⁵ that governs the exchange of angular momentum is (Murray & Dermott 1999):

\begin{eqnarray} \mathcal {R}_j &=& \frac{1}{2}B_{jj} \, i_j^2 + B_{jj-1} \, i_j\,i_{j-1} \cos (\Omega _j - \Omega _{j-1}) \nonumber \\ &&+ B_{jj+1} \, i_j\,i_{j+1} \cos (\Omega _j - \Omega _{j+1}), \end{eqnarray}

(21)

where i is orbital inclination, Ω is the longitude of ascending node, and B's are interaction coefficients that depend only on the semimajor axis ratios and masses.

Although Keplerian elements i and Ω do not constitute a set of canonically conjugated variables, the quantities

\begin{eqnarray} p=i\,\sin (\Omega )\quad q=i\,\cos (\Omega ) \end{eqnarray}

(22)

do, where p is interpreted as the coordinate and q is the momentum (with |$\mathcal {R}_j$| acting as the Hamiltonian). Note that physically, p and q represent a measure of the angular momentum deficit of wire j in the |$\hat{z}$|-direction. Collecting these variables into a single complex coordinate akin to equation (3):

\begin{eqnarray} \eta =\frac{q+\imath \,p}{\sqrt{2}}=\frac{i}{\sqrt{2}}\exp (\imath \, \Omega), \end{eqnarray}

(23)

we rewrite equation (21) as follows:

\begin{eqnarray} \mathcal {R}_j &=& B_{jj} \, \eta _j \,\eta ^*_j + B_{jj-1} \, \big (\eta _j \,\eta ^*_{j-1}+\eta _{j-1} \,\eta ^*_{j} \big ) \nonumber \\ &&+\, B_{jj+1} \, \big (\eta _j \,\eta ^*_{j+1}+\eta _{j+1} \,\eta ^*_{j} \big ). \end{eqnarray}

(24)

This expression is clearly reminiscent of Hamiltonian (4). However, in order to draw further analogy between equations (24) and (4), we must evaluate the relationship between the coefficients B_jj, B_jj − 1, and B_jj + 1. Under the assumption that m_j ≪ M, these quantities are expressed as follows (Murray & Dermott 1999):

\begin{eqnarray} B_{jj-1}&=& \frac{n_j}{4}\frac{m_{j-1}}{M}\,\alpha \,\tilde{b}_{3/2}^{(1)}\lbrace \alpha \rbrace \nonumber \\ B_{jj+1}&=& \frac{n_j}{4}\frac{m_{j+1}}{M}\,\alpha ^2\,\tilde{b}_{3/2}^{(1)} \lbrace \alpha \rbrace \nonumber \\ B_{jj}&=& -\big ( B_{jj-1}+B_{jj+1} \big ), \end{eqnarray}

(25)

where |$\tilde{b}_{3/2}^{(1)}\lbrace \alpha \rbrace$| is the Laplace coefficient of the first kind. Recalling from equation (14) that α ≈ 1 − β, expressions (25) imply B_jj − 1 ≈ B_jj + 1. In particular, upon defining the quantity

\begin{eqnarray} \mathcal {B}&=& \frac{n_j}{4}\frac{m_{j}}{M}\,\frac{\alpha }{2}\,\tilde{b}_{3/2}^{(1)} \lbrace \alpha \rbrace \bigg [ \frac{m_{j-1 }}{m_j}+\frac{m_{j+1}}{m_j}\frac{1}{1+\beta } \bigg ] \nonumber \\ &=& \frac{n_j}{4}\frac{m_{j}}{M}\,\alpha \,\tilde{b}_{3/2}^{(1)} \lbrace \alpha \rbrace \bigg [ \frac{2+2\,\beta +\beta ^2}{2\,(1+\beta )^{3/2}} \bigg ], \end{eqnarray}

(26)

we have:

\begin{eqnarray} \vert B_{jj-1} - \mathcal {B} \vert &=& | B_{jj+1} - \mathcal {B} | \approx \frac{\beta }{4} + \mathcal {O}(\beta ^2) \ll 1 \nonumber \\ B_{jj} &=& - 2\,\mathcal {B}. \end{eqnarray}

(27)

At this point, the advantage of choosing the specific form of the surface density profile (16) becomes evident. Since |$m_j\propto \sqrt{a^3}$| and |$n_j\propto 1/\sqrt{a^3}$|⁠, the two dependencies cancel, rendering |$\mathcal {B}$| constant throughout the disc. Thus, to an excellent approximation,

\begin{eqnarray} \frac{d\eta _j}{dt} = \imath \frac{\mathrm{\partial} \mathcal {R}_j}{\mathrm{\partial} \eta _j^*}\approx \imath \, \mathcal {B}\, \big (\eta _{j-1}-2\eta _j+\eta _{j+1} \big ). \end{eqnarray}

(28)

We note that under a different choice for the radial dependence of the surface density profile, the quantity |$\mathcal {B}$| would become an explicit function of the semimajor axis.

Because the disc annuli are spaced geometrically in a, they are equidistant in ρ (equation 15), with

\begin{eqnarray} \delta \rho =\rho _{j+1}-\rho _{j}=\rho _{j}-\rho _{j-1}=\log (1+\beta ). \end{eqnarray}

(29)

Hence, employing the finite-difference approximation (equation 6) and recalling that the semimajor axes of the wires are secularly invariant, the continuum limit of equation (28) takes the form:

\begin{eqnarray} \frac{\mathrm{\partial} \eta }{\mathrm{\partial} t} = \imath \, \mathcal {B}\,\big (\log (1+\beta ) \big )^2\,\frac{\mathrm{\partial} ^2\eta }{\mathrm{\partial} \rho ^2}. \end{eqnarray}

(30)

Due to the proximity of the neighbouring wires to each other, it is sensible to evaluate equation (26) in the limit where α → 1. However, Laplace coefficients are notoriously singular at α = 1 (this is simply a re-statement of the fact that the gravitational potential becomes infinite at null separations). This mathematical inconvenience is easily circumvented by accounting for the inherent velocity dispersion of particles within the disc and softening the Laplace coefficient by the disc aspect ratio (Touma 2002; Hahn 2003):

\begin{eqnarray} \tilde{b}_{\ell }^{(j)} \lbrace \alpha \rbrace =\frac{2}{\pi }\int _{0}^{\pi }\frac{\cos (j\,\psi )}{\big (1-2\,\alpha \cos (\psi )+\alpha ^2+\beta ^2 \big )^{\ell }}\, \text{d}\psi . \end{eqnarray}

(31)

These softened Laplace coefficients are well-behaved in the α → 1 limit, and the quantity |$\tilde{b}_{3/2}^{(1)}$| can now be expressed in terms of standard elliptic integrals. For the time being, it suffices to evaluate this function at α = 1/(1 + β) to obtain:

\begin{eqnarray} \alpha \,\tilde{b}_{3/2}^{(1)} \lbrace \alpha \rbrace \approx \frac{1}{\pi \,\beta ^2} + \mathcal {O}(\beta ^{-1}). \end{eqnarray}

(32)

Substituting this result back into equation (26), multiplying both sides by ı ω_i, while expanding to leading order in β, and noting that log (1 + β) ≈ β, we obtain the potential-free Schrödinger's equation of inclination dynamics within the disc

\begin{eqnarray} \imath \, \omega _i \, \frac{\mathrm{\partial} \eta }{\mathrm{\partial} t} = - \omega _i^2 \,\frac{\mathrm{\partial} ^2\eta }{\mathrm{\partial} \rho ^2}, \end{eqnarray}

(33)

where

\begin{eqnarray} \omega _i = \frac{n}{4\pi }\frac{m}{M} = \frac{\beta \,\Sigma _0 \sqrt{\mathcal {G}M a_0}}{2 M}. \end{eqnarray}

(34)

Thus, the mathematical description of secular angular momentum exchange within self-gravitating discs parallels that of a quantum particle confined to an infinite square potential well, with boundaries that extend from ρ = 0 to |$\rho =\mathcal {L}$|⁠. Note further that in this formalism, the fundamental frequency ω_i – a quantity related to the angular momentum budget of the disc – takes the place of ℏ⁶.

From intuitive grounds, one may speculate that because we have limited ourselves to only considering near-neighbour interactions in the treatment outlined above, ω_i must substantially underestimate the true wave propagation frequency within the disc. This is indeed the case, and the limitations arising from this approximation will be addressed quantitatively in Section 3.2.5. In the meantime, however, it is advantageous to temporarily ignore this problem and examine the character of the normal modes dictated by the above equation.

3.2.2 Boundary conditions

The specific character of the solution to equation (33) is determined by the imposed boundary conditions. In the well-known problem of a quantum-mechanical infinite square well, it is appropriate to employ the Dirichlet boundary conditions, η = 0 at ρ = 0 and |$\rho =\mathcal {L}$|⁠, since the wavefunction must vanish at the boundaries. On the contrary, for the problem at hand, there is no reason to require the orbital inclination to approach a particular value at the margins of the disc, and instead the boundary conditions must be deduced from the behaviour of the discrete system at the disc's inner and outer edges.

Let us consider the outer edge first. The disturbing function for the wire indexed by j = N is

\begin{eqnarray} \mathcal {R}_N = - \mathcal {B} \, \eta _N \,\eta ^*_N + \mathcal {B} \, \big (\eta _N \,\eta ^*_{N-1}+\eta _{N-1} \,\eta ^*_{N} \big ), \end{eqnarray}

(35)

where we have recycled the same approximations from the preceding sub-section. The resulting equation of motion has the form:

\begin{eqnarray} \frac{\text{d}\eta _N}{\text{d}t} = \imath \frac{\mathrm{\partial} \mathcal {R}_N}{\mathrm{\partial} \eta _N^*} = - \imath \, \mathcal {B} \big ( \eta _N - \eta _{N-1} \big ). \end{eqnarray}

(36)

The RHS of this expression is a finite difference approximation for the first derivative of η with respect to ρ. Correspondingly, recalling that |$\mathcal {B} \, (\delta \rho )\approx \omega _i/\beta$|⁠, we obtain the boundary condition at |$\rho =\mathcal {L}$|⁠:

\begin{eqnarray} \imath \frac{\mathrm{\partial} \eta }{\mathrm{\partial} \rho } = - \frac{\beta }{\omega _i} \frac{\mathrm{\partial} \eta }{\mathrm{\partial} t}. \end{eqnarray}

(37)

An identical procedure can be carried out for the inner boundary by change of indexes. Specifically, the appropriate disturbing function is

\begin{eqnarray} \mathcal {R}_1 = - \mathcal {B} \, \eta _1 \,\eta ^*_1 + \mathcal {B} \, \big (\eta _1 \,\eta ^*_{2}+\eta _{2} \,\eta ^*_{1} \big ). \end{eqnarray}

(38)

However, because the order of the wires is now reversed, the derivative appears with the opposite sign. The boundary condition at ρ = 0 is thus

\begin{eqnarray} \imath \frac{\mathrm{\partial} \eta }{\mathrm{\partial} \rho } = \frac{\beta }{\omega _i} \frac{\mathrm{\partial} \eta }{\mathrm{\partial} t}. \end{eqnarray}

(39)

With these specifications in place, we can now proceed to write down the solution to the governing equation (33).

3.2.3 Solution

Let us begin by noting that the potential-free Schrödinger equa-tion (33) itself is a diffusion equation in imaginary time. Physically, this means that the ‘diffusion’ process must conserve the phase-space volume of the underlying distribution function (i.e. the time evolution must be unitary; Stone 2000). This is perfectly sensible since in Section (2) we derived Schrödinger's equation from Hamilton's equations, which are themselves rooted in Liouville's theorem. Accordingly, rather than describing decay as its solution, equation (33) must instead be satisfied by standing waves – i.e. normal modes of a disc.

Let ℓ denote the index of a mode (characterized by frequency ω_ℓ). By separation of variables

\begin{eqnarray} \eta _{\ell }=c_\ell \,\exp (-\imath \,\omega _\ell \, t)\,\mathcal {I}_\ell , \end{eqnarray}

(40)

where c_ℓ is a constant, we obtain the familiar harmonic oscillator equation

\begin{eqnarray} \omega _\ell \, \mathcal {I}_\ell +\omega _i\frac{\mathrm{\partial} ^2\mathcal {I}_\ell }{\mathrm{\partial} \rho ^2}=0, \end{eqnarray}

(41)

subject to the boundary conditions

\begin{eqnarray} \frac{\mathrm{\partial} \mathcal {I}_\ell }{\mathrm{\partial} \rho } = -\beta \bigg (\frac{\omega _\ell }{\omega _i} \bigg ) \mathcal {I}_\ell \ \bigg |_{\rho =0}\quad \frac{\mathrm{\partial} \mathcal {I}_\ell }{\mathrm{\partial} \rho } = \beta \bigg (\frac{\omega _\ell }{\omega _i} \bigg ) \mathcal {I}_\ell \ \bigg |_{\rho =\mathcal {L}}. \end{eqnarray}

(42)

Note that for long-wavelength modes, where ω_ℓ does not exceed the fundamental frequency ω_i by a large margin, expressions (42) are well approximated by the rudimentary Neumann boundary conditions |$\mathrm{\partial} \mathcal {I}_\ell /\mathrm{\partial} \rho = 0$| at |$\rho =0,\mathcal {L}$|⁠, characteristic of a string with free ends.

The solution to equation (41) is

\begin{eqnarray} \mathcal {I}_\ell = \cos \bigg ( \sqrt{\frac{\omega _\ell }{\omega _i}} \rho \bigg ) - \beta \sqrt{\frac{\omega _\ell }{\omega _i}} \sin \bigg ( \sqrt{\frac{\omega _\ell }{\omega _i}} \rho \bigg ), \end{eqnarray}

(43)

with the quantization of ω_ℓ/ω_i stemming from the condition

\begin{eqnarray} \bigg [\frac{\beta \,\omega _\ell }{2\,\omega _i} -\frac{1}{2\, \beta } \bigg ]\sin \bigg ( \sqrt{\frac{\omega _\ell }{\omega _i}} \mathcal {L} \bigg )= \sqrt{\frac{\omega _\ell }{\omega _i}} \cos \bigg ( \sqrt{\frac{\omega _\ell }{\omega _i}} \mathcal {L} \bigg ). \end{eqnarray}

(44)

Although this expression does not admit a simple solution, in the β → 0 limit of a razor-thin disc, it simplifies to |$\sin \left(\sqrt{\omega _\ell /\omega _i}\,\mathcal {L}\right)=0$|⁠, which is trivially solvable. Adopting the razor-thin limit as a leading-order approximation, it is straightforward to derive a correction to the quantization condition:

\begin{eqnarray} \frac{\omega _\ell }{\omega _i}=\bigg (\frac{\ell \,\pi }{\mathcal {L}} \bigg )^2\bigg (1-4\,\frac{\beta }{\mathcal {L}} \bigg )+\mathcal {O}(\beta ^2). \end{eqnarray}

(45)

Fig. 3 depicts stationary states |$\mathcal {I}_\ell$| with ℓ = 0, 1, ..., 5.

Figure 3.

Normal modes of a razor-thin (β → 0) disc. Six low-frequency modes are shown, and are labeled by the appropriate index ℓ = 0, 1, ...5. Stationary modes of a geometrically thin disc with β = 0.01 appear nearly indistinguishable from those depicted in this figure.

Collecting the above results, the normal inclination modes of a self-gravitating disc are expressed as follows:

\begin{eqnarray} \eta _{\ell }&=& c_\ell \,\exp \bigg [-\imath \, \bigg (\frac{\ell \,\pi }{\mathcal {L}} \bigg )^2\bigg (1-4\,\frac{\beta }{\mathcal {L}} \bigg )\,\omega _i\, t \bigg ] \nonumber \\ &&\times \bigg [ \cos \bigg ( \frac{\ell \,\pi \,\rho }{\mathcal {L}}\sqrt{1-4\,\frac{\beta }{\mathcal {L}} } \bigg ) - \beta \,\frac{\ell \,\pi }{\mathcal {L}} \nonumber \\ &&\times \sqrt{1-4\,\frac{\beta }{\mathcal {L}}}\, \sin \bigg ( \frac{\ell \,\pi \,\rho }{\mathcal {L}}\sqrt{1-4\,\frac{\beta }{\mathcal {L}}} \bigg ) \bigg ]. \end{eqnarray}

(46)

We note that this solution corresponds to a special case of constant ω_i, which is in turn facilitated by a particular choice of the surface density profile (equation 16). Employing the β → 0 approximation once again, we obtain the simple expression:

\begin{eqnarray} \eta _{\ell }\approx c_\ell \,\exp \bigg [-\imath \, \bigg (\frac{\ell \,\pi }{\mathcal {L}} \bigg )^2\,\omega _i\, t \bigg ] \cos \bigg ( \frac{\ell \,\pi \,\rho }{\mathcal {L}}\bigg ). \end{eqnarray}

(47)

Within the framework of this solution, the coefficients c_ℓ are determined from Fourier decomposition of the initial conditions. Upon determination of these quantities, superposition of the eigenstates (47) fully describes the secular evolution of a self-gravitating disc. Physically, these eigenstates describe stationary nodal bending waves, with regression frequencies given by ω_ℓ. A number of low-frequency modes are shown in Fig. 4, as they appear in physical space. Note that unlike the quantum infinite square well, where the ground state of the wavefunction corresponds to ℓ = 1, the lowest index allowed within the context of our formalism is ℓ = 0, which describes a uniformly inclined, static disc.

Figure 4.

Six low-frequency inclination normal modes of a razor-thin disc, rendered in physical space. Note that owing to the definition of the logarithmic coordinate ρ, large-scale variations associated with higher mode indexes are concentrated towards small orbital radii. Here, mode amplitudes of c_ℓ = π/10 are assumed throughout.

3.2.4 Beyond nearest neighbours

The preceding analysis was carried out entirely within the limited framework of the nearest-neighbour interactions. Let us now quantify the validity of this assumption. We begin by extending the range of interactions, such that wire j can now interact with neighbours up to index j ± 2. The relevant disturbing function of the discrete system is then:

\begin{eqnarray} \mathcal {R}_j &=& B_{jj} \, \eta _j \,\eta ^*_j + B_{jj-1} \, \big (\eta _j \,\eta ^*_{j-1}+\eta _{j-1} \,\eta ^*_{j} \big ) \nonumber \\ &&+\, B_{jj+1} \, \big (\eta _j \,\eta ^*_{j+1}+\eta _{j+1} \,\eta ^*_{j} \big ) + B_{jj-2} \, \big (\eta _j \,\eta ^*_{j-2} \nonumber \\ &&+\, \eta _{j-2} \,\eta ^*_{j} \big ) + B_{jj+2} \, \big (\eta _j \,\eta ^*_{j+2}+\eta _{j+2} \,\eta ^*_{j} \big ). \end{eqnarray}

(48)

For the sake of this demonstration, let us assert that because β ≪ 1, all of the annuli in question are still in sufficient proximity to one another for the α → 1 limit to apply. Then, we can crudely assume that most of the variation among the coefficients B will stem from evaluation of the Laplace coefficient at different values of α. Thus, setting α = 1/(1 + β)^ν, where ν is an integer that indexes the interaction length (i.e. j:j ± ν coupling), we have

\begin{eqnarray} \alpha \,\tilde{b}_{3/2}^{(1)} \approx \frac{2}{(1+\nu ^2)\,\pi \,\beta ^2} + \mathcal {O}(\beta ^{-1}), \end{eqnarray}

(49)

which yields B_jj ± 2 ≈ (2/5) B_jj ± 1.

Following up on the same procedure as before, we deduce the equation of motion for wire j:

\begin{eqnarray} \frac{\text{d}\eta _j}{\text{d}t} &\approx & \imath \,\bigg ( \frac{n_j}{4\pi }\frac{m_j}{M}\frac{1}{\beta ^2} \bigg )\bigg [ -\bigg (2+\frac{4}{5} \bigg )\eta _j+\eta _{j+1}+\eta _{j-1} \nonumber \\ &&+\, \frac{2}{5}\left(\eta _{j-2} + \eta _{j+2} \right) \bigg ]. \end{eqnarray}

(50)

The new terms give rise to a higher order derivative, such that the continuum limit takes the form:

\begin{eqnarray} \frac{\mathrm{\partial} \eta }{\mathrm{\partial} t} \approx \imath \,\omega _i\,\bigg [ \frac{13}{5} \frac{\mathrm{\partial} ^2\eta }{\mathrm{\partial} \rho ^2}+\beta ^2\,\frac{2}{5}\frac{\mathrm{\partial} ^4\eta }{\mathrm{\partial} \rho ^4} \bigg ]. \end{eqnarray}

(51)

Continuing this procedure to j ± 3 yields

\begin{eqnarray} \frac{\mathrm{\partial} \eta }{\mathrm{\partial} t} \approx \imath \,\omega _i\,\bigg [\frac{22}{5} \frac{\mathrm{\partial} ^2\eta }{\mathrm{\partial} \rho ^2}+\beta ^2\,\frac{8}{5}\frac{\mathrm{\partial} ^4\eta }{\mathrm{\partial} \rho ^4} +\beta ^4\,\frac{1}{5}\frac{\mathrm{\partial} ^6\eta }{\mathrm{\partial} \rho ^6} \bigg ] \end{eqnarray}

(52)

and so on.

The general form of the above expression is reminiscent of boundary layer analysis of fluid mechanics (see e.g. section IV of Landau & Lifshitz 1959), and elucidates that as long as the characteristic wavelength of interest exceeds the gravitational softening length by a significant margin (which is an implicit assumption of the continuum limit), the contributions due to higher order derivatives can be safely neglected. In other words, the functional form of the eigenstates derived within the framework of the nearest neighbour interactions represents a perfectly adequate approximation to the true long-wavelength secular modes of a globally coupled system. Intuitively, the validity this approximation within the context of low wavenumber perturbations can be understood as a consequence of the fact that if 2π/k ≫ β, the complex field η does not change dramatically on scales comparable to the softening length, meaning that its relevant features are encapsulated within its curvature to a good approximation.

3.2.5 Collective effects

Although the proceeding analysis demonstrates that eigenmodes deduced from Schrödinger's equation (33) apply even when global effects are taken into account (we will check this assertion more thoroughly below), it also illuminates a key shortcoming of the local treatment of the dynamics. Namely, equations (51) and (52) clearly suggest that the eigenfrequencies obtained within the context of the nearest-neighbour approximation severely underestimate the true evolution rates of the disc's normal modes, since each incremental increase in interaction range significantly boosts the effective value of ω_i. Accordingly, let us now alleviate this limitation and compare the resulting formulae with a complete Lagrange–Laplace model of a disc.

The specific goals of the following calculations are two-fold. First, by initializing the full Lagrange–Laplace solution in eigenstates derived from the Schrödinger's equation and examining the resulting temporal evolution, we will determine how close these functions are to the true normal modes of the fully coupled system. That is, if the coherence of the initial condition is retained to a good approximation as time marches forward, the initial state is an eigenfunction of the system. On the other hand, if the inclinations oscillate with a large amplitude, then the assumed initial condition cannot be considered a stationary mode of the disc. Secondly, we will obtain the nodal regression rates corresponding to eigenfunctions (47), fully accounting for long-range gravitational interactions. For succinctness, we will obtain the analytic expression for the evolution rates first.

In order to properly describe collective effects within the system, let us return to a discrete representation of the dynamics of wire j, accounting for the potential of the full disc. The appropriate form for the disturbing function reads (Murray & Dermott 1999):

\begin{eqnarray} \mathcal {R}_j = \frac{1}{2}B_{jj} \, (p_j^2+q_j^2) + \sum _{k=1,j\ne k}^{N} B_{jk} (p_{j}\,p_{k}+q_{j}\,q_{k}). \end{eqnarray}

(53)

Recalling the definitions of p and q from equations (22), the expression for the longitude of ascending node is:

\begin{eqnarray} \Omega =\arctan \bigg (\frac{p}{q}\bigg ). \end{eqnarray}

(54)

Since we are interested in the the evolution rate of a specific mode, we can readily assume that a common phase (modulo sign) is shared by all constituent annuli within the disc, and its rate of change is independent of its value. Thus, without loss of generality we can set p = 0 everywhere, which yields a simplified expression for the derivative:

\begin{eqnarray} \frac{\text{d}\Omega _j}{\text{d}t}\bigg |_{p=0} &=& \frac{1}{q_j}\frac{\text{d}p_j}{\text{d}t}=\frac{1}{q_j}\frac{\mathrm{\partial} \mathcal {R}_j}{\mathrm{\partial} q_j} \nonumber \\ &=& \frac{1}{q_j}\Bigg (B_{jj}q_j + \sum _{k=1,j\ne k}^{N} B_{jk} q_{k} \Bigg ). \end{eqnarray}

(55)

In similarity with equations (25), the interactions coefficients take the form:

\begin{eqnarray} B_{jj} &=& -\frac{n_j}{4}\sum _{k=1,j\ne k}^{N}\frac{m_{k}}{M}\,\alpha _{jk}\bar{\alpha }_{jk}\,\tilde{b}_{3/2}^{(1)}\lbrace \alpha _{jk} \rbrace \nonumber \\ B_{jk} &=& \frac{n_j}{4}\frac{m_{k}}{M}\,\alpha _{jk}\bar{\alpha }_{jk}\,\tilde{b}_{3/2}^{(1)} \lbrace \alpha _{jk} \rbrace \end{eqnarray}

(56)

where we have adopted the notation of Murray & Dermott (1999) to set α_jk = (a_j/a_k), |$\bar{\alpha }_{jk}=(a_j/a_k)$| if j < k (external perturbation) and α_jk = (a_j/a_k), |$\bar{\alpha }_{jk}=1$| if j > k (internal perturbation). To this end, note that since we are no longer working within the confines of nearest neighbour interactions, expression (32) does not apply as a valid estimate of |$\tilde{b}_{3/2}^{(1)}$|⁠, meaning that Laplace coefficients must be evaluated at different values of α_jk explicitly⁷.

Substituting the β → 0 solution (47) for q, i.e.

\begin{eqnarray} q^\ell _j = \cos \bigg ( \frac{\ell \pi }{\mathcal {L}} \log \bigg (\frac{a_j}{a_0} \bigg ) \bigg ), \end{eqnarray}

(57)

and taking the continuous limit, the scaled precession rate at a given semimajor axis |$\hat{a}$| is given by

\begin{eqnarray} &&\bigg (\frac{\text{d}\hat{\Omega }}{\text{d}t} \bigg )_\ell = -\!\frac{\hat{n}}{4} \!\bigg ( \int _{a_{\rm {in}}}^{\hat{a}(1-\beta /2)} \!\frac{\bar{\Sigma } \bar{a}}{M} \!\bigg ( \frac{\bar{a}}{\hat{a}} \bigg ) \tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \frac{\bar{a}}{\hat{a}} \bigg \rbrace \text{d}\bar{a}\nonumber \\ &&- \int _{a_{\rm {in}}}^{\hat{a}(1-\beta /2)} \!\frac{\bar{\Sigma } \bar{a}}{M} \!\bigg ( \frac{\bar{a}}{\hat{a}} \bigg ) \tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \frac{\bar{a}}{\hat{a}} \!\bigg \rbrace \frac{\cos \bigg (\frac{\ell \pi }{\mathcal {L}} \log \bigg (\frac{\bar{a}}{a_0} \bigg )\bigg )}{\cos \bigg (\frac{\ell \pi }{\mathcal {L}} \log \!\bigg (\frac{\hat{a}}{a_0} \bigg )\bigg )} \text{d}\bar{a}\nonumber \\ &&- \int _{\hat{a}(1+\beta /2)}^{a_{\rm {out}}} \!\frac{\bar{\Sigma } \bar{a}}{M} \!\bigg ( \frac{\hat{a}}{\bar{a}} \bigg )^2 \!\tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \frac{\bar{a}}{\hat{a}} \!\bigg \rbrace \frac{\cos \bigg (\frac{\ell \pi }{\mathcal {L}} \log \!\bigg (\frac{\bar{a}}{a_0} \bigg )\bigg )}{\cos \!\bigg (\frac{\ell \pi }{\mathcal {L}} \log \!\bigg (\frac{\hat{a}}{a_0} \bigg )\bigg )} \text{d}\bar{a}\nonumber \\ &&+\int _{\hat{a}(1+\beta /2)}^{a_{\rm {out}}} \!\frac{\bar{\Sigma } \bar{a}}{M} \!\bigg ( \frac{\hat{a}}{\bar{a}} \bigg )^2 \tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \!\frac{\bar{a}}{\hat{a}} \!\bigg \rbrace \text{d}\bar{a}\bigg ), \end{eqnarray}

(58)

where the Laplace coefficients (equation 31) are evaluated setting α to the semimajor axis ratio inside the curly brackets, and the hatted and barred quantities are evaluated at |$\hat{a}$| and |$\bar{a}$|⁠, respectively. Assuming perfect rigidity of the mode under consideration, we assert that the true evolution rate is given by the average of the regression rates of the constituent annuli of the disc, weighted by their angular momentum deficit:

\begin{eqnarray} \bar{\omega }_\ell &=& \int _{a_{\rm {in}}}^{a_{\rm {out}}} \bigg (\frac{\text{d}\hat{\Omega }}{\text{d}t} \bigg )_\ell \frac{\hat{\Sigma } \hat{a}}{M} \sqrt{\mathcal {G}M \hat{a}} \cos ^2\bigg (\frac{\ell \pi }{\mathcal {L}} \log \bigg ( \frac{\hat{a}}{a_0} \bigg ) \bigg ) \text{d}\hat{a}\nonumber \\ &&\times \bigg ( \int _{a_{\rm {in}}}^{a_{\rm {out}}} \frac{\hat{\Sigma } \hat{a}}{M} \sqrt{\mathcal {G}M \hat{a}} \cos ^2\bigg (\frac{\ell \pi }{\mathcal {L}} \log \bigg ( \frac{\hat{a}}{a_0} \bigg ) \bigg ) \text{d}\hat{a}\bigg )^{-1}. \end{eqnarray}

(59)

In order to evaluate the validity of the above equations quantitatively, we compare the derived analytic theory to the full Lagrange–Laplace solution of a disc comprised of N = 100 discrete annuli⁸. Specifically, adopting a softening parameter of β = 10⁻³ as in (Kocsis & Tremaine 2011) and following the recipe outlined in Section 3.1, we initialize the disc in a ‘pure’ state given by equation (47). Obtaining the solution via the standard approach of matrix diagonalization (see e.g. chapter 7 of Morbidelli 2002), we examine the extent to which the disc remains in the prescribed stationary state, and compare the global evolution rate to that dictated by equation (59).

Fig. 5 shows the comparison between the discrete Lagrange–Laplace model of the disc and the analytical theory described above for ℓ = 0, 1, ..., 5. Specifically, the left-hand panel of the figure depicts the oscillations amplitude of inclination as a function of semimajor axis of the full Lagrange–Laplace system (red dots) over many secular precession time-scales. The right panel depicts temporal evolution of the node, wherein values corresponding to all orbital radii are plotted together as a function of time.

Figure 5.

Comparison between a full Lagrange-Laplace model of a self-gravitating disc (red dots) and analytical theory based upon the continuum treatment presented herein (purple lines). The left panels depict stationary ℓ = 0, 1, ..., 5 inclination modes of the disc as a function of semimajor axis, and the right panels show the corresponding rates of nodal regression. Over secular time-scales, the discrete solution experiences low-amplitude oscillations around its assumed initial conditions, but nevertheless remains approximately confined to the derived eigenstates of the system. In this example, we adopted a gravitational softening parameter of β = 10⁻³, and considered a disc composed of N = 100 wires. The disc was taken to orbit a M = 1 M_⊙ star, comprise |$m_{\rm {\text{disc}}} = 1\,{\rm M}_{\oplus}$| in total, and extend from a_in = 1 AU to a_out = (1 + β)^N ≈ 1.1 AU. Clearly, the agreement between analytical theory and the discrete model is satisfactory, especially for low-frequency modes.

Clearly, the eigenstates derived from Schrödinger's equation provide an adequate approximation to the global behaviour of the simulated system, since the full solution never drifts away from its initial condition (purple curves on the left-hand panel), and instead only experiences low-amplitude oscillations around the derived eigenfunctions. Moreover, the analytically computed evolution rates (purple lines on the right panel) also match the discrete solution well. Accordingly, we conclude that a linear superposition of normal modes of the form:

\begin{eqnarray} \eta _{\ell }\approx c_\ell \,\exp \bigg [\imath \, \bar{\omega }_\ell \, t \bigg ] \cos \bigg ( \frac{\ell \,\pi \,\rho }{\mathcal {L}}\bigg ). \end{eqnarray}

(60)

provides a simple and easily computable avenue towards evaluating the evolution of geometrically thin, self-gravitating discs on secular time-scales.

4 PERTURBED discS

In the last section, we analysed the dynamics of isolated self-gravitating discs, in the continuum limit of the Lagrange–Laplace secular theory. While examples of such cloistered systems surely exist in nature, real astrophysical discs often reside in phenomenologically rich, dynamic environments, and can be subject to substantial external (or internal) perturbations. Extending the formalism developed above to account for a select class of such perturbations is the primary goal of this section. For definitiveness, here we restrict ourselves to considerations of external gravitational forcing, although it is understood that other extrinsic effects (e.g. radiative stripping, turbulent infall of material, etc.) can also influence astrophysical discs’ long-term evolution. As in the proceeding section, we concentrate on strictly secular perturbations.

As already discussed above, secular interactions generically arise as a consequence of the orbit-averaged gravitational field of bound companions. Within the context of circumstellar discs, such companions can be binary stars (Batygin 2012; Lai 2014), or massive planets (Matsakos & Königl 2017; Nesvold et al. 2017). Under appropriate conditions, the ambient potential of a stellar birth cluster can also be modelled in this manner. In the Galactic centre, perturbations due to a molecular torus residing at ∼1.5–7 pc (Christopher, Scoville & Stolovy 2005) as well as the gravitational effects of a putative intermediate-mass black hole residing outside the stellar disc (Yu et al. 2007) can be treated within the secular framework. Given the well-known limitations of Lagrange–Laplace theory, we do not aim to provide a complete description of secular dynamics that covers every imaginable regime. Instead, here we focus on deriving a quantitative measure of the disc's tendency towards deformation in face of external excitations.

4.1 Secular forcing

Consider a perturbing companion of mass m^΄, residing on an arbitrarily inclined orbit with eccentricity e^΄ and semimajor axis a^΄ ≫ a_out. Envision that the angular momentum of the companion greatly exceeds the angular momentum budget of the disc, such that the back-reaction of the disc upon the companion can be neglected. Under these assumptions, we may orient the coordinate system to coincide with the plane of the companion's orbit (such that i^΄ = 0) and expand the gravitational potential in powers of the semimajor axis ratio⁹, (a/a^΄) (Kaula 1962). To quadrupole order¹⁰, the secular disturbing function associated with the orbit-averaged gravitational potential of the companion has the form (Mardling 2010):

\begin{eqnarray} \mathcal {R}&=&\frac{n}{4} \frac{m^{\prime }}{M}\bigg (\frac{a}{a^{\prime }}\bigg )^3\frac{1}{\sqrt{1-e^{\prime 2}}^3} \bigg [\bigg (1+\frac{3}{2}e^2 \bigg )\bigg (\frac{3\,\cos ^2(i)-2}{2} \bigg )\nonumber \\ &&+\, \frac{15}{4}e^2\,\sin ^2(i)\,\cos (2(\varpi -\Omega )) \bigg ]. \end{eqnarray}

(61)

The harmonic term that appears on the second line of equation (61) governs the Kozai-Lidov effect (Lidov 1962; Kozai 1962), a flavour of dynamical evolution that can manifest as coupled large-scale oscillations of eccentricity and inclination. A number of recent works (Martin et al. 2014; Fu et al. 2015; Nesvold et al. 2017; Zanazzi & Lai 2017) have studied the Kozai–Lidov effect pertaining to astrophysical discs, and have shown that it can indeed operate under specific conditions. Here, we assume that these conditions are not satisfied. Practically, this assumption can be translated to a restriction on the mutual disc-companion inclination (namely |$i<\arccos (\sqrt{3/5})$|⁠). However, we also note that even in discs whose inclination exceeds this critical value, the Kozai effect is typically suppressed due to self-induced regression of the argument of perihelion (Batygin et al. 2011; Xiang-Gruess & Papaloizou 2014). In such instances, this term can be readily averaged away and discarded from the disturbing function.

Defining the scaled Poincar|$\acute{\rm {e}}$| action-angle variables (Morbidelli 2002)

\begin{eqnarray} \tilde{Z}=\sqrt{1-e^2}\big (1-\cos (i) \big ) \approx 1-\cos (i)\quad \tilde{z}=-\Omega , \end{eqnarray}

(62)

and assuming that e ≪ 1 as before, the disturbing function (61) becomes

\begin{eqnarray} \mathcal {R}=\frac{n}{4}\frac{m^{\prime }}{M}\bigg (\frac{a}{a^{\prime }}\bigg )^3\bigg ( \frac{1-3\,\tilde{Z}+3\,\tilde{Z}^2/2}{\sqrt{1-e^{\prime 2}}^3} \bigg ). \end{eqnarray}

(63)

The resulting nodal regression rate is given by

\begin{eqnarray} \frac{\text{d}\Omega }{\text{d}t}=\frac{\mathrm{\partial} \mathcal {R}}{\mathrm{\partial} \tilde{Z}}\approx -\frac{3}{4}\frac{m^{\prime }}{M} \bigg (\frac{a}{a^{\prime }}\bigg )^3\frac{n}{\sqrt{1-e^{\prime 2}}^3}, \end{eqnarray}

(64)

where we have taken the mutual inclination to be small enough for the approximation cos (i) ≈ 1 to hold.

Equation (64) describes phase rotation, the rate of which is independent of η. Thus, carrying out the analysis outlined in the previous section, it is immediately evident that this effect can be incorporated into Schrödinger's equation as a potential term. Recalling the definition of the logarithmic coordinate, ρ (equation 15), we have:

\begin{eqnarray} \imath \, \omega _i\,\frac{\mathrm{\partial} \eta }{\mathrm{\partial} t} = - \omega _i^2\, \frac{\mathrm{\partial} ^2\eta }{\mathrm{\partial} \rho ^2} +\omega _i\,\omega ^{\prime }\,\exp (3\rho /2)\,\eta . \end{eqnarray}

(65)

Maintaining the definitions outlined in Section (3.1), in the above expression, we have scaled the perturbation constant by the parameters relevant to the inner edge of the disc:

\begin{eqnarray} \omega ^{\prime }=\frac{3}{4} \frac{m^{\prime }}{M} \sqrt{\frac{\mathcal {G}\,M}{a_{\rm {in}}^3}} \bigg ( \frac{a_{\rm {in}}}{a^{\prime }} \bigg )^3 \frac{1}{\sqrt{1-e^{\prime 2}}^3}. \end{eqnarray}

(66)

In the context of quantum mechanics, a Schrödinger's equation with an exponential potential constitutes a rudimentary model for molecular interactions (Atabek et al. 1982; Amore & Fernández 2008). Therefore, equation (65) implies that gravitational perturbations exerted upon a quasi-Keplerian astrophysical disc can be understood within a framework that is closely related to quantum scattering theory.

Employing the standard separation of variables (40) once again, we obtain the following equation for the stationary states:

\begin{eqnarray} -\!\omega ^{\prime }\exp (3\rho /2)\,\mathcal {I}_\ell +\omega _\ell \, \mathcal {I}_\ell +\omega _i\,\frac{\mathrm{\partial} ^2\mathcal {I}_\ell }{\mathrm{\partial} \rho ^2}=0. \end{eqnarray}

(67)

A solution to equation (67) was first identified¹¹ by Ma (1946), under the condition that |$\mathcal {I}_\ell$| must vanish at the origin (see also Bargmann 1949). Although a useful starting point, Ma's solution is not directly applicable to secular disc dynamics. Instead, the boundary conditions relevant to the problem at hand are obtained using the procedure outlined in Section 3.2.2:

\begin{eqnarray} \frac{\mathrm{\partial} \mathcal {I}_\ell }{\mathrm{\partial} \rho } &=& -\beta \bigg (\frac{\omega _\ell }{\omega _i} - \frac{\omega ^{\prime }}{\omega _i} \bigg )\mathcal {I}_\ell \ \bigg |_{\rho =0} \nonumber \\ \frac{\mathrm{\partial} \mathcal {I}_\ell }{\mathrm{\partial} \rho } &=& \beta \bigg (\frac{\omega _\ell }{\omega _i} - \frac{\omega ^{\prime }}{\omega _i}\exp (3\mathcal {L}/2) \bigg ) \mathcal {I}_\ell \ \bigg |_{\rho =\mathcal {L}}. \end{eqnarray}

(68)

Together, equations (67) and (68) admit an analytic solution, which is expressed in terms of the following quantity

\begin{eqnarray} \Lambda &=& \Gamma \bigg [1+\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}} \bigg ] \nonumber \\ &&\times\, \Bigg (\mathcal {X}_{\frac{4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }\,\exp (3\rho /2)}{\omega _i}} \bigg ) \nonumber \\ &&-\,\Bigg [\mathcal {X}_{-\frac{4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }\,\exp (3\rho /2)}{\omega _i}} \bigg ) \nonumber \\ &&\times\, \bigg (2\,\beta \,(\omega _\ell -\omega ^{\prime })\, \mathcal {X}_{\frac{4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) \nonumber \\ &&+\,\sqrt{\omega _i\,\omega ^{\prime }}\, \mathcal {X}_{\frac{4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}-1} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) \nonumber \\ &&+\,\mathcal {X}_{\frac{4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}+1} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) \bigg ) \Bigg ] \nonumber \\ &&\times\, \Bigg [ 2\,\beta \,(\omega _\ell -\omega ^{\prime })\,\mathcal {X}_{-\frac{4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) \nonumber \\ &&+\,\sqrt{\omega _i\,\omega ^{\prime }}\,\bigg ( \mathcal {X}_{\frac{-4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}-1} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) \nonumber \\ &&+\, \mathcal {X}_{1-\frac{4\,\imath }{3} \sqrt{\frac{\omega _\ell }{\omega _i}}} \bigg ( \frac{4}{3} \sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) \bigg ) \Bigg ]^{-1}\Bigg ), \end{eqnarray}

(69)

where Γ is the Gamma function and |$\mathcal {X}$| is the modified Bessel function of the first kind. Suitably, the expression for |$\mathcal {I}_\ell$| reads:

\begin{eqnarray} \mathcal {I}_\ell =\Lambda \,\exp \bigg [-\imath \arctan \bigg ( \frac{\Im (\Lambda )}{\Re (\Lambda )} \bigg ) \bigg ]. \end{eqnarray}

(70)

It is noteworthy that even though the expression for Λ itself is comprised of Bessel functions with imaginary indexes, they are summed together in such a way as to render |$\mathcal {I}_\ell$| purely real (Bocher 1892).

Unfortunately, the quantization condition that defines the frequencies ω_ℓ for equation (67) is exceedingly cumbersome. Thus, here we restrict ourselves to reporting its approximate form, which is attained by expanding the full expression in β to zeroth order:

\begin{eqnarray} &&\Bigg [\mathcal {X}_{-1-\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }\,\exp (3\mathcal {L}/2)}{\omega _i}} \bigg ) \nonumber \\ &&\times\, \Bigg (\mathcal {X}_{\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}-1}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) +\mathcal {X}_{\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}+1}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg )\Bigg ) \nonumber \\ &&+\,\mathcal {X}_{1-\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }\,\exp (3\mathcal {L}/2)}{\omega _i}} \bigg ) \nonumber \\ &&\times\, \Bigg (\mathcal {X}_{\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}+1}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg ) +\mathcal {X}_{\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}-1}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg )\Bigg ) \nonumber \\ &&-\,\Bigg ( \mathcal {X}_{\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}-1}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }\,\exp (3\mathcal {L}/2)}{\omega _i}} \bigg ) \nonumber \\ &&+\,\mathcal {X}_{\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}+1}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }\,\exp (3\mathcal {L}/2)}{\omega _i}} \bigg ) \Bigg ) \nonumber \\ &&\times\, \Bigg (\mathcal {X}_{-\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}-1}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg )+\mathcal {X}_{1-\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}}}\bigg (\frac{4}{3}\sqrt{\frac{\omega ^{\prime }}{\omega _i}} \bigg )\Bigg )\Bigg ]\nonumber \\ &&\times\, \Gamma \bigg [1-\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}} \bigg ]\, \Gamma \bigg [1+\frac{4\,\imath }{3}\sqrt{\frac{\omega _\ell }{\omega _i}} \bigg ] =0. \end{eqnarray}

(71)

We note that this condition is exact in the limit of a razor-thin disc, and for (ω^΄/ω_i) of order unity or less, the zeroes of this function are well approximated by |$\omega _\ell /\omega _i\simeq (\ell \,\pi /\mathcal {L})^2$|⁠.

Recall from Section 3 that the ground state of an unperturbed disc has ℓ = 0, which corresponds to η = const. everywhere. In presence of external forcing, this simple solution no longer holds, as the nodal regression induced by the companion can lead to a prominent distortion of the lowest energy state. Importantly, a controlling parameter that determines the extent of this distortion is |$\upsilon =(\omega ^{\prime }/\omega _i)\exp (3\mathcal {L}/2)$|⁠. Physically, this ratio corresponds to the relative magnitudes of the differential precession induced by the companion and nodal regression facilitated by self-gravity.

Forgoing collective effects for the moment, Fig. 6 depicts the lowest frequency mode (as dictated by equation 71) for various proportions of the aforementioned parameter υ. As may be intuitively expected, the lowest energy eigenfunctions become progressively deformed as the external forcing strength is increased. Moreover, it is worth noting that the extent of deviation away from the unperturbed solution becomes comparable to the amplitude of the unperturbed solution when υ ∼ 1. We will revisit this notion again below.

Figure 6.

Ground state of a perturbed self-gravitating disc, at various secular forcing strengths. As the magnitude of differential precession induced by the external perturber is increased (which corresponds to an increasing υ), the ground state of the disc becomes progressively more deformed. Accordingly, a flat, uniformly inclined system is replaced by a globally warped disc. Deformation amplitude of order the unperturbed mode amplitude is attained for |$\omega \sim \omega ^{\prime }\exp (3\mathcal {L}/2)$|⁠. This condition motivates the criterion for the gravitational rigidity of astrophysical discs presented in Section 4.2.

4.2 Gravitational rigidity

A key question that can be addressed using the framework derived in the previous subsection is: ‘under what conditions, will a disc develop significant structure due to external perturbations?’ However, to answer this question quantitatively, we must re-examine the extent to which the ground state of the perturbed system is deformed by external forcing, accounting for collective effects within the disc.

Fig. 6 shows that to a fair approximation, the functional form of the ground-state of a secularly forced disc can be represented as a superposition of ℓ = 0 and ℓ = 1 states of the unperturbed system:

\begin{eqnarray} \mathcal {I}_0\approx c_0 \bigg [1+\epsilon \,\cos \bigg ( \frac{2\,\pi \,\rho }{\mathcal {L}} \bigg ) \bigg ]. \end{eqnarray}

(72)

In this approximation, the deviation away from strict coplanarity, ε, fully encapsulates the strength of external perturbation. The key is thus to compute its magnitude directly from the parameters of the system.

Following the discussion outlined in Section 3.2.5, we recall that the evolution rate of the unperturbed ℓ = 1 mode, which we denote as |$\bar{\omega }_1$|⁠, is given by equations (58) and (59). Qualitatively, this frequency can be interpreted as the rate of angular momentum redistribution within the disc, facilitated by self-gravity. The resulting value should be compared with the rate of differential nodal regression induced upon the system by the external companion:

\begin{eqnarray} \bar{\omega }^{\prime }=\omega ^{\prime }\bigg (\bigg (\frac{a_{\rm {out}}}{a_{\rm {in}}}\bigg )^{3/2}-1 \bigg ). \end{eqnarray}

(73)

Guided by the findings of the previous subsection, we argue that the gravitational rigidity of the system is characterized entirely by the relative importance of these two effects, and assert that the dimensionless parameter ε is given by their ratio. Explicitly, in terms of physical parameters of the system, the expression reads:

\begin{eqnarray} \epsilon &=& \frac{\bar{\omega }^{\prime }}{\bar{\omega }_1} \!=\! \!\Bigg [ \int _{a_{\rm {in}}}^{a_{\rm {out}}} \!\frac{\mathcal {G}M}{4\hat{a}} \!\frac{\Sigma _0 \sqrt{a_0 \hat{a}}}{M} \cos ^2\bigg (\frac{\pi }{\mathcal {L}} \log \bigg ( \frac{\hat{a}}{a_0} \bigg ) \bigg ) \nonumber \\ &&\times\, \!\Bigg ( \int _{a_{\rm {in}}}^{\hat{a}(1-\beta /2)} \!\frac{\Sigma _0 \sqrt{a_0 \bar{a}}}{M} \!\bigg ( \frac{\bar{a}}{\hat{a}} \bigg ) \!\tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \frac{\bar{a}}{\hat{a}} \!\bigg \rbrace \text{d}\bar{a}\nonumber \\ &&-\, \int _{a_{\rm {in}}}^{\hat{a}(1-\beta /2)} \!\frac{\Sigma _0 \sqrt{a_0 \bar{a}} \bar{a}}{M} \!\bigg ( \frac{\bar{a}}{\hat{a}} \bigg ) \!\tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \frac{\bar{a}}{\hat{a}} \bigg \rbrace \frac{\cos \bigg (\frac{\pi }{\mathcal {L}} \log \bigg (\frac{\bar{a}}{a_0} \bigg )\bigg )}{\cos \bigg (\frac{ \pi }{\mathcal {L}} \log \!\bigg (\frac{\hat{a}}{a_0} \bigg )\bigg )} \text{d}\bar{a}\nonumber \\ &&-\, \int _{\hat{a}(1+\beta /2)}^{a_{\rm {out}}} \!\frac{\Sigma _0 \sqrt{a_0 \bar{a}}}{M} \!\bigg ( \frac{\hat{a}}{\bar{a}} \bigg )^2 \!\tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \frac{\hat{a}}{\bar{a}} \bigg \rbrace \!\frac{\cos \bigg (\frac{ \pi }{\mathcal {L}} \log \bigg (\frac{\bar{a}}{a_0} \bigg )\bigg )}{\cos \bigg (\frac{ \pi }{\mathcal {L}} \log \!\bigg (\frac{\hat{a}}{a_0} \bigg )\Bigg )} \text{d}\bar{a}\nonumber \\ &&+\, \int _{\hat{a}(1+\beta /2)}^{a_{\rm {out}}} \!\frac{\Sigma _0 \sqrt{a_0 \bar{a}}}{M} \!\bigg ( \frac{\hat{a}}{\bar{a}} \bigg )^2 \tilde{b}_{3/2}^{(1)} \!\bigg \lbrace \frac{\hat{a}}{\bar{a}} \!\bigg \rbrace \text{d}\bar{a}\Bigg )\text{d}\hat{a}\Bigg ] \nonumber \\ &&\times\, \Bigg [\frac{3}{4} \frac{m^{\prime }}{M} \sqrt{\frac{\mathcal {G}M}{a_{\rm {in}}^3}} \!\bigg ( \frac{a_{\rm {in}}}{a^{\prime }} \bigg )^3 \!\frac{1}{\sqrt{1-e^{\prime 2}}^3} \!\bigg (\bigg (\frac{a_{\rm {out}}}{a_{\rm {in}}}\bigg )^{3/2}-1 \bigg ) \nonumber \\ &&\times\, \!\int _{a_{\rm {in}}}^{a_{\rm {out}}} \!\frac{\Sigma _0 \!\sqrt{a_0 \hat{a}}}{M} \!\sqrt{\mathcal {G}M \hat{a}} \!\cos ^2\bigg (\frac{\pi }{\mathcal {L}} \log \!\bigg ( \!\frac{\hat{a}}{a_0} \bigg ) \!\bigg ) \text{d}\hat{a}\!\Bigg ]^{-1} \ll 1. \end{eqnarray}

(74)

We remind the reader that encoded in this equation (as well as in the form of the solution 69 and 70 itself) is the assumption of a |$\Sigma \propto 1/\sqrt{a}$| surface density profile, meaning that for a different radial dependence of Σ, the rigidity criterion would take on a quantitatively distinct form.

In order to check that equations (72) and (74) truly represent an adequate approximation to the lowest energy state of a secularly perturbed system, we once again employ the discrete Lagrange–Laplace model of the disc shown in Fig. 5. Subjecting the disc to secular perturbations arising from a m^΄ = 10⁻² M perturber residing on a circular orbit at (a_in/a^΄) = 1/3, we tune the disc mass to correspond to ε = 0.01, 0.1, 0.3, and 1, contrasting our analytic results with the full model at each iteration. A variant of Fig. 5 depicting this comparison is shown in Fig. 7. As is made evident by the essentially flat nature of the inclination ground state depicted on the left-hand panels of Fig. 7, and the coherence of ascending nodes shown on the right-hand panels, discs characterized by ε ≪ 1 maintain near-perfect gravitational rigidity, and show only minimal deformation in face of external forcing.

Figure 7.

Ground states of an externally perturbed self-gravitating disc, corresponding to different values of the deformation parameter, ε. As in Fig. 5, we consider the dynamics of a narrow disc with β = 10⁻³, composed of N = 100 massive wires. In the depicted calculations, however, this disc is perturbed by a m = 10⁻²M companion residing at a^΄ = 3 au on a circular orbit. In each panel, the disc is initialized in a state given by equation (72), and its mass tuned via equation (74) to give the desired value of ε. As expected from theoretical arguments, the gravitational rigidity of the system is ensured for ε ≲ 0.1. When this condition is satisfied, the disc's longitude of ascending node executes a slow regression with respect to the plane defined by the companion's orbit. For reference, the characteristic value of the Toomre's Q within the disc is quoted on the right-hand panel for each realization.

Qualitatively speaking, the capacity of a disc to retain its unperturbed form when ε ≪ 1 stems from nothing other than adiabatic invariance of the system's actions (Henrard 1993). Suitably, ε itself represents a readily computable measure of the extent to which the emergent quasi-integrals are conserved. Accordingly, the propensity of self-gravitating near-Keplerian discs to deform under extrinsic forcing is fully encapsulated within this parameter. Just as Toomre's Q ≳ 1 criterion corresponds to a disc that is stable against self-gravitational fragmentation, a perturbed disc with ε ≪ 1 is stable against external secular excitation. In other words, systems with ε substantially smaller than unity will behave as gravitationally rigid bodies.

5 DISCUSSION

Although the time-dependent Schrödinger equation is often thought of as a mathematical description that is reserved for quantum mechanics alone, its non-linear counterpart (equation 10) appears in numerous instances of classical physics. Examples of such contexts include small-amplitude (fluid) gravity waves, Langmuir plasma waves, and non-linear optics (Agrawal 2013; Bellan 2013). In this work, we have shown that the linear Schrödinger equation is also keenly relevant to the long-term evolution of astrophysical discs and constitutes a continuum description of the secular dynamics of stable self-gravitating systems. Remarkably, this result stems from the well-known Lagrange–Laplace perturbation theory of celestial mechanics, the origins of which were established over two centuries ago.

As a perturbative model, the Schrödinger equation lends easy access to orbital evolution that unfolds on time-scales that greatly exceed the longest Keplerian period of the system, but are nevertheless shorter than its physical lifetime. While the resulting characterization of long-term angular momentum transfer is not readily attainable by other (non-secular) methods, it is of key importance to a judicious interpretation of the observed structure of astrophysical discs. Simultaneously, we emphasize that because of the orbit-averaged framework within which the theory was constructed, it is not designed to capture the full richness of dynamical phenomena that may ensue within self-gravitating quasi-Keplerian systems¹².

By virtue of being exactly solvable, Schrödinger's equation provides an enthralling reduction of the full solution of the gravitational N-body problem. Within the context of this description, the eigenstates derived from Schrödinger's equation translate to normal modes of quasi-Keplerian discs, and the physical interpretation of radially travelling inclination pulses corresponds to propagation of nodal bending waves (Binney & Tremaine 1987). The derived formalism further allows us to formulate a measure of the gravitational rigidity of near-Keplerian systems, subject to extrinsic secular forcing. The resulting ε ≪ 1 stability criterion (equation 74) quantifies the response of self-gravitating discs to external perturbations, and compliments Toomre's Q ≳ 1 stability criterion for self-gravitational fragmentation.

In light of our model's simplicity, it is important to keep in mind that Schrödinger's equation only provides an adequate description for the angular momentum exchange within self-gravitating discs in a specific parameter regime (i.e. sufficiently low eccentricities and inclinations, small aspect ratio, etc.). That is to say that the reduction of 6N non-linear ordinary differential equations to a single linear partial differential equation necessarily entails some approximations. As a result, the outlined theory cannot serve as a general replacement for more complex numerical simulations (Touma et al. 2009). Instead, our analytical model can be meaningfully used to provide qualitative context for numerical results.

There exist numerous ways in which our model can be extended. As a first step, it is straightforward to generalize the derived formalism to account for arbitrary radial dependence of the surface density profile, and to lift the assumption of constant aspect ratio within the disc. While of considerable practical use, these generalizations will endow the fundamental frequency ω_i with a dependence upon the semimajor axes, compromising the applicability of simple solutions (such as that outlined in equation 47) to the governing Schrödinger-like equation.

Another curious extension of our model lies in incorporation of non-linearity. As already mentioned in Section 2, adding non-linear action terms to the governing Hamiltonian yields the non-linear variant of Schrödinger's equation in the continuum limit. Importantly, such terms arise naturally in the secular disturbing function at fourth order in eccentricities and inclination (Murray & Dermott 1999), and can yield significant coupling among the two degrees of freedom (see e.g. Batygin et al. 2015). Simultaneously, it is worth noting that adding non-linear terms to the governing Hamiltonian turns the description of self-gravitating discs into a chain of anharmonically coupled oscillators – a system closely related to the Fermi–Pasta–Ulam lattice (see Ford 1992 for a review). As a consequence, it is reasonable to speculate that the resulting evolution will be satisfied by breather solutions, and that Fermi–Pasta–Ulam recurrence can occur in self-gravitating discs.

Because we have considered purely gravitational coupling in this work, the obtained formulae are only applicable to particle discs, formally speaking. Although the strictly gravitational picture can in some cases serve as a good approximation of a hydrodynamic disc¹³ (see e.g. Fragner & Nelson 2010; Batygin 2012), a comprehensive theory for fluid disc evolution must include an account for internal forces. To this end, Ogilvie (2001, 2006) has demonstrated that the dynamics of a fluid disc subject exclusively to pressure forces (i.e. with negligible viscosity and self-gravity) can be understood within the framework of a non-linear variant of Schrödinger's equation (see also Barker & Ogilvie 2016). Correspondingly, fusing such formalism with the results developed herein may provide a viable way to extend our theory to simultaneously account for gravitational and hydrodynamic effects.

We conclude this work with a brief discussion of perturbations exerted upon astrophysical discs. As discussed in Section 4, the flavour of interactions considered in this manuscript represents only a subset of the full range of possibilities. For example, carrying out the expansion of the secular disturbing function (equation 61) to higher order (Naoz et al. 2013) in the semimajor axis ratio can introduce purely real forcing terms that are independent of η into the governing equations. Along similar lines of reasoning, inclination damping can be modelled as a purely real diffusion term within this framework. On the contrary, stochastic effects, such as those arising from passing stars or turbulent infall of disc material can be simulated by adding appropriately modulated noise to the exterior boundary condition (Spalding et al. 2014; Li & Adams 2015). Taken together, incorporation of aforementioned effects constitutes a intriguing avenue towards further development of the model and a more complete characterization of the long-term evolution of self-gravitating discs.

Acknowledgements

I am grateful to G. Laughlin, F. Adams, D. Stevenson, A. Morbidelli, C. Spalding, S. Michalakis, A. Kitaev, S. Tremaine, and J. Touma for illuminating discussions. Additionally, I would like to thank the anonymous referee for providing a thorough and insightful report that has led to a considerable improvement of the manuscript, as well as the David and Lucile Packard Foundation for their generous support.

Footnotes

1

Intriguingly, owing to their mildly collisional nature, Saturn's rings reside in between the ideal fluid and particle parameter regimes.

2

Well-known examples of extrasolar debris discs are found around Vega, Fomalhaut, β Pictoris, and ε Eridani.

3

Intuitively, softening the gravitational potential of a point-mass by a length h can be thought of as spreading the mass of the object over a Plummer sphere of radius h.

4

Coupling terms arise at fourth order in the perturbation series, and are responsible for secular chaos within the Solar system (Laskar 1994; Batygin et al. 2015).

5

The unscaled disturbing function is a measure of energy, like the Hamiltonian. Reduced by a characteristic angular momentum |$m\sqrt{\mathcal {G}\,M\,a}$|⁠, the scaled disturbing function becomes a measure of frequency (and has units of inverse time). For the remainder of the paper, we will refer to the scaled disturbing function as simply the disturbing function.

6

We remind the reader that assuming a functional form for the surface density profile other than that given by equation (16) would render ω_i a function of a.

7

For computational ease, it is useful to express |$\tilde{b}_{3/2}^{(1)}$| in terms of standard elliptic integrals (Hahn 2003).

8

Despite being analytic in nature, the actual solution of a Lagrange–Laplace system with the aid of computer algebra becomes computationally taxing for N substantially greater than ∼100.

9

Unlike the ‘literal’ expansion of the disturbing function employed within the framework of the Lagrange–Laplace theory (which assumes small eccentricities and inclinations while placing no restrictions upon the semimajor axis ratio), expansion of the potential in terms of the semimajor axis ratio assumes that (a/a^΄) ≪ 1 but places no restrictions on the orbital eccentricities and inclinations.

10

For an extensive exploration of secular dynamics governed by higher order terms, see e.g. Naoz et al. (2013); Li, Naoz & Holman (2014).

11

Through a change of variables, equation (67) can be turned into Bessel's equation.

12

For example, spiral density waves that cause the central object to be displaced from the centre of mass and rely upon the resulting indirect potential to further amplify their magnitude (Adams et al. 1989) are not captured in our model.

13

Generically, it is reasonable to neglect internal pressure forces if |$\sqrt{\mathcal {G}\,\Sigma \,a} \gg c_{\rm {s}}$|⁠, where c_s is the speed of sound (Tremaine 2001).