Nonlinear resonant torus oscillations as a model of Keplerian disc warp dynamics

Callum W. Fairbairn1 and Gordon I. Ogilvie,1

1Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences,

Wilberforce Road, Cambridge CB3 0WA, UK E-mail: cwf29@cam.ac.ukE-mail: gio10@cam.ac.uk (Accepted XXX. Received YYY; in original form ZZZ)

Abstract

Observations of distorted discs have highlighted the ubiquity of warps in a variety of astrophysical contexts. This has been complemented by theoretical efforts to understand the dynamics of warp evolution. Despite significant efforts to understand the dynamics of warped discs, previous work fails to address arguably the most prevalent regime – nonlinear warps in Keplerian discs for which there is a resonance between the orbital, epicyclic and vertical oscillation frequencies. In this work, we implement a novel nonlinear ring model, developed recently by Fairbairn and Ogilvie, as a framework for understanding such resonant warp dynamics. Here we uncover two distinct nonlinear regimes as the warp amplitude is increased. Initially we find a smooth modulation theory which describes warp evolution in terms of the averaged Lagrangian of the oscillatory vertical motions of the disc. This hints towards the possibility of connecting previous warp theory under a generalised secular framework. Upon the warp amplitude exceeding a critical value, which scales as the square root of the aspect-ratio of our ring, the disc enters into a bouncing regime with extreme vertical compressions twice per orbit. We develop an impulsive theory which predicts special retrograde and prograde precessing warped solutions, which are identified numerically using our full equation set. Such solutions emphasise the essential activation of nonlinear vertical oscillations within the disc and may have important implications for energy and warp dissipation. Future work should search for this behaviour in detailed numerical studies of the internal flow structure of warped discs.

keywords:

hydrodynamics – waves – accretion discs ††

1 Introduction

1.1 Astrophysical motivation

The traditional model for astrophysical discs assumes the simplest coplanar configuration with fluid streamlines on circular orbits. However, there has been growing interest in the behaviour of these systems when they become distorted by a warp. This introduces a radial variation in the inclination of the circular streamlines which might drastically alter the disc dynamics. Indeed, there is an ever expanding host of observational evidence for warped discs in a variety of contexts, which demands an improved theoretical understanding.

Warped discs have been indirectly inferred from the long period luminosity variations of ‘superorbital’ X-ray binary systems where a precessing warped structure periodically obscures light from a central source (e.g. Katz, 1973; Kotze & Charles, 2012). In a similar vein, intensity deficits in the outer regions of protoplanetary discs may be explained by shadows cast by an inner tilted precessing disc (e.g. Debes et al., 2017; Muro-Arena, G. A. et al., 2020). Comparison of radiative models with observed shadows have suggested even more extreme inclination variations in transition discs, where large radial gaps divide the inner and outer regions (e.g. Marino et al., 2015; Pinilla et al., 2015; Stolker et al., 2016; Benisty et al., 2017; Casassus et al., 2018). In some systems these distinct rings are thought to form by disc tearing and breaking, as found in several numerical simulations. Nixon & King (2012) find that Lense-Thirring torque around a spinning black hole can induce disc breaking whilst Facchini et al. (2013) find breaking of a circumbinary disc when it is sufficiently tilted with respect to the plane of the binary. Radiative post-processing of such structures produces images capable of explaining observed precessing shadows (Facchini et al., 2017). More recently, there has been an observation of the spectacular triple star system GW Orionis wherein gravitational effects may have torn the disc into independently precessing rings (Kraus et al., 2020).

These indirect cases have been complemented by direct observations of maser emission lines tracing warped galactic midplanes, as for the spiral galaxy NGC 4258 (M106) (Miyoshi et al., 1995). More recently, the Atacama Large Millimeter/submillimeter Array (ALMA) has measured dust emission in young protostellar discs with misaligned inner and outer regions (Sakai et al., 2019). ALMA has also traced gas kinematics through CO and HCO+ molecular line emission which is consistent with warped inner regions (Rosenfeld et al., 2012; Loomis et al., 2017).

1.2 Warped disc theory

In order to understand this host of observational phenomena, we require theoretical models for the evolution of warped discs. Much of the mathematical language underpinning these was laid down by the work of Petterson (1977a, b) and Hatchett et al. (1981) wherein the warp is described as a series of nested, interacting rings. Understanding the evolution is then a question of determining the time dependence of the inclination of each ring. Petterson (1977a) included a viscous torque between the rings which naturally led to the diffusion of warp on a viscous timescale. However, Papaloizou & Pringle (1983) showed that this simple model neglected the internal flow dynamics established by the warp itself, which enhance the angular momentum transport and accelerate the warp evolution. They found that the evolution is diffusive (but faster than the viscous timescale) when α>H/R, where α is the Shakura-Sunyaev viscosity parameter and H/R is the angular semi-thickness of the disc. Later, Papaloizou & Lin (1995) and Lubow & Ogilvie (2000) investigated the nearly inviscid regime for which α<H/R. Here the linearised evolution takes the form of a non-dispersive bending wave in Keplerian discs and a dispersive bending wave when the degeneracy between the epicylic and vertical frequencies is broken.

All these models focus on linear warps, but of course it is crucial to extend this understanding into the nonlinear regime where there are observational consequences. Ogilvie (1999) improved on the efforts of Pringle (1992) and developed a self consistent, fully nonlinear model of diffusion in Keplerian discs and bending waves in non-Keplerian discs. This theory has been shown to agree well with numerical simulations of warps (Lodato & Price, 2010). Despite such success, this model is unable to describe arguably the most important case – inviscid Keplerian discs where the epicyclic motion is resonantly driven by the warping geometry. Ogilvie (2006) attempted to explore this missing regime by performing a weakly non-linear analysis of Keplerian bending waves. However, the strongly nonlinear case still lacks a complete theory and requires further attention.

1.3 Outline of this paper

In order to address this problem we previously introduced a novel ring model, capable of describing the fully nonlinear hydrodynamic oscillations of an ideal, non-self gravitating torus (Fairbairn & Ogilvie 2021, hereafter Paper I). We found that small amplitude tilting oscillations in this local model could be identified with global linear bending waves. Indeed, our shearing box formulation effectively captures the evolution of a warp as we zoom in on a localised patch of the disc. In this picture, ring oscillations over the fast orbital timescale track the azimuthal changes in the disc geometry as the shearing box moves around the orbit. Thus tilting motions are associated with streamlines on inclined orbits and hence warps about the midplane. In this paper we aim to advance this theory and examine the nonlinear extension of bending modes. We begin by summarising the derivation and interpretation of the ring model equations in section 2. We then motivate our analytical progress by performing some numerical runs in section 3. We will find that as the initialised warp amplitude is increased, two distinct nonlinear regimes arise. We will tackle the first in section 4 by using an averaged Lagrangian method which describes a smooth modulation of the warp amplitude and phase. This behaviour drastically changes beyond some critical warp amplitude, at which point the disc enters into an extreme bouncing regime. To this end we develop an impulsive bouncing theory in sections 5 and 6 which predicts a family of highly compressive, warped solutions. These analytical predictions are confirmed within the full ring model equation set in section 7 before we discuss the implications for astrophysical discs and warp theory in section 8.

2 Summary of ring model

2.1 Model assumptions

In Paper I we constructed a ring model for oscillating tori which will prove a useful framework in our current study. In this section we will briefly revisit the key assumptions and the resulting equations. Following the standard shearing box construction (e.g. Hill, 1878; Hawley et al., 1995), we expand the ideal hydrodynamic equations about a local circular reference orbit at r₀ with angular velocity @boldsymbol_0=(r_0)@boldsymbolz, assuming an axisymmetric potential Φ(r,z). This orbit has an attached, co-rotating coordinate system (x,y,z) which is defined by x=(r-r₀), y=r₀(ϕ-Ω₀ t) and z=z, such that x, y and z are the radial, azimuthal and vertical directions respectively. This leads to the usual shearing box equations

D@boldsymbolu+2@boldsymbol_0×@boldsymbolu=-_t-(1)/(ρ) p,

(1)

where

D=_t+@boldsymbolu·

(2)

is the Lagrangian derivative, @boldsymbolu is the velocity, p is the pressure and ρ is the density. The tidal potential is expanded as

Φₜ=-Ω₀ S₀ x²+1/2ν₀² z²,

(3)

where S₀=-(rdΩ/dr)₀ is the orbital shear rate and ν₀²=(∂_zz Φ)₀ is the square of the vertical oscillation frequency of a test particle perturbed from its circular orbit. Similarly, inertial restorative forces cause a natural radial oscillation about this orbit which is defined by the epicyclic frequency κ₀ given by,

κ₀²=2Ω₀(2Ω₀-S₀).

(4)

Henceforth we will drop the subscript on the orbital velocity, shear rate, and vertical/epicyclic frequencies in order to simplify our notation. We restrict our basic model to an isentropic energy equation with adiabatic index γ but allow for compressibility. Density ρ and pressure p are then governed by

Dρ=-ρΔ, Dp=-γ pΔ,

(5) (6)

where

Δ=·@boldsymbolu

(7)

is the velocity divergence. We look for axisymmetric dynamical solutions for which the density and pressure are described by a common materially invariant function f(x,z,t). This allows us to perform the separation of variables

ρ=ρ̂(t)ρ̃(f), p=p̂(t)p̃(f).

(8) (9)

We then enforce linear flow fields, which capture the lowest order global motions supported by the tori, such that

uᵢ=Aᵢⱼ xⱼ,

(10)

where Aᵢⱼ is a time-dependent, square flow matrix. Since this linear flow maps ellipses to ellipses, the materially conserved function f should be a quadratic function of the coordinates such that

f=C-1/2Sᵢⱼ xᵢ xⱼ,

(11)

where C is some constant and Sᵢⱼ(t) is a time dependent, positive-definite shape matrix with Sᵢ₂=S₂ᵢ=0 in the y-independent case. Thus contours of equal density and pressure trace out elliptical contours, which are described by the time evolution of the shape matrix. Imposing the material conservation condition for f at all points in space requires that

dₜ Sᵢⱼ+Sᵢₖ Aₖⱼ+Sⱼₖ Aₖᵢ=0,

(12)

which gives three independent ODEs for S₁₁, S₁₃ and S₃₃. Meanwhile, the flow matrix evolution is deduced by inserting our assumptions into the equation of motion and gathering terms linear in each spatial coordinate. We must also impose dp̃/df=ρ̃ so that the pressure gradient term is compatible with this linear form in the coordinates. This gives rise to six ODEs,

dₜ A₁₁+A₁₁²+A₁₃ A₃₁-2Ω A₂₁=2Ω S+T̂S₁₁, dₜ A₁₃+A₁₁ A₁₃+A₁₃ A₃₃-2Ω A₂₃=T̂S₁₃, dₜ A₂₁+A₂₁ A₁₁+A₂₃ A₃₁+2Ω A₁₁=0, dₜ A₂₃+A₂₁ A₁₃+A₂₃ A₃₃+2Ω A₁₃=0,

dₜ A₃₁+A₃₁ A₁₁+A₃₃ A₃₁=T̂S₁₃, dₜ A₃₃+A₃₁ A₁₃+A₃₃²=-ν²+T̂S₃₃,

(13) (14) (15) (16) (17) (18)

where T̂(t)=p̂/ρ̂ is a characteristic temperature. This evolves according to

dₜT̂=-(γ-1)T̂Δ,

(19)

where Δ=A₁₁+A₃₃ is the velocity divergence.

2.2 Lagrangian perspective

This model may alternatively be reformulated from a Lagrangian perspective. We construct a material mapping of points from an arbitrary, stationary reference state, denoted by @boldsymbolx_0=(x_0,y_0,z_0), to the dynamical state @boldsymbolx by means of the linear transformation

@boldsymbolx_0@boldsymbolx: x_i=J_ijx_0,j,

(20)

where Jᵢⱼ is the time dependent Jacobian matrix. For the assumed axisymmetric setup J₁₂=J₃₂=0 and J₂₂=1, so the 6 remaining independent components describe the linear flow field uᵢ=J̇ᵢⱼ xⱼ. We load mass in the reference state such that the materially conserved density and pressure contours lie on circles with radius LR. Here, R=√(2(C-f)) is a dimensionless radius measured in units of the characteristic length L, which arises when taking the second mass weighted moment of the reference distribution. Our separation of variables then becomes

ρ_0(@boldsymbolx_0)=ρ_0ρ(R(@boldsymbolx_0)), and p_0(@boldsymbolx_0)=p_0p(R(@boldsymbolx_0)),

(21)

where ρ₀ and p₀ denote the density and pressure in the reference state whilst ρ̂₀ and p̂₀ are characteristic density and pressure factors.

In Paper I we outline the construction of a Lagrangian composed of the kinetic, rotational, internal and potential energies

L=1/2(J̇₁₁²+J̇₁₃²+J̇₂₁²+J̇₂₃²+J̇₃₁²+J̇₃₃²)-(T̂₀)/((γ-1)Jᵞ⁻¹ L²) +Ω S(J₁₁²+J₁₃²)-1/2ν²(J₃₁²+J₃₃²)+2Ω(J₁₁ J̇₂₁+J₁₃ J̇₂₃),

(22)

where J=det(Jᵢⱼ)=J₁₁ J₃₃-J₁₃ J₃₁ is proportional to the area of the elliptical cross-section of the ring and T̂₀=p̂₀/ρ̂₀ is a characteristic temperature. The usual Euler-Lagrange equations then give the dynamical equations

J_11=2J_21+2 SJ_11+T_0J^γL^2J_33, J_13=2J_23+2 SJ_13-T_0J^γL^2J_31, J_21=-2J_11, J_23=-2J_13,

J̈₃₁=-ν² J₃₁-(T̂₀)/(Jᵞ L²)J₁₃, J̈₃₃=-ν² J₃₃+(T̂₀)/(Jᵞ L²)J₁₁.

(23) (24) (25) (26) (27) (28)

