Nonlinear dynamics of hydrodynamic tori as a model of oscillations and bending waves in astrophysical discs

Understanding oscillations and waves in astrophysical fluid bodies helps to elucidate their observed variability and the underlying physical mechanisms. Indeed, global oscillations and bending modes of accretion discs or tori may be relevant to quasi-periodicity and warped structures around compact objects. While most studies rely on linear theory, observationally significant, nonlinear dynamics is still poorly understood, especially in Keplerian discs for which resonances typically demand a separate treatment. In this work we introduce a novel analytical model which exactly solves the ideal, compressible fluid equations for a non-self-gravitating elliptical cylinder within a local shearing sheet. The aspect ratio of the ring is an adjustable parameter, allowing a continuum of models ranging from a torus of circular cross-section to a thin ring. We restrict attention to flow fields which are a linear function of the coordinates, capturing the lowest order global motions and reducing the dynamics to a set of coupled ordinary differential equations (ODEs). This system acts as a framework for exploring a rich range of hydrodynamic phenomena in both the large amplitude and Keplerian regimes. We demonstrate the connection between tilting tori and warped discs within this model, showing that the linear modes of the ring correspond to oppositely precessing global bending modes. These are further confirmed within a numerical grid based simulation. Crucially, the ODE system developed here allows for a more tractable investigation of nonlinear dynamics. This will be demonstrated in a subsequent paper which evidences mode coupling between warping and vertical motions in thin tilted rings.


Astrophysical motivation
Discs appear ubiquitously in the zoo of astrophysical phenomena and occur in a range of systems. From protoplanetary discs to thick accretion flows around black holes, discs are now regarded to be as dynamically important as the central host objects themselves. These structures are inherently dynamic and support a range of oscillations and waves which may be responsible for observational signatures. The key restorative effect which modifies these modes in a non-inertial frame is the Coriolis force, inherited from the rotational motion of the fluid. Whilst a local treatment can categorise a host of wave regimes (e.g. Carroll et al. 1985;Kato 2001), observationally significant signatures require global coherent modes of large amplitude (Okazaki et al. 1987).
Several authors have captured such global modes by confining oscillations within tori of limited radial extent. These structures ★ E-mail: cwf29@cam.ac.uk † E-mail: gio10@cam.ac.uk are thought to form around black holes and neutron stars as puffed up hot accretion flows (Frank et al. 2002). The initial theoretical expositions (e.g. Abramowicz et al. 1978;Blaes 1985) garnered further astrophysical interest in light of the Rossi X-Ray Timing Explorer (RXTE) observations. Quasi-periodic oscillations (QPOs) in accreting black holes and neutron stars have often been given interpretations that involve oscillating or precessing rings (e.g. Nowak et al. 1997;Stella & Vietri 1999). High-frequency QPOs sometimes exhibit frequency ratios that have been modelled as the resonances occurring between orbital, vertical and radial oscillations (e.g. Abramowicz & Kluźniak 2001;Rezzolla et al. 2003;Fragile et al. 2016;de Avellar et al. 2018).
Other oscillating structures have also been identified. In particular warped discs, wherein the orbital plane varies with radius, may explain a range of puzzling observations. Warps can be interpreted as global disc modes which break the mid-plane symmetry of the system and occur whenever a misalignment is present in a system. For example, the seminal work of Bardeen & Petterson (1975) studied the warped structure due to the misaligned black hole spin and angular momentum axis of the disc. Since then, a range of other scenarios have been investigated. When the magnetic dipole of a central source is tilted with respect to the disc plane it can induce a warp (Lai 1999). Alternatively, companion stars or planets introduce gravitational torques. These may excite tilt in discs via an inclination instability as investigated by Lubow (1992). This misalignment then allows for the torquing of fluid rings, driving differential precession which manifests as a warp (Papaloizou & Terquem 1995;Lubow & Ogilvie 2000).
These theoretical efforts have been supported by several key observations. The first indirect observation of a warped disc came from the long period luminosity variations in the Her X-1 X-ray binary source (Katz 1973). The line of sight flux is modulated periodically as light from the source is attenuated by the precessing tilted outer edges. More examples of these so called 'superorbital' X-ray binary systems have since been identified by Kotze & Charles (2012). Direct observations were then provided by maser emission which traces the warped structure of discs. The active galaxy NGC4258 (M106) is one such prominent example. This intermediate spiral galaxy exhibits a strong emission feature tracking a warped accretion disc around a central black hole (Miyoshi et al. 1995).
More recently, interferometric techniques employed by the Atacama Large Millimeter/submillimeter Array (ALMA) have heralded a new era of disc observations. Sakai et al. (2019) reported dust continuum observations of a young protostellar disc with misaligned inner and outer discs, possibly due to anisotropic gas accretion or a misaligned magnetic field. Other protoplanetary discs have proved fertile ground for warp hunting. The Hubble Space Telescope images of TW Hya have captured a shadow moving around the outer disc regions at a rate of 22.7 • yr −1 -faster than any feature possibly advected with the flow. This suggests an inner tilted precessing disc is blocking light from the central T Tauri star (Debes et al. 2017). ALMA detection of CO molecular line emission also allows for a kinematic study of the gas flow which is effectively modelled with a warped inner structure (Rosenfeld et al. 2012). These warps inferred from shadows are even more exaggerated in transition discs where large radial gaps divide the inner and outer regions and they present extreme inclination differences. The comparison of the HD 142527 gapped protoplanetary disc with parametric radiative transfer models allowed Marino et al. (2015) to deduce a large relative inclination of 70 • between the inner and outer discs, which may be caused by a companion planet. Similar results have been found for other systems including the DoAr 44 T Tauri transition disc (Casassus et al. 2018). SPHERE+IRDIS observations reveal clear azimuthal dips in reflected infrared polarised intensity which are effectively modelled by an inner disc tilted at ∼ 30 • . In some systems these distinct rings are thought to form by disc tearing and breaking (Nixon & King 2012). Indeed, there has been a recent observation of the spectacular triple star system GW Orionis wherein gravitational effects have torn the disc into independently precessing rings (Kraus et al. 2020).
The growing host of observational tilted rings and warps emphasises the need for a comprehensive theoretical understanding of these systems. Whilst there is a body of theory investigating linear warps in discs, the nonlinear dynamics is still poorly understood, particularly in the most important Keplerian regime. The resonances between the vertical and radial motions make this regime trickier to understand but all the more important for the dynamical evolution.

