A non-linear theory of vertical resonances in accretion discs

An important and widely neglected aspect of the interaction between an accretion disc and a massive companion with a coplanar orbit is the vertical component of the tidal force. As shown by Lubow, the response of the disc to vertical forcing is resonant at certain radii, at which a localized torque is exerted, and from which a compressive wave (p mode) may be emitted. Although these vertical resonances are weaker than the corresponding Lindblad resonances, the m=2 inner vertical resonance in a binary star is typically located within the tidal truncation radius of a circumstellar disc. In this paper I develop a general theory of vertical resonances, allowing for non-linearity of the response, and dissipation by radiative damping and turbulent viscosity. The problem is reduced to a universal, non-linear ordinary differential equation with two real parameters. Solutions of the complex non-linear Airy equation are presented to illustrate the non-linear saturation of the resonance and the effects of dissipation. It is argued that the m=2 inner vertical resonance is unlikely to truncate the disc in cataclysmic variable stars, but contributes to angular momentum transport and produces a potentially observable non-axisymmetric structure.


Horizontal and vertical tidal forces
Accretion discs in binary stars and protoplanetary systems are subject to periodic tidal forcing from the orbiting companion objects. Similar effects can occur in planetary ring systems, where the planetary satellites provide the perturbing forces, and in galactic discs. Even in the simple case of a companion with a circular orbit that is coplanar with the disc, the effects of the tidal force are manifold.
Consider a system of two spherical bodies of masses M1 and M2 in a bound orbit. Consider a third, test body orbiting about M1. In a non-rotating coordinate system centred on M1, the force per unit mass experienced by the test body is where d(t) is the position vector of the companion M2 with respect to M1, and the third term is the fictitious force arising from the acceleration of the origin of the coordinate system. Inasmuch as the second and third terms perturb the Keplerian motion of the test body about M1, they represent the tidal force in this system. Consider a thin accretion disc around M1, and let (R, φ, z) be cylindrical polar coordinates such that z = 0 is the midplane of the disc. Suppose that the companion has a circular, coplanar orbit of radius d and angular velocity ω. Then the components of the tidal acceleration are f tφ = −GM2[R 2 + d 2 − 2dR cos(φ − ωt) + z 2 ] −3/2 d sin(φ − ωt) + GM2d −2 sin(φ − ωt), The horizontal components of the tidal acceleration are nearly independent of z within the disc and cause a general non-axisymmetric distortion of its streamlines and surface density distribution (Papaloizou & Pringle 1977). Indeed, the streamlines of the disc correspond closely to a family of simple periodic orbits of the restricted three-body problem found by Paczyński (1977). The tidal distortion is enhanced in the neighbourhood of Lindblad resonances, radii at which the forcing frequency resonates with the natural epicyclic oscillations of the disc. Usually at a Lindblad resonance, a non-axisymmetric wave will be launched and will then propagate some way radially through the disc before dissipating and transferring its angular momentum to the disc (Goldreich & Tremaine 1979). The resonant torque exerted between the companion and the disc may be calculated from a standard formula that applies under quite general circumstances in linear theory, for example if the response is dominated by viscosity and no wave is emitted (Meyer-Vernet & Sicardy 1987), or if the vertical structure of the disc is taken fully into account (Lubow & Ogilvie 1998).
Tidal torques usually limit the outer radius of a circumstellar disc (or the inner radius of a circumbinary disc). Except in the case of a very low-mass companion, the Lindblad resonances are so strong that they are excluded from the disc. In a binary star with a mass ratio q = M2/M1 of order unity, a circumstellar disc is tidally truncated well inside the innermost (m = 2) Lindblad resonance, which in any case lies outside the Roche lobe (Paczyński 1977;Papaloizou & Pringle 1977). For smaller mass ratios typical of giant planets orbiting stars, the Lindblad resonances play a more direct role in truncating the disc (Lin & Papaloizou 1986).
Another horizontal tidal effect is the eccentric interaction between the companion and the disc. At eccentric Lindblad resonances such as the 3:1 resonance, a local eccentric instability occurs, which may be able to sustain a large-scale eccentric distortion of the disc (Lubow 1991). The eccentric corotational resonances act to damp eccentricity, however.
Studies of the vertical component of the tidal force have mostly been concerned with the case of a companion with an inclined orbit. In that case, the vertical tidal acceleration is nearly independent of z within the disc. As a result, nonaxisymmetric bending waves are launched at vertical resonances where the forcing frequency resonates with the natural vertical oscillations of the disc (Shu, Cuzzi & Lissauer 1983). In a Keplerian disc the Lindblad and vertical resonances coincide.
For a companion with a coplanar orbit, however, the vertical tidal acceleration is proportional to z and vanishes on the mid-plane. Other than those concerned with bending waves, most analytical and numerical studies have relied on a twodimensional treatment of the disc, and the vertical component of the tidal force has been widely neglected. However, Lubow (1981) analysed the response of a vertically isothermal disc to vertical tidal forcing by a companion with a circular, coplanar orbit. The general response consists of a non-axisymmetric vertical expansion and contraction of the disc, characterized by a vertical velocity proportional to z. This 'breathing' motion corresponds to a compressive mode of the disc, the n = 1 p mode in the notation of Lubow & Pringle (1993), or the p e 1 mode in the classification of Ogilvie (1998). At particular radii where the forcing frequency resonates with the natural frequency of this mode, a non-axisymmetric p-mode wave may be launched and a resonant torque exerted. These are also vertical resonances, but are distinct from those associated with bending waves because a different mode is involved. Lubow's vertical resonances are further from the companion than the corresponding Lindblad resonances and, in a thin disc, the resonant torques are much weaker. However, Lubow (1981) pointed out that, in a binary star with a mass ratio of order unity, the m = 2 inner vertical resonance typically lies within the standard tidal truncation radius, and the associated torque can compete with the viscous torque in the disc.
Another consequence of the vertical tidal force is a local tilt instability, closely related to the eccentric instability (Lubow 1992). The tilt instability occurs at inclination resonances that are essentially coincident with the eccentric Lindblad resonances, but the tilt instability is much weaker.