The conservation of angular momentum gives rise to the integrability of equations (25) and (26) which allows us to reduce this system to 4 second-order, coupled ODEs,

J̈₁₁+κ² J₁₁=2C_z+(T̂₀)/(Jᵞ L²)J₃₃, J̈₁₃+κ² J₁₃=2Cₓ-(T̂₀)/(Jᵞ L²)J₃₁, J̈₃₁+ν² J₃₁=-(T̂₀)/(Jᵞ L²)J₁₃, J̈₃₃+ν² J₃₃=(T̂₀)/(Jᵞ L²)J₁₁,

(29) (30) (31) (32)

where Cₓ and C_z represent the constants arising from circulation conservation and we have made use of the definition of the epicyclic frequency κ to eliminate the shear rate S.

2.3 Physical interpretation and connection with warping

It is worth emphasising the physical intuition behind these variables. J₁₁ and J₃₃ describe a radial and vertical stretching of the ring respectively, capable of capturing breathing motions. J₁₃ corresponds to the vertical shear of horizontal flows whilst J₃₁ describes the ring tilting as one moves radially outwards (for helpful visualisations refer to Paper I). The Lagrangian form of the equations clearly elucidates the oscillatory structure underlying the ring system. The left-hand side terms correspond to free harmonic oscillators, whilst on the right-hand side, matters are complicated by the pressure terms which couple the oscillators together.

As discussed in Paper I, the off-diagonal Jacobian elements act to break the midplane symmetry of the elliptical rings and can be identified with bending modes. Indeed, these tilting motions, as observed within the shearing box orbital frame, can be reinterpreted in a global, non-rotating reference frame as a series of nested circular orbits with a radially dependent inclination. This tilting of streamlines may be thought of as an m=1 azimuthal mode and hence associated with a warped structure (Ogilvie & Latter, 2013). To illustrate this, imagine setting up a line of test particles on circular orbits with a radial, linear variation in inclination about the reference shearing box orbital plane. In a Keplerian potential these orbits are closed and describe a fixed, warped annulus. However, when viewed from the rotating shearing box frame, the line of test particles rock up and down, simply tracking the geometry of the tilted annulus.

More generally, the inclusion of pressure in a gaseous disc couples these particle orbits and may introduce some precession of the streamlines. As the global structure slowly rotates, the oscillation period in the shearing box frame will depart from the orbital period. These two perspectives are connected by Doppler shifting the m=1 warping mode such that

ωₚ=Ω-ω,

(33)

where ωₚ is the precessional frequency in the global frame and ω is the frequency of the tilting mode in the local model. Thus ωₚ<0 and ωₚ>0 correspond to retrograde and prograde precessing warped structures respectively.

In order to quantitatively connect this with global warped theory, we will introduce a local measure of the warp amplitude. Consider a test particle on an inclined orbit such that it undergoes vertical oscillations in the local model according to z= (Ze^-i t) where Z is a complex amplitude. Then the magnitude of Z is proportional to the orbit inclination whilst the argument is related to the longitude of ascending node. Thus we may express this quantity in terms of the classic complex tilt variable W=lₓ+il_y such that Z=-r₀ W. Here @boldsymboll is the unit tilt vector, pointing normal to the circular orbits of the test particle, which clearly encapsulates the amplitude and phase of the z motion about the reference plane. The value of Z is extracted from the local model via

$$Z=\left(z+\frac{i\dot{z}}{\Omega}\right)e^{i\Omega t}.$$

(34)

Consider a set of particles along the midplane of the ring z₀=0 labelled by reference coordinate x₀, such that z=J₃₁ x₀ and

$$Z=\left(J_{31}+\frac{i}{\Omega}\dot{J}_{31}\right)x_{0}e^{i\Omega t}.$$

(35)

Then along this line x=J₁₁ x₀ such that the warp amplitude, defined as the gradient ψ≡dZ/dx, is given by

$$\psi=\frac{1}{J_{11}}\left(J_{31}+\frac{i}{\Omega}\dot{J}_{31}\right)e^{i\Omega t}.$$

(36)

This connection between the local tilting modes and the global warped perspective is crucial for understanding the solutions derived later.

3 Numerical motivation

The simplicity of the linear harmonic form presented by the left hand side of equations (29) – (32) makes them an attractive framework to explore the nonlinear effects introduced by the right hand side pressure terms. The large number of degrees of freedom and significant nonlinearity introduced by the pressure couplings means it is instructive to first numerically solve this system of ODEs. This will reveal a rich range of dynamical behaviour. Using an implicit Runge-Kutta integrator we test a range of tilted initial conditions which break the midplane symmetry of the ring. As we increase the amplitude of the tilt and depart further from equilibrium, we identify two distinct nonlinear warping regimes which will motivate our analysis in subsequent sections.

3.1 Tilting setup

As demanded by the gap in the current warped disc theory, we will focus on the resonant regime for which the epicyclic, vertical and orbital frequencies are all equal with κ=ν=Ω. Without loss of generality we can choose our units such that Ω=1 and L=1. We assume a typical adiabatic index γ=5/3 and set up a thin equilibrium ring with J₁₁=100 and J₃₃=1 such that the aspect ratio is given by ε=J₃₃/J₁₁=0.01. This choice ensures that the length scale of the warp is much longer that the disc scale-height. As described in Paper I, the vertical equilibrium is established via the hydrostatic balance described by equation (32) which sets the value of T̂₀=ε Jᵞ. The finite width of the ring then incurs a radial pressure gradient which is balanced by an enhanced shear. This manifests as a reduced value of the Bjerknes circulation constant C_z=(1/2)J₁₁(1-ε²), as the local shear flow vorticity component counteracts the global rotational vorticity.

In Paper I, we investigated linear tilting modes by slightly perturbing this equilibrium ring and found close correspondence with linear bending-wave theory. We now gain a foothold on the transition to nonlinear tilting dynamics by releasing the ring from increasing tilt angles θₜ. We simply rotate the equilibrium ring so that θₜ corresponds to the angle between the ellipse’s major axis and reference plane measured in radians. Releasing from this rotated state presents a general configuration which naturally engages the warping motions of interest. In order to interpret the change in the dynamics as the amplitude is increased we will examine the warp amplitude ψ as defined in equation (36). This is a useful diagnostic for understanding the tilting and precession of the ring as the warped structure evolves. Linear bending waves generally trace out elliptical paths in a polar plot of ψ. As the amplitude of the tilting perturbation increases we expect the nonlinearities to significantly distort this picture, as we shall soon see.

Figure 1: The polar plots of the warp amplitude ψ tracked over 200 orbital periods for different initial tilt angles. The complex nature of this variable means that the plots display the evolution of both the magnitude and the phase of the warp. These govern the linear rate at which fluid streamlines are tilted from the reference plane when moving radially, and the global orientation of this tilting, respectively. Upper panel: θₜ=0.01 is a small tilt. The elliptical track is indicative of the linear bending mode regime. Middle panel: θₜ=0.14 is a moderate tilt. The elliptical track is smoothly distorted as amplitude and phase of tilt and shear oscillators are modulated on secular timescales. Lower panel: θₜ=0.15 is a critical tilt. The smooth track suddenly changes behaviour as the vertical oscillator is resonantly driven into a bouncing regime and feedback onto the warp occurs impulsively.

3.2 Transition to nonlinearity

For small amplitude θₜ, the J₁₃ and J₃₁ oscillators exhibit a beating pattern as both the in-phase and anti-phased linear tilting modes are excited by a general initial condition. This corresponds to elliptical tracks traced out by the warp amplitude as seen in Paper I. Here we observe that for a tilt angle of θ=0.01, ψ traces a squashed elliptical track as seen in the upper panel of Fig. 1. This path is indicative of the secular precession of the tilted ring structure over many orbital timescales. As the initial tilt amplitude is increased, the system smoothly extends into the nonlinear regime. Whilst the shear and tilt oscillators remain largely harmonic in their behaviour, the vertical oscillator J₃₃ is driven to nonlinear amplitudes and becomes dynamically important. This nonlinearity feeds back onto the warp, driving a slow modulation of the phase and amplitude of the tilt and shear oscillators. This distorts the linear warp amplitude elliptical tracks into more interesting configurations, as shown for the θₜ=0.14 run in the middle panel of Fig. 1. The Jacobian variables for this run are plotted in the upper four panels of Fig. 2.

3.3 Critical onset of bouncing regime

When θ_t 0.15 a dynamically distinct behaviour arises. Note that this critical angle generally depends on the system parameters i.e. ε and γ. The lower panel in Fig. 1 plots ψ when the ring is released from this initial tilt and shows a rapid, possibly chaotic evolution. To gain further insight, the individual Jacobian components are plotted in the lower four panels of Fig. 2. Comparison with the θₜ=0.14 run in the upper four panels emphasises a drastically different behaviour. This demonstrates a resonant coupling between the tilting motions associated with the J₃₁ and J₁₃ components and the breathing motions indicated by the J₃₃ component. We see that the initial beating envelopes of the tilt and shear terms are disrupted as their combined effect drives a growth in the breathing motion. This is shown by the large amplitude, compressive bumps in the lower right panel. The system enters into a quasi-periodic regime with strong mode coupling between the vertical breathing and warping oscillations. The driving of such extreme breathing modes has also been separately recognised in the periodically forced scale heights associated with elliptical fluid flows (Ogilvie & Barker, 2014).

Whilst in the smooth nonlinear regime the breathing and warping modes remain largely disconnected reservoirs of energy, in this compressive nonlinear phase the pressure couplings facilitate a large energy exchange flowing back and forth between these motions. This is visualised clearly in Fig. 3 which shows how the energy is partitioned between the different modes over time. The red line plots the energy terms in the Lagrangian corresponding to the warping motions i.e. the kinetic and potential energies involving J₁₃ and J₃₁. Meanwhile, the blue line plots the kinetic and potential energies of the J₃₃ oscillator plus the contribution from the internal energy. It is natural to combine the internal energy with the breathing mode since only compressive motions can heat the ring. Indeed, in linear theory the tilting modes are incompressible and internal energy is conserved. We see that both lines are essentially symmetric about the average energy, denoted by the black dashed line. A large dip in warping energy is balanced by an increase in breathing energy and vice versa. This emphasises the mode coupling channel which is clearly active. Furthermore, we note the red and blue lines appear to vary in a step like manner. This is not an artefact of numerical resolution but in fact a key part of the dynamical behaviour. Each step coincides with a compression of the ring where J₃₃ is squashed. At these discrete times, the cross sectional area of the ring is small and the determinant value J is minimised. It is at these instances that the pressure coupling terms on the right hand side of equations (29) – (32) dominate and allow for an impulsive exchange of energy. It is this impulsive coupling mechanism that will motivate our analytical progress in this regime in the following sections.

Figure 2: The numerically integrated solution to equations (29) – (32) for resonant runs with κ=ν=Ω=1 and γ=5/3. The ring is initialised with aspect ratio ε=0.01 and then rotated from this equilibrium state by θₜ. This measures the angle between the major axis of the ellipse and the x-axis. The J₁₃ and J₃₁ panels encapsulate the shear and tilting motions respectively whilst the J₃₃ panel captures the compressive breathing motions. Upper four panels: θₜ=0.14 corresponding to the middle panel of Fig. 1 shows the smooth modulation regime. Lower four panels: θₜ=0.15 corresponding to the lower panel of Fig. 1 shows extreme vertical bouncing motions.

Figure 3: Comparison of the energy partitioning between the breathing and warping motions for the run with θₜ=0.15. The red line plots the tilting energy contribution appearing in the Lagrangian, Eₜᵢₗₜ=1/2(J̇₁₃²+J̇₃₁²+J₁₃²+J₃₁²)-2Cₓ J₁₃. The blue line plots the vertical breathing contributions E_breathe=1/2(J̇₃₃²+J₃₃²)+(T̂₀)/((γ-1)Jᵞ⁻¹). The dashed black line plots the average between these two energies (L_tilt+L_breath)/2.

4 Smooth modulation theory

4.1 Asymptotic scalings

We first confront the smooth nonlinear regime where we expect a secular modulation of the linear oscillatory solutions. For a thin ring we anticipate that the radial breathing motions are not dynamically important and so we ignore equation (23) and set J₁₁ to be constant. For the equilibrium ring with small aspect ratio ε we have the characteristic temperature T̂₀=ε Jᵞ. The scale invariance of ideal hydrodynamics means we are free to adopt a reference state with area of order unity, so we take the determinant J∼ O(1) and T̂₀=ε which directly introduces a small parameter into the governing equations. To facilitate this area scaling we will take the radial and vertical deformations to be J₁₁ ∼ 0(ε⁻¹⁄²) and J₃₃ ∼ O(ε¹⁄²) respectively. We are interested in exploring nonlinear warps, so adopt the scalings J₁₃ ∼ J₃₁ ∼ O(1) such that they contribute at leading order to J. Inserting these scalings into the reduced Lagrangian for the tilt, shear and vertical oscillators is

$$L=\underbrace{\frac{1}{2}\left(\dot{J}_{13}^{2}+\dot{J}_{31}^{2}-J_{13}^{2}-J_{31}^{2}\right)}_{L_{0}\sim O(1)}+\underbrace{\frac{1}{2}\left(\dot{J}_{33}^{2}-J_{33}^{2}-\frac{\epsilon}{(\gamma-1)J^{\gamma-1}}\right)}_{\epsilon L_{1}\sim O(\epsilon)}.$$

(37)

We see that the Lagrangian is split into a leading order component which just describes harmonic motion of the tilt and shear. At higher order we see the contribution from the vertical oscillator kinetic, potential and internal energies. Notably the tilt and shear are coupled to the vertical oscillator through the internal energy term and will drive a slow modulation of the harmonic motion phase and amplitude over longer timescales. To capture the fast harmonic motion and the slow evolution owing to the nonlinearities we introduce the multiple timescales expansion