Plan of this paper
The main purpose of this paper is to introduce a novel ring model which offers a theoretical framework for advancing our understanding of oscillations in tori and warped discs. In particular, this model will allow for nonlinear and resonant phenomena to be investigated, extending the previous body of theory. Whilst thick tori and thin warped disc regimes seem to probe very different physical settings, by varying the aspect ratio of the torus, we expect many of the dynamical results to be smoothly related.
Here we will focus on the development of the ring model and establish its validity by connecting it with previous linear theory. We begin by localising our fluid equations about a reference circular orbit and then seek axisymmetric, linear flow field solutions which represent the lowest order global modes for an oscillating torus. This reduces to solving a set of coupled ordinary differential equations (ODEs) which exactly describe the flow and shape of the ring. We will also present an alternative formulation of the model using a Lagrangian construction which may prove more amenable to future nonlinear analysis. In order to confirm the correspondence with warped disc theory we will examine the small amplitude modes of the torus and compare these with linear bending waves. These may be understood as global rotating structures exhibiting a hierarchy of precessional timescales. To support these findings, we will briefly present a numerical setup using a grid based solver which is capable of capturing these linear modes and will highlight some of the difficulties in modelling torus oscillations. Finally, we will discuss the connection of our work to previous theoretical efforts. In a subsequent paper, we will apply this model more rigorously to examine the extreme nonlinear dynamics supported by the ring. We find that highly nonlinear mode coupling can occur as energy is interchanged between warping and vertical bouncing motions. These mixed modes merit further attention, in particular with relation to highly warped disc and precessing tilted rings in broken disc systems.

Localised fluid equations
We will begin by introducing a fiducial frame which allows for a localisation of the fluid equations. Following the standard shearing box construction (e.g. Hill 1878; Hawley et al. 1995), we expand our equations about a local circular reference orbit at 0 with angular velocity Ω 0 = Ω( 0 ), assuming an axisymmetric potential Φ( , ). This orbit has an attached, co-rotating coordinate system ( , , ) which is defined by = ( − 0 ), = 0 ( − Ω 0 ) and = , such that , and are the radial, azimuthal and vertical directions respectively. Within the local expansion, for which the size of the domain is much less than 0 , the coordinate system is effectively Cartesian and curvature effects are neglected. The equations of motion are then found to be is the Lagrangian derivative, is the velocity, is the pressure and is the density. The local expansion of the tidal potential, Φ t , is given by where 0 = −( Ω/ ) 0 is the orbital shear rate and 2 0 = ( Φ) 0 is the vertical oscillation frequency squared of a test particle perturbed from its circular orbit. The first term in Φ t is a reservoir of shear flow energy and describes the destabilising centrifugal contribution arising from the non-inertial acceleration of the orbiting frame. Meanwhile, the second term provides the restorative vertical gravitational term as test particles oscillate harmonically about the mid-plane. A test particle will similarly oscillate radially about the reference orbit at the natural epicyclic frequency 0 , thanks to the restorative Coriolis forces appearing in the equation of motion. This is defined by For any spherically symmetric potential, 0 = Ω 0 . Furthermore, for a Keplerian disc in a Newtonian point-mass potential, 0 = (3/2)Ω 0 and so 0 = 0 = Ω 0 . This resonance can lead to distinct dynamical behaviour of interest and should be treated carefully. Our model will initially explore the simplest case of an ideal, adiabatic fluid so the dynamical equations are complemented by the mass continuity and energy equations where is the velocity divergence and is the adiabatic index (assumed constant). Together these govern the evolution of density and pressure.

A simple equilibrium solution
We can motivate our general ring solutions by first examining a simple family of equilibrium configurations for the flow. Indeed, the simplest steady state solution of these equations consists of the shear flow = − 0 , which describes the circular orbit of fluid particles without any radial pressure gradient support, whilst hydrostatic balance holds in the vertical direction. We can now look for solutions where the shear deviates from this trivial orbital rate and instead consider zonal flows where a geostrophic balance holds between the Coriolis force and radial pressure gradients. Let us introduce this zonal flow to be = − where is some constant defining the linear shear flow. We look for equilibrium solutions in which and depend only on ( , ). Equations (1)-(3) are then balanced when which is equivalent to It then follows that both and are functions of only, with the differential relation = − . If > 0 , the contours of are similar concentric ellipses and an equilibrium exists in the form of a ring with an elliptical cross-section. This occurs when the shear flow exceeds the orbital shear rate. The aspect ratio of this ring, , is then defined by the ratio of the radial width to the vertical thickness and is calculated to be Tuning of then allows for a ring of arbitrary aspect ratio. As the the zonal flow increases, the counteracting pressure gradients must become steeper and the ring becomes more radially confined. Of course, for the choice = 0 , then = 0 and the ring returns to the limiting radially extended case.