Aims of the paper
In this paper I develop a non-linear theory of Lubow's vertical resonances. I generalize the analysis of Lubow (1981) to nonisothermal discs allowing for non-linearity, viscosity and radiative damping. By reducing the problem to a universal differential equation and studying its solutions, I aim to demonstrate the effects of non-linear saturation and dissipation on the resonant response of differentially rotating discs to periodic forcing in general.
In a related paper (Ogilvie 2001, hereafter Paper I) I have presented an analysis of the non-linear tidal distortion of a thin, three-dimensional accretion disc by a binary companion on a circular orbit. I showed that the resulting distortion can plausibly explain the two-armed non-axisymmetric features seen in the Doppler tomograms of IP Peg and other dwarf novae in outburst (e.g. Steeghs 2001). Although the m = 2 inner vertical resonance plays an important role in regulating the amplitude and phase of the tidal distortion, the method of analysis in Paper I does not allow for the emission of a wave from the vertical resonance. One of the aims of the present paper is to justify that limitation by showing that, under typical conditions, vertical resonances are broadened and damped by non-linearity and dissipation to the extent that the emission of a wave may reasonably be neglected.
Before embarking on the detailed analysis, it may be helpful to look ahead to the equation to be derived, In this equation, x represents the radial distance from the resonant orbit and y the complex amplitude of the response, both in dimensionless terms. The linear terms −y ′′ + xy represent the ability of the disc to support freely propagating waves in the region x < 0. In the present case these waves correspond to the compressive p e 1 mode. The point x = 0 is the location of the resonance, which is also the turning point of free waves. The constant term a on the right-hand side represents the tidal forcing in the neighbourhood of the resonance. The inhomogeneous Airy equation −y ′′ + xy = a is familiar in studies of resonant wave excitation in differentially rotating discs and features in the analysis of Lubow (1981). The new, non-linear term |y| 2 y derives from couplings between the resonant mode and non-resonant modes. It leads to new effects such as the non-linear broadening and saturation of the resonance. The correct derivation of this term accounts for much of the complexity of the analysis in the present paper. The linear term iby represents the dissipation of vertical motions by shear and bulk viscosity and by radiative damping. It leads to attenuation of the waves and broadening of the resonance.
The remainder of this paper is organized as follows. In Section 2 I recall the coordinate system and the basic equations used in Paper I. In Section 3 I expand the equations in the neighbourhood of the vertical resonance and solve them in a systematic manner to derive the complex non-linear Airy equation. In Section 4 I discuss the properties of the equation and present numerical solutions for a range of parameter values. The astrophysical consequences of the analysis are explored in Section 5.