$$J_{11} =\epsilon^{-1/2},\qquad J_{13} =J_{13,0}(t,T)+\epsilon J_{13,1}(t,T)+O(\epsilon^{2}),\qquad J_{31} =J_{31,0}(t,T)+\epsilon J_{31,1}(t,T)+O(\epsilon^{2}),\qquad J_{33} =\epsilon^{1/2}J_{33,0}(t,T)+O(\epsilon^{3/2}),$$

(38) (39) (40) (41)

where T=ε t is a slow timescale treated as an independent parameter. Thus the full time derivatives become

$$d_{t}=\partial_{t}+\epsilon\partial_{T},\quad d_{tt}=\partial_{tt}+2\epsilon\partial_{tT}+\epsilon^{2}\partial_{TT}.$$

(42)

Inserting this expansion into the dynamical equations (30)–(32) yields a hierarchy of equations ordered in powers of ε. One should note that the scale invariance of the equations of motion means that we are in fact free to stretch the results provided the underlying aspect ratio is preserved. This scale invariance may be parameterised relative to the width of the ring J₁₁ so the dynamics is similar if we re-scale variables such that J₁₃ and J₃₁ are of order O(ε¹⁄² J₁₁), whilst J₃₃ ∼ O(ε J₁₁). This is important to remember later on when comparing our theory to general numerical runs where the scaling of the elliptical area measure J is not necessarily of order unity.

4.2 Modulation equations

As anticipated, at leading order O(ε) we have

$$\partial_{tt}J_{13,0}+J_{13,0}=0\quad\text{and}\quad\partial_{tt}J_{31,0}+J_{31,0}=0.$$

(43)

These have harmonic solutions

$$J_{13,0}=\Re\left[A(T)e^{-it}\right]\quad\text{and}\quad J_{31,0}=\Re\left[B(T)e^{-it}\right],$$

(44)

where $\Re$ denotes the extraction of the real part and A(T) and B(T) are complex amplitudes encoding the slow modulation of oscillator amplitude and phase. At order O(ε¹⁄²) we obtain the leading order equation for the vertical oscillator

$$\partial_{tt}J_{33,0}+J_{33,0}=\frac{1}{H^{\gamma}},$$

(45)

where H≡J₃₃,₀-J₁₃,₀ J₃₁,₀. This equation may be tackled by changing variables in favour of H such that

$$\partial_{tt}H+H-\frac{1}{H^{\gamma}}=-\partial_{tt}(J_{13,0}J_{31,0})-J_{13,0}J_{31,0},$$

(46)

showing that the compressional motion is driven by the product of the tilt and shear. The right-hand side of this equation is periodic with frequency 2 (i.e. twice the orbital frequency). Periodic solutions with frequency 2 are possible for a certain range of forcing amplitudes, as we shall discuss in Section 4.4 below. We assume here that the solutions are indeed periodic in t, rather than the more general quasi-periodic solutions that include a free oscillation as well as the forced one. In the meantime we will expand to next order in the aspect ratio hierarchy so at O(ε¹) we have

$$\partial_{tt}J_{13,1}+J_{13,1} =-2\partial_{tT}J_{13,0}-\frac{J_{31,0}}{H^{\gamma}}\equiv f_{13,1},\qquad \partial_{tt}J_{31,1}+J_{31,1} =-2\partial_{tT}J_{31,0}-\frac{J_{13,0}}{H^{\gamma}}\equiv f_{31,1}.$$

(47) (48)

On the left hand side we see the linear operator ∂ₜₜ+1 which yields complementary harmonic solutions. We will also denote the right hand side forcing terms as f₁₃,₁ and f₃₁,₁. A necessary condition for periodic solutions requires that the forcing on the right hand side contains no e^± it resonant Fourier components. This is equivalent to the Fredholm solvability conditions

$$\left\langle f_{13,1}e^{\pm it}\right\rangle=\left\langle f_{31,1}e^{\pm it}\right\rangle=0,$$

(49)

where

$$\left\langle\cdot\right\rangle\equiv\frac{1}{2\pi}\int_{0}^{2\pi}\cdot\,dt$$

(50)

denotes averaging over the fast orbital timescale. Evaluating these conditions yields

$$-i\partial_{T}A=-\left\langle\frac{J_{31,0}}{H^{\gamma}}e^{it}\right\rangle\quad\text{and}\quad-i\partial_{T}B=-\left\langle\frac{J_{13,0}}{H^{\gamma}}e^{it}\right\rangle$$

(51)

which gives the evolution of the complex amplitudes over secular timescales based on the fast averaging of the lower order equations.

4.3 Averaged Lagrangian formulation

Since we are dealing with ideal hydrodynamics as derived from a variational principle, we anticipate that the averaged terms can in fact be related to the averaged Lagrangian. This idea was first introduced by Whitham (1965) with application to wave trains propagating through a slowly varying background and has applications in a wide variety of contexts. Returning to the forced vertical oscillator described by equation (45) we see this can be derived from a Lagrangian

$$L_{10}=\frac{1}{2}(\dot{J}_{33,0}^{2}-J_{33,0}^{2})-\frac{H^{-(\gamma-1)}}{\gamma-1},$$

(52)

which is the leading order O(ε¹) contribution from the vertical part of the full Lagrangian. The Lagrangian explicitly depends on time and the forcing parameters A and B through the product of tilt and shear oscillations appearing in H, and implicitly through the dependence of the J₃₃,₀ solution as forced by the warp. These complex amplitudes encode two degrees of freedom each, encapsulating the amplitude and phase, so the complex conjugated quantities A̅ and B̅ may also be treated as independent quantities. Thus consider L₁₀=L₁₀(A,A̅,B,B̅) and first compute

$$\frac{\partial L_{10}}{\partial\bar{A}}=-\frac{H^{-\gamma}}{2}e^{it}J_{31,0}+H^{-\gamma}\frac{\partial J_{33,0}}{\partial\bar{A}}+\dot{J}_{33,0}\frac{\partial\dot{J}_{33}}{\partial\bar{A}}-J_{33,0}\frac{\partial J_{33,0}}{\partial\bar{A}}.$$

(53)

Using equation (45) to replace H⁻ᵞ in the second right-hand side term and averaging over the orbital period yields

$$\frac{\partial}{\partial\bar{A}}\left\langle L_{10}\right\rangle=-\left\langle\frac{H^{-\gamma}}{2}e^{it}J_{31,0}\right\rangle+\left\langle\dot{J}_{33,0}\frac{\partial\dot{J}_{33,0}}{\partial\bar{A}}+\ddot{J}_{33,0}\frac{\partial J_{33,0}}{\partial\bar{A}}\right\rangle.$$

(54)

Integrating by parts shows that

$$\left\langle\dot{J}_{33,0}\frac{\partial\dot{J}_{33,0}}{\partial\bar{A}}+\ddot{J}_{33,0}\frac{\partial J_{33,0}}{\partial\bar{A}}\right\rangle=0,$$

(55)

and so we have

$$\frac{\partial}{\partial\bar{A}}\left\langle L_{10}\right\rangle=-\left\langle\frac{H^{-\gamma}}{2}e^{it}J_{31,0}\right\rangle.$$

(56)

Similarly we find that

$$\frac{\partial}{\partial\bar{B}}\left\langle L_{10}\right\rangle=-\left\langle\frac{H^{-\gamma}}{2}e^{it}J_{13,0}\right\rangle.$$

(57)

These can be inserted into the solvability conditions given by equation (51), which then read

$$\partial_{T}A =2i\frac{\partial}{\partial\bar{A}}\left\langle L_{10}\right\rangle\qquad \partial_{T}B =2i\frac{\partial}{\partial\bar{B}}\left\langle L_{10}\right\rangle.$$

(58) (59)

These may be identified as the Euler-Lagrange equations for the orbital period time-averaged Lagrangian at leading order

$${} \left\langle L\right\rangle =\frac{1}{2\pi}\int_{0}^{2\pi}L_{0}(A,B,\partial_{T}A,\partial_{T}B;c.c.)+\epsilon L_{1}(A,B;c.c.)\,dt\qquad {}=0+\epsilon\left[\frac{i}{4}(\bar{A}\partial_{T}A-A\partial_{T}\bar{A}+\bar{B}\partial_{T}B-B\partial_{T}\bar{B})+\left\langle L_{10}\right\rangle\right]+O(\epsilon^{2}),$$

(60)

where c.c. denotes the complex conjugated variables. The variational principle for minimising the action ∫⟨L⟩dT with respect to the generalised coordinates (A,B,∂_T A,∂_T B;c.c.), recovers equations (58) and (59). Alternatively we may identify ⟨H⟩=-⟨L₁₀⟩ as the Hamiltonian governing the secular evolution of the warp. Indeed, the Legendre transform of the averaged Lagrangian can be written as

$$\left\langle\mathcal{H}\right\rangle=\partial_{T}A\frac{\partial\langle L\rangle}{\partial(\partial_{T}A)}+\partial_{T}B\frac{\partial\langle L\rangle}{\partial(\partial_{T}B)}+c.c.-\langle L\rangle.$$

(61)

Since the Lagrangian is linear in the ‘velocity’ coordinates, all terms cancel apart from ⟨L₁₀⟩ which is reversed in sign. In this case, the modulation equations formally have the structure of the complex Hamilton’s equations

$$iz=\frac{\partial\langle H\rangle}{\partial\bar{z}},$$

(62)

where the canonical variable is identified as z=A/√(2) or B/√(2) for equations (58) and (59) respectively.

4.4 Forced vertical oscillator

In order to ground this formalism, it remains to determine the evolution of the forced vertical oscillator at leading order J₃₃,₀ so we can compute the averaged Lagrangian ⟨L₁₀⟩. Recall, the dynamics of the forced vertical oscillator is given by equation (46), where the forcing term on the right-hand side is given as a product of the tilt and shear

$$f_{ts}\equiv J_{13,0}J_{31,0}=\Re\left[Ae^{-it}\right]\Re\left[Be^{-it}\right].$$

(63)

Expanding this forcing yields

$$f_{ts}=\frac{1}{2}\Re\left\lbrace |A||B|e^{-i[2t-\arg(A)-\arg(B)]}+A\bar{B}\right\rbrace ,$$

(64)

where |·| and (·) denote the modulus and argument respectively. We are free to choose the time origin since equation (46) has no explicit temporal dependence. Taking t t-(B) gives the forcing form

$$f_{ts}=\frac{1}{2}\Re\left(Z_{1}e^{-2it}+Z_{1}\right),$$

(65)

where Z₁=AB̅. The net forcing on the right hand side then becomes

$$F\equiv-\partial_{tt}f_{ts}-f_{ts}=\Re\left(\frac{3}{2}Z_{1}e^{-2it}-\frac{1}{2}Z_{1}\right).$$

(66)

The solution for H is therefore only dependent on the value of Z₁ so the modulation equations become

$$\partial_{T}A=2i\frac{\partial\left\langle L_{10}\right\rangle}{\partial\bar{Z_{1}}}B\quad\text{and}\quad\partial_{T}B=2i\frac{\partial\left\langle L_{10}\right\rangle}{\partial Z_{1}}A,$$

(67)

where ⟨L₁₀⟩ is now regarded as a function of Z₁ and its complex conjugate. In order to find a solution for small Z₁, we first perform a weakly nonlinear analysis and find a series expansion solution for H. Taking Z₁=δ Z₁,₀, with δ≪1 and Z₁,₀ ∼ O(1), we expand the vertical oscillator equation in terms of

$$H=H_{0}(t)+\delta H_{1}(t)+\delta^{2}H_{2}(t)+\cdots,$$

(68)

which again generates a hierarchy of equations. At leading order we recover the unforced vertical oscillator

$$\partial_{tt}H_{0}+H_{0}-H_{0}^{-\gamma}=0.$$

(69)

In order to conform with the 2π periodic boundary conditions required by the solvability conditions discussed previously, we will set the free oscillation to zero and adopt the equilibrium value H₀=1. At the nᵗʰ order expansion we observe the general form

$$\left[\partial_{tt}+\omega_{n}(\gamma)^{2}\right]H_{n}(t)=F_{n}(\gamma,Z_{1},t,H_{i\lt n}),$$

(70)

where the forcing term on the right-hand side, Fₙ, depends on the lower order solutions. Again, we ignore the complementary solution so as to avoid quasiperiodic solutions and simply extract the forced oscillation at each order. This weakly nonlinear solution can be evaluated to arbitrary order and used to calculate the Lagrangian given by equation (52). In terms of the H variable this may be written as

$$L_{10}=\frac{1}{2}\left(\partial_{t}H+\partial_{t}f_{ts}\right)^{2}-\frac{1}{2}\left(H+f_{ts}\right)^{2}-\frac{H}{\gamma-1}\left(\partial_{tt}H+H-F\right).$$

(71)

Computing the average then gives

$${} \left\langle L_{10}\right\rangle= -\frac{1+\gamma}{2(\gamma-1)}-\frac{1}{4}(Z_{1}+\bar{Z_{1}})\qquad {}-\gamma\frac{\left[(\gamma-3)Z_{1}^{2}+(\gamma-3)\bar{Z_{1}}^{2}-4(\gamma+3)Z_{1}\bar{Z_{1}}\right]}{32(\gamma-3)(\gamma+1)}$$

(72)

to second order in Z₁. As a preliminary check on this result we can test the linear limit. The modulation equations (67) may be combined into the oscillator equation

$$\partial_{TT}A=-2i\left|\frac{\partial^{2}\langle L_{10}\rangle}{\partial T\partial Z_{1}}\right|B-4\left|\frac{\partial\left\langle L_{10}\right\rangle}{\partial Z_{1}}\right|^{2}A.$$

(73)

Inserting our weakly nonlinear averaged Lagrangian and retaining terms at linear order yields oscillatory solutions A∝exp(iωₚ T) with precessional frequency ωₚ=± 1/2. This matches onto the linear bending modes with tilting frequency ω=1±ε/2 in the local frame, as found previously in Paper I. Retaining higher order contributions allows us to extend this result for weakly nonlinear forcing warps. More generally, for larger amplitude oscillations we must solve for the forced vertical motions numerically. We will demonstrate this semi-analytical procedure in section 4.5.