Dynamical ideal ring solutions
These elliptical ring equilibria can be generalised to dynamically evolving tori. We will again look for solutions which are locally axisymmetric and independent of . We assume that the spatial density and pressure can be described by a common function ( , , ) which is a material invariant advected by the fluid motion, i.e.
Furthermore, to accommodate the changes in density and pressure as the fluid undergoes compressive motions, we write the density and pressure using the separation We now enforce the key assumption of this framework, positing that the fluid motion is described by a linear flow field in the Cartesian coordinates such that where is a time-dependent, square flow matrix. Since the flow is independent of we also recognise that 2 = 0. This allows us to write the divergence as which depends only on time and corresponds to a homogeneous expansion or contraction of the fluid. With these assumptions the mass continuity and energy equations are satisfied provided where = / . Such a linear velocity field will map ellipses into ellipses. Hence we adopt our materially conserved function to be a time-dependent quadratic function of the Cartesian coordinates, where is some constant and ( ) is a time dependent, positivedefinite shape matrix with 2 = 2 = 0 in the -independent case. Thus contours of equal density and pressure lie on ellipses akin to the equilibrium discussed in Section 2.2. The shape of these elliptical contours is described by the evolution of which is found by imposing the material conservation condition (15). Applying the material derivative to and using the linear flow field assumption (18) yields This can only be satisfied for all points in space provided which reduces to 11 + 2( 11 11 + 13 31 ) = 0, 13 + 11 13 + 33 31 + 13 Δ = 0, 33 + 2( 13 13 + 33 33 ) = 0.
In order to determine the evolution of the flow matrix, we must return to the equations of motion (1)-(3). Our assumptions ensure that all the terms are manifestly linear in the Cartesian coordinates, except possibly those involving the pressure gradient which will have the form Since ∇ is linear in the coordinates, we can obtain a consistent solution provided˜/ ∝˜. Without loss of generality, we can define the functions such that whereby the proportionality has been absorbed into the time dependent termsˆandˆ. There are many choices available which satisfy this condition and link the density and pressure. An oft used family of solutions is found by assuming a polytropic relationship with index such that In the case that = ∞ we recover a spatially isothermal ring with a Gaussian distribution of pressure and density. Otherwise, a finite value of gives a definite ellipse boundary where density and pressure drop to zero and match onto the surrounding vacuum. Indeed, the limiting case = 0 describes a homogeneous ring. However, it should be noted that the dynamical evolution of the system doesn't depend on the specific choice of relationship and only requires that equation (29) be satisfied. Using this condition in combination with the previous assumptions, the equation of motion is now satisfied provided 11 + 2 11 + 13 31 − 2Ω 21 = 2Ω +ˆ1 1 , 13 + 11 13 + 13 33 − 2Ω 23 =ˆ1 3 , 21 + 21 11 + 23 31 + 2Ω 11 = 0, 23 + 21 13 + 23 33 + 2Ω 13 = 0, where we have dropped the subscript on the vertical, epicyclic and orbital frequencies and also the orbital shear rate for ease of notation. Hereˆ( ) =ˆ/ˆis a characteristic temperature, which evolves according tô which can be shown using equations (20) and (21). Together, the evolutionary equations (25)-(27) for the shape matrix, (31)-(36) for the flow matrix and (37) for the characteristic temperature, constitute a closed system of 10 first-order, coupled, non-linear ODEs. These govern the dynamics of the ideal ring and despite the various assumptions imposed, we have in fact introduced no approximation. Indeed, these equations present an exact nonlinear solution framework to investigate the largest scale global modes of the tori.

Defining the reference state
We can alternatively construct our model from a Lagrangian perspective whereby we track the evolution of individual particles. This will naturally elucidate the conservation properties of the flow and emphasise the underlying oscillatory behaviour of individual fluid particles. Consider a reference state denoted by the Lagrangian coordinates ( 0 , 0 , 0 ). In general this does not have to be an equilibrium state but is instead an arbitrary, fixed configuration which simply acts to label fluid particles. In order to describe the linear flow fields, we map these fluid particles to their dynamical state by means of a linear transformation denoted where is the time dependent Jacobian matrix. For the assumed axisymmetric setup this will consist of 6 components describing the flow, (39) Using this formalism, the fluid velocity is simply = J 0 so comparison with equation (18) yields the relation between the flow matrix and the Jacobian elements, Recall that we require elliptical contours of density, as assumed in the Eulerian model. Thus, without loss of generality, we load mass in the reference state such that particles with constant lie on circular contours. When these circular fluid rings are acted upon by the linear transformation J, they will form ellipses with constant . This may be formulated mathematically as with Here I is the identity matrix and we have introduced as the characteristic radius of the reference state which is defined by the second mass weighted moment, Here, is the total mass per unit length and denotes integration over the 2D mass distribution. Therefore ( ) defines a dimensionless radius of the circular contours in the reference state. Such a mapping is visualised in Fig. 1. This parameterises the spatial part of the density and pressure separation such that where 0 and 0 denote the density and pressure in the reference state whilstˆ0 andˆ0 are characteristic density and pressure factors. Note that we have recast the functional dependence of˜and˜from to via = − (1/2) 2 . The action of the Lagrangian mapping then scales area elements by the determinant of the Jacobian matrix, such that the density and pressure transform as Furthermore, the homogeneous compression/expansion of the ring is related to this Jacobian area scaling by This result can be inserted into equation (37) and integrated to yield the evolution forˆin terms of the Jacobian elements, where comparison with equations (44) and (46) showsˆ0 =ˆ0/ˆ0. Evidently, this integration constant arises as an expression of the conservation of entropy on each fluid particle in the ring.

