Eigenstates of Quasi-Keplerian Self-Gravitating Particle Discs
MMNRAS , 1–9 (2021)
Eigenstates of Quasi-Keplerian Self-Gravitating Particle Discs
Walker Melton, , ★ Konstantin Batygin, Deparment of Physics, Harvard University, 17 Oxford St., Cambridge, MA 02138 Division of Physics, Mathematics, and Astronomy, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA Division of Geological and Planetary Sciences, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA
Feb. 2021
ABSTRACT
Although quasi-Keplerian discs are among the most common astrophysical structures, compu-tation of secular angular momentum transport within them routinely presents a considerablepractical challenge. In this work, we investigate the secular small-inclination dynamics of arazor-thin particle disc as the continuum limit of a discrete Lagrange-Laplace secular per-turbative theory and explore the analogy between the ensuing secular evolution – includingnon-local couplings of self-gravitating discs – and quantum mechanics. We find the ‘quantum’Hamiltonian that describes the time evolution of the system and demonstrate the existence of aconserved inner product. The lowest-frequency normal modes are numerically approximatedby performing a Wick rotation on the equations of motion. These modes are used to quantifythe accuracy of a much simpler local-coupling model, revealing that it predicts the shape ofthe normal modes to a high degree of accuracy, especially in narrow annuli, even though itfails to predict their eigenfrequencies.
Key words: protoplanetary discs – Galaxy: nucleus – methods: analytical
Astrophysical discs are ubiquitous, and their long-term evolutionunder self-gravity constitutes one of the classic problems of dynam-ical astronomy. Having a generic formation channel characterizedby energy dissipation at constant angular momentum (Goldreich& Tremaine 1982), astrophysical discs exhibit significant diversity.Some discs, including active galactic nuclei and young protoplane-tary discs, are primarily composed of gaseous hydrogen and helium.In these fluid discs, self-gravity, viscosity, and pressure support allplay significant roles. Other systems, including debris discs anddiscs of stars around supermassive black holes, have no appreciableinternal pressure and undergo purely gravity-dominated evolution.Moreover, when the mass of the system is dominated by the cen-tral body, the trajectories of the individual particles resemble theclosed elliptical orbits of a binary system, rendering the dynamicsquasi-Keplerian.Although astrophysical discs are often envisioned to be pla-nar and perfectly circular, observations suggest that in reality, theycan possess considerable large-scale structure. One common fea-ture present in discs spanning a broad range of scales is warping- a process routinely attributed to an interplay between externalperturbations and internal angular momentum transfer (Binney &Tremaine 2008). For example, the origin of the warped stellar discat the center of the Milky Way has been ascribed to a mix of self-gravitational effects and torques from nearby star clusters (Kocsis &Tremaine 2011). Similarly, large-scale warps of young circumstellar ★ E-mail: [email protected] discs are often imagined to arise due to both self-interactions withinthe disc and external perturbations (Nesvold et al. 2016). Gravita-tional torqueing of protoplanetary discs may even be responsiblefor planetary spin-orbit misalignments (Bate et al. 2010; Spalding& Batygin 2014).In these and other quasi-Keplerian discs, self-gravitationalcoupling leads to long-period variations in the orbital elements.Although qualitatively simple, obtaining quantitative descriptionsof the long-term consequences of this process can present a consid-erable practical challenge. In this work, we attack one aspect of thisproblem and consider angular momentum transport on timescalesmuch longer than the orbital period but considerably shorter thanthe lifespan of the system. In particular, we focus on the secularinclination evolution of a dynamically cold particle disc.Adopting an orbit-averaged framework, we investigate the sec-ular evolution as the continuum limit of a discrete model of discdynamics, assuming that the degree of warping is small. To this end,Batygin (2018) has recently demonstrated that the dynamics cap-tured by the Lagrange-Laplace perturbative treatment of disc evo-lution mirror the structure of the free-particle Schrödinger equationwhen coupling between distant parts of the disc is neglected. In thislocal-coupling model, the dynamics are equivalent to a particle-in-a-box with Neumann boundary conditions; in terms of this analogy,propagating waves in the disc formed from a superposition of nor-mal modes would correspond with nodal bending waves (Binney &Tremaine 2008).The local-coupling model presents a compellingly simple pic-ture of the secular small-inclination evolution of a self-gravitatingrazor thin particle disc, but more work is needed to understand its © a r X i v : . [ a s t r o - ph . E P ] F e b W. Melton and K. Batygin limits of applicability. Notably, the secular eigen-frequencies pre-dicted by the naive local-coupling model are significantly lower thanthose seen in simulations including non-local coupling (Batygin2018). Additionally, because asymmetric gravitational couplingsbetween narrow rings in the disc are neglected, the local-couplingmodel does not in general conserve angular momentum. Despitethese shortcomings, numerical simulations indicate that they areclose to the true normal modes for thin annuli, such as those thatmay occur during pebble accretion in protoplanetary discs (Mor-bidelli, A. 2020). Even for wider discs, the equations of motion arealmost diagonal in the local-coupling normal mode basis.In this paper, we extend the local-coupling Schrödinger modelto include non-local couplings and show that a conserved innerproduct exists on the space of solutions to the equations of motion.By exploiting this, we can extract close approximations of the truesmall-inclination normal modes. We compare these normal modesto those predicted by the local-coupling model to understand itslimits of applicability, finding that the local-coupling normal modesprovide an excellent approximation to the true normal modes innarrow discs. The remainder of the paper is structured as follows.We first introduce the secular theory of a discretized self-gravitatingdisc and take the limit as the discretization becomes infinitely fineto obtain an effective model for the long-wavelength excitations ofa continuous disc in section 2. We then derive a conserved innerproduct on the space of solutions of the equations of motion. Thefirst four nontrivial lowest-frequency, or low-lying, normal modes ofthe disc are extracted by formally rotating time in the complex planeand these modes are compared to the local-coupling normal modesand the conventional Lagrange-Laplace normal modes found bydiagonalizing the equations of motion. Our results are summarizedin section 4.
To understand secular angular momentum transport in razor-thinparticle discs, we adopt an orbit-averaged Lagrange-Laplace per-turbative treatment. We partition the disc into a finite set of rings,interacting through a lowest-order expansion of the secular perturb-ing function (see chapters 6 and 7 of Murray & Dermott 1999), andtake the continuum limit by letting the spacing of the rings approachzero while holding fixed the extent of the disc. We then derive theequations of motion in the continuum limit and find a conservedinner product on the set of solutions.
We assume that the mass of the central body dominates the short-term dynamics of the system so that the orbits of individual particlesappear Keplerian over an orbital period, with a mean motion that iswell-approximated by Kepler’s third law 𝑛 = √︂ G 𝑀𝑎 , (1)where 𝑀 is the mass of the central body and 𝑎 is the semimajoraxis. Additionally, we assume that orbital eccentricities and relativeinclinations are small throughout the disc. This ensures that the ve-locity dispersion – which sets the vertical thickness of the system –is small compared to the Keplerian velocity at a fixed semi-majoraxis 𝑎 . Accordingly, we consider the aspect ratio of the disc 𝛽 = ℎ / 𝑎 as a small parameter and assume that it is constant throughout thedisc. Finally, we assume that the disc extends from 𝑎 in to 𝑎 out , anddefine L = log 𝑎 out / 𝑎 in . Realistic systems have L on the order ofa few; L exceeds 10 only in exceptional circumstances (Batygin2018). Location within the disc is parametrized by 𝜌 = log 𝑎 / 𝑎 ,where 𝑎 is taken to be 𝑎 in by convention. The inclination and lon-gitude of ascending node are parametrized by a complex coordinate 𝜂 such that 𝜂 = 𝑖 exp ( 𝜄 Ω )√ , (2)where 𝑖 is the inclination and Ω is the longitude of ascending node. 𝜂 𝑗 refers to the complex number parametrizing the inclination andlongitude of ascending node of the 𝑗 -th ring of the discretized disc, 𝑚 𝑗 the mass , 𝑛 𝑗 the mean motion, and 𝑎 𝑗 the semimajor axis. Aftertaking the continuum limit, we parametrize position within the ringby the variable 𝜌 = log 𝑎 / 𝑎 in , so that 𝜌 extends from 𝜌 = 𝜌 = L .In our analysis of the local-coupling model, we assume thatthe surface density profile is described by the power law: Σ ( 𝑎 ) = Σ ( 𝑎 in / 𝑎 ) / . (3)We focus on this case to make contact with the numerical workperformed in Batygin (2018), but the formalism described in thispaper applies more generally, and, when applicable, we will describethe generalization to arbitrary surface density profile.Moreover, by virtue of being orbit-averaged, we assume thatthe system is gravitationally stable, implying that 𝑚 disc (cid:47) 𝛽𝑀 / To correct the local-coupling model and understand its successesand deficiencies, a similar method must be developed taking non-local coupling into account. Accordingly, we will partition the discinto a finite set of rings, indexed by 𝑖 = , . . . , 𝑁 , and will eventuallytake 𝑁 → ∞ . In the ring-partitioned disc system, the disturbingfunction including all couplings is (Batygin 2018) R 𝑗 = 𝐵 𝑗 𝑗 𝜂 ∗ 𝑗 𝜂 𝑗 + ∑︁ 𝑖 ≠ 𝑗 𝐵 𝑗𝑖 ( 𝜂 𝑗 𝜂 ∗ 𝑖 + 𝜂 ∗ 𝑗 𝜂 𝑖 ) 𝐵 𝑗 𝑗 = − ∑︁ 𝑖 ≠ 𝑗 𝐵 𝑗𝑖 𝐵 𝑗𝑘 = 𝑛 𝑗 𝑚 𝑘 𝑀 𝑎 𝑘 𝑎 𝑗 ˜ 𝑏 ( ) / ( 𝑎 𝑘 / 𝑎 𝑗 ) 𝑘 ≤ 𝑗 𝑎 𝑗 𝑎 𝑘 ˜ 𝑏 ( ) / ( 𝑎 𝑗 / 𝑎 𝑘 ) 𝑘 > 𝑗 , (4)where ˜ 𝑏 ( ) / ( 𝛼 ) is a Laplace coefficient (Binney & Tremaine 2008).To prevent ˜ 𝑏 ( ) / ( 𝛼 ) from diverging as 𝛼 →
1, we soften the Laplacecoefficient by the disc aspect ratio 𝛽 (Hahn 2003):˜ 𝑏 ( 𝑗 ) 𝑘 ( 𝛼 ) = 𝜋 ∫ 𝜋 cos 𝑗𝜓 ( − 𝛼 cos 𝜓 + 𝛼 + 𝛽 ) 𝑘 . (5) This is a standard technique to account for the finite vertical extent ofthe disc; if we were performing an area integral 𝑑𝜌𝑑𝑧 , the result would befinite. When we pass to a description ignoring the ˆ 𝑧 axis, we must includethis softening parameter to avoid unphysical singularities.MNRAS , 1–9 (2021) igenstates of Quasi-Keplerian Discs Because 𝐵 𝑗𝑘 ∝ 𝑚 𝑘 for 𝑗 ≠ 𝑘 , we can take the continuum limit ofthe disturbing function to obtain R( 𝜌 ) = ∫ 𝐶 ( 𝜌, 𝜌 (cid:48) ) (cid:2) 𝜂 ( 𝜌 (cid:48) ) ∗ 𝜂 ( 𝜌 ) + 𝜂 ( 𝜌 (cid:48) ) 𝜂 ( 𝜌 ) ∗ (cid:3) 𝑑𝑚 (cid:48) − (cid:18)∫ 𝐶 ( 𝜌, 𝜌 (cid:48) ) 𝑑𝑚 (cid:48) (cid:19) 𝜂 ( 𝜌 ) ∗ 𝜂 ( 𝜌 ) 𝐶 ( 𝜌, 𝜌 (cid:48) ) = 𝑛 ( 𝜌 ) 𝑀 𝑒 𝜌 (cid:48) − 𝜌 ˜ 𝑏 ( ) / ( 𝑒 𝜌 (cid:48) − 𝜌 ) 𝜌 (cid:48) ≤ 𝜌𝑒 𝜌 − 𝜌 (cid:48) ˜ 𝑏 ( ) / ( 𝑒 𝜌 − 𝜌 (cid:48) ) 𝜌 (cid:48) > 𝜌 . (6)Applying the equation of motion 𝜄 𝜕𝜂𝜕𝑡 = − 𝜕 R( 𝜌 )/ 𝜕𝜂 ( 𝜌 ) ∗ , wefind that 𝜄 𝜕𝜂 ( 𝜌 ) 𝜕𝑡 = ∫ 𝐶 ( 𝜌, 𝜌 (cid:48) )( 𝜂 ( 𝜌 ) − 𝜂 ( 𝜌 (cid:48) )) 𝑑𝑚 (cid:48) = ∫ 𝜋𝑎 𝑖𝑛 𝑒 𝜌 (cid:48) Σ ( 𝑎 𝑖𝑛 𝑒 𝜌 (cid:48) ) 𝐶 ( 𝜌, 𝜌 (cid:48) )( 𝜂 ( 𝜌 ) − 𝜂 ( 𝜌 (cid:48) )) 𝑑𝜌 (cid:48) = ∫ 𝐵 ( 𝜌, 𝜌 (cid:48) )( 𝜂 ( 𝜌 ) − 𝜂 ( 𝜌 (cid:48) )) 𝑑𝜌 (cid:48) . (7)This can be written in the form of a general time-dependentSchrödinger equation with Hamiltonian H , in which H 𝜂 = ∫ 𝐵 ( 𝜌, 𝜌 (cid:48) )( 𝜂 ( 𝜌 ) − 𝜂 ( 𝜌 (cid:48) )) 𝑑𝜌 (cid:48) . (8)If we expand 𝜂 ( 𝜌 (cid:48) ) in a power series around 𝜌 , we find that H = − ∞ ∑︁ 𝑘 = 𝑐 𝑘 ( 𝜌 ) 𝜕 𝑘 𝜕𝜌 𝑘 𝑐 𝑘 ( 𝜌 ) = 𝑘 ! ∫ 𝐵 ( 𝜌, 𝜌 (cid:48) )( 𝜌 (cid:48) − 𝜌 ) 𝑘 𝑑𝜌 (cid:48) . (9)We should stress that we have so far only assumed that the relativeinclination is perturbatively small throughout the disc and that therelevant dynamics are secular. We have not assumed anything aboutthe coupling between different parts of the disc; hence, the equationof motion presented above is equivalent to the Lagrange-Laplaceequations of motion written in terms of the complex variable 𝜂 .We obtain the local-coupling equation by keeping only the termproportional to 𝜕 / 𝜕𝜌 and averaging over the spatial dependenceof 𝑐 ( 𝜌 ) : 𝜄 𝜕𝜂𝜕𝑡 = H 𝐿𝐶 𝜂 = − 𝑐 𝜕 𝜂𝜕𝜌 𝑐 = L ∫ L 𝑐 ( 𝜌 ) 𝑑𝜌. (10)In the interior of the disc, 𝑐 𝑘 is small for odd 𝑘 due to theapproximate symmetric coupling, while 𝑐 𝑘 remains significant foreven 𝑘 , so that the dominant contribution to the equations of motionhas 𝑘 =
2. While 𝑐 𝑘 tends to drop in magnitude as 𝑘 increases,it is worth noting that the higher-order contributions can becomesignificant for states with large wavenumbers. To extend the analogy between the small-inclination dynamics of aself-gravitating particle disc and the structure of quantum mechan-ics, we need not only the emergence of a general time-dependentSchrödinger equation but also the presence of a conserved innerproduct. Many of the tools of quantum mechanics rely on the ex-istence of a set of orthogonal eigenmodes with real oscillation fre-quencies, which is guaranteed by the conservation of some inner product (cid:104) 𝜓 | 𝜙 (cid:105) . This will be true if the Hamiltonian is Hermitian withrespect to the given inner product. Correspondingly, in this sectionwe first show that there exists an inner product such that (cid:104) 𝜂 | 𝜂 (cid:105) isconserved by requiring that angular momentum is conserved; then,we show that (cid:104) 𝜂 | 𝜈 (cid:105) must be conserved as well for any solutions tothe equations of motion 𝜂 and 𝜈 .Consider the angular momentum of a small ring of mass 𝑑𝑚 and semi-major axis 𝑎 . Its ˆ 𝑧 angular momentum deficit is, in thesmall inclination approximation, √ 𝐺 𝑀𝑎𝑖 𝑑𝑚 / 𝑧 -angular momentumdeficit is ∫ √ 𝐺 𝑀𝑎𝑖 𝑑𝑚 ∝ ∫ 𝜋𝑎 Σ ( 𝑎 )√ 𝑎𝑖 𝑑𝑎 ∝ ∫ 𝑎 Σ ( 𝑎 )√ 𝑎𝑖 𝑑𝑎 = (cid:104) 𝜂 | 𝜂 (cid:105) . (11)For the choice of disc parameters above, Σ ( 𝑎 )√ 𝑎 is constant, andthe inner product is given by (dropping an insignificant constant) (cid:104) 𝜂 | 𝜂 (cid:105) = ∫ 𝑎𝜂 ∗ 𝜂𝑑𝑎. (12)For these dynamics to mirror completely the structure of quan-tum mechanics, it must be that (cid:104) 𝜂 | 𝜈 (cid:105) is conserved for any twosolutions 𝜂 and 𝜈 . Indeed, this is the case for stable dynamics suchas those described by the Lagrange-Laplace model. Because the‘matrix elements’ 𝐵 ( 𝜌, 𝜌 (cid:48) ) of H are real, (H 𝜂 ) ∗ = H [ 𝜂 ∗ ] . Thus,if 𝜂 is a normal mode with frequency 𝜔 , 𝜂 ∗ must be a normalmode with frequency 𝜔 ∗ . For stable dynamics, 𝜔 = 𝜔 ∗ so thatthe amplitude of normal modes does not grow with time.. Hence, 𝜂 ∗ + 𝜂 is also a normal mode, and the normal modes can be takento be real functions of the coordinate 𝜌 . Suppose that 𝜂 𝑎 and 𝜂 𝑏 are two real and distinct eigenfunctions. Then, (cid:104) 𝜂 𝑎 + 𝜂 𝑏 | 𝜂 𝑎 + 𝜂 𝑏 (cid:105) is conserved, which, by the conservation of (cid:104) 𝜂 | 𝜂 (cid:105) , implies that (cid:104) 𝜂 𝑎 | 𝜂 𝑏 (cid:105) = − (cid:104) 𝜂 𝑏 | 𝜂 𝑎 (cid:105) = − (cid:104) 𝜂 𝑎 | 𝜂 𝑏 (cid:105) ∗ . Therefore, (cid:104) 𝜂 𝑎 | 𝜂 𝑏 (cid:105) = 𝛾𝜄 for 𝛾 ∈ R . (cid:104) 𝜂 𝑎 | 𝜂 𝑏 (cid:105) must also be real, however, so 𝛾𝜄 ∈ R . The only wayboth conditions can be satisfied is if (cid:104) 𝜂 𝑎 | 𝜂 𝑏 (cid:105) = 𝜂 = (cid:205) 𝑎 𝛼 𝑎 𝜂 𝑎 , 𝜈 = (cid:205) 𝑎 𝛽 𝑎 𝜂 𝑎 . Then, (cid:104) 𝜂 | 𝜈 (cid:105) = ∑︁ 𝑎,𝑏 𝛼 ∗ 𝑎 𝛽 𝑏 (cid:104) 𝜂 𝑎 | 𝜂 𝑏 (cid:105) = ∑︁ 𝑎 𝛼 ∗ 𝑎 𝛽 𝑎 (cid:104) 𝜂 𝑎 | 𝜂 𝑎 (cid:105) . (13)Because this quantity is independent of time, the inner product isconserved, provided that the model conserves angular momentum.We will provide a direct proof of the conservation of this innerproduct in Appendix A.For a general surface density profile Σ ( 𝑎 ) , it is easy to extendthe analysis above to show that the inner product (cid:104) 𝜂 | 𝜈 (cid:105) ∝ ∫ 𝜂 ∗ ( 𝑎 ) 𝜈 ( 𝑎 ) 𝑎 / Σ ( 𝑎 ) 𝑑𝑎 (14)is conserved. This completes the analogy with quantum mechanics:not only do the dynamical equations take the form of a generalizedtime-dependent Schrödinger equation, the space of solutions of theequations of motion possesses a conserved inner product. This im-plies that the Hamiltonian is Hermitian with respect to the innerproduct above.It is worth noting that Equation 14 is not the same inner productas that preserved by the local coupling model. The local couplingmodel preserves the inner product (cid:104) 𝜂 | 𝜈 (cid:105) ∝ ∫ L 𝜂 ∗ ( 𝜌 ) 𝜈 ( 𝜌 ) 𝑑𝜌 MNRAS , 1–9 (2021)
W. Melton and K. Batygin
Figure 1.
A quantitative measure of the non-conservation of angular mo-mentum wihtin the framework of the local coupling model. The quantity 𝑟 ( L) is shown for a disc with an inverse square root density profile from L = log 𝑎 out / 𝑎 in = L =
3. Note that the degree of angular momentumnon-conservation is only small for
L (cid:28)
1, implying that the local model isapplicable to annuli but not extended discs. which is the sum of the angular momentum deficit of each infinites-imal ring with a nontrivial weighting. The source of this differenceis the assumption that the coupling was symmetric when derivingthe local coupling model. As the true normal modes of the discare orthogonal with respect to the inner product given by Equa-tion 14 rather than that preserved by the local coupling model,this inner product is not preserved by the true, non-local dynam-ics. As a benchmark of the failure of the local coupling model toconserve angular momentum, we compute a measure of the varia-tion of angular momentum if the local coupling model was correct.Since we focus on the long-wavelength normal modes of the disc,we examine the failure of angular momentum to be conserved instates dominated by the longest-wavelength normal modes. If thelocal coupling model was correct, we could consider a system withstate 𝜂 ( 𝜌, 𝑡 ) = 𝑐 ( 𝜌 ) 𝑒 − 𝜄𝜔 𝑡 + 𝑐 ( 𝜌 ) 𝑒 − 𝜄𝜔 𝑡 , letting 𝑐 ℓ ( 𝜌 ) be thelocal coupling modes. Then, to linear order, the angular momentumof this state is (cid:104) 𝜂 | 𝜂 (cid:105) = (cid:104) 𝑐 | 𝑐 (cid:105) + (cid:104) 𝑐 | 𝑐 (cid:105) + (cid:104) 𝑐 | 𝑐 (cid:105) 𝑒 𝜄 ( 𝜔 − 𝜔 ) 𝑡 +(cid:104) 𝑐 | 𝑐 (cid:105) 𝑒 𝜄 ( 𝜔 − 𝜔 ) 𝑡 = (cid:104) 𝑐 | 𝑐 (cid:105) + (cid:104) 𝑐 | 𝑐 (cid:105) + (cid:104) 𝑐 | 𝑐 (cid:105) cos ( 𝜔 − 𝜔 ) 𝑡 ,so that the width of the variation of angular momentum is 2 (cid:104) 𝑐 | 𝑐 (cid:105) .We thus define 𝑟 (L) = | (cid:104) 𝑐 | 𝑐 (cid:105) | √︁ (cid:104) 𝑐 | 𝑐 (cid:105) (cid:104) 𝑐 | 𝑐 (cid:105) (15)which provides a measure of the extent of variation around theaverage angular momentum that is insensitive to the normalizationof the inner product and the normalization of the solutions 𝑐 , 𝑐 .For a disc with an inverse square root surface density profile, 𝑟 (L) = L √︁ − L /( 𝜋 + L ) coth L 𝜋 + L (16)which is plotted in Figure (1). The error is small (on the order of5%) for L (cid:47) .
1, so the local coupling model may be expected tohold for narrow annuli, but becomes progressively inaccurate as L increases.s The non-local equations of motion found in the previous sectionare not amenable to analytical techniques in the same way as the local model. Variational methods can be applied, but without a clearunderstanding of the space of states described by these dynamics,we risk finding false solutions. Additionally, these results may besensitive to the choice of family of trial functions. An approach thatextracts the normal modes without requiring assumptions abouttheir form would be far more useful.Such a method can be found by rotating time in the com-plex plane, turning the oscillatory solutions of the true equations ofmotion into exponentially damped solutions. By applying the imag-inary time equations of motion and the orthogonality rules derivedabove, we can successively extract the low-lying normal modes ofthe disc. In this section, the imaginary time formalism is describedand extended to the equations of motion for a disc composed of alarge but finite number of rings.
In real time, the equations of motion are 𝜄 𝜕𝜂𝜕𝜏 = H 𝜂𝜂 ( 𝜌, 𝑡 ) = 𝑒 − 𝜄 H 𝑡 𝜂 ( 𝜌, ) . (17)The solutions take the familiar form 𝜂 ( 𝜌, 𝑡 ) = (cid:205) 𝑘 𝑐 𝑘 𝜂 𝑘 ( 𝜌 ) 𝑒 − 𝜄𝜔 𝑘 𝑡 .When we perform a Wick rotation, sending 𝑡 → − 𝜄𝜏 , the equationsof motion take the form 𝜕𝜂𝜕𝜏 = −H 𝜂𝜂 ( 𝜌, 𝜏 ) = 𝑒 −H 𝜏 𝜂 ( 𝜌, ) . (18)These equations have solutions of the form 𝜂 ( 𝜌, 𝜏 ) = (cid:205) 𝑘 𝑐 𝑘 𝜂 𝑘 ( 𝜌 ) 𝑒 − 𝜔 𝑘 𝜏 where 𝜂 𝑘 ( 𝜌 ) is the shape of the 𝑘 -th normalmode. These solutions are exponentially damped, with higher-frequency normal modes decaying more quickly than lower-frequency modes. At large 𝜏 , 𝜂 ( 𝜌, 𝜏 ) will be dominated by thenormal mode of lowest frequency whose coefficient in the modeexpansion above is non-zero.By choosing 𝜂 ( 𝜌, ) such that (cid:104) 𝜂 ( )| 𝜂 𝑘 (cid:105) = 𝑘 < ℓ and (cid:104) 𝜂 ( )| 𝜂 ℓ (cid:105) ≠ 𝜂 ( 𝜌, 𝜏 ) will be dominated by 𝜂 ℓ ( 𝜌 ) for large 𝜏 .Because the ground state 𝜂 ( 𝜌 ) = 𝜂 ( 𝜌 ) , the state of the disc will be described by a finite setof complex numbers 𝜂 𝑗 , with 𝑗 = . . . 𝑁 . The Hamiltonian is alinear function of the space of states with coefficients given by thediscrete counterpart of equation 4.3, and the inner product is (cid:104) 𝜂 | 𝜈 (cid:105) = ∑︁ 𝑗 𝑚 𝑗 √ 𝑎 𝑗 𝜂 ∗ 𝑗 𝜈 𝑗 (19)so that (cid:104) 𝜂 | 𝜈 (cid:105) is proportional to the ˆ 𝑧 angular momentum deficit. Wewill use the equations of motion for this discretized model to extractthe normal modes of the disc. MNRAS , 1–9 (2021) igenstates of Quasi-Keplerian Discs Table 1.
The periods and best-fit wavenumbers of the first numericallyextracted modes. i 𝑘 𝑇 𝑖 = 𝜋𝜔 𝑖 (years)1 24.1748 39472 61.2425 16703 90.6545 10414 120.424 756 Σ ∝ 𝑎 − / disc A distinct prediction made by the local coupling model outlined inBatygin (2018) and recalled in section 2.2 is that the normal modesof the disc have the form: 𝜂 𝑘 ( 𝜌 ) = cos 𝑘𝜋𝜌 L = cos (cid:18) 𝑘𝜋 L log 𝑎𝑎 (cid:19) . (20)Correspondingly, in order to quantify the accuracy of the local-coupling model, we simulated a trial disc numerically. For defini-tiveness, we adopted the same example as that considered by Batygin(2018). The disc had total mass 1 𝑀 ⊕ , orbited a star of 1 𝑀 (cid:12) , andextended from 𝑎 = 𝑎 = . Σ ( 𝑎 ) ∝ 𝑎 − / as in section2.1. The aspect ratio of the disc was 𝛽 = ℎ / 𝑎 = − . In the simu-lation, the disc was partitioned into a series of 𝑁 =
100 rings, eachof which coupled gravitationally to every other ring.The imaginary time equations of motion 𝜕𝜂 / 𝜕𝜏 = −H 𝜂 weresolved using a conventional 4th-order Runge-Kutta integration.The initial conditions were chosen so that 𝜂 ( 𝜌, ) was close tocos 𝑘𝜋𝜌 /L for 𝑘 = , , ,
4, with contributions from the lower-frequency normal modes removed using Gram-Schmidt orthogo-nalization. To prevent the desired state from becoming dampedbelow machine precision, the Hamiltonian was shifted by a termproportional to the identity so that the frequency of the desiredmode was close to 0. The calculations were run for several thousandyears of simulation time, by which point higher-frequency modeshad damped to the point that they were negligible compared to themode of interest. After the state was extracted from the simula-tion, Gram-Schmidt orthogonalization was applied again to removeresidual contributions from low-lying modes that grew exponen-tially. The modes derived were then used as the input for a real-timesimulation; upon completing this exercise, no significant evolutionof the magnitude of the inclination was observed, indicating thatthese modes are indeed very close to the true normal modes of thesystem. The results of these simulations are presented in Figure 2.The periods of these normal modes were then estimated fromthe real time simulations; the results are presented in Table 1 alongwith best fit wavenumbers 𝑘 . We reiterate that these normal modeswere not derived by assuming a specific form of the normal modesbeyond the initial condition 𝜂 ( 𝜏 = ) = cos 𝑘𝜋𝜌 /L , and the resultsare generally insensitive to the initial choice as long as the correctorthogonalization procedures are followed. Hence, we use thesenormal modes to test the accuracy of the local-coupling normalmode approximations, which given by Equation 20.While this may seem a numerically tedious way of extractingthe normal modes, convergence to the true normal mode is fast asthe difference to the true normal modes falls exponentially in time.As such, when it is computationally simpler to simulate a systemover a small number of periods than to numerically diagonalize Table 2.
The periods in years from explicit diagonalization and the imaginarytime simulation i 𝑇 𝑖 (Imaginary Time) 𝑇 𝑖 (Explicit Diagonalization)1 3947 39472 1670 16713 1041 10414 756 757 the Lagrange-Laplace equations of motion (as becomes progres-sively the case with increasing 𝑁 ), this will provide a useful way ofnumerically extracting low-lying modes of the disc. Now that we have used the Wick-rotated equations of motion tonumerically evaluate the normal modes, we compare to the resultsfound by diagonalizing the Lagrange-Laplace equations of motion.Specifically, we will show that the normal modes for the 𝑁 = 𝑑𝜂 𝑗 𝑑𝑡 = − 𝜄𝐵 𝑗𝑘 𝜂 𝑘 (21) 𝐵 𝑗𝑘 was explicitly diagonalized and its eigenvalues and eigenvectorsexplicitly computed. The periods predicted by explicit diagonaliza-tion of 𝐴 𝑖 𝑗 are compared to those predicted from the imaginary-time simulations in Table 2. The agreement between the predictionsderived using the imaginary time formalism are clearly excellentapproximations to those derived by explicit diagonalization of thenormal modes.The normal modes derived from the imaginary time simulationand from diagionalizing the equations of motion, normalized sothat 𝜂 ( ) =
1, are shown in Figure 3 below. Differences betweenthe exact modes of the discrete system and those derived from theimaginary time procedure were on the order of 10 − to 10 − asa fraction of the exact mode, with particularly strong agreementfor the two lowest normal modes. Better agreement could be foundby running the imaginary time simulation for a longer period oftime. We conclude that the modes derived using the imaginarytime simulations are in excellent agreement with the modes foundby explicitly diagonalizing the equations of motion, validating ourmethod, that the Gram-Schmidt orthogonalization is sufficientlyreliable, and indicating that we can use them to quantify the accuracyof the local coupling modes. Although the curves delineated in Figure 2 appear close to theprediction of the local model (Eq. 20), their correspondence is notexact. But just how inexact are they? Having strong estimates of thetrue normal modes of the disc, we can now numerically measurethe accuracy of these local-coupling modes.Let us consider a few different measures that are sensitive todifferent aspects of the local-coupling model. The first, based on theFourier series expansion of the true normal mode, is sensitive to boththe sinusoidal nature of the mode as well as the boundary conditionson the disc. Others, based on the accuracy of sinusoidal fits to thetrue normal mode, are insensitive to the boundary conditions.
MNRAS000
MNRAS000 , 1–9 (2021)
W. Melton and K. Batygin
Figure 2.
The first four non-constant normal modes of the local coupling model (right) and the corresponding numerically extracted modes for the fully coupledmodel with the same disc parameters (left). The red line is the lowest-frequency non-constant mode ( 𝑘 = 𝑘 = 𝑘 = 𝑘 = 𝑁 =
100 equally spaced, interacting wires. The aspect ratio is taken to be 𝛽 = .
001 throughout. The cumulative disc mass isequal to that of the Earth and the mass of the central body is equal to that of the sun.
Figure 3.
The first four normal modes as extraced by the imaginary time formalism (red) and direct diagonalization (blue). The two are in excellent agreement.MNRAS , 1–9 (2021) igenstates of Quasi-Keplerian Discs If we assume that 𝜂 ℓ ( 𝜌 ) is the true normal mode most closelyapproximated by cos ℓ𝜋𝜌 /L , we can expand 𝜂 ℓ ( 𝜌 ) as a Fourierseries: 𝜂 ℓ ( 𝜌 ) = 𝜂 + ∑︁ 𝑘 𝑐 𝑘 cos 𝑘𝜋𝜌 /L + ∑︁ 𝑞 𝑠 𝑞 sin 𝑞𝜋𝜌 /L . (22)Then, if 𝜂 ℓ is well approximated by cos ℓ𝜋𝜌 /L , we expect 𝑐 ℓ (cid:29) 𝑐 ℓ (cid:48) , 𝑠 𝑘 in which ℓ (cid:48) , 𝑘 ∈ N and ℓ (cid:48) ≠ ℓ . As 𝜂 ℓ becomes more andmore like the local-coupling mode, 𝑐 ℓ will approach unity for aproperly normalized state and 𝑐 ℓ (cid:48) , 𝑠 𝑘 will approach 0. As such, thenumber 𝑓 ℓ = | 𝑐 ℓ | (cid:205) 𝑘 | 𝑐 𝑘 | + | 𝑠 𝑘 | = | 𝑐 ℓ | L ∫ 𝑑𝜌 | 𝜂 ℓ ( 𝜌 )| (23)provides a reasonable measure of the similarity of the ℓ -th true normal mode to the local-coupling normal mode. As ∫ 𝑑𝜌 cos 𝑘𝜋𝜌 /L cos 𝑞𝜋𝜌 /L = L 𝛿 𝑘𝑞 / 𝑐 ℓ = L ∫ L 𝑑𝜌𝜂 ℓ ( 𝜌 ) cos ℓ𝜋𝜌 /L . (24)Hence, 𝑓 ℓ = L (cid:16)∫ 𝑑𝜌𝜂 ℓ ( 𝜌 ) cos ℓ𝜋𝜌 /L (cid:17) ∫ L 𝜂 ℓ ( 𝜌 ) . (25)The values of 𝑓 ℓ for ℓ = , , , 𝑓 𝑘 appears to monotonicallydecrease as 𝑘 increases, which may indicate that the local-couplingnormal mode provides a more accurate model of large-wavelengthexcitations than small-wavelength excitations, although our methodsof extracting the normal modes become infeasible for investigatingthese high-frequency modes. Because the coefficient 𝑓 𝑘 represents only the contribution from co-sine terms to the power spectrum, it is sensitive to assumptions aboutthe boundary conditions on the disc; in the local-coupling model(Batygin 2018), these are taken to be the Neumann boundary con-ditions 𝜕𝜂 / 𝜕𝜌 = 𝜌 = , L . However, the local-coupling modelmay still provide an accurate picture of the normal modes even ifthe boundary conditions differ from these Neumann conditions. Itis therefore important to develop a measure of the accuracy of thelocal-coupling model that is sensitive only to the sinusoidal natureof the normal modes and not to the specific choice of boundaryconditions on those modes.One measure of the accuracy of the fit is the root-mean-squareof the residuals when a sinusoid is fit to a properly normalized state.To carry out this analysis, a standard least-squares algorithm wasused to fit a function of the form 𝐴 sin 𝑘𝑥 + 𝐵 cos 𝑘𝑥 to each of the firstfour numerically extracted normal modes, each normalized so that 𝜂 ( 𝜌 = ) =
1. We should expect that the residuals are proportionalto the magnitude of the state, so this is a reasonable measure ofthe size of deviations from harmonicity that is insensitive to thenormalization of the state. This method was applied to the first four
Table 3.
The Accuracy of the Local-Coupling Normal Modes for a ThinAnnulus ℓ 𝑓 ℓ RMS 𝑅 Table 4.
The variation of the accuracy of local-coupling modes with discsize for 𝑘 = 𝑓 RMS Residual R
50 0.995754 0.00847134 0.99992100 0.994026 0.0170737 0.99967150 0.990185 0.0246727 0.999286300 0.972407 0.0728 0.994461000 0.777252 0.2818 0.993884 non-constant normal modes numerically extracted from the trialdisc. The results are presented in the middle column of Table 3.Another reasonable measure of the accuracy of the fit is the ad-justed 𝑅 value. The results of these measurements are reported inthe final column of Table 3. Just as in the power-spectrum measureof the accuracy of the local-coupling model, these measures indi-cate that the accuracy of the local-coupling model decreases withincreasing wavenumber. More investigation is needed for ℓ ≥ Solutions of the form (cid:205) 𝑘 𝑐 𝑘 𝑒 𝜄𝜔𝑡 cos 𝑘𝜋𝜌 /L do not conserve an-gular momentum if 𝑐 𝑘 ≠ 𝑐𝛿 𝑘,𝑘 for some 𝑐, 𝑘 . For the trial discsimulated, the disc was rather narrow, with 𝑎 out / 𝑎 in = .
1. For sucha narrow disc, the local-coupling modes are close to orthogonal,with only small variations of angular momentum over the courseof a period. As the disc becomes wider, however, this failure ofangular momentum conservation can become more significant. It isthus reasonable to question whether the excellent matches betweenthe local-coupling and general-coupling modes found here are anartifact of the narrow-disc limit.To investigate this possibility, the first normal modes of discswith the same density profile were simulated for different values ofthe width of the disc. More specifically, this was done by keep-ing the separation between adjacent wires fixed, but increasingor decreasing the number of rings of the partition. The resultsfor 𝑁 = , , , , 𝑛 𝑏 ≥ ℓ = ℓ = , ℓ = MNRAS000
1. For sucha narrow disc, the local-coupling modes are close to orthogonal,with only small variations of angular momentum over the courseof a period. As the disc becomes wider, however, this failure ofangular momentum conservation can become more significant. It isthus reasonable to question whether the excellent matches betweenthe local-coupling and general-coupling modes found here are anartifact of the narrow-disc limit.To investigate this possibility, the first normal modes of discswith the same density profile were simulated for different values ofthe width of the disc. More specifically, this was done by keep-ing the separation between adjacent wires fixed, but increasingor decreasing the number of rings of the partition. The resultsfor 𝑁 = , , , , 𝑛 𝑏 ≥ ℓ = ℓ = , ℓ = MNRAS000 , 1–9 (2021)
W. Melton and K. Batygin
Figure 4.
The first normal mode of a disc with 𝑁 = 𝑎 + cos ( 𝑘 𝑥 + 𝜙 ) (blue). decreases as the wavenumber increases. Although local couplingeigenfunctions still offer adequate results for disc widths of thesame order the mean semi-major axis, their degradation with in-creasing L suggests that for Δ 𝑎 / 𝑎 significantly in excess of unity,these expressions become highly approximate. Based upon the re-sults summarized in Table 3, we further expect that the agreementwill gradually worsen with decreasing wavelength. Despite this, thelocal-coupling basis remains useful as the Hamiltonian is almostdiagonal in the local-coupling normal mode basis. Even this levelof agreement may be surprising, however, as the local-couplingmodel is a drastic simplification of the notoriously rich dynamicsof self-gravitating discs. The analogy between the secular small-inclination evolution ofa self-gravitating razor-thin particle disc and the time-dependentSchrödinger equation of quantum mechanics was extended beyondthe local-coupling model developed in Batygin (2018). While thetime-dependent Schrodinger equation 𝜄 (cid:164) 𝜓 = H 𝜓 is most often seenin quantum mechanical regimes, it can appear more generally insituations in which the dynamics are linear and conservative. Insuch a case, the map 𝑡 ↦→ 𝑈 ( 𝑡, 𝑡 ) where 𝑈 ( 𝑡 ) is the generator of aunitary representation of the Lie group R . The Hamiltonian, then,is the generator H = 𝜄𝑑𝑈 ( 𝑡 )/ 𝑑𝑡 (Stillwell 2008).In this work, we have demonstrated that the stable small-inclination dynamics of a self-gravitating disc form a unitary repre-sentation of the time-translation Lie group R . This structure and thedeeper analogy with quantum mechanics provided by the existenceof the conserved inner product (cid:104)·|·(cid:105) give us analytical and numer-ical tools to investigate the stable, secular inclination dynamics ofparticle discs, both when the parameters of the disc are constant andwhen they vary on a secular timescale.By rotating time in the complex plane, we extracted the normalmodes of the disc. This rotation turned the conservative, oscillatorydynamics of the real-time equations of motion into a purely dissi-pative dynamical system in which higher-frequency normal modesdecay away. While explicit diagonalization is also feasible for thediscs discussed in this paper, the imaginary-time formalism pro-vides a numerical method that may prove useful for systems whereexplicit diagonalization is computationally disadvantageous.After using this formalism to find assumption-free approxima- tions for the first four non-constant normal modes, Fourier analysisand least-squares fitting were used to estimate their similarity tothe purely sinusoidal normal modes predicted by the local-couplingmodel. While there were indications that the agreement betweenthe two sets of normal modes decreases with increasing wavenum-ber and with the width of the disc, we found that the modes areindeed very close to the sinusoidal modes predicted by the localcoupling model for narrow annuli, indicating that there is an inter-esting region of parameter space where the local-coupling normalmodes, which have an explicit analytical form in the continuumlimit, well-approximate the true normal modes of the disc. Evenin wider discs, the equations of motion for the disc are almost di-agonal when written in the local-coupling mode basis, which candramatically simplify numerical calculations.This quantum-mechanics-like formalism also provides a con-venient tool for estimating the response of discs to external forces.If the perturbing system preserves the rotational symmetry of thedisc, the perturbing Hamiltonian generated by the external systemwill be Hermitian, and quantum perturbation theory can be used toestimate how the frequencies and shapes of the normal modes willchange and provide a measure of the rigidity of the disc to externalperturbations. This perturbative technique is justified by the exis-tence of a ‘quantum’ Hamiltonian that preserves an inner product,as has been shown in this paper.This work clarifies the applicability of the local couplingmodel. In thin, narrow discs whose dynamics are dominated byself-gravitational effects, , the local-coupling normal modes pro-vide an excellent approximation to the true normal modes of thedisc. These systems include, for example, narrow annuli of dustformed during planet formation in protoplanetary discs. Moreover,the Terrestrial planets of the solar system themselves are routinelyenvisioned to have formed from a narrow annulus of solid debris(Hansen 2009). By carrying out the analysis of this work, we havealso lifted one of the chief limitations of the approach adopted by(Batygin 2018): the assumption that the dynamics are effectivelylocal. While this breaks the theoretical tractability of the normalmodes, it provides a self-consistent method of finding the normalmodes and their frequency of evolution, as well as ways to computeangular momentum and the rate of change of angular momentum.Nevertheless, there are several important limitations to our updatedmodel. First, the dynamics considered in this paper occur only on thesecular timescale; this model fails to capture dynamics that occuron timescales comparable to the orbital period or comparable to thelifespan of the system. Correspondingly, the orbit-averaged natureof the dynamics discussed in this paper preclude it from describingsome interesting behavior.Second, while this paper usefully extends the local-couplingSchrödinger model of the secular small-inclination dynamics of arazor-thin self-gravitating particle disc, significant work remains.Most prominently, we assumed that the dynamics were generatedonly through purely gravitational forces, and that the disc is a parti-cle disc. While neglecting pressure forces is a good approximationif √G Σ 𝑎 is much greater than the speed of sound (Tremaine 2001),any comprehensive theory needs to account for these internal ef-fects. Other authors have argued that discs evolving purely throughpressure forces, neglecting viscosity and self-gravity, can be de-scribed by a non-linear Schrödinger equation. Incorporating thisformalism with the theory contained in this paper may allow us toextend this work to fluid pressure-supported discs. Indeed, such adevelopment constitutes an interesting avenue for future work. MNRAS , 1–9 (2021) igenstates of Quasi-Keplerian Discs ACKNOWLEDGEMENTS
W. M. is grateful to the Caltech Summer Undergraduate ResearchFellowship program, during which this work was initiated. K.B. isgrateful to the David and Lucile Packard Foundation and the AlfredP. Sloan Foundation for their generous support. Additionally, wewould like to thank the anonymous referees for providing usefulinsights that have led to an improvement of the manuscript.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Armitage P. J., 2009, Astrophysics of Planet Formation. Cambridge Univer-sity Press, doi:10.1017/CBO9780511802225Bate M., Lodato G., Pringle J., 2010, MNRAS, 401, 1505Batygin K., 2018, MNRAS, 475, 5070Binney J., Tremaine S., 2008, Galactic Dynamics. Princeton University PressGoldreich P., Tremaine S., 1982, ARA&A, 20, 249Hahn J. M., 2003, ApJ, 595, 531Hansen B. M. S., 2009, The Astrophysical Journal, 703, 1131Kocsis B., Tremaine S., 2011, MNRAS, 412, 187Laskar J., 1997, A&A, 317, L75Morbidelli, A. 2020, A&A, 638, A1Murray C. D., Dermott S. F., 1999, Solar System Dynamics, 1 edn. Cam-bridge University PressNesvold E., Naoz S., Vican L., Farr W. M., 2016, ApJ, 826Spalding C., Batygin K., 2014, ApJ, 790Stillwell J., 2008, Naive Lie Theory. SpringerStrocchi F., 1966, Rev. Mod. Phys., 38Tremaine S., 2001, AJ, 121, 1776
APPENDIX A: A DIRECT PROOF OF THECONSERVATION OF THE INNER PRODUCT
In this appendix, we present a proof that the inner product (cid:104) 𝜂 | 𝜈 (cid:105) defined in Equation 14 is conserved without assuming that angularmomentum ( (cid:104) 𝜂 | 𝜂 (cid:105) ) is conserved. Earlier, we took the continuumlimit of the disturbing function using a logarithmic parametriza-tion of radial distance 𝜌 = log 𝑎 / 𝑎 in . Once we take the continuumlimit, we can choose to reparametrize the problem in terms of an-other coordinate 𝜒 , such that the semimajor axis of the infinitesimaldisc parametrized by 𝜒 is 𝑎 ( 𝜒 ) . In terms of 𝜒 , the Hamiltoniandescribing the coupling becomes H 𝜂 = ∫ 𝐶 ( 𝑎 ( 𝜒 ) , 𝑎 ( 𝜒 (cid:48) )( 𝜂 ( 𝜒 ) − 𝜂 ( 𝜒 (cid:48) )) 𝜋𝑎 ( 𝜒 (cid:48) ) Σ ( 𝜒 (cid:48) ) 𝑑𝑎𝑑𝜒 (cid:48) 𝑑𝜒 (cid:48) , (A1)with coupling between the points 𝜒 and 𝜒 (cid:48) being described by thefunction 𝐵 ( 𝜒, 𝜒 (cid:48) ) = 𝜋𝑎 ( 𝜒 (cid:48) ) Σ ( 𝑎 ( 𝜒 (cid:48) )) 𝑑𝑎𝑑𝜒 ( 𝜒 (cid:48) ) 𝐶 ( 𝑎 ( 𝜒 ) , 𝑎 ( 𝜒 (cid:48) )) . (A2)Suppose that we choose 𝑎 ( 𝜒 ) such that 𝐵 ( 𝜒, 𝜒 (cid:48) ) = 𝐵 ( 𝜒 (cid:48) , 𝜒 ) . Rear-ranging this equation, we find that this implies that 𝑑𝑎𝑑𝜒 ( 𝜒 ) = 𝑎 ( 𝜒 (cid:48) ) Σ ( 𝑎 ( 𝜒 (cid:48) )) 𝐶 ( 𝜒, 𝜒 (cid:48) ) 𝑎 ( 𝜒 ) Σ ( 𝑎 ( 𝜒 )) 𝐶 ( 𝜒 (cid:48) , 𝜒 ) 𝑑𝑎𝑑𝜒 ( 𝜒 (cid:48) ) . (A3) Because 𝐶 ( 𝑎, 𝑎 (cid:48) ) 𝐶 ( 𝑎 (cid:48) , 𝑎 (cid:48)(cid:48) ) 𝐶 ( 𝑎 (cid:48)(cid:48) , 𝑎 ) = 𝐶 ( 𝑎 (cid:48) , 𝑎 ) 𝐶 ( 𝑎 (cid:48)(cid:48) , 𝑎 (cid:48) ) 𝐶 ( 𝑎, 𝑎 (cid:48)(cid:48) ) , this equation provides a consistent definition of 𝑑𝑎 / 𝑑𝜒 . With thisparametrization, we can find a classical Hamiltonian that describesthe full system: H 𝑐 = ∬ 𝐵 ( 𝜒, 𝜒 (cid:48) )( 𝜂 ( 𝜒 ) 𝜂 ( 𝜒 (cid:48) ) ∗ ) 𝑑𝜒𝑑𝜒 (cid:48) − ∫ (cid:18)∫ 𝐵 ( 𝜒, 𝜒 (cid:48) ) 𝑑𝜒 (cid:48) (cid:19) 𝜂 ( 𝜒 ) 𝜂 ( 𝜒 ) ∗ 𝑑𝜒. (A4)Following the derivation above with this parametrization leads to aHermitian Hamiltonian with respect to the inner product (cid:104) 𝜂 | 𝜈 (cid:105) = ∫ 𝜂 ∗ 𝜈𝑑𝜒 (Strocchi 1966). Incorporating formula (14)above, we have that the conserved inner product is (cid:104) 𝜂 | 𝜈 (cid:105) ∝ ∫ 𝜂 ∗ ( 𝑎 ) 𝜈 ( 𝑎 ) 𝑎 Σ ( 𝑎 ) 𝐶 ( 𝑎 , 𝑎 ) 𝑎 Σ ( 𝑎 ) 𝐶 ( 𝑎, 𝑎 ) 𝑑𝑎 / 𝑑𝜒 ( 𝑎 ) 𝑑𝑎. (A5)For our choice of disc parameters, we have that (cid:104) 𝜂 | 𝜈 (cid:105) = ∫ 𝜂 ∗ 𝜈𝑎𝑑𝑎 ascan be directly computed from the above expression. For arbitrarysurface density profile Σ ( 𝑎 ) , the inner product derived here takesthe form given by Equation 14. We have not assumed that (cid:104) 𝜂 | 𝜂 (cid:105) isconserved; hence, we have proved that this model conserves angularmomentum. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000