4.5 Precessing warp solutions

With this semi-analytical modulation theory in hand, we will look for a pair of special solutions which correspond to the nonlinear extension of the normal bending modes. In line with the equipartition of tilt and shear found for linear bending waves, we restrict attention to complex amplitudes for which the magnitudes are equal and perfectly in-phase or anti-phased. In this case B=± A and thus Z₁=±|A|² is a real quantity, where the positive/negative sign describes in/anti-phase tilt and shear. Thus we can restrict attention to the averaged Lagrangian along the real line for which we denote Z₁=X. This may be computed numerically using a shooting code which converges to the 2π periodic solutions for H as shown in Fig. 4. Periodic solutions are found for all X<0, which correspond to anti-phased tilt and shear forcing. Meanwhile, the solution terminates in a saddle node bifurcation for sufficiently large X>0 (as previously noted by Ogilvie & Latter (2013) in the case γ=1), whereupon this theory breaks down. Observe that the weakly nonlinear solution for H, computed to second order in equation (4.4), is plotted as the red dashed line and provides a good fit for small X.

Figure 4: Upper panel: ⟨L₁₀⟩ is calculated by averaging equation (52) over the orbital timescale, having numerically identified the periodic solutions for H forced by real Z₁ according to equation (46). The gradient of this is then plotted in the lower panel and is related to the precessional frequency as per equation (77). The weakly-nonlinear expansion obtained up to O(Z₁²) in equation (4.4) is then over-plotted as dashed red lines, which give the leading order contribution of the nonlinearity. Note the solutions terminate at X∼ 0.4 at which point this smooth modulation theory will break down. Meanwhile the solutions may be continued indefinitely for X<0.

Making use of Wirtinger complex differentiation,

$$\frac{\partial}{\partial Z_{1}}=\frac{1}{2}\left(\frac{\partial}{\partial X}-i\frac{\partial}{\partial Y}\right)\quad\text{and}\quad\frac{\partial}{\partial\bar{Z_{1}}}=\frac{1}{2}\left(\frac{\partial}{\partial X}+i\frac{\partial}{\partial Y}\right),$$

(74)

and inserting B=± A, the amplitude modulation equations become

$$\partial_{T}A=\pm i\frac{\partial\left\langle L_{10}\right\rangle}{\partial X}A$$

(75)

when evaluated along the real Z₁ line. Note that Y derivatives disappear as ⟨L₁₀⟩ possesses reflectional symmetry about the X-axis. This must be the case since the forcing function fₜₛ obtained upon conjugating Z₁ is the identical up to a shift in phase. Therefore the periodic solutions for H and hence the averaged Lagrangian must be the same for Z_1Z_1. Examining the form of the modulation equation (75) we see it corresponds to a rotation of the complex amplitude whilst the magnitude remains constant. Both oscillator amplitudes rotate at equal rates (since they are described by identical equations) and so the forcing product Z₁ remains constant and on the real axis:

$$\frac{\partial Z_{1}}{\partial T}=\pm\frac{\partial|A|^{2}}{\partial T}=\pm\left(A\frac{\partial\bar{A}}{\partial T}+\bar{A}\frac{\partial A}{\partial T}\right)=0.$$

(76)

We seek oscillatory solutions of the form eⁱωₚ T such that

$$\omega_{p}=\pm\frac{\partial\left\langle L_{10}\right\rangle}{\partial X}$$

(77)

where the - solution corresponds to the anti-phased solutions with X<0 and the + solutions correspond to the in phase solutions with X>0. Reconstructing the tilting oscillator motion,

$$J_{31,0}=\sqrt{|X|}\Re[e^{-i(1-\epsilon\omega_{p})t}],$$

(78)

aids the interpretation of the result. The frequency observed in the local model is ω=1-εωₚ. Numerically we see from Fig. 4 that ∂_X ⟨L₁₀⟩<0 so for the in-phase motions ωₚ<0. Thus the local frequency is enhanced whilst the period

$$T=\frac{2\pi}{1-\epsilon\omega_{p}}$$

(79)

is reduced. This recovers the retrograde precession expected for the nonlinear extension of the in-phase bending modes. Similarly, for the anti-phase tilt and shear, ωₚ>0 and the oscillation period is less than the orbital period. This may be interpreted as prograde precession of the warped torus structure from a non-rotating global frame. We will return to these solutions in section 7 where we will verify this theory against the full equation set.

5 Bouncing regime

Upon reaching a critical warping amplitude, the smooth modulation theory of section 4 will break down. Indeed, Fig. 4 shows that the averaged Lagrangian solution terminates past a certain forcing amplitude. Furthermore, in section 3, we numerically identified a qualitatively distinct behaviour where the nonlinear vertical oscillator resonantly grows to large amplitudes and becomes extremely compressive. In this section we will develop a separate analytical theory for understanding this regime. We will begin by focusing our attention on the vertical oscillator forced by the warp, which is the defining feature of this bouncing regime, before incorporating the feedback self-consistently onto the tilt and shear.

5.1 Bouncing vertical oscillator

Initially we will ignore the dynamical evolution of the warp, as described by the tilt and shear equations for J₃₁ and J₁₃ respectively. Instead we treat the warp as being fixed and look for the response of the vertical oscillator as described by equation (32) for J₃₃. This approach is similar to that taken by Ogilvie & Latter (2013), to which we uncover a close mathematical correspondence. In order to draw a formal comparison with their analysis, we motivate a coordinate transformation which essentially subtracts the tilting motion and isolates the compressive behaviour. We take H=J/J₁₁, which can be interpreted as a measure of the disc thickness since J is proportional to the cross-sectional area and J₁₁ approximates the width of the ring. We also fix the the tilt and shear coordinates to oscillate harmonically with some arbitrary phase relationship. Thus, J_13=[Aexp(it)] and J_31=[Bexp(it)] where we redefine the complex amplitudes A=aexp(iθ_A) and B=bexp(iθ_B). We will assume that the radial extent of the ring, described by J₁₁, is held constant. Indeed, we see in Figs. 2 and 3 that the J₁₁ oscillator evolves independently from the mode coupling phenomenon so we will ignore equation (29). Inserting these transformations into equation (32) yields

$$\ddot{H}+H-\frac{\hat{T}_{0}}{J_{11}^{\gamma-1}}\frac{1}{H^{\gamma}}=\frac{1}{J_{11}}\Re\left(\frac{3}{2}Z_{2}e^{2it}-\frac{1}{2}Z_{2}\right),$$

(80)

where the forcing product Z₂ is defined as AB̅. Here we are allowing for an arbitrary choice of J₁₁ as opposed to the convenient scaling chosen previously in equation (38). This re-scaled definition is simply related to that introduced in our modulation theory by a multiplicative factor,

$$Z_{2}=J_{11}^{2}\epsilon Z_{1}.$$

(81)

Clearly when J₁₁=ε⁻¹⁄² as before, we recover the equality between the two definitions. The left-hand side of equation (80) represents a free oscillator where the harmonic trajectory is interrupted by the pressure based anharmonic restoring force as the ring is compressed. The right-hand side is a forcing term with a strength proportional to the product of the shear and tilt magnitudes Z₂. As we have already seen in section 4, equation (80) once again emphasises the generic effect of warped geometries forcing vertical motions.

5.1.1 Free vertical oscillation period

We will now concentrate on the properties of the free non-linear vertical oscillator. To this end we set Z₂=0 and work with

$$\ddot{H}+H-\frac{\hat{T}_{0}}{J_{11}^{\gamma-1}}\frac{1}{H^{\gamma}}=0.$$

(82)

This can be derived from a conserved energy Hamiltonian H_f composed of the sum of kinetic, potential and internal energies respectively

$$\mathcal{H}_{f}=\frac{1}{2}\dot{H}^{2}+\frac{1}{2}H^{2}+\frac{\hat{T}_{0}}{J_{11}^{\gamma-1}}\frac{H^{1-\gamma}}{\gamma-1}.$$

(83)

Equation (82) clearly permits an equilibrium at H=(T̂₀/J₁₁ᵞ⁻¹)ᵞ⁺¹ and a linear perturbation then yields a natural oscillation frequency of √(γ+1), as expected from our previous analysis of breathing modes in Paper I. As the amplitude increases into the non-linear regime, we can qualitatively see that the frequency tends monotonically towards 2. Indeed, in this case H behaves predominantly as a harmonic oscillator in a quadratic potential with an impulsive pressure reversal acting when H→ 0 which rectifies the motion. As the amplitude becomes ever larger, the harmonic motion dominates for the majority of the trajectory. Thus we expect the period of the free oscillator to be T∼π to leading order with a small correction due to the phase shift incurred by pressure. Using the Hamiltonian energy function we can construct an integral for the period as follows:

$$T=2\int_{H_{min}}^{H_{max}}\frac{dH}{\sqrt{2\mathcal{H}_{f}-\frac{2\hat{T}_{0}H^{1-\gamma}}{J_{11}^{\gamma-1}(\gamma-1)}-H^{2}}},$$

(84)

where

$$H_{min}\sim\left[\frac{\hat{T}_{0}}{(\gamma-1)J_{11}^{\gamma-1}\mathcal{H}_{f}}\right]^{\frac{1}{\gamma-1}}\quad\text{and}\quad H_{max}\sim\sqrt{2\mathcal{H}_{f}}$$

(85)

denote the minimum and maximum turning points of the vertical oscillator to leading order. Unfortunately this integral cannot be analytically evaluated except in the special case γ=3, for which we find a period of exactly π. For other values of γ we turn to a range splitting technique which allows us to construct an asymptotic expression in the limit of large amplitude oscillations. This involves approximating the integrand in three distinct intervals and then matching them together such that the errors are subdominant. The leading order deviation from period π is found to be

$$T-\pi= \varpi(\gamma,H_{max})=c(\gamma)H_{max}^{-1-\gamma},\gamma\lt 2,\quad \eta(\gamma,H_{min})=d(\gamma)H_{min}^{(1+\gamma)/2},\gamma\gt 2,$$

(86)

with γ dependent coefficients

$$c(\gamma)=\frac{4\sqrt{\pi}\hat{T}_{0}}{(\gamma-2)(\gamma-1)J_{11}^{\gamma-1}}\frac{\Gamma(2-\gamma/2)}{\Gamma(1/2-\gamma/2)},\qquad d(\gamma)=-\sqrt{\frac{2\pi(\gamma-1)J_{11}^{\gamma-1}}{\hat{T}_{0}}}\frac{\Gamma(1+1/(1-\gamma))}{\Gamma(1/2+1/(1-\gamma))}.$$

(87) (88)

where Γ denote gamma functions. The key point here is that the phase delay has two separate asymptotic limits set by the value of γ. When γ<2 the ring is more compressible and the period offset is attributed to the cumulative extended effects of pressure over the trajectory. Meanwhile when γ>2 the ring is less compressible and the pressure effects are localised near the minimum turning point. For γ<3 both c(γ) and d(γ) are greater than 0 so the period is slightly greater than π. For the special integrable case γ=3, d(3)=0 as expected and the period is exactly π. For γ greater than this, d(γ) is negative and the period is slightly less than π. In the upcoming sections we will assume a typical γ<2 and hence adopt a period offset $\varpi$ from bounce to bounce.

5.1.2 Impulsively forced vertical oscillation

We now extend this analysis to the case where this large amplitude bouncing mode is forced by a fixed warp, with non-zero Z₂, as described by equation (80). In fact, we can recast this equation using an intuitive coordinate transformation which reinterprets this forcing term as a localised bouncing off an oscillating boundary. Indeed, if we write

h(t)≡H(t)+f(t)=J₃₃,

(89)

where f(t)=J₁₃ J₃₁/J₁₁, and substitute into (80) we recover

$$\ddot{h}+h=\frac{\hat{T}_{0}}{J_{11}^{\gamma-1}}H^{-\gamma}.$$

(90)

This is simply equation (32) in disguise, which may seem a rather circular procedure. However the purpose of introducing this coordinate transformation lies in the helpful physical reinterpretation of the problem. We can view h(t) as the extension of a mass on a spring from an equilibrium position. This wants to undergo harmonic motion according to Hooke’s law until the motion is interrupted by an oscillating wall at position f(t). Thus H(t) is the distance between the mass and the wall as visualised in Fig. 5.

Figure 5: The forced vertical oscillator described by equation (90) may be reinterpreted as a harmonically oscillating mass bouncing off a moving wall. The equilibrium position of the spring-mass system is denoted by the dashed line x=0. The wall oscillates about x=0 according to f(t) whilst the mass is a distance H(t) from this wall. The spring then has a total extension from its equilibrium h(t)=f(t)+H(t) which incurs a restoring force according to Hooke’s law.

When H(t) tends to zero from above the relative velocity between the mass and the wall will reverse in an ideal, elastic bounce. This is very similar to the problem investigated by Holmes (1982) and Luo & Han (1996) with regards to a ball bouncing off an oscillating table, where the motion is reduced to a discrete mapping from bounce to bounce. We proceed similarly by neglecting the pressure contribution from H⁻ᵞ in between bounces. Instead, we assume that it acts impulsively to reverse the direction of motion upon each elastic collision with the wall. Furthermore, it incurs a small phase delay (γ,H_max), in accordance with our asymptotic investigation of the free vertical oscillator as presented in equation (86). These assumptions are valid provided the amplitude and velocity of the mass motion are much larger than the wall position f and velocity ḟ at the time of impact. In this case, the bouncing period only slightly departs from π and thus the phase relationship with respect to the oscillating wall evolves slowly. Assume that the nᵗʰ bounce occurs at tₙ=(n-1)π+ϕₙ, where ϕₙ is a phase offset which evolves slowly. Just after the bounce we have

$$H_{n}=0,\quad\dot{H}_{n}\equiv v_{n}\gt 0,$$

(91)

and the wall has position and velocity given by

$$f_{n}=\frac{ab}{2J_{11}}[\cos(\theta_{A}+\theta_{B}+2\phi_{n})+\cos(\theta_{A}-\theta_{B})],\qquad \dot{f}_{n}=-\frac{ab}{J_{11}}\sin(\theta_{A}+\theta_{B}+2\phi_{n}).$$