Identifying the Lagrangian
Understanding the evolution of the ring now reduces to finding the time evolution of the Jacobian transformation elements, . We proceed by constructing a Lagrangian which encapsulates the relevant physics. Within the local model this has the form where the first term represents the kinetic energy, the second term arises from the non-inertial effects of the rotating frame, is the internal energy and Φ t is the tidal potential. With invariance along , we take to represent the mass per unit length and integrate over the total cross-sectional area. Each term must now be reformulated in terms of the Jacobian elements which act as the generalised coordinates of the model. The kinetic energy term is expanded to We can express the velocity in terms of the Jacobian elements by switching to the reference state, so for example, is the mass weighted covariance measure of the reference coordinate moments. The circular symmetry of the reference state allows us to identify the second density moments as Inserting this into (51) and repeating for the remaining terms yields L kinetic = 2 2 2 11 + 2 13 + 2 21 + 2 23 + 2 31 + 2 33 .
One can proceed similarly for the rotational and tidal contributions which respectively give Finally we examine the internal energy contribution. Using the ideal equation of state, the energy density integral becomes Inserting the separation ansatz (17) and transforming to Lagrangian coordinates gives where the area scaling incurs a factor of . In order to progress we recall that the spatial pressure term is constant on circular Lagrangian contours as per equation (44), so˜=˜( ). Thus it is convenient to switch to polar Lagrangian coordinates such that where is the usual polar angle in the Lagrangian reference state and so Using condition (29) allows us to relate˜and˜viã Integrating equation (61) by parts, inserting the expression above and converting back to an integral over mass gives The integral can now be computed from the standard covariance results (53) and the ratioˆ/ˆis given by equation (48) such that Finally, summing all contributions (54), (55), (56) and (63) and cancelling the 2 factor gives the total Lagrangian, L = 1 2 ( 2 11 + 2 13 + 2 21 + 2 23 + 2 31 + 2 33 ) −ˆ0 ( − 1) −1 2 + Ω ( 2 11 + 2 13 ) − 1 2 2 ( 2 31 + 2 33 ) + 2Ω( 11 21 + 13 23 ).
We can now derive the dynamical evolution by means of the usual Euler-Lagrange equations This gives Formally this is a coupled system of 6 non-linear, second order differential equations which govern the evolution of the flow. These can in fact also be derived as a specific sub-case of the general Lagrangian framework for the affine motion of discs as developed by Ogilvie (2018) and discussed in appendix A.

Conserved circulation integrability
Equations (68) and (69) are clearly integrable, allowing for the reduction in order of the system of equations. This naturally arises as a consequence of conservation laws which restrict the degrees of freedom of this system. The conservation of entropy has already introduced the constantˆ0 into equation (48). Now, we will make use of the conservation of circulation to interpret the integrability of the 21 and 23 terms. The Poincaré-Bjerknes circulation theorem extends Kelvin's result to the rotating shearing box reference frame such that where with the line integral taken around the material loop C( ). Since material contours are fixed in the Lagrangian reference frame it is convenient to switch to this perspective. Using the transformation = 0, we have that For any material contour C 0 , the integral can be evaluated to give a constant matrix which encapsulates the relation between the Jacobian coordinates through (74). This matrix is evaluated using Stokes' Theorem which gives where, is the anti-symmetric Levi-Civita symbol and A 0 is the vector area of any open surface bounded by C 0 . We are free to choose the material contours however we like, so in particular we may take 0,1 ≡ and 0,2 = 0,3 = 0 so that the contour lies entirely in the − plane and has a conserved circulation denoted Γ . Inserting this into (74) gives where −Γ / ≡ is obviously the constant identified by integrating (69). Equally we can choose 0,3 ≡ and 0,1 = 0,2 = 0 so the contour lies entirely in the − plane and the conserved circulation is denoted Γ . In this case and Γ / ≡ is the constant arising from the integration of (68). Using this integrability, we may reduce equations (66) -(71) to where the shear rate has been eliminated in favour of the epicylic frequency. These equations clearly highlight the harmonic oscillatory structure inherited from the simple test particle motion. These are then coupled together by means of the pressure term on the right hand side.

LINEAR MODES
Having developed the ring model, we may now examine the simplest linear dynamics. To this end we will explore the linear modes supported by the system. Note that the linear transformations permitted by equation (18) leave the centre of mass of the ring fixed at the origin of the shearing box. Thus we are neglecting transformations which translate the whole ring vertically or radially. Such vertical oscillations at the vertical frequency, , would simply correspond to a global inclination of the ring with respect to the chosen shearing box orbital plane. Meanwhile, radial oscillations at the epicylic frequency would describe an eccentric offset of the torus. Here we focus instead on the physically interesting breathing modes which engage the thermodynamics of the ring and tilting modes which rock back and forth about the centre of the shearing box.

Eulerian perspective
In accordance with Section 2.2 we can identify the equilibrium solutions about which to perform the perturbation analysis. We take to be constant with 13 = 0, such that the ring has an aspect ratio = √︁ 11 / 33 . This shape is supported by a geostrophic balance where the flow matrix vanishes apart from the 21 component. In this case equilibrium is achieved provided −2Ω( + 21 ) =ˆ1 1 and 2 =ˆ3 3 .
We can now introduce small perturbations about this equilibrium ring in the form = , + , = , + andˆ=ˆ+ˆ , where the subscript denotes the equilibrium background quantity.
Inserting into (25) -(27) and (31)-(37) and dropping non-linear terms yields the matrix equation where X = [ 11 , 13 , 33 , 11 , 13 , 21 , 23 , 31 , 33 ,ˆ ], and is a sparse 10 × 10 matrix which depends on the background equilibrium. We look for oscillatory solutions with time dependence exp[ ], which results in the eigenvalue problem X = X. Solving this gives a 10th order dispersion relation for the 10 associated eigenmodes. Alternatively a more enlightening approach is to identify the separate mode families by carefully motivating the form of the initial perturbations.
These modes generally couple the vertical and radial directions as indicated by the presence of both natural frequencies and . Furthermore, the dependence on the adiabatic index signifies that the modes are compressive in nature. Indeed, as we let → ∞ in the incompressible limit, the frequencies become infinite as the sound speed diverges. In the opposite isothermal regime, = 1 and the radial and vertical modes decouple. The radial modes take on the form of an inertial-acoustic oscillation whilst the vertical mode reduces to the usual = √ 2 as expected for perturbed isothermal discs (Ogilvie & Latter 2013). We also note the dependence on the aspect ratio of the ring. As → 0 the ring becomes extended and the radial pressure gradients are negligible. Thus the radial oscillation breathing mode is no longer driven by the vertical motion and the frequency reduces to .