Orbital coordinates
Consider a binary system with two stars in a circular orbit. The stars are approximated as spherical masses M1 and M2. Let d be the binary separation and the binary orbital frequency. Consider the binary frame rotating about the centre of mass with angular velocity ω, and let (R, φ, z) be cylindrical polar coordinates such that star 1 (about which the disc orbits) is at the origin, while star 2 has fixed coordinates (d, 0, 0). In Paper I a system of non-orthogonal orbital coordinates (λ, φ, z) was introduced instead of cylindrical polar coordinates. These are based on Paczyński's orbits for the restricted three-body problem, and naturally account for the principal tidal distortion of the disc. The coordinate λ is a quasi-radial coordinate that labels the orbits and differs only slightly from R, reducing exactly to R in the limit M2 → 0. In this paper I again consider the case of a prograde circumstellar disc. However, the generalization to retrograde and circumbinary discs is straightforward.
The standard notation of tensor calculus is adopted, with g ab being the metric tensor, J = det(g ab ) the Jacobian of the coordinate system and η abc the Levi-Civita alternating tensor. Expressions for the metric coefficients and connection components for the orbital coordinate system may be found in Paper I, Section 2.

Basic equations
The basic equations governing a fluid disc in three dimensions are the same as in Paper I. However, the effect of a non-zero relaxation time of the turbulent stress will not be considered here, as it was found to introduce algebraic complications without affecting the final result of this paper significantly. Accordingly, a purely viscous model of the turbulent stress is adopted here.
The equation of mass conservation is where ρ is the density and u a the velocity relative to the rotating frame. The equation of motion is Here ω a is the angular velocity of the rotating frame (only ω z = ω being non-zero), p is the pressure, is the part of the effective potential not due to star 2, and is the tidal potential due to star 2. The turbulent stress tensor T ab is assumed to be given by where µ is the effective shear viscosity and µ b the effective bulk viscosity. The energy equation is where γ is the adiabatic exponent and F a the radiative energy flux, given in the Rosseland approximation for an optically thick medium by where σ is the Stefan-Boltzmann constant, T the temperature and κ the opacity. The equation of state of an ideal gas, is adopted, where k is Boltzmann's constant, µm the mean molecular weight and mH the mass of the hydrogen atom. The opacity is assumed to be of the power-law form where Cκ is a constant. This includes the important cases of Thomson scattering opacity (x = y = 0) and Kramers opacity (x = 1, y = −7/2). The effective viscosity coefficients are assumed to be given by an alpha parametrization. The precise form of the alpha prescription relevant to a tidally distorted disc is to some extent debatable. It will be convenient to adopt the form where α and α b are the dimensionless shear and bulk viscosity parameters. In the limit of a circular disc in the absence of star 2, this prescription reduces to the usual one, µ = αp/Ω, etc., whereΩ is the orbital angular velocity in the inertial frame. For simplicity, the five dimensionless parameters of the disc, (α, α b , γ, x, y) will be assumed to be constant throughout this paper. The other parameter in the problem is the mass ratio q = M2/M1.