(92) (93)

Thus the position and velocity of the mass are

$$h_{n}=f_{n},\quad\dot{h}_{n}=\dot{H}_{n}+\dot{f}_{n}=v_{n}+\dot{f}_{n}.$$

(94)

Between bounces we assume purely harmonic motion governed by ḧ+h=0. The initial conditions (94) then determine the trajectory

$$h=f_{n}\cos(t-t_{n})+(v_{n}+\dot{f}_{n})\sin{(t-t_{n})}.$$

(95)

However, we wish to capture the retarding effect of pressure so we incorporate the phase offset $\varpi$ taken from our asymptotic analysis of the free non-linear vertical oscillator as follows:

$$h =f_{n}\cos(t-t_{n}-\varpi_{n})+(v_{n}+\dot{f}_{n})\sin(t-t_{n}-\varpi_{n}),\qquad \dot{h} =-f_{n}\sin(t-t_{n}-\varpi_{n})+(v_{n}+\dot{f}_{n})\cos(t-t_{n}-\varpi_{n}).$$

(96) (97)

The next bounce occurs at tₙ₊₁=nπ+ϕₙ₊₁. Substituting this into the above expressions allows us to relate successive bounces as

$$f_{n+1}= -f_{n}\cos(\phi_{n+1}-\phi_{n}-\varpi_{n})\qquad {}-(v_{n}+\dot{f}_{n})\sin(\phi_{n+1}-\phi_{n}-\varpi_{n}),\qquad v_{n+1}= -f_{n}\sin(\phi_{n+1}-\phi_{n}-\varpi_{n})\qquad {}+(v_{n}+\dot{f}_{n})\cos(\phi_{n+1}-\phi_{n}-\varpi_{n})+\dot{f}_{n+1}.$$

(98) (99)

In the large amplitude limit with vₙ≫|fₙ|,|ḟₙ|, the terms in the equation (5.1.2) can only be consistently balanced provided |ϕ_n+1-ϕ_n-_n| 1. Since $\varpi$ is a small phase correction, this in turn ensures |ϕₙ₊₁-ϕₙ|≪1. This agrees with our expectation that the phase evolves slowly from bounce to bounce. With this assumption, these equations can be simplified to leading order giving the recursive update scheme

$$\phi_{n+1}-\phi_{n} =\varpi_{n}-\frac{2f_{n}}{v_{n}},\qquad v_{n+1}-v_{n} =2\dot{f}_{n}.$$

(100) (101)

The ϕ update is composed of two parts – the $\varpi$ contribution from pressure and also the effect of the oscillating impact position. Recognising that the amplitude of the oscillating mass is approximately equal to the impact velocity with the wall, we can then express the phase delay as _n(v_n)=c(γ)v_n^-(γ+1) in accordance with (86).

5.1.3 Hamiltonian structure

Considering the variable updates from bounce to bounce are small, we may take the continuous ODE analogue of these discrete mappings to be

$$\frac{d\phi}{dn}=\varpi-\frac{2f}{v},\qquad \frac{dv}{dn}=2\dot{f}.$$

(102) (103)

As we might anticipate for an ideal system, these equations possess an autonomous symplectic structure. This is best seen by the change of variables I=1/2v², which is the classical action of a harmonic oscillator. The Hamiltonian is then found to be

$$\mathcal{H}=\frac{c(\gamma)}{(1-\gamma)}(2I)^{(1-\gamma)/2}-2(2I)^{1/2}f(\phi),$$

(104)

with the canonical equations of motion

$$\frac{d\phi}{dn}=\frac{\partial\mathcal{H}}{\partial I},\quad\frac{dI}{dn}=-\frac{\partial\mathcal{H}}{\partial\phi}.$$

(105)

Note that when the warped forcing is absent, the Hamiltonian is independent of the phase angle and hence the action is invariant whilst the phase advances uniformly. This simply corresponds to the free harmonic oscillator with constant amplitude and phase delay from bounce to bounce. More generally for non-zero forcing, the contours of the Hamiltonian trace out the trajectories in phase space. An example of this structure is shown in Fig. 6 for the particular choice θ_A=θ_B=0, which corresponds to the tilt and shear oscillators being in-phase. Here we set the value of c(γ;J₁₁,T̂₀) for γ=5/3, which also depends on the scaling parameters chosen for the ellipse. As per our numerical runs in section 3 we choose J₁₁=100 and the value of T̂₀ so the associated equilibrium ring has aspect ratio ε=0.01. Note that the structure is π periodic since the phase variable ϕ is measured modulo the rectified harmonic period of π. The red dashed line plots the hetero-clinic separatrix structure emanating from the unstable saddle point located at

{0,[c(γ)J₁₁/(2ab)]¹⁄ᵞ} .

(106)

This delimits a circulating solution from a resonantly growing solution which becomes phase locked as the bounce amplitude tends to infinity.

Figure 6: ϕ-v phase plane portrait for the specific choice θ_A=θ_B=0, γ=5/3, J₁₁=100 and T̂₀=100ᵞ⁻¹. The coloured solid lines denote contour levels of the Hamiltonian H, along which the system evolves from bounce to bounce. The red dashed lines denote the separatrix curves originating from the critical unstable saddle points and delimit the resonant, phase locking behaviour from the periodic behaviour.

This resonant phase locking is observed in numerical solutions of equation (80) and helps elucidate the physical mechanism responsible for the growth of compressive vertical motions. When the phase delay incurred by the pressure retardation is sufficiently counteracted by the changing phase relationship with the wall, energy is constructively input into the breathing mode over many cycles. This increases its amplitude and reduces the rate of future phase evolution, further locking it into a resonant relationship. This distinct behaviour for sufficiently large warps hints towards the existence of a critical warping amplitude above which the oscillator will be driven into the bouncing regime as we found in our numerical experiments in section 3.

5.1.4 Periodic structure using shooting method

Figure 7: The thick line tracks periodic solutions of equation (80) with the phase variables set to be θ_A=θ_B and for γ=5/3, J₁₁=100 and T̂₀=100ᵞ⁻¹. The x-axis plots the value of ab which is a measure of the tilt and shear forcing amplitude. The y-axis plots the maximum value v=Ḣ for each periodic solution. The branch colour denotes the maximum magnitude of the eigenvalues obtained from the monodromy matrix associated with a Floquet analysis about each periodic solution. Values above 1 indicate instability. The dashed red line denotes the analytical prediction for the location of the saddle point in the large bounce regime according to equation (106).

The analytical progress made in the previous section accurately describes the forced vertical oscillator in the extreme bouncing regime. However, in the case of low amplitude oscillations, not far from the equilibrium of the disc, we might expect our approximations to break down. Indeed, previous work by Ogilvie & Latter (2013) found stable periodic solutions for the simple laminar flows in a warped disc, provided the enforced warp amplitude is sufficiently low. In order match the high and low amplitude regimes, we use a shooting scheme to identify the existence of periodic solutions as the imposed warp is varied through the tilt and shear product ab. The shooting code implemented solves equation (80), with γ=5/3, J₁₁=100 and T̂₀ such that the equilibrium ring has aspect ratio ε=0.01. We use a typical Runge-Kutta integrator with adaptive step-size and then minimise residuals at the boundary in accordance with Levenberg–Marquardt least squares optimisation (Dednam & Botha, 2014).

This efficiently converges onto 2π periodic solutions which are plotted in Fig. 7. The tilt and shear oscillators are set with the phase relationship θ_A=θ_B. The y-axis denotes the maximum value of v=Ḣ which is the appropriate amplitude measure for the periodic solutions. Meanwhile, the x-axis describes the forcing product ab. When ab>0 the tilt and shear are in phase, whilst when ab<0 the tilt and shear are in anti-phase. For each periodic solution we also perform a Floquet stability analysis. We calculate the monodromy matrix and extract the eigenvalues with maximum magnitude. Since we are expanding about a periodic solution, there always exists an eigenvalue equal to 1 which corresponds to a perturbation tangential to the periodic solution. If there exists an eigenvalue with absolute magnitude greater than 1 (i.e. outwith the complex unit circle), the periodic solution is unstable. The value of the periodic solutions are coloured according to the maximum magnitude eigenvalue, with purple denoting the stable solution baseline with eigenvalue equal to 1. Finally, the saddle point location, predicted by the high amplitude bouncing theory, is plotted as the red dashed line for comparison.

In the high amplitude limit, we do indeed converge to the unstable saddle-point solutions with an in-phase forcing ab>0. For large values of v the analytically predicted solution tends asymptotically towards our numerical findings. As we move down this branch towards the kink, the numerical shooting code deviates from our prediction as the impulsive approximation breaks down. The saddle point exhibits a peak unstable growth rate before the kink turns over and enters the stable lower branch. This can be continued indefinitely towards large negative values of ab. When ab becomes less than 0, this is equivalent to the tilt and shear becoming π out of phase. This in turn changes the phase relationship with the driven vertical oscillator so the velocity amplitude becomes negative. Whilst the upper branch describes a saddle point, the numerically identified lower branch represents a stable centre. As the forcing warp is increased towards the turning point, these two points converge and eventually collide in a saddle node bifurcation at ab∼ 40. This behaviour is consistent with the termination of solutions found previously for the forced vertical oscillator within the context of our modulation theory in section 4.4. In Fig. 4 the solution branch ends abruptly at Z₁=0.4. When this is re-scaled by J₁₁² ε=100 to account for our arbitrary choice of ring width we find agreement with Z₂=ab=40. This bifurcation point sets a critical warping amplitude beyond which no periodic solutions can be found for in phase tilt and shear. Instead, trajectories are carried up along the steep contours as shown in Fig. 6 tending asymptotically to the fixed resonant phase relationship ϕ=±π/2.

These results agree with the findings of Ogilvie & Latter (2013). They also find that for a sufficiently large positive warped forcing, the π periodic solutions terminate. This offers a mechanism for which a system with no initial vertical motion can be driven to large amplitudes, provided the warp amplitude lies beyond this critical turning point. This is what we see in Fig. 2 where, for moderate warp, the the vertical motion becomes highly activated. As growth continues, the feedback of the vertical motion onto the warp will become important and the enforced warp assumption will also break down. We will address this via a self-consistent coupling of the warp to the vertical bouncing in the next section.

5.2 Feedback onto the warp

The previous analysis assumes that the warp is fixed, with the tilt and shear oscillating sinusoidally at the orbital frequency. We found that this leads to resonant growth if the phase becomes locked and energy continues to be injected into the bouncing motions. In reality, total energy is conserved and energy flowing into one mode must be coupled with energy leaving another, as seen in the motivating plots of Fig. 3. Indeed, we must consider the back-reaction onto the warp which is then allowed to evolve.

Let us consider the case that the breathing mode has entered into the highly compressive non-linear regime. We have seen that the effect of pressure can be treated as an impulsive forcing which reverses the direction of the bouncing mass. This also gives us reason to believe that the pressure terms in the equations (30) and (31) also enter as time localised impulsive forces. Indeed, in this regime we expect the tilt and shear oscillators to undergo linear harmonic motion which is periodically kicked, causing an instantaneous change in their amplitude and phase from bounce to bounce. By connecting the piece-wise harmonic intervals between the nᵗʰ and (n+1)ᵗʰ bounces, we can create an iterable mapping for the evolution of the system. We see from equation (90) that the impulsive forcing on the right hand side T̂₀/J₁₁ᵞ⁻¹ Hᵞ provides this Dirac delta forcing. It reverses the impact velocity vₙ₊₁ of the mass relative to the wall such that

$$\frac{\hat{T}_{0}}{J_{11}^{\gamma-1}}\frac{1}{H^{\gamma}}=\frac{\hat{T}_{0}J_{11}}{J^{\gamma}}=2v_{n+1}\delta(t-t_{n+1}),$$

(107)

where δ(t) is the Dirac delta function. Thus our tilt and shear oscillator equations have the form

$$\ddot{J}_{13}+J_{13} =-\frac{2v_{n+1}J_{31}}{J_{11}}\delta(t-t_{n+1}),\qquad \ddot{J}_{31}+J_{31} =-\frac{2v_{n+1}J_{13}}{J_{11}}\delta(t-t_{n+1}).$$

(108) (109)

where we have neglected (the often small constant) Cₓ. These have the general form of harmonic oscillators undergoing impulsive kicks as described by the equation

$$\ddot{x}(t)+\omega^{2}x(t)=F_{0}\delta(t-t_{n+1}),$$

(110)

where F₀ is the momentum impulse, such that integration over the equation gives an instantaneous change in velocity Δẋ=F₀. This problem is completed by furnishing it with the initial conditions x(0)=x₀ and ẋ(0)=v₀. This equation has wide reaching physical applications and has been studied extensively with application to both classical and quantum problems. The solution is easily found by converting it to an algebraic equation via the Laplace transform and then inverting back to the original variable domain. We find the solution to be

$$x(t)=x_{0}\cos(\omega t)+\frac{v_{0}}{\omega}\sin(\omega t)+\frac{F_{0}}{\omega}u(t-t_{n+1})\sin[\omega(t-t_{n+1})],$$

(111)

where u is the unit-step function. Clearly the amplitude and phase of the oscillator are modified after the impact. We will find it convenient to describe this in terms of a complex amplitude χ=χ_R+iχ_I=|χ|exp(iθ_χ) such that

$$x(t) =\Re[\chi e^{i\omega t}]\qquad {}=\frac{1}{2}\left(\chi e^{i\omega t}+\bar{\chi}e^{-i\omega t}\right)\qquad {}=\chi_{R}\cos(\omega t)-\chi_{I}\sin(\omega t)\qquad {}=|\chi|\cos(\omega t+\theta_{\chi})$$

(112) (113) (114) (115)

By comparing the sine and cosine coefficients before and after the bounce we find the complex amplitude mapping