Zero frequency modes
The above analysis identifies 4 eigenmodes out of the 6 available for the 6th order system. The remaining 2 modes correspond to the zero frequency perturbations which essentially shift the initial equilibrium to a neighbouring equilibrium ring. The fact these zero frequency modes appear within the breathing mode analysis is to be expected since the equilibrium configurations all possess even − parity so are reached by symmetric perturbations. This shift to a neighbouring equilibrium is constrained by the conservation laws which are associated with each zero frequency mode. The conservation of entropy ensures that 1 2 11 ,11 whilst the conservation of angular momentum enforces 21 2Ω + ,21 − 11 2 11 = 0.

Tilting modes
The remaining 4 modes associated with the full 10th order system are now identified by introducing perturbations which break the reflectional symmetry. Retaining the off-diagonal terms and setting 11 = 33 = 21 = 11 = 33 =ˆ = 0 yields a 4th order eigenvalue problem for the remaining variables. Eliminating 23 and 13 in favour of 13 and 31 gives the eigenvalue problem as The associated frequencies are then given by These tilting modes break the symmetry of the ring as the major and minor axes are now inclined with respect to the mid-plane. The appearance of and again signifies the strong coupling of the vertical and radial oscillations. However, now the absence of implies that they are incompressible. This is indeed the case as 11 = 33 = 0 ensures that the divergence Δ = 0, in accordance with equation (19). For small aspect ratios the approximate solutions are which shows a distinction between the dominant radial epicyclic shearing mode of the ring as opposed to the dominant vertical tilting mode of the ring. However, in the Keplerian case = and there is resonant degeneracy between the epicyclic frequencies so instead When viewed from the local model, these modes appear to rock back and forth periodically about the centre of the shearing box. However an alternative perspective is to Doppler shift these solutions back into the global reference frame. Now the dynamical tilting can be interpreted as the passage of the global ring structure seen by an observer moving around on the fiducial orbit. For example, a stationary global torus, with some azimuthally dependent tilt about this orbit, will appear to oscillate within the local model. This is visualised in figure 2. It is worth re-emphasising the key but subtle distinction between a local tilt of the elliptical cross-section about the centre of the shearing box, versus a global inclination of the whole ring structure which is simply a trivial rotation within a spherically symmetric potential. In this paper a tilting motion refers to the former, non-trivial mode and is viewed as a warping of the disc about the shearing box orbital plane, as will be explored further in section 5.

Bending wave theory
The linear tilting modes identified in Section 4.1.3 are analogous to warped structures in discs. Indeed, warps have previously been examined within the context of a local model by Ogilvie & Latter (2013). They note that a fixed, global warp appears as a uniform rocking of the disc within the local model. Similarly in this ring model, the tilting torus will further illuminate the behaviour of warps. To this end we will formally illustrate the correspondence between our tilting modes and linear warps, as described by the bending wave theory of Lubow & Ogilvie (2000). Let us restrict our attention to the thin disc, inviscid regime without nodal precession such that = Ω. Within this framework, the warp structure in the global inertial frame is governed by a pair of partial differential equations, where Σ and are the vertically integrated density and pressure respectively. l and G are then the horizontal components of the unit tilt vector and the internal torque, which describe the evolution of the warp. The term ≡ (Ω 2 − 2 )/2Ω encapsulates the detuning from Keplerian resonance which results in an apsidal precession. Indeed, when the epicyclic frequency only differs slightly from the orbital frequency, ≈ Ω − to leading order. Now, in order to compare with our local model we must expand these equations about the reference radius 0 : These equations can be combined into a single, second order PDE for the tilt vector. This has two independent components so is amenable to a complex description, where = + is the complex tilt variable. The vertically integrated density and pressure structure is now set in correspondence with the ring model. We perform a separation of variables as before such that =ˆ˜where˜= ∫˜ (  99) and Note once again the ratioˆ/ˆ=ˆ. Differentiating (99), inserting (29) and using the definition of given by equation (22) yields Now, recall equation (83) for the hydrostatic equilibrium 33ˆ= 2 = Ω 2 . Thus, 11ˆ= 2 Ω 2 . Using this result and inserting into equation (98) We now seek oscillatory tilting solutions which are linear in the coordinate . These will take the form = 0 exp[ ] . Inserting this ansatz yields the frequencies, Physically these oscillations correspond to a rotation of the tilt vector according to where 0 is some arbitrary phase offset. Thus we interpret as a precessional frequency -as the warped global structure slowly rotates, an observer fixed at a point in the ring will see it rock back and forth. Clearly > 0 and < 0 result in prograde and retrograde precession respectively. In the limit → 0 the ring reduces to a thin disc with = 0 or = . Thus it permits a stationary warped structure or a shearing mode which is driven by the radial apsidal precession. In the Keplerian resonant regime with = 0 we see that = ± Ω/2 so the precession is driven entirely by the radial pressure gradients set up due to the finite extent of the ring.

Identification with tilting ring modes
Indeed we can formally identify these bending modes with the linear tilt solutions identified in Section 4.1.3. The bending wave analysis is performed about some local radius in a non-rotating, inertial frame where the warped structure corresponds to a global pattern with azimuthal wavenumber = 1. In order to compare with the frequency observed in the orbiting frame, we must perform a Doppler shift which changes the frequency by 1 times Ω such that The linear bending wave theory is formally valid in the case of thin discs with small deviations from Keplerian resonance, where 1 parameterises the deviation from resonance. Thus and are both small. We now focus on two regimes in turn. First consider the non-resonant regime for which / 1. In this case Taking the (+) prograde solution and Doppler shifting yields This will match onto the result found for the non-resonant tilting mode given by equation (92) provided the ( 2 ) terms are negligible. Taking ∼ and balancing the order of the last two terms yields the range of agreement 2/3 < < 1. Meanwhile the (-) retrograde mode gives These match onto the results found for the non-resonant tilting modes in (92). Similarly, we can examine the resonant case where = 0 and = ± Ω/2. The Doppler shifted frequency is then which agrees with the resonant ring tilting modes found in (93). This agreement is visualised in the upper panel of figure 3 which compares 2 for the Doppler shifted bending wave theory (dashed lines) and the ring model (solid lines) as is varied whilst = 0.01 is held fixed. We see that the upper frequency branch corresponding to the retrograde precession matches well for all epicyclic frequencies. However, the lower branches corresponding to the prograde precession begin to deviate noticeably when falls below 2/3.