Thin disc
Analytical progress can be made only if the disc is thin so that certain approximations may be used. The systematic way of introducing these approximations is through an asymptotic analysis. Let a system of units be adopted such that the radius of the disc and the orbital frequency are O(1). Let the small parameter ǫ be a characteristic value of the angular semi-thickness H/R of the disc. To resolve the vertical structure of the disc, introduce a stretched vertical coordinate ζ according to so that ζ = O(1) within the disc. Any slow time-dependence of the unperturbed disc, which occurs on a characteristic viscous time-scale α −1 ǫ −2 Ω −1 , is consistently neglected. The desired solution is stationary in the binary frame.

Orbital motion
denote the total potential in the binary plane. The equations governing the shape R(λ, φ) of Paczyński's orbits and their orbital angular velocity Ω(λ, φ) can be expressed in the form (Paper I, Section 5.1) where R φ = (∂R/∂φ) λ , etc. These also correspond to the horizontal components of the equation of motion of the disc at leading order in ǫ. The divergence of the orbital motion is The primary potential can be expanded in a Taylor series about the mid-plane, The tidal potential can be expanded similarly, but is assumed formally to be smaller by a factor O(ǫ) than the primary potential, i.e.
t , and Paczyński's orbits can be considered as deformed circular orbits such that the orbital quantities possess expansions of the form where Let denote the angular velocity in the inertial frame. An explicit calculation of R (1) , Ω (1) and ∆ (1) is deferred until Section 3.8. At present it is sufficient to note that they possess Fourier expansions of the form where R (1)(m) and Ω (1)(m) are real owing to the symmetry of the orbits, while ∆ (1)(m) is imaginary.
The assumption that the tidal potential scales as O(ǫ) is a formal device that will turn out to provide a critical level of non-linearity in the response. Similarly, it is expedient to choose viscosity parameters that scale formally as O(ǫ 2/3 ), i.e.
This will turn out to provide a critical level of dissipation in the response.

Expansion about the resonant orbit
Let λ = λ * be the resonant orbit, to be identified subsequently, and let a subscript * generally denote a quantity evaluated on the resonant orbit. It is known from Lubow (1981) that the characteristic radial extent of the resonance is (H 2 R) 1/3 , as for a Lindblad resonance. One therefore introduces a scaled quasi-radial coordinate ξ according to so that ξ = O(1) in the resonant region. Any function of λ may then be expanded in a Taylor series about the resonant orbit. Thus where etc., where Φ (2) * =Ω 2 * and Φ (2) 2/3 = 2Ω * Ω 2/3 .

Fluid variables
A steady solution of the basic equations is sought. After careful analysis, the following expansion scheme for the fluid variables in the resonant region has been found to be appropriate and self-consistent.
Here s is a positive parameter, which drops out of the analysis, although formally one requires s = (10 − 6y)/(6 + 3x) in order to balance powers of ǫ in the opacity law. Other components of F a are smaller and will not be required. According to this scheme, the relative amplitude of the perturbation in the resonant region is O(ǫ 1/3 ). This is the magnitude of the fractional perturbations of density and pressure, and of the Mach number of the vertical velocity perturbation. The relative amplitude of the perturbation is greater than the relative magnitude of the tidal potential, O(ǫ), because of the effect of the resonance. The horizontal velocity perturbations are smaller than the vertical one by a factor O(ǫ 1/3 ).

Equations to be solved
The above expansions for the fluid variables are now substituted into the basic equations, which are expanded only as far as is required for the analysis that follows.

Unperturbed vertical structure
From equation (55), one has the important relation which expresses the vertical hydrostatic equilibrium of the unperturbed disc. The final two terms in equation (51) may be taken to balance, since these represent the unperturbed thermal equilibrium of the disc, i.e.