$$\chi_{R,n+1} =\chi_{R,n}-\frac{F_{0}}{\omega}\sin(\omega t_{n+1}),\qquad \chi_{I,n+1} =\chi_{I,n}-\frac{F_{0}}{\omega}\cos(\omega t_{n+1}),$$

(116) (117)

which is equivalent to

$$\chi_{n+1}=\chi_{n}+\frac{F_{0}}{\omega}e^{-i(\omega t_{n+1}+\pi/2)}.$$

(118)

We now apply this method to the equations for J₁₃ and J₃₁. Inserting the relevant Dirac-delta forcing coefficients leads to the iterative scheme

$$A_{n+1} =A_{n}-\frac{2v_{n}}{J_{11}}\Re(B_{n}e^{it_{n+1}})e^{-i(t_{n+1}+\pi/2)},\qquad B_{n+1} =B_{n}-\frac{2v_{n}}{J_{11}}\Re(A_{n}e^{it_{n+1}})e^{-i(t_{n+1}+\pi/2)}.$$

(119) (120)

The forcing coefficient is proportional to vₙ which is set by the bouncing vertical oscillator. Thus in order to close this discrete set we must couple it to the mappings of vₙ and ϕₙ derived previously in equations (100) and (101). As the tilt and shear evolve from bounce to bounce, this in turn modifies the forcing on the breathing motions in accordance with

$$f(t_{n}) =\frac{a_{n}b_{n}}{J_{11}}\cos(t_{n}+\theta_{A})\cos(t_{n}+\theta_{B}).$$

(121)

6 Resonant centres for coupled system

Having developed an impulsive theory for the bouncing regime, we will now seek special resonant solutions for which the amplitude of the oscillators is constant and the phase relationship between them remains fixed. Physically, these describe large amplitude, globally precessing warped solutions with extreme compressions twice per orbit.

6.1 Hamiltonian structure and resonant centres

We have derived a self-consistent set of discrete equations mapping the non-linear breathing and warping motions at each compression of the ring. This is still highly coupled and requires some further simplification to gain more dynamical insight. We again take the continuous limit (as done previously for the non-linear vertical oscillator) and study the resulting system of ODEs. By writing A=aeⁱθ_A and B=beⁱθ_B we can split up the real and imaginary parts of (119) and (120) to generate evolutionary equations for the amplitudes and phases:

(da)/(dn)=(vb)/(J₁₁)[ sin(θ_A-θ_B)+ sin(θ_A+θ_B+2ϕ)], (db)/(dn)=(va)/(J₁₁)[- sin(θ_A-θ_B)+ sin(θ_A+θ_B+2ϕ)], a(dθ_A)/(dn)=(vb)/(J₁₁)[cos(θ_A-θ_B)+cos(θ_A+θ_B+2ϕ)], b(dθ_B)/(dn)=(va)/(J₁₁)[cos(θ_A-θ_B)+cos(θ_A+θ_B+2ϕ)].

(122) (123) (124) (125)

Since the motions of all three oscillators are essentially harmonic between bounces, we anticipate that simple action-angle coordinates will further elucidate the structure of our equation set. We naturally adopt θ_A, θ_B and θ_C≡-ϕ as our angles whilst the energy of each oscillator gives the actions I_A=a²/2, I_B=b²/2 and I_C=v²/2. Note that taking the negative of ϕ makes sense as an angle since an increase in ϕ represents a delay to the bounce time. This corresponds to a negative shift in the phase angle of a rectified harmonic oscillator. Using these transformations leads to a set of 6 equations:

(dθ_A)/(dn)=(√(2))/(J₁₁)√((I_B I_C)/(I_A))[cos(θ_A-θ_B)+cos(θ_A+θ_B-2θ_C)], (dI_A)/(dn)=(2√(2))/(J₁₁)√(I_A I_B I_C)[ sin(θ_A-θ_B)+ sin(θ_A+θ_B-2θ_C)], (dθ_B)/(dn)=(√(2))/(J₁₁)√((I_A I_C)/(I_B))[cos(θ_A-θ_B)+cos(θ_A+θ_B-2θ_C)], (dI_B)/(dn)=(2√(2))/(J₁₁)√(I_A I_B I_C)[- sin(θ_A-θ_B)+ sin(θ_A+θ_B-2θ_C)],

$$\frac{d\theta_{C}}{dn}=\frac{\sqrt{2}}{J_{11}}\sqrt{\frac{I_{A}I_{B}}{I_{C}}}\left[\cos(\theta_{A}-\theta_{B})+\cos(\theta_{A}+\theta_{B}-2\theta_{C})\right]-\varpi,\qquad \frac{dI_{C}}{dn}=-\frac{4\sqrt{2}}{J_{11}}\sqrt{I_{A}I_{B}I_{C}}\sin(\theta_{A}+\theta_{B}-2\theta_{C}).$$

(126) (127) (128) (129) (130) (131)

These possess a symplectic structure amenable to a Hamiltonian formalism. The appropriate Hamiltonian is found to be

$$\mathcal{H}=\frac{2\sqrt{2}}{J_{11}}\sqrt{I_{A}I_{B}I_{C}}[\cos(\theta_{A}-\theta_{B})+\cos(\theta_{A}+\theta_{B}-2\theta_{C})]+\Delta(I_{C})$$

(132)

where Δ(I_C) is defined such that dΔ/dI_C=-(γ,I_C) and characterises the retarding phase offset from the vertical oscillator. Hamilton’s equations are then given by

$$\frac{d\theta_{i}}{dn}=\frac{\partial H}{\partial I_{i}},\quad\frac{dI_{i}}{dn}=-\frac{\partial H}{\partial\theta_{i}}.$$

(133)

This nicely extends the Hamiltonian structure for the forced vertical oscillator found previously in section 5.1.3, which is recovered by fixing the action and angle variables corresponding to J₁₃ and J₃₁. Notice also the inherent symmetry in the Hamiltonian upon a constant translation in the angles θᵢ →θᵢ+δθ. This canonical transformation is facilitated by the arbitrariness of setting the phase origin and, via Noether’s theorem, is generated by the conserved total action Iₜ=I_A+I_B+I_C. This is reminiscent of the energy conservation as seen in Fig. 3 where now the action variables are a proxy for the energy contained in the different modes. The dynamical evolution in phase space allows for an action interchange between modes but constrains trajectories to lie on contours of conserved total action.

Phase locking now occurs if the resonant angle combinations ξ₁≡θ_A-θ_B and ξ₂≡θ_A+θ_B-2θ_C are librating around a fixed centre rather than circulating. These fixed points are located by solving ξ̇₁=ξ̇₂=İᵢ=0. Resonance requires that ξ₁=lπ and ξ₂=mπ where l,m∈ℤ. The freedom to choose the phase relationship permits two separate resonant centres; an upper and a lower branch corresponding to the plus and minus sign respectively in cosξ₁+cosξ₂=± 2. The upper branch requires θ_A=θ_B so the tilt and shear are in phase. Meanwhile the lower branch requires that θ_A=θ_B+π so the tilt and shear are exactly out of phase. θ_C is either 0 or π such that the bounce point of the compressive breathing mode coincides with the points of maximal tilt and shear. This allows us to solve for the action centres

$$I_{0}(I_{C})\equiv I_{A}=I_{B}=I_{C}\pm\frac{\varpi(I_{C})\sqrt{I_{C}}J_{11}}{2\sqrt{2}},$$

(134)

where the common function I₀ equals the equipartition between I_A and I_B. Again here, the plus family of solutions correspond to the in-phase tilt and shear, whilst the minus branch correspond to the out of phase solution. We see a continuous family of resonant centres parameterised by the breathing action I_C.

Figure 8: A plot of the upper (red line) and lower (blue line) resonant branches corresponding to the centres identified in equation (134). The bifurcation points in the dynamical behaviour are noted with black dots and the asymptotic limit of I₀=I_C is plotted as a dashed black line.

6.2 Exploring the branch structure

This branch structure is plotted in Fig. 8 for γ=5/3, J₁₁=100 and T̂₀=J₁₁ᵞ⁻¹. The x-axis corresponds to the purely vertical breathing modes with only the vertical I_C action excited. Numerical evaluation of Hamilton’s equations (126) – (131) suggests that this is a stable periodic solution until reaching a pitchfork bifurcation at the point

$$I_{3,bif}=\left[\frac{c(\gamma)J_{11}}{4}2^{-\gamma/2}\right]^{\frac{2}{\gamma+2}},$$

(135)

where the lower blue branch intercepts the I_C axis. At this point, the breathing mode becomes unstable to warping motions and a stable mixed-mode branch is spawned. The tilt and shear are out of phase and their action is slightly less than that of the vertical oscillator as indicated by the dashed line in Fig. 8. Formally this bifurcation point arises as a parametric instability of the tilt and shear oscillators which we will now demonstrate. Akin to the classic example of a swing being pumped at twice its natural frequency, the breathing mode pumps the tilt and shear as we increase its amplitude and the frequency becomes sufficiently close to 2. Diagonalising equations (30) and (31) we have

$$\ddot{q}_{1}+\left(1+\frac{\hat{T}_{0}}{J_{11}^{\gamma}H^{\gamma}}\right)q_{1}=0\quad\text{and}\quad\ddot{q}_{2}+\left(1-\frac{\hat{T}_{0}}{J_{11}^{\gamma}H^{\gamma}}\right)q_{2}=0,$$

(136)

where we have used the change of basis

$$\begin{pmatrix}J_{13}\\ J_{31}\end{pmatrix}=\begin{pmatrix}1&1\\ 1&-1\end{pmatrix}\begin{pmatrix}q_{1}\\ q_{2}\end{pmatrix}.$$

| | (137) |

Thus the q₁ component corresponds to the in-phase tilt and shear contribution and q₂ the anti-phase component. Since it is the lower, anti-phased branch which is spawned from the vertical breathing mode x-axis in Fig. 8, we proceed to look for parametric instability in the q₂ equation. As before, we treat the forcing by the large amplitude breathing mode impulsively so T̂₀/(J₁₁ H)ᵞ=2vᵢₘₚ δ(t-tᵢₘₚ)/J₁₁ where vᵢₘₚ and tᵢₘₚ denote the impact velocity and time respectively. Then we can write the impulsive system of differential equations as

$$\dot{\begin{pmatrix}q\\ p\end{pmatrix}}=\begin{pmatrix}0&1\\ -1&0\end{pmatrix}\begin{pmatrix}q\\ p\end{pmatrix}\quad\text{for}\quad t\neq t_{imp},$$

| | (138) |

$$\begin{pmatrix}\Delta q\\ \Delta p\end{pmatrix}=\begin{pmatrix}0&0\\ \frac{2v_{imp}}{J_{11}}&0\end{pmatrix}\begin{pmatrix}q\\ p\end{pmatrix}\quad\text{for}\quad t=t_{imp},$$

| | (139) |

where q=q₂ and p=q_2. Δ p and Δ q denote the discrete jump in the quantities at the impact times of the breathing mode tᵢₘₚ. This form is now amenable to the impulsive Floquet theory developed by Bainov & Simeonov (1993). Similar to the usual continuous Floquet analysis, we construct the monodromy matrix M which captures the evolution of the system over one bounce period T,

$$M=\begin{pmatrix}\cos T&\sin T\\ \frac{2v_{imp}}{J_{11}}\cos T-\sin T&\frac{2v_{imp}}{J_{11}}\sin T+\cos T\end{pmatrix}.$$

(140)

The eigenvalues of M correspond to Floquet multipliers μᵢ which determine the stability of the trivial solution (q,p)=(0,0). Stability requires that |μᵢ|≤1 for all i. It should also be noted from the determinant of M that μ₁ μ₂=1. The characteristic equation for M gives

$$\mu_{1,2}=\cos{T}+\frac{v_{imp}}{J_{11}}\sin{T}\pm\sqrt{(\cos{T}+\frac{v_{imp}}{J_{11}}\sin{T})^{2}-1}.$$

(141)

If |cosT+(vᵢₘₚ)/(J₁₁) sinT|<1 then the Floquet multipliers are complex conjugates. Since μ₁ μ₂=1 this requires |μᵢ|=1 and the trivial solution is stable. In contrast, if |cosT+(vᵢₘₚ)/(J₁₁) sinT|>1 then the multipliers are real and distinct. Therefore one must be greater than 1 and the trivial solution is unstable to parametric growth. Let us insert the bouncing period T=π+ into this criterion such that

$$\cos{\varpi}+\frac{v_{imp}}{J_{11}}\sin{\varpi}\gt 1.$$

(142)

Expanding terms to O (^2) in the small phase delay allows us to deduce the instability criterion v_imp/J_11-/2> 0. As before, the phase delay can be written asymptotically in accordance with equation (86) as =c(γ)v_imp^-1-γ, so rearrangement yields the critical value

$$v_{imp,bif}=\left[\frac{c(\gamma)J_{11}}{2}\right]^{\frac{1}{2+\gamma}}.$$

(143)

This agrees exactly with the bifurcation point identified as the x-intercept of the lower resonant branch in equation (135) and demonstrates the underlying parametric mechanism.

We are also able to deduce the stability of the non-trivial resonant branches themselves. Consider the equations for ξ₁, ξ₂, I_A, I_B and I_C. Linearising about the fixed resonant solutions, parameterised by I_C, yields a 5× 5 Jacobian matrix which encapsulates the stability as we move along the branches. The solution is stable provided no eigenvalues have a real component. We find that the lower branch is stable for all I_C above the parametric bifurcation point from whence it originates. Meanwhile the upper branch shows a transition from an unstable to a stable region. In Fig. 9 we plot the maximum real part of the eigenvalues λᵢ, corresponding to the Jacobian computed about the upper branch. We see that the branch is unstable when I_C<9.05 which corresponds to the region left of the black dot as plotted on the upper branch of Fig. 8.

Figure 9: The maximum real eigenvalue for the Jacobian of our action-angle equations evaluated about the upper resonant branch for the case J₁₁=100 and γ=5/3. This is parameterised by the vertical mode action I_C. We see instability for I_C<9.05 and stability for I_C>9.05.