Precession of resonant tilting modes
We have seen that when the epicyclic and vertical frequencies match, the resonant response leads to an enhanced precession of the ring. Indeed, the local resonant tilting mode frequency departs from the orbital rate as ∼ ( ) whilst the non-resonant case goes as ∼ ( 2 ). This highlights the importance of resonances in driving distinct dynamics.
Whilst we have considered the precessional behaviour of pure modes in the above analysis, it is important to understand the long term evolution for general initial conditions. It is useful here to proceed using the Lagrangian linear theory. Similar to the Eulerian analysis discussed in Section 4 we perturb the terms which break the midplane symmetry of the equilibrium. This leads to a simple eigenvalue problem for which the general solution is given by Here, = 0 relates to the perturbative constant of integration and 1,2 = bending wave theory. The constants are then constrained using the initial conditions (0) = ,0 and (0) = ,0 . In order to extract the long timescale precessional behaviour from this local oscillatory solution, we require some way of isolating the global orientation of the ring. This is motivated by thinking about a test particle on an inclined orbit. Locally the vertical motion can be described by = ( − Ω ) where is a complex amplitude and we have assumed a spherically symmetric potential for which = Ω. The magnitude of is proportional to the orbit inclination whilst the argument is related to the longitude of the ascending node. Thus, we may express this quantity in terms of the classic complex tilt variable = + by means of = − 0 . The evolution of then tracks the global evolution of the horizontal angular momentum vector. Ogilvie & Latter (2013) construct the local analogue of the horizontal angular momentum to be where e and e are the radial and azimuthal unit vectors in the local model. Upon normalising against the total angular velocity for a circular orbit, 2 0 Ω, and multiplying by exp ( Ω ) to counteract the local frame rotation, we find In the absence of any torques which drive the evolution of this orbit, angular momentum is conserved and is a constant -essentially describing the fixed amplitude and phase of a harmonic oscillator in terms of its displacement and velocity at any given time.
More generally in our fluid continuum, we will regard the mass weighted average of as our measure of the tilt orientation of the ring. Our Lagrangian reference state is symmetric about the midplane so a vertical averaging extracts the line 0 = 0. Inserting our Lagrangian variables, elements on this line then have the complex amplitude The local measure of warp is then given by the gradient of this line The prefactor 33 / is simply constant in linear theory, corresponding to the width of the ring which scales the magnitude of the tilt. Let us instead focus on the evolution of the reference warp term We now insert the general solution for 31 and the resonant tilting mode frequencies expanded to first order 1,2 = 1± /2. Expanding terms and gathering the leading order contributions shows that This describes a linear transformation of a circular path, so in general the global warp quantity 0 evolves on elliptical tracks. Only when the system is initialised in one of the normal modes, with equal amplitude 31 and 13 in perfect phase or anti-phase, do we see circles. This describes the precession of resonant tilting modes on a timescale ∼ ( −1 ) for arbitrary initial conditions. An example of a general precession track in shown in Fig. 4. The angular coordinate measures the phase of the tilt and thus tracks the changing orientation of the longitude of the ascending node. The radial coordinate then measures the amplitude of the warp which varies as a consequence of the beating phenomenon between the two slightly detuned normal modes. For a ring with aspect ratio = 0.01, a full precessional rotation occurs after 200 orbital timescales. If we expand the resonant frequencies to second order in the aspect ratio, such that 1,2 = 1 ± /2 − 2 /8, this introduces a long timescale prograde bias to the structure. Indeed, we can show that 0 ↦ → 0 exp[ 2 /8] such that the elliptical tracks rotate anticlockwise over time. This longer timescale is visualised in Fig. 4 as the later time yellow tracks have advanced prograde to the earlier purple tracks. This hierarchy of timescales will be important in understanding the evolution of tilted rings in Keplerian systems and the resonant evolution of warped discs as will be explored further in a companion paper. ( 2 ) these elliptical tracks precess in the prograde direction, as indicated by the purple ellipses at earlier times rotating counter-clockwise to the yellow ellipses at later times.

NUMERICAL EXCITATION OF RING MODEL MODES
Having analytically constructed our model and identified the linear modes, we will now demonstrate the excitation of these within a grid based code. Such torus oscillations are difficult to simulate given that the finite extent of the structure leads to zero density regions which are not well described on a fixed grid. However, here we are in fact able to find good agreement between our analytical ODE model and a full hydro solution found using the finite-volume code PLUTO.
PLUTO naturally lends itself to our local ring model as we can invoke the native shearing box module. Furthermore, the FARGO scheme is implemented, wherein the background Keplerian shear flow is subtracted and the code solves for the residual velocity components. This is particularly useful for thin rings where the underlying Keplerian shear flow can become large and would severely limit the integration time-step.