Form of the resonant mode
Equations (47), (50) and (56) may be combined to give where is a linear operator. Normally L would be a non-singular operator and equation (67) would have no non-trivial solution. However, the resonant orbit is, by definition, one on which L is a singular operator possessing a null eigenfunction of the form such that Lw = 0. This requires which is the condition given by Lubow (1981) for a vertical resonance for azimuthal wavenumber m (m > 0 always). Accordingly, the general solution of equation (67) is where f (ξ) is a dimensionless function to be determined. The corresponding solutions of equations (47) and (50) are One also finds This motion represents a simple 'breathing mode' of the disc. In the notation of Ogilvie (1998) it is the p e 1 mode, i.e. the first p mode of even symmetry about the mid-plane.

Solvability condition for the linear operator
At higher orders one naturally obtains equations of the form where X is unknown and F is known. Now L is a singular operator, and it is important to understand the conditions under which such an equation is soluble. Since L is self-adjoint, the solvability condition is easily obtained in the form where w is the null eigenfunction (69), and the integration is over the full azimuthal and vertical extent of the disc.

Horizontal velocities
Equations (53) and (54) may be combined to give The solution is where equation (70) has been used.

Relation to physical units
When the asymptotic scalings are removed, the dimensionless parameters can be expressed directly in terms of the physical variables evaluated at the location of the resonance. The ordering parameter ǫ then cancels out. However, the integral involved in C1 introduces the dimensionless parameterǫ defined by This is a more specific measure of H/R, and depends on the vertical structure of the disc. For a vertically isothermal disc, where H is the isothermal scale height. For a vertically polytropic disc with polytropic index 3/2, where H is the true semi-thickness. For a vertically polytropic disc with polytropic index 3,ǫ 2 = (1/44)(5γ 2 − 8γ + 12)(H/R) 2 . A quantity deriving from C2 is the rate of detuning of the vertical resonance, Another quantity having the dimensions of frequency-squared is the forcing combination appearing in C5, One then finds a = γ 2 (γ + 1)(3 − γ) 12 This shows why the tidal potential was chosen to scale formally as O(ǫ), in order to have a = O(1). Similarly, The scaling between x and the physical length λ is given by

Evaluation of the forcing terms
There are three contributions to the forcing coefficient Ψ. The direct vertical forcing involves two terms: Φ (2)(m) 1 , which is the perturbation of the primary potential due to the tidal distortion of the orbit, and Φ (3)(m) t * , which is the tidal potential itself. The term proportional to ∆ (m) 1 is an indirect forcing of the p mode resulting from the divergence of the tidally distorted orbital motion.
In the following it is assumed that the resonance is an inner vertical resonance with m ≥ 2. The Fourier components of the tidal potential can be expressed in terms of the Laplace coefficients where β = λ/d. Thus The location of the resonance is c 0000 RAS, MNRAS 000, 000-000 The orbital equations (19) and (20) yield the following linear equations for the perturbed orbital quantities: Thus and Evaluating these quantities at the resonance, one finds and therefore Here the condition for an inner vertical resonance, mΩ * = (γ + 1) 1/2Ω * has been used, as well as the fact that the tidal potential satisfies Laplace's equation, One also finds

Illustrative parameters
The following parameter values might be appropriate for a system such as IP Peg: γ = 5/3, m = 2 and q = 0.5. For the purposes of illustration, suppose that α b = 0, x = 1 and y = −7/2. Then one finds For these parameters, radiative damping provides a somewhat greater contribution to b than does viscous damping. A more non-linear and less dissipative example occurs when γ = 6/5, m = 2 and q = 1:

Physical interpretation
The non-linear ordinary differential equation (97) describes the response of the disc in the neighbourhood of the vertical resonance. All the fluid variables can be determined from the solution of this equation. A detailed study of the solutions is therefore called for.
tidal forcing mode m mode 0 mode 2m Figure 1. Mode couplings producing the cubic non-linearity. Quadratic non-linear terms in the basic equations, such as u z ∂zu z , couple together the different azimuthal mode numbers. The combination of any two modes m 1 and m 2 (shown by arrows converging on a dot) forces modes |m 1 ± m 2 | (each shown by an arrow leaving a dot). Although all modes are forced at similar levels by the tidal potential, the resonant mode m responds with anomalously large amplitude. Its non-linear self-interaction excites modes 0 and 2m at larger amplitudes than they would acquire by direct forcing, and their own interactions with mode m feed back on mode m. The net result is that mode m affects itself through a cubic non-linearity.
Equation (97) does not appear to have been studied before. It may be conveniently described as the complex non-linear Airy equation. The five terms may be interpreted as follows. The linear terms −y ′′ + xy represent the ability of the disc to support freely propagating waves in the region x < 0. In the present case these waves correspond to the compressive p e 1 mode. The point x = 0 is the location of the resonance, which is also the turning point of free waves. The non-linear term |y| 2 y derives from mode couplings. The linear term iby represents the dissipation of vertical motions by shear and bulk viscosity and by radiative damping. The term a is the tidal forcing.
It is important to realize that the cubic form of the non-linear term is not a low-order truncation of a more complicated expression, but is the exact form of the non-linearity within the consistent ordering scheme described above. The cubic term arises as follows. For simplicity, consider only the vertical velocity u z , which is proportional to z. In the equation for u z there is a non-linear term u z ∂zu z , also proportional to z, that couples together different azimuthal modes. With the scalings adopted in this paper, the vertical tidal acceleration −∂zΦt is O(ǫ 2 ). All azimuthal modes respond to this forcing with u z = O(ǫ 2 ), except the resonant mode m, which responds with u z = O(ǫ 4/3 ). The quadratic self-interaction of mode m forces modes 0 and 2m at O(ǫ 5/3 ), i.e. at lower order than the tidal forcing, and they respond with u z = O(ǫ 5/3 ). The interaction of these modes with mode m forces modes m and 3m at O(ǫ 2 ), i.e. at the same order as the tidal forcing. This modifies the leading-order resonant mode, leading to a cubic non-linearity (Fig. 1).

Special limits
When a ≪ 1, the response is y = O(a) and the cubic term may be neglected. One then has − d 2 y dx 2 + xy + iby = a.
In the absence of dissipation, b = 0 and one obtains an inhomogeneous Airy equation, which is well known in studies of resonant wave excitation in discs. The general solution is where c1 and c2 are constants, Ai and Bi are the Airy functions, and Gi is the inhomogeneous Airy function defined by, e.g., Abramowitz & Stegun (1965). As x → +∞, both Gi(x) and Ai(x) decay monotonically, while Bi(x) diverges exponentially. Therefore c2 = 0 for an acceptable solution. For x < 0, both Gi(x) and Ai(x) are wavelike. The desired solution is which has the asymptotic form as x → −∞, and represents a purely outgoing wave emitted from the resonance and travelling to the left. For any other choice of c1, the solution would contain an incoming wave component that is unphysical unless the outgoing wave is able to reach an edge of the disc, reflect and return to the site of launching with a non-negligible amplitude.
When dissipation is present so that b > 0, the emitted wave is damped. Formally, the solution may be obtained simply by replacing x by x + ib. This effectively moves the resonance off the real axis.