Moreover we find that the eigenvectors associated with the unstable growth correspond to equal in-phase perturbations of tilt and shear. i.e. those which maintain I_A=I_B and θ_A=θ_B. We will make use of this fact and restrict our attention to the equal amplitude in-phase tilt and shear. This reduces our system of equations to

$$\frac{d\theta_{A}}{dn} =\frac{\sqrt{2}}{J_{11}}\sqrt{I_{C}}\left[1+\cos{(2\theta_{A}-2\theta_{C})}\right],\qquad \frac{dI_{A+B}}{dn} =\frac{2\sqrt{2}}{J_{11}}I_{A+B}\sqrt{I_{C}}\sin(2\theta_{A}-2\theta_{C}),\qquad \frac{d\theta_{C}}{dn} =\frac{\sqrt{2}}{2J_{11}}\frac{I_{A+B}}{\sqrt{I_{C}}}\left[1+\cos{(2\theta_{A}-2\theta_{C})}\right]-\varpi(I_{C}),\qquad \frac{dI_{C}}{dn} =-\frac{2\sqrt{2}}{J_{11}}I_{A+B}\sqrt{I_{C}}\sin(2\theta_{A}-2\theta_{C}),$$

(144) (145) (146) (147)

where I_A+B=I_A+I_B and Iₜ=I_A+B+I_C is still clearly a conserved quantity. This can be derived from the reduced Hamiltonian

$$\mathcal{H}=\frac{\sqrt{2}}{J_{11}}I_{A+B}\sqrt{I_{C}}\left[1+\cos{(2\theta_{A}-2\theta_{C})}\right]+\Delta(I_{C}).$$

(148)

We can simplify this if there exists a canonical transformation which invokes the conserved total action as one of our momenta. Indeed, the point transformation ϕ₁=θ_A-θ_C and ϕ₂=θ_C with conjugate momenta P₁=I_A+B and P₂=I_A+B+I_C=Iₜ yields the simplified Hamiltonian

$$\mathcal{H} =\frac{\sqrt{2}}{J_{11}}P_{1}\sqrt{P_{2}-P_{1}}\left[1+\cos{(2\phi_{1})}\right]+\Delta(P_{2}-P_{1})\qquad {}=\frac{\sqrt{2}}{J_{11}}I_{A+B}\sqrt{I_{t}-I_{A+B}}\left[1+\cos{(2\phi_{1})}\right]+\Delta(I_{t}-I_{A+B}).$$

(149)

Figure 10: Top Left: The resonant branch structure as per Fig. 8 is overlaid with green dashed lines corresponding to I_A+B=Iₜ-I_C for three choices of Iₜ={40.0,59.77,80.0}. Positions where the green dashed line intersects the upper branch result in fixed points, denoted by black dots. The remaining panels show contours of the Hamiltonian (6.2) for the different values of Iₜ. We see that as Iₜ is increased and the green line crosses the upper branch, a saddle-node bifurcation spawns a centre and saddle-point, plotted as black dots in the bottom panels.

Clearly the absence of ϕ₂ ensures Iₜ is conserved. Now we can examine slices in phase space for a choice of constant Iₜ and visualise the reduced two dimensional structure. Here trajectories of I_A+B and ϕ₁ are traced out by contours of the Hamiltonian. Examining the evolution of the phase portrait as we vary the constant Iₜ, helps us gain further insight into the upper resonant branch. In Fig. 10 we plot the phase portrait for three different values of Iₜ={40.0,59.77,80.0}. The conservation of Iₜ ensures that the trajectories must follow tracks where I_A+B=Iₜ-I_C. These are plotted as the dashed green lines in the upper left panel. The value of Iₜ sets the intercept of these lines with the I_C axis and as we increase Iₜ they are translated upwards. For Iₜ=40 the green line never intersects the red upper branch and so the phase portrait has no fixed points. When Iₜ is set such that the green line just touches the upper branch this results in a saddle node bifurcation, spawning two resonant centres. By combining the conservation of action constraint with the upper branch equation, we find that this bifurcation point occurs at

$$I_{c,saddle}=\left[\frac{12}{c(\gamma)2^{-\gamma/2}J_{11}\gamma}\right]^{-\frac{2}{2+\gamma}}.$$

(150)

Increasing Iₜ beyond this shows that the two fixed points diverge as the green dashed line intersects the upper branch in two locations. The point to the left of the saddle-node is an unstable saddle whilst the point to the right is a stable centre. This elucidates the stability structure discussed previously for the upper branch. Indeed, inputting γ=5/3 and J₁₁=100 into equation (150) yields I_[c,saddle]=9.05, which agrees with the critical value of I_C separating the stable and unstable regime as seen in Fig. 9.

An equivalent analysis can be performed for the lower branch. Now the full set of action-angle equations are reduced by restricting our attention to the case I_A=I_B and θ_B=θ_A+π for which we only permit out of phase tilt and shear motions. This is the correct simplification since we found it is the out of phase tilt and shear mode which is susceptible to the parametric instability. We again reduce the dimensionality of the original system and perform the same canonical transformation as before. This is then described by the Hamiltonian

$$\mathcal{H} =-\frac{\sqrt{2}}{J_{11}}P_{1}\sqrt{P_{2}-P_{1}}\left[1+\cos{(2\phi_{1})}\right]+\Delta(P_{2}-P_{1})\qquad {}=-\frac{\sqrt{2}}{J_{11}}I_{A+B}\sqrt{I_{t}-I_{A+B}}\left[1+\cos{(2\phi_{1})}\right]+\Delta(I_{t}-I_{A+B}).$$

(151)

Figure 11: Top Left: The resonant branch structure as per Fig. 8 is overlaid with green dashed lines corresponding to I_A+B=Iₜ-I_C for three choices of Iₜ={5.0,12.48,25.0}. Positions where the green dashed line intersects the lower branch result in fixed points, denoted by black dots. The remaining panels show contours of the Hamiltonian (6.2) for the different values of Iₜ. We see that as Iₜ is increased and the green line crosses the parametric instability threshold for the vertical breathing mode, a pitchfork bifurcation spawns two saddle-points and a centre, plotted as black dots in the bottom panels.

Again we can visualise this for slices through constant Iₜ as shown in Fig. 11. The three choices of Iₜ={5.0,12.48,25.0} correspond to the three green dashed lines in the upper left panel, along which I_A, I_B and I_C are constrained to move. When Iₜ=1.0 the green line intersects the I_C axis before the onset of parametric instability and we see that the purely vertical oscillator is stable. However when Iₜ=12.5 the system crosses the bifurcation point defined by equation (143). Beyond this, we see the formation of unstable saddle-points along the I_A+B=0 axis which once again emphasises the instability of purely vertical breathing modes here. The intercept of the green dashed line with the lower branch yields the stable centre for the anti-phased mixed mode. Of course this agrees with the stable behaviour for the lower branch, found earlier using linear perturbation techniques.

6.3 Interpreting these solutions

Having developed a thorough understanding of our resonant equilibria it is important to relate these back to their physical interpretation. Returning to the more intuitive Jacobian coordinates, these stable action branches correspond to constant amplitude oscillatory solutions for the variables J₁₃, J₃₁ and J₃₃. The process of mapping from bounce to bounce of the non-linear vertical oscillator effectively removes the harmonic motion in between. In this sense, our technique effectively captures the slow timescale associated with amplitude and phase evolution.

Both resonant branches predict a highly non-linear family of bouncing modes with excited shearing and warp. Since I_A=I_B, the J₁₃ and J₃₁ oscillators are excited with equal amplitude, akin to the equipartition seen in the linear modes for which κ=ν. Whilst the resonant angles are constant, θ_A, θ_B and θ_C each advance at a steady rate,

$$\frac{d\theta_{A,B,C}}{dn}=\pm\frac{2\sqrt{2I_{C}}}{J_{11}},$$

(152)

where the plus and minus signs correspond to the upper and lower branches respectively. The progression of tilt, shear and bouncing phase results in a precession of the modes when viewed from a global reference frame. To see this, consider a global ring with azimuthal variation in tilting and thickness. If this torus is fixed in space, an orbiting fluid parcel would see the periodic structure pass by at the orbital frequency with a fixed phase set by the azimuthal origin. However, if the structure rotates and the azimuthal origin evolves, the orbiting observer would see this time dependent phase manifest as a modification to the periodic frequency. This precessional frequency is then defined by

$$\omega_{p}\equiv-\frac{d\theta}{dt}=-\frac{d\theta/dn}{\pi-d\theta/dn}.$$

(153)

For the upper branch we have retrograde precession ωₚ<0 and for the lower branch prograde precession ωₚ>0. Within the local model these predict periodic solutions with angular frequency ω=Ω-ωₚ, so when the global torus rotates with the orbit, the local frequency decreases. Meanwhile if the ring rotates against the orbit, the local frequency increases. Ogilvie & Latter (2013) previously showed that discs with a fixed global warping geometry permit periodic solutions provided the epicyclic and vertical frequencies are sufficiently detuned or a viscosity is introduced to temper the resonant flows. Here however, we see the Keplerian resonance drives a precession of the ring which acts as an effective detuning from the orbital frequency.

7 Numerical verification

The smooth modulation theory developed in section 4 and the bouncing theory developed in sections 5 and 6 may now be tested by returning to our full equation set (29) – (32) and numerically finding the periodic solutions. We select the same parameters as described in the setup of section 3.1 which we will now reiterate. Of course we are examining the resonant case with κ=ν=Ω and adopt units so Ω=L=1. We take γ=5/3 and choose the characteristic temperature and C_z circulation constant so that the equilibrium ring has J₁₁=100 and J₃₃=1, corresponding to an aspect ratio of ε=0.01.

We proceed with the same shooting scheme previously used to identify the periodic solutions for the forced vertical oscillator in section 5.1.4. However, as our analysis has shown, the feedback of the vertical oscillator onto the warp results in a phase modulation of the tilt and shear oscillations. These may be interpreted as precessing modes with a period which now deviates from the orbital timescale. Thus our shooting code is generalised to incorporate the period as a parameter which should also be determined. Furthermore, our theory predicts that the periodic solutions correspond to the nonlinear extension of bending waves for which the tilt and shear are in equipartition with phase relationship 0 or π. We use this to inform our initial guesses in the shooting method.

We converge to the periodic branch structure which is plotted in the left panel of Fig. 12. The solid lines mark the periodic solutions found, whilst the dashed lines correspond to the solution branches predicted from our theory. The shooting method solutions are identified in terms of the Lagrangian variables, Jᵢⱼ, which are then approximately converted into action variables by identifying the maximum values of 1/2J̇₁₃²=1/2J̇₃₁² as the warping action I₀=I_A=I_B which is plotted along the y-axis. Then the value of J₃₃ is extracted at times for which the forcing product J₁₃ J₃₁=0, such that 1/2J₃₃² is our proxy for the vertical action variable I_C which is plotted along the x-axis. For each identified solution we perform a Floquet stability analysis, as per the method described in section 5.1.4. The maximum eigenvalue from the computed monodromy matrix determines the colour along the branches, with values greater than 1 (departing from purple) indicating instability. In the right hand panel we plot the period of these solutions against the warping amplitude as the solid black lines. Again these are compared with the analytical theory predictions which are plotted as dashed and dotted lines.

We number the qualitatively distinct branches (i)–(iv), and show typical solutions for each regime in the rows of Fig. 13. In branches (i) and (ii) we see the anti-phased and in-phase smooth nonlinear branches respectively. These stem from the equilibrium configuration for which there is no tilt or shear and a constant value of J₃₃=1 and J₁₁=100 such that the aspect ratio of the thin base state is ε=0.01. As we expect from the continuation of the averaged Lagrangian for X<0 in section 4.5, the solutions for the anti-phased tilt and shear may be continued indefinitely to large warp amplitudes. Indeed, our smooth modulation theory agrees very well as indicated by the over-plotted blue dashed line. This plots the value J₃₃ inherited from the periodic solutions found for equation (46) at times for which the the forcing product J₁₃,₀ J₃₁,₀=0 (setting a consistent phase relationship with the warp as compared with the choice described above). The right panel of Fig. 12 shows that the period for this branch is slightly greater than the orbital period and agrees very well with the precessional frequency offset as deduced from the gradient of the average Lagrangian, as per equation (79). Branch (ii) meanwhile shows a more interesting behaviour. The red dashed line from the modulation theory agrees very well with the identified periodic structures for low to intermediate warp amplitudes. There is also good agreement for the predicted period within this range, which is slightly less than the orbital period as expected for the extension of the in-phase bending modes. However, the red dashed line eventually terminates at the saddle-node bifurcation, as seen for the computed average Lagrangian at some critical in-phase forcing X>0 – see Fig. 4. Beyond this point the modulation theory breaks down and we expect some different behaviour to arise.

Here, the periodic solutions begin to deviate from our modulation theory and bend round onto branch (iv). Now the vertical oscillator action begins to grow rapidly as it enters into the extreme bouncing regime. The red dashed line, showing the predicted in-phase bouncing centres as described by equation (134), converges to the periodic solutions as the bounce amplitude increases. Note, the periodic solution space identified avoids the unstable portion of the upper bouncing branch since the vertical action is in fact too low here and the bouncing approximations break down. Instead there is a smooth transition connecting onto the modulation theory. The period of these in-phase bouncing solutions also shows a dramatic change in behaviour as the retrograde detuning from the orbital rate becomes more pronounced. At large warp amplitudes (and hence bouncing amplitudes), the analytical period predictions deduced from equations (152) and (153) agree well.

Along the I_C axis of Fig. 12 we see the non-linear vertical mode with no tilt and shear activation. As discussed in section 6.2, this undergoes parametric instability and spawns the lower anti-phased bouncing branch as labelled by (iii). Beyond this point the departure from purple colouration emphasises the instability of the pure bouncing mode with no warp activation. The lower bouncing branch incurs both growing tilt/shear and extreme bouncing motions as predicted from the lower branch of equation (134). This analytical result is over-plotted as a dashed blue line which agrees remarkably well and nicely intersects the parametric instability threshold along the x-axis. This correspondence with theory is further confirmed in the period plot where there is almost perfect overlap between the dashed blue line and black line in branch (iii). We see that the bouncing solution incurs a large prograde departure from the orbital frequency as T becomes longer for larger warp amplitudes.