Setting up an equilibrium
In order to avoid problems associated with zero density regions, we set up an isothermal equilibrium ring with a Gaussian density profile. This is as valid as any other choice since we found that the ring model modes are insensitive to the details of the spatial pressure and density structure set by˜and˜. Note however, that the dynamical evolution takes place under an ideal adiabatic energy equation with = 5/3, which does affect the linear breathing modes (see equation (87)).
We will investigate the resonant regime for which the central potential sets = = Ω and = (3/2)Ω. We must then carefully set up an equilibrium compatible with the shearing periodic boundary conditions. This requires that the zonal flow reduces to the Keplerian value at the boundaries. Thus the equilibrium constructed in Section 2.2 is now modified according to the piece-wise velocity field, where is the total width of the numerical domain in and 0 is an arbitrary position which delimits the transition between the linear zonal flow and Keplerian flow. This quadratic modification matches onto the Keplerian flow at the boundary. Although this breaks the linear flow field assumption used to construct the ring model, the dynamics of the modes should be unaffected provided this correction is switched on at suitably large 0 where the density is negligible. Modification of the flow field is accompanied by a change in the radial profiles of density and pressure, as determined by geostrophic balance. For an isothermal equilibrium with = 2 , where 2 is the isothermal sound speed, we have The density only exhibits a slight deviation from the Gaussian in the tails such that where ( ) = exp[−Ω 2 2 /2 2 ] is set by the vertical hydrostatic balance and is the characteristic density value at the centre of the ring. The particular equilibrium used in our experiment is visualised in figure 5. This corresponds to a disc with aspect ratio = 0.14, requiring a flow coefficient = 1.51. For ideal hydrodynamics we are also free to choose our non-dimensional units such that = Ω = 1. The domain is taken to have dimensions ∈ [−30, 30] and ∈ [−5, 5] such that = 60. The piece-wise delimiter 0 is taken to be 25, so the density has dropped off sufficiently, as shown by the red dashed lines in figure 5.

Exciting linear modes
This equilibrium is stable when implemented in PLUTO and acts as a suitable background which can then be perturbed. We would like to excite both tilting and breathing modes and compare with the solution obtained by integrating our ring model ODEs derived in Section 2. To this end we employ an implicit Runge-Kutta Radau solver and numerically integrate the corresponding initial value problem. We will consider two different perturbations. Taking the initial value of 31 = 0.01 excites tilting mode motions in the ring as the mid-plane symmetry is broken. Meanwhile, taking 33 = 0.01 initialises the ring with vertical velocity corresponding to a breathing mode.
In order to compare the subsequent evolution of the ring between the PLUTO and ODE schemes, we require some diagnostic which is sensitive to the small oscillations. We will examine the mass weighted covariance moments which are defined by and characterise the shape of the ellipse. Indeed, one can show that this measure is intrinsically related to the shape matrix components as follows. Consider and convert to Lagrangian coordinates: = 2 11 0 0 + 2 13 0 0 = 2 ( 2 11 + 2 13 ).
Combining these relations allows us to write These are plotted as the dashed red lines in Fig. 6. Meanwhile the numerical covariance results, plotted as solid black lines, are calculated by numerical integration of the mass distribution across the grid. In the left hand panel we see the diagonal covariance moment which tracks the dominant motion for the tilting mode run. We see almost perfect overlap between the ODE and PLUTO runs over the course of 10 orbital periods. The beating phenomenon seen indicates that both tilting modes are excited with the beating interval indicative of the frequency difference between the modes. Indeed, in this near resonant case with a relatively thin ring equation (93) predicts Δ ∼ which corresponds to a beating interval of / = 1/ . For our chosen ring this corresponds to / ≈ 7.1 which agrees well with the pattern observed. Similarly, in the right hand panel we examine the vertical covariance moment and again observe good agreement between PLUTO and the ODE comparison. In this case, the kick to 33 only excites the vertical breathing mode and a single frequency is present.

CONNECTION WITH PREVIOUS THEORY
In this work we have developed a novel ring model for investigating hydrodynamical phenomena. The development of this model, in part, has been inspired by the historical treatment of equilibrium fluid figures. Indeed, the elliptical cylinders constructed in this paper are reminiscent of the classical work on homogeneous ellipsoids, as summarised by Chandrasekhar (1969). Whilst this classical work incorporates self-gravity, which is easier to do in the case of a homogeneous, incompressible fluid, we have so far neglected this effect. Many astrophysical discs are non-self-gravitating, such that their dynamics is dominated by the tidal gravitational potential and the fluid pressure, as in our model. By comparing the pressure term and destabilising self-gravity term in the 2D dispersion relation for linear density waves we find that the mass ratio of the ring to the central star should satisfy / * << ( / ) 2 . However, young massive protoplanetary discs are susceptible to the gravitational instability and in future work we will extend our model to include these neglected terms. Indeed, the modular nature of the Lagrangian formalism allows for additional physics to easily be incorporated.
Roche and Darwin ellipsoids which include tidal effects have since been explored in a host of astrophysical applications. Carter & Luminet (1983) developed an affine model for a star being tidally disrupted within the tidal radius of a black hole. Akin to our Jacobian matrix, they restrict the allowed dynamical degrees of freedom to components of a deformation matrix which then allows for a Lagrangian formalism. This affine model crucially permits compressible deformations of the star, which becomes squashed in a 'pancake' phase as the star passes through pericentre. Subsequently, the dynamical model of viscous tidal disruption developed by Sridhar & Tremaine (1992) also bears formal similarity to our oscillating ring. They examined an ellipsoidal planetisimal described by a time dependent shape matrix in a locally expanded tidal potential. By assuming a flow field linear in the Cartesian coordinates, they reduce the fluid equations of motion and Poisson's equation to a set of coupled differential equations. Whilst their study assumes an incompressible, uniform density ellipsoid, the work of Goodman et al. (1987) models blobby fluid 'planets' as polytropic ellipsoids and finds the equilibrium structure corresponds to linear flow fields.
In fact, their paper is addressing the nonlinear outcome of the nonaxisymmetric breakup of tori due to the Papaloizou-Pringle instability (Papaloizou & Pringle 1983) which may lead one to question the validity of our axisymmetric ring model. However, this instability is suppressed for steep angular momentum gradients (Papaloizou & Pringle 1985) for which our torus has a small aspect ratio. Thus the thin ring regime, in which we are most interested for the application to warped discs, is resistant to such breakups.
More recently, Ogilvie (2018) has developed an an affine model of the dynamics of thin discs, which incorporates vertical degrees of freedom wherein the evolution is described at each radial location by the time-dependent affine transformation of a fluid column. We can show that our ring model is a special case of this more general theory. Specifically the model must make the assumption that the deformation gradient tensor is constant in each fluid column. However, in our model this approximation is exact as the contraction and expansion of the ring is homogeneous with having no spatial dependence. A detailed comparison between the our model and this previous work is presented in Appendix A.
Our analysis of linear modes in Section 4 agrees with the lowest order global modes of slender tori identified by the analysis of Blaes et al. (2006) in the non-relativistic limit. The so called X and inertial modes correspond to our incompressible tilt and shear modes whilst the + and breathing modes correspond to our radial and vertical compressive modes respectively. A more general perturbation of our equilibrium fluid ring yields an eigenvalue problem which permits higher order polynomial eigenfunctions. These are of course excluded by our linear flow field assumption.