Angular momentum flux and tidal torque
One can multiply equation (97) by y * and take the imaginary part to deduce the conservation law In the absence of dissipation and forcing, there is a strictly conserved quantity which is proportional to the angular momentum flux of the wave. The term a Im(−y) is proportional to the tidal torque density, and the term −b|y| 2 represents dissipative losses. The total tidal torque exerted on the disc is therefore proportional to In the linear, inviscid limit (a ≪ 1, b = 0) this quantity is equal to Define the tidal torque parameter In linear theory fT = 1. When b = 0 the torque goes entirely into launching the wave. When dissipation is included, it is still true that fT = 1, because the contour integral is unchanged when the function is shifted parallel to the imaginary axis. The torque then goes partly into launching the wave and partly directly into changing the angular momentum of the disc. As the wave attenuates, its angular momentum is also transferred to the disc. In non-linear theory fT = 1 in general. It can be shown that the tidal torque, in physical terms, is equal to where all quantities are to be evaluated at the resonance λ = λ * . This is closely related to the standard formula for the torque exerted at a Lindblad resonance (Goldreich & Tremaine 1979) except that (i) vertical rather than horizontal forcing is involved, (ii) the second vertical moment of the density, rather than the surface density, appears, and (iii) the denominator involves the rate of detuning of the p-mode resonance condition rather than the epicyclic resonance. The factor fT also accounts for a non-linear correction to the torque. In the case fT = 1 and for an isothermal disc with γ = 1, this formula agrees with equation (40) of Lubow (1981).

Numerical method of solution
In the non-linear case, equation (97) must be solved numerically. In principle, it might be solved as a two-point boundary-value problem with suitable conditions at large positive and large negative x. However, no satisfactory method was found to solve the non-linear equation by this method owing to the existence of unwanted growing solutions. Instead, equation (97) has been solved by converting it into a time-dependent partial differential equation, The steady solutions of this equation are identical to the solutions of equation (97), and can be found numerically by solving the time-dependent equation as an initial-value problem starting from a suitable initial condition such as y = 0. In fact, this is not a purely formal device. It is possible to extend the analysis of this paper to permit the solution to vary slowly in time in the binary frame. This can be done by allowing the perturbations to depend on a slow time coordinate ǫ 2/3 t in addition to ξ, φ and ζ. The resulting amplitude equation is precisely equation (134), with a suitably scaled dimensionless time variable. Therefore time-dependent solutions of equation (134) that approach a steady state genuinely describe the development of a steady wave pattern in an initially undisturbed disc, and also suggest the stability of the steady solution. Equation (134) is a dispersive, non-linear wave equation related to the non-linear Schrödinger equation and the Ginzburg-Landau equation. The term xy breaks the translational invariance, providing a turning point for waves at x = 0. The term iby with b > 0 ensures that all solutions decay in the absence of the forcing term a.
Equation (134) has been solved by discretizing it in space on a regular grid of 5000 points on the interval −100 < x < 100. A second-order, centred finite difference approximation of the second derivative was used. The desired solution is evanescent for large positive x and has the form of an outgoing wave at large negative x. In order to select this solution and to prevent wave reflection as far as possible, the computational domain was extended to −110 < x < 110 with an additional critical damping effect applied in the added intervals. The boundary condition y = 0 was then applied at x = ±110. The resulting ordinary differential equations in the temporal domain were solved using a fifth-order Runge-Kutta method with adaptive step-size. Starting from the initial condition y = 0, the solution was advanced until a steady solution was obtained. The torque parameter was then evaluated by integrating the solution over −100 < x < 100.
Some illustrative solutions of the complex non-linear Airy equation are shown in Figs 2 and 3. For each solution, the estimated value of the torque parameter fT is quoted. It should be borne in mind that, for those solutions in which a wave of significant amplitude reaches the edge of the computational domain at x = −100, the integrated torque continues to fluctuate spatially at the level of a few per cent, and so the value of fT depends to this extent on exactly where the boundary is drawn.
In Fig. 2 the effect of increasing the dissipation parameter b on solutions in the linear regime a ≪ 1 is investigated. For b < ∼ 1 the most noticeable effect is the attenuation of the emitted wave. For b > ∼ 1 the emitted wave is almost completely absent and the resonant response is a broad, smooth distortion. The estimated torque parameters are consistent with the analytical expectation that fT = 1 when a ≪ 1, for any value of b.
In Fig. 3 the effect of increasing the amplitude parameter a is investigated. The most obvious effect is a broadening of the resonance, in which the first wavelength of the emitted wave becomes distorted into a much wider feature. The amplitude of the emitted wave increases initially more rapidly than a, and fT achieves values in excess of unity. For larger a, however, fT declines to values less than unity, perhaps indicating a saturation of the resonance.
When either a or b is large, the resonance is substantially broadened and the scaling assumptions used to derive the complex non-linear Airy equation become inapplicable. Under these conditions, the vertical resonance can be treated, as in Paper I, as a part of the smooth tidal distortion of the disc.