Figure 12: Left panel: Branches of periodic solutions are plotted as thick lines, coloured according to the maximum eigenvalue of the monodromy matrix associated with each periodic solution. The x-axis plots the initial value of J₃₃²/2 which is a measure of the vertical action I_C. The y-axis plots the maximum value of J̇₁₃²/2=J̇₃₁²/2 which measures the warping action I₀. Four branches of qualitatively different solutions are labelled (i)–(iv). Branches (i) and (ii) represent the smoothly modulated anti-phased and in-phase tilt and shear solutions respectively, stemming from the equilibrium at J₃₃=1. The analytical predictions are over-plotted as dotted blue and red lines. Branches (iii) and (iv) then represent the extreme bouncing regime. The over-plotted blue and red dashed lines show the correspondence with theory encapsulated in the branch equations (134). Right panel: The converged period normalised against the orbital value is plotted as a black line against the I₀ warping action. The qualitative regimes (i)–(iv) are identified corresponding to the different solution branches in the left panel. The dotted red and blue lines over-plotting branches (i) and (ii) show the predicted period according to our smooth modulation theory. Similarly the dashed red and blue lines overlying branches (iii) and (iv) plot the period predicted from our impulsive bouncing theory. Note the in-phase branches fall below the orbital period whilst the anti-phased solutions are longer than the orbital period.

Figure 13: Example periodic solutions from the four key branches identified in Fig. 12 are labelled (i)–(iv). The left columns plot the J₁₃ shear (black) and J₃₁ tilt (red) variables across one period. The right hand panels then plot the corresponding evolution of J₃₃.

8 Discussion

These numerical results confirm our smooth modulation theory and the connection to the predicted bouncing regime. In both cases the key effect is the feedback of the vertical oscillator onto the warp which has not been taken into account in previous work. Here we see that the period of the solutions deviates from the orbital value in order to circumvent the Keplerian resonance for which κ=ν=Ω.

The periodic solution branches found using our local model may be reinterpreted as large-scale precessing structures when viewed from a non-rotating, global reference frame. In Paper I we saw that by Doppler shifting the linear tilting modes of our ring model into the non-rotating frame, they may be interpreted as global bending waves. Essentially the orbital time within the local model can be mapped onto the azimuthal coordinate as the shearing box performs its orbit. The ring evolution over the orbital timescale simply corresponds to the azimuthal variation in the geometry of the disc as elucidated in section 2.3.

We might then think of our tilting ring as a model which approximately zooms in on a local patch of a globally warped disc. Thus we expect the qualitative solution families found in this paper to be applicable to a radially extended, globally warped Keplerian disc. The linear tilting modes extend into branches (i) and (ii) where the smooth modulation theory applies. The anti-phased solutions have prograde precession whilst the in-phase solutions exhibit a retrograde precession. We have focused on finding special periodic solutions for which the amplitude of the tilt and shear are constant. However, we might speculate that the combination of general tilt and shear initialisations, as plotted in the middle panel of Fig. 1 and the upper four panels of Fig. 2 for example, might be some modified superposition of these nonlinear precessing modes. We see that the retrograde precession dominates over the prograde precession as the warp amplitude increases and branch (ii) bends away from the orbital period in the right panel of Fig. 12. Hence we can expect a typical retrograde bias for warped structures.

Crucially we found that this behaviour breaks down as the tilt and shear grow to sufficient amplitudes. For γ=5/3 we found that the smooth modulation theory breaks down for Z_[1,c]=0.4, where the solutions for the forced vertical oscillator terminate in a saddle node bifurcation. More generally, Z₁ can can be connected with the global warping amplitude by the following scaling argument. Consider the radial tilting of the reference midplane line z₀=0 and the shearing of the vertical axis x₀=0. If the vertical and horizontal displacements from equilibrium are denoted by ξ_z and ξₓ respectively, the associated gradients are ∂ξ_z/∂ x∼ J₃₁/J₁₁=ψ (where ψ is the warp amplitude) and ∂ξₓ/∂ z∼ J₁₃/J₃₃. For a Keplerian bending wave, equipartition of tilt and shear energy demands that the displacements are of the same order. Identifying the typical warp length scale λ as the width of our ring and taking the scale height H, the characteristic tilt and shear displacements balance provided ∂ξₓ/∂ z∼(ψλ)/H. Noting the relation Z₂ ∼ J₁₃ J₃₁ ∼ J₁₁ J₃₃ Z₁ as described by equation (81) we see that

$$Z_{1}\sim\frac{J_{13}J_{31}}{J_{11}J_{33}}\sim\frac{\partial\xi_{z}}{\partial x}\frac{\partial\xi_{x}}{\partial z}\sim\psi^{2}\frac{\lambda}{H},$$

(154)

so the critical warp amplitude scales as ψ_c ∼√(H/λ), i.e. as √(H/r) in the case of a global warp (λ∼ r). Whilst we have used ε=0.01 throughout the course of this paper to emphasise that the warping length scale is much longer than the disc scale-height, our results also extend through to thicker discs with ε=0.1 where we have verified the critical warp scaling law above. Beyond this value of warp, we would expect extreme vertical bouncing motions to be activated in the warped disc. This would correspond to the transition towards branch (iv), where the global warped geometry indicated by the oscillating J₃₁ component is now accompanied by extreme compression of J₃₃ twice per orbit, as seen in Fig. 13. The disc would present locations which are extremely thin, whilst other regions are vertically extended. This may lead to observational signatures sensitive to enhanced density. Furthermore, puffed up regions or sufficient warp amplitudes may obscure light from a central source and cast shadows as found in the various observations discussed in section 1.1.

The solution families predicted here are found using ideal hydrodynamics where we have no dissipation, despite the extreme compressive behaviour. However, by incorporating some viscosity prescription and a more general energy equation we might expect that the compressive motions would lead to a significant damping of the warp. Indeed, similar ‘nozzle-like’ compressive structures occur in eccentric disc models of tidal disruption events (TDEs) wherein bouncing modes are forced periodically as gravity is enhanced at pericenter (Ogilvie & Barker, 2014; Lynch & Ogilvie, 2020). These motions may release significant amounts of energy as the gas is compressed at closest approach (Zanazzi & Ogilvie, 2020; Ryu et al., 2021). Furthermore, global warped disc simulations performed by Sorathia et al. (2013) exhibit an enhanced damping of the warp. This is not explained in their paper but might be attributed to the conversion of warp action to extreme vertical motions, via the nonlinear mode coupling, which is then damped due to the bulk artificial viscosity. Future numerical work should examine if these extreme phenomena are in fact present in the simulations and then establish observational consequences.

In fact the prediction of a critical warp amplitude in our work is reminiscent of the recent quest to understand ring breaking phenomena which are believed to occur in sufficiently warped discs (e.g. Nixon & King, 2012; Doǧan et al., 2018). The interplay and connection between our critical warp with this previous work is unclear and merits future investigation. Indeed, a variety of other effects might modify our solution families, including the parametric instability proposed by Gammie et al. (2000). This has been shown to be active in global disc simulations by Deng et al. (2021) and may present an enhanced turbulent viscosity affecting the evolution of our periodic modes. In the future, we propose setting up numerical simulations which target the internal flow structure of warped discs and test how robust they are in the presence of more general physics.

9 Conclusions

In this paper we have performed an extensive nonlinear analysis of the local ring model equations derived in Paper I and uncovered two distinct regimes relevant to the nonlinear dynamics of warped Keplerian discs. We find the extension of the linear bending modes at larger warp amplitudes is well described using an asymptotic averaged Lagrangian theory whereby the amplitude and phase of the warp smoothly vary over a long timescale. However, beyond some critical warp (which scales as the aspect ratio of the ring or disc), the in-phase product of tilt and shear motions resonantly force the vertical oscillation to large amplitudes. The disc becomes extremely compressed and feeds back impulsively onto the warp. We have identified periodic solutions using a variety of careful approximations which have then been confirmed within the full equation set. These local modes map onto globally precessing warped structures with compressions and expansions twice per orbit. These regions could manifest observationally as regions of enhanced emission or by casting shadows to outer regions of the disc. Although we have analytically extracted special solutions, we expect these compressive motions to be present in more general setups, as evidenced in our motivating numerical experiments. This may have profound consequences for the evolution of warped discs as such compressions might lead to an enhanced dissipation of energy and warp in Keplerian systems. This demands attention in future numerical simulations, with detailed analysis of the flow structure as the warp amplitude is varied.

Acknowledgements

The authors would like to thank the anonymous reviewer for their helpful comments and suggestions. This research was supported by an STFC studentship and STFC grants ST/P000673/1 and ST/T00049X/1.

Data Availability

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

References

• Bainov & Simeonov (1993)

  Bainov D., Simeonov P., 1993, Impulsive Differential Equations: Periodic Solutions and Applications. Monographs and Surveys in Pure and Applied Mathematics, Longman, Essex

• Benisty et al. (2017)

  Benisty M., et al., 2017, A&A, 597, A42

• Casassus et al. (2018)

  Casassus S., et al., 2018, MNRAS, 477, 5104

• Debes et al. (2017)

  Debes J. H., et al., 2017, ApJ, 835, 205

• Dednam & Botha (2014)

  Dednam W., Botha A. E., 2014, Engineering with Computers, 31, 749–762

• Deng et al. (2021)

  Deng H., Ogilvie G. I., Mayer L., 2021, MNRAS, 500, 4248

• Doǧan et al. (2018)

  Doǧan S., Nixon C. J., King A. R., Pringle J. E., 2018, MNRAS, 476, 1519

• Facchini et al. (2013)

  Facchini S., Lodato G., Price D. J., 2013, MNRAS, 433, 2142

• Facchini et al. (2017)

  Facchini S., Juhász A., Lodato G., 2017, MNRAS, 473, 4459

• Fairbairn & Ogilvie (2021)

  Fairbairn C. W., Ogilvie G. I., 2021, MNRAS

• Gammie et al. (2000)

  Gammie C. F., Goodman J., Ogilvie G. I., 2000, MNRAS, 318, 1005

• Hatchett et al. (1981)

  Hatchett S. P., Begelman M. C., Sarazin C. L., 1981, ApJ, 247, 677

• Hawley et al. (1995)

  Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742

• Hill (1878)

  Hill G. W., 1878, American Journal of Mathematics, 1, 5

• Holmes (1982)

  Holmes P., 1982, Journal of Sound and Vibration, 84, 173

• Katz (1973)

  Katz J. I., 1973, Nature Physical Science, 246, 87

• Kotze & Charles (2012)

  Kotze M. M., Charles P. A., 2012, MNRAS, 420, 1575

• Kraus et al. (2020)

  Kraus S., et al., 2020, Science, 369, 1233

• Lodato & Price (2010)

  Lodato G., Price D. J., 2010, MNRAS, 405, 1212

• Loomis et al. (2017)

  Loomis R. A., Öberg K. I., Andrews S. M., MacGregor M. A., 2017, ApJ, 840, 23

• Lubow & Ogilvie (2000)

  Lubow S. H., Ogilvie G. I., 2000, ApJ, 538, 326

• Luo & Han (1996)

  Luo A. C. J., Han R. P. S., 1996, Nonlinear Dynamics, 10, 1

• Lynch & Ogilvie (2020)

  Lynch E. M., Ogilvie G. I., 2020, MNRAS, 500, 4110

• Marino et al. (2015)

  Marino S., Perez S., Casassus S., 2015, ApJ, 798, L44

• Miyoshi et al. (1995)

  Miyoshi M., Moran J., Herrnstein J., Greenhill L., Nakai N., Diamond P., Inoue M., 1995, Nature, 373, 127

• Muro-Arena, G. A. et al. (2020)

  Muro-Arena, G. A. et al., 2020, A&A, 635, A121

• Nixon & King (2012)

  Nixon C. J., King A. R., 2012, MNRAS, 421, 1201

• Ogilvie (1999)

  Ogilvie G. I., 1999, MNRAS, 304, 557

• Ogilvie (2006)

  Ogilvie G. I., 2006, MNRAS, 365, 977

• Ogilvie & Barker (2014)

  Ogilvie G. I., Barker A. J., 2014, MNRAS, 445, 2621

• Ogilvie & Latter (2013)

  Ogilvie G. I., Latter H. N., 2013, MNRAS, 433, 2403–2419

• Papaloizou & Lin (1995)

  Papaloizou J. C. B., Lin D. N. C., 1995, ApJ, 438, 841

• Papaloizou & Pringle (1983)

  Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181

• Petterson (1977a)

  Petterson J. A., 1977a, ApJ, 214, 550

• Petterson (1977b)

  Petterson J. A., 1977b, ApJ, 216, 827

• Pinilla et al. (2015)

  Pinilla P., et al., 2015, A&A, 584, L4

• Pringle (1992)

  Pringle J. E., 1992, MNRAS, 258, 811

• Rosenfeld et al. (2012)

  Rosenfeld K. A., et al., 2012, The Astrophysical Journal, 757, 129

• Ryu et al. (2021)

  Ryu T., Krolik J., Piran T., 2021, arXiv e-prints, p. arXiv:2105.09434

• Sakai et al. (2019)

  Sakai N., Hanawa T., Zhang Y., Higuchi A. E., Ohashi S., Oya Y., Yamamoto S., 2019, Nature, 565, 206

• Sorathia et al. (2013)

  Sorathia K. A., Krolik J. H., Hawley J. F., 2013, ApJ, 768, 133

• Stolker et al. (2016)

  Stolker T., et al., 2016, A&A, 595, A113

• Whitham (1965)

  Whitham G. B., 1965, Journal of Fluid Mechanics, 22, 273–283

• Zanazzi & Ogilvie (2020)

  Zanazzi J. J., Ogilvie G. I., 2020, MNRAS, 499, 5562