CONCLUSION
In this paper we have presented a novel analytical framework for exploring a range of hydrodynamic phenomenon, with particular application to torus oscillations and warped disc theory. We have constructed a local ring model within a shearing sheet approximation, whereby variation of the aspect ratio can probe both the thick torus and thin ring regime. This is introduced from an Eulerian perspective and then using a Lagrangian formalism. We have calculated the linear modes and drawn a close correspondence with linear bending wave theory. The tilting ring within the local model can be interpreted as a precessing inclined disc as seen from a global view, with both prograde and retrograde motion depending on whether the tilt and shear are in phase or anti-phased. The existence of these modes is also identified within a grid based numerical simulation, which highlights the difficulties of treating boundary conditions. Instead our model allows for a dynamical evolution via the comparatively simple task of solving a set of coupled ordinary differential equations. In a future paper we will exemplify the key features of this model by exploring the extremely nonlinear dynamics of a thin torus. This will reveal a strong mode coupling between warping motions and vertical bouncing of the disc. This phenomenon requires further attention in misaligned disc systems with regards to future numerical experiments and observational implications.

DATA AVAILABILITY
Data used in this paper is available from the authors upon reasonable request.

APPENDIX A: COMPARISON WITH AFFINE MODEL
Our ring model presents a special case of the affine model for warped discs as developed by Ogilvie (2018), hereafter OL18. Indeed, our Lagrangian formalism can be derived from this theory upon introducing some simplifying assumptions about the reference state and the nature of the linear Jacobian mapping. In this Section we will compare the two models and identify the necessary restrictions to obtain equations (66) -(71).

A1 Affine mapping
OL18 uses a Lagrangian framework to construct a model for warped discs. This involves mapping material points from a reference state to a dynamical state, under the action of a time dependent transformation, akin to the procedure in Section 3. Whereas we restrict attention to linear transformations, OL18 permits more general affine mappings. In this framework the reference state, described by coordinates ( 0 , 0 ), consists of vertical fluid columns as visualised in Fig. A1. Each fluid column is in hydrostatic equilibrium and centred on 0 = 0. The distance along each column is parameterised by the scaled coordinate where 0 ( 0 ) is the scale-height for the fluid column labelled by the coordinate 0 . The affine mapping then consists of a translation, followed by stretching and rotation of the columns. This is described by in line with equation (19) in OL18. Here X = ( , , ) denotes the position vector of the column centres in the dynamical state Figure A1. Left: The reference state consists of vertical fluid columns with the height chosen to represent circular density contours. The centre of these columns is denoted by black dots lying on the midplane. Right: Under the linear Jacobian transformation, each of these columns undergoes a translation which stretches and tilts the midplane line, followed by a rotation and stretching of the columns. This results in generally elliptical density contours.
and H = ( , , ) is a vector scaling and rotating each fluid column. The resulting velocity is then given by where v = X/ and w = H/ . Thus each column has six degrees of freedom (X, H), which can be identified with our six evolutionary equations for the Jacobian coordinates. In our model, we focus on the dynamics of the largest-scale motions, which correspond to purely linear transformations of the reference state under the Jacobian mapping (39). By assuming linear transformations, we ensure that the deformations of the columns are homogeneous and that they cannot clash by tilting into one another. The midplane line through the centre of each column is then mapped to a straight line through the origin as described by X = ( 11 , 21 , 31 ) 0 + (0, 1, 0) 0 .
The mapping of each column for varying 0 at fixed 0 then gives the correspondence with the remaining Jacobian coordinates H = 0 ( 13 , 23 , 33 ). (A5)

A2 Density and pressure structure
Our assumptions concerning the elliptical ring profiles will also impose conditions on the nature of the pressure and density structure entering the model of OL18. Indeed, our reference state consists of circular pressure and density contours, so the fluid column properties must be chosen to reproduce this. It should also be noted that the framework developed by OL18 considers the reference state to be in vertical hydrostatic equilibrium. We see that this is compatible with equation (62) By vertically integrating the exemplar density and pressure relationships proposed for our ring model (e.g. equation (30)) we can obtain the homogeneous, polytropic and isothermal distributions for Σ 0 and 0 which will automatically satisfy conditions (A6) and (A8).

A3 Equations of motion
With the correspondence between model quantities made, all that remains is to compute the equations of motion. The Lagrangian presented by OL18 equation (68) is slightly modified to account for the local shearing sheet perspective as used in the ring model. Indeed, as per our Lagrangian integral (49), we must include the Coriolis terms which arise from expanding about an orbiting reference frame. The Lagrangian density then becomes Then the dynamics is governed by the Euler-Lagrange field equa- which give the six equations of motion 2 2 = 2Ω Here = − / is the projected scale-height quantity equivalent to The dynamical state density and pressure are related to the prescribed reference state quantities by means of the Jacobian determinants Then inserting equations (A4) -(A5) and making use of the relation (A8), reduces equations (A12) -(A17) to our ring model equations (66) -(71) for the choiceˆ0 = 2 2 . This paper has been typeset from a T E X/L A T E X file prepared by the author.