DISCUSSION
In this paper I have developed a general theory of vertical resonances, first analysed by Lubow (1981), which are an important aspect of the interaction between an accretion disc and a massive companion with a coplanar orbit. When Lubow's analysis is generalized to allow for non-linearity of the response, and dissipation by radiative damping and turbulent viscosity, the problem reduces to a universal, non-linear ordinary differential equation with two real parameters.
Numerical solutions of a time-dependent version of the complex non-linear Airy equation describe the process by which a steady and stable pattern of non-axisymmetric distortion is established in an initially undisturbed disc. For small values of the dissipation parameter b, a p-mode wave is launched at the resonance and propagates radially through the disc, carrying angular momentum. For larger values of b the wave is attenuated or not launched at all, and the resonant response is broadened into a smooth feature. The effect of increasing the amplitude parameter a is to broaden the resonance and ultimately to saturate it. The total tidal torque exerted at the resonance is independent of b in the linear regime a ≪ 1, but can be either greater or less than the linear prediction when non-linearity is important.
If the magnitude of the resonant tidal torque exceeds that of the local viscous torque, the disc may be truncated at the resonance (Lubow 1981). The viscous torque in the absence of tidal distortions is and a comparison with equation (133) shows that truncation occurs if α < αtrunc, where αtrunc = fT mΨ 2 3DΩ 2 .
However, if α is so small that the launched wave does not damp significantly in the vicinity of the resonance, the tidal torque is not transmitted to the disc locally and the process of truncation cannot occur in a straightforward manner. In the derivation of the complex non-linear Airy equation it was assumed that α = O(ǫ 2/3 ) while Ψ = O(ǫ). Therefore the process of tidal truncation is not described by that equation, but occurs in a different parameter regime. Some illustrative values of αtrunc for the m = 2 inner vertical resonance are given in Table 1, where it is assumed that fT = 1. It appears that truncation by this resonance is unlikely in cataclysmic variable discs, even in quiescence, and especially if γ ≈ 5/3. Although Lubow (1981) argued that truncation is likely, his estimates were based on the case γ = 1, for which the resonant radius is maximal. This results in a much larger torque than is obtained when γ = 5/3, because of the proximity to the companion. However, for the same reason, the resonance is less likely to lie inside the disc in the case γ = 1.
Much more probable is that a non-destructive, two-armed structure is formed, as described by the complex non-linear Airy equation. In Paper I I have argued that this effect, in the limit in which the resonant response is broadened and damped by non-linearity and dissipation, could explain the non-axisymmetric structures seen in Doppler tomograms of dwarf novae in  outburst. Indeed, estimates of the parameters for a system such as IP Peg (Section 3.9) suggest that the resonance is strongly damped and mildly non-linear. It is likely that the techniques used in this paper will be useful in other circumstances where spatial resonances occur in wave-bearing media subject to periodic forcing. The equation derived has a concise form and its terms have a clear physical interpretation, suggesting that it may be generic to some extent. In any case it provides a simple mathematical model in which the effects of non-linearity and dissipation on the resonant launching of waves can be studied.

ACKNOWLEDGMENTS
I thank Jon Dawes, Steve Lubow and Jim Pringle for helpful discussions. I acknowledge the support of the Royal Society through a University Research Fellowship.