Radiation Field

The general method of Chapter 1 is now applied to the electromagnetic field. In this chapter we disregard entirely the effect of the sources and study only the quantum-mechanical formalism of the field far removed from the sources, or the radiation field. The field equations are especially simple if they are expressed in terms of the vector potential subject to a subsidiary condition (Sections 2-1, 2-2, and 2-3). This condition introduces a complication in the quantum-mechanical formalism which was not discussed in the general theory of Chapter 1. Apart from this, the general theory can be readily applied to derive the correct expressions for the commutation rules and the momentum operators (Sections 2-4 and 2-5).

structure and the influence it has on the possibility of instability to warping.
We derive a pair of equations describing the evolution of a small tilt as a function of radius in the small amplitude regime that applies to both the diffusive and bending wave regimes. We also study the non linear vertical response of the disc numerically using an analogous one dimensional slab model. For global warps, we find that in order for the disc vertical structure to respond as a quasi uniform shift or tilt, as has been assumed in previous work, the product of the ratio of the external radiation momentum flux to the local disc mid plane pressure, where it is absorbed, with the disc aspect ratio should be significantly less than unity. Namely, this quantity should be of the order of or smaller than the ratio of the disc gas density corresponding to the layer intercepting radiation to the mid plane density, λ ≪ 1.
When this condition is not satisfied the disc surface tends to adjust so that the local normal becomes perpendicular to the radiation propagation direction. In this case dynamical quantities determined by the disc twist and warp tend to oscillate with a large characteristic period T * ∼ λ −1 T K , where T K is some 'typical' orbital period of a gas element in the disc. The possibility of warping instability then becomes significantly reduced.

INTRODUCTION
A significant number of X-ray binaries exhibit long-term periodicities on time-scales of 10-100 d. Examples are Her X-1, SS 433 and LMC X-4 see e.g. Clarkson et al. 2003. Precession and warping of a tilted accretion disc has been proposed as an explanation (Katz 1973;Petterson 1975) which has subsequently been found to have observational support (Clarkson et al. 2003). The effect of radiation pressure on a twisted tilted accretion disc was first considered by Petterson 1977a,b who noted that when radiation from the central source is absorbed at the disc surface and re-emitted, a potentially important torque will result.
Iping & Petterson 1990 later suggested that such torques determined the shape of the disc and its precession rate. Pringle 1996 subsequently showed that an initially axisymmetric thin disc could be unstable to warping as a result of interaction with an external radiation field (see e.g. Maloney et a. 1998;Ogilvie &Dubus 2001 andFoulkes et al. 2006 for later developments).
The analysis considered the disc to behave as a collection of rings, which could interact by transferring angular momentum through the action of viscosity, but otherwise behaved as if they were rigid. In particular, the vertical displacement or tilt was assumed to be essentially uniform and independent of the vertical coordinate. Thus the pressure applied at the surface at optical depth unity is assumed to be effectively communicated through the disc vertical structure so as to give a near uniform response. In this paper we extend the theoretical treatment of the interaction of a disc with an external radiation field to take account of possible significant departures of the tilt or displacement from uniformity in the vertical direction. One of our objectives is to determine the conditions under which the assumption of a uniform response is valid, and then to estimate some of the consequences when they are not satisfied, including potential additional dissipation resulting from non linear effects such as the production of shock waves. The latter is done using a one dimensional slab analogue model with the required high resolution in the vertical direction.
Although we focus on the effects of a surface pressure induced by an external radiation ⋆ E-mail:pbi20@cam.ac.uk (PBI) J.C.B.Papaloizou@damtp.cam.ac.uk (JCBP) field, very similar considerations are likely to apply when warps are induced by a surface pressure resulting from interaction with an external wind (e.g. Quillen, 2001) or surface forces produced by the interaction of the disc with an external magnetic field originating in the central star ( e.g. Pfeiffer & Lai, 2004).
We begin by considering the relevant issues using simple physical arguments. Following Pringle 1996 we consider a thin disc immersed in an external radiation field with mid plane initially coinciding with a Cartesian (x, y) plane. It is supposed that the upper surface is parallel to the external rays so that there is initially no interaction with the radiation. The upper surface is then given a vertical elevation h(r, φ), where we now use polar coordinates (r, φ). As a result, a disc element with surface area dA absorbs momentum from the radiation field at a ratė F = F 0 dA r ∂(h/r) ∂r . (1) Here F 0 is the momentum flux at radius r associated with the external radiation field. The factor in brackets gives the angle through which the local normal is rotated as a result of the perturbation (see Appendix A for more details). Assuming the absorbed momentum is reradiated isotropically above the disc by the surface layers at an optical depth of unity, there will be an applied pressure there of magnitude P = 2F 0 3 r ∂(h/r) ∂r . ( This external pressure, when applied to a complete elementary ring, produces a net torque per unit length of magnitude dT dr = 2πF 0 3 where the complex vector l has components equal to (i, 1, 0) in the Cartesian coordinate system. Here, in performing the azimuthal integration, we take into account that the azimuthal dependence of h is through a factor of the form exp(−iφ) and work with the radial amplitude from now on.
To find the consequent evolution of the disc, one requires the component of the angular momentum per unit length perpendicular to its unperturbed direction. This is given by with Ω and Σ being the near Keplerian local disc angular velocity and the disc surface density, respectively. h = dζρh/Σ where ζ is a vertical coordinate and ρ is the gas density.
For an elementary ring the condition that the rate of change of angular momentum equal the applied torque gives a tilt evolution equation of the form Assuming the vertical response is uniform, so that we can set h = h ≡ rW, we obtain a description of warp evolution equivalent to Pringle 1996 when effects due to viscosity and bending wave propagation are neglected (see sections 5.1 and 6.5 below). This description indicates the possibility of instability to radiation pressure warping.
However, here we stress the fact that the averaged elevation h , enclosed in angled brackets, applies to the bulk of the inertia of the disc and can therefore be shown to be close to the elevation of the mid plane. On the other hand h as used in the torque formula applies to the disc surface elevation. An important aspect of this paper is to distinguish these two elevations and investigate under what conditions they can be taken to be equal. This would be possible if the vertical structure responds as a rigid body to the external pressure forcing and we find the conditions for this to occur. The general requirement is found to be that the density at the surface of the disc where the pressure is applied should not be too small.
It is possible to estimate in a simple way when the vertical displacement response of the disc to the external pressure forcing becomes non uniform. Let us suppose that the external pressure P is applied at the surface layer where the density ρ = ρ * . The induced vertical displacement h, now should be considered to be also a function of the vertical coordinate ζ.
Assuming a linear response, which should be appropriate for sufficiently small elevations If the vertical extent of the disc is ζ 0 , and h varies on a radial scale comparable to r, the characteristic change in h induced over the vertical thickness is easily estimated to be where we use the approximate relation c s = Ωζ 0 .
For a uniform response, one requires that ǫδ 2 ≪ λ. This is equivalent to the requirement that the product of the ratio of the external radiation momentum flux to the local disc pressure with the disc aspect ratio should be significantly less than unity. We recall here that because of the geometrical configuration, the momentum flux locally reradiated by the disc is in general much less than the external momentum flux at that location.
In this paper we extend the treatment of disc warping induced by the action of an external pressure to include the effects of the response of the vertical structure particularly when this is non uniform as is the case when the above condition is not satisfied. In that case we find that the disc dynamics enters a different regime. This is such that the exposed surface tends to align so as to reduce the momentum absorption and hence the applied pressure. In the extreme limit of this regime the upper surface acts as if it is in contact with a rigid wall with the tendency to radiation warping instability tending to vanish. In this limit quantities characterising the disc twist and warp tend to oscillate at typical frequency ∼ λΩ.
To illustrate these effects we derive a description of the one dimensional evolution of the disc inclination in radius and time that incorporates the effects of the vertical structure response, radiation torques and which applies both to the high viscosity regime, when the evolution is diffusive, and to the low viscosity regime when the evolution is wavelike. In all cases, the efficacy of surface radiation pressure driven instabilities is found to be reduced once the model parameters are such that the vertical response is significantly non uniform.
We go on to estimate conditions for the response to be nonlinear and investigate the development of shock waves in the response using a one dimensional slab analogue model which has the same characteristic behaviour of linear perturbations as the full disc model.
We find that such shocks potentially provide an important dissipation mechanism.
The plan of this paper is as follows. In section 2 we describe some aspects of the thin disc model used. In section 3 we go on to introduce the twisted coordinate system used together with the notation convention, giving the basic equations in section 4. By integrating over the vertical direction we use these to derive a single equation governing the dynamics of a twisted disc in section 5. When the disc behaves like a set of rigid rings for which the vertical response is uniform, this equation can be used to give a complete specification of the warp evolution. We note that this provides an extension of previous formulations to be able to consider the case when warps propagate as waves rather than diffuse radially (e.g. Nelson & Papaloizou 1999).
In this Paper we use the twisted coordinate system formalism first introduced by Petterson 1977aPetterson , 1978 where the dynamical equations take the most simple form. When this formalism is adopted and an accurate description of all components of the equations of motion is needed as in the problem on hand we show that, in general, there is an ambiguity in choice of twisted coordinates with a set of these describing the same physical situation.
As discussed in section 5 a choice of a most appropriate twisted coordinate system can be motivated by the condition that perturbations of all dynamical quantities determined by the disc twist and warp are small. The transformation law between different twisted coordinates corresponding to the same physical situation is derived in appendix C for an inviscid disc.
In order to obtain a complete description of the warp evolution when the vertical response is non uniform, we begin by obtaining a complete solution of the vertical problem for a polytropic model in section 6.1. This is used to obtain a pair of equations governing the radial warp evolution. We also indicate how the results can be extended to apply to more general models. We confirm the condition for the disc response to be like that of a series of rigid rings as λ ≫ ǫ(ζ 0 /r) 2 .
We go on to perform a linear stability analysis of the radial evolution equations adopting a WKB approach in section 7 obtaining a maximum potential growth/decay rate in section 7.3. In section 8 we discuss and confirm the analogy between the response of the disc and the linear and non linear dynamics of a vertically stratified one dimensional slab. We consider the development of shock waves and the formation of a rarefied hot atmosphere in section 8.3 giving a crude estimate of the warp dissipation rate in section 8.4. Finally in section 9 we discuss and summarize our results.

A THIN DISC MODEL
To begin we briefly summarize some of the properties of an unperturbed axisymmetric thin accretion disc model which are used in the discussion below.
However, the physical basis of our results does not depend on the detailed properties of this model which has been used for the purpose of making specific calculations. In order to make our treatment as simple as possible we use a highly simplified version of the Shakura & Sunyaev 1973 α-disc model. We assume a polytropic equation of state for the disc gas where P and ρ are the gas pressure and density, r is the radial coordinate and γ is the assumed constant specific heat ratio. It is assumed that the "polytropic constant" K is, in general, a function of the radial coordinate, r, but does not depend on the local vertical coordinate ζ, see equations (14), (15) below for a definition of the coordinate system used in our study. The equilibrium variation of density and pressure with height follows from the equation of state (8) together with the equation of hydrostatic equilibrium Here c s is the adiabatic sound speed, Ω = (GM)/r 3 is the Keplerian angular velocity and M is the mass of a central object.
Integration of equation (9) gives where ρ c and P c are the mid plane values of the density and pressure, respectively and ζ 0 is the disc semi-thickness. These quantities are related to each other through As shown below, many of our equations take an especially simple form for the case γ = 5/3, accordingly we adopt this value below.
It is assumed below that the standard Navier-Stokes equations are valid as the dynamical equations of our problem with the viscosity law set up according to the standard prescription (Shakura & Sunyaev 1973 ) where ν is the kinematic viscosity. For simplicity we assume below that α is a small constant parameter and consider only leading terms in corresponding expansion series.
The radiation intercepted by the disc is considered to result in an external pressure applied to a surface near to the free surface of the disc. This surface is formally defined through the condition ρ = ρ * (τ = 1) where τ is the optical depth. It is assumed hereafter that the ratio of the density at this surface to the midplane density ρ * /ρ c ≪ 1. Thus, we neglect the effect of heating of the upper disc layers by incident radiation. Also, for simplicity it is assumed that the radiation does not affect the structure of the unperturbed disc. This can be achieved by setting the disc semi-thickness to be proportional to the radial distance where δ ≪ 1 is the disc opening angle which is constant.

COORDINATE SYSTEM AND NOTATION CONVENTION
The dynamical equations governing a twisted disc take their simplest form in the so-called twisted coordinate system first introduced by Petterson 1977a. Since this system has been discussed extensively elsewhere (for details see e.g . Petterson 1977a;1978), in this paper we give only some basic definitions and relations important for the future discussion.
In order to define the twisted coordinates let us consider a fixed Cartesian coordinate system (x 1 , x 2 , x 3 ) with origin at the central star and (x 1 , x 2 ) plane coinciding with the midplane of the unperturbed disc. We then perform two elementary rotations of the coordinate axes parametrised by the two Euler angles γ and β. The first rotation is of the (x 1 , x 2 ) axes through the angle γ keeping the x 3 axis fixed. As a result of this rotation the x 1 and x 2 axes change their directions. The second rotation is of the (x 2 , x 3 ) axes in the plane perpendicular to the new fixed direction of the x 1 axis and is through the angle β. The cylindrical coordinates (r, ψ, ζ) defined with respect to the rotated axes now determine the twisted coordinate system.
Formally, in the linear approximation that the inclination β is small, the transformation to the twisted coordinate system (r, ψ, ζ) from the initial Cartesian coordinates (x 1 , x 2 , x 3 ) is given by where B is the rotation matrix and β(t, r) and γ(t, r) are, in general, functions of r and the time t.
Instead of the original Euler angles it is useful to introduce the combinations ψ 1 = β cos γ, ψ 2 = β sin γ. Further, we define the angle φ = ψ + γ which is used to replace the angle ψ. We are thus able to consider the case β(r, t) = 0 for which the transformation (15) is degenerate. The rotation matrix giving the transformation from the twisted coordinates back to the Cartesian coordinates is given by transpose of B, B T . The cylindrical coordinates (r c , φ c , z c ) associated with the original Cartesian system are defined through: x 1 = r c cos φ c , x 2 = r c sin φ c and x 3 ≡ z c . From equations (14) and (15) we obtain and From equations (16) and (17) it follows that when ζ = 0 we have φ c = φ and r c = r. In addition, the difference φ c − φ is also first order in the small angle β. Thus we can set φ c = φ in equations for quantities that are first order in β.
All vectors and tensors appearing in our calculations are projected onto a local orthonormal basisē r ,ē ψ , andē ζ (Petterson 1977a(Petterson , 1978 where the bar is used to denote a vector quantity from now on. The explicit expressions for the above basis vectors are not important for our purposes and can be found elsewhere, see e.g. Petterson 1977a, 1978, Demianski & Ivanov 1997 hereafter DI. We note, however, that the basis can be uniquely defined by the requirement thatē ζ is the coordinate vector such thatē ζ · ∇ ≡ ∂ ∂ζ for ζ → 0. All other vectors can be found from the orthonormality conditions together with the requirement that the coordinate axes form a right handed system.

State variables
Since the coordinate lines of the twisted coordinate system are not orthogonal and both the coordinate lines and directions of the basis vectors depend on time and on the radial coordinate, differentiation of quantities is best performed with help of appropriate connection coefficients. Their explicit form can be found in Petterson 1977aPetterson , 1978 and DI. Also, we note that the projections of the velocity vectorv onto the basis vectors are not in general proportional to the time derivatives of the appropriate coordinates of a particular fluid element, see Petterson 1978. Accordingly, we give relations between these components and the time derivative of the appropriate gas element coordinates where this is explicitly needed.
Writing them in the twisted coordinate system, in the approximation that the inclination angle is small enough that a linearization procedure in which quadratic and higher powers of β can be neglected, the Navier-Stokes equations and the continuity equation then lead to: 1) a set of equations for unperturbed quantities Q 0 formally coinciding with the equations describing them for flat disc accretion.
2) A set of equations for quantities Q 1 , which are non axisymmetric perturbations to Q 0 induced by deformation of the disc shape. The dependence of the latter quantities on φ is harmonic: It is convenient to represent these perturbations as well as the angles β and γ using the complex quantities, Q and W, which are defined through Here we note that we adopt the convention that quantities in bold represent complex perturbations. This convention applies also to the components of a vector with vectors themselves denoted with an over bar. When real as in the perturbed equations of motion (23)-(29) given below, these are not in bold. Note also that for convenience we use ρ 1 and ξ ζ to denote the complex density perturbation and vertical component of the Lagrangian displacement below, and these are not in bold.

BASIC EQUATIONS
As we mentioned above we write each of the variables entering the continuity and the Navier-Stokes equations as the sum of an unperturbed part, formally coinciding with a solution of the equations for axisymmetric thin disc accretion, and a perturbation. Accordingly, we have for the components of velocity and the gas density respectively. But as the only non zero component of the unperturbed velocity is in the azimuthal direction and corresponds to Keplerian rotation which we deal with explicitly, we shall drop the subscript 1 from the velocity perturbation.
The unperturbed gas density ρ 0 is given by equation (10). We assume that only circular Keplerian motion of the gas is present in the unperturbed state and accordingly have v φ 0 = 1 Strictly speaking, in order to satisfy boundary conditions at the disc surface we can assume a harmonic dependence of some of these quantities only over a limited domain of the angle φ, see Appendix A for details. However, as follows from the results provided there, we can formally proceed assuming such quantities depend harmonically on φ over the appropriate domain and then extend to the whole φ domain using the rule discussed in Appendix A. rΩ = GM r . Thus we neglect the radial component of the unperturbed velocity v r 0 present in the thin disc mainly due to action of viscosity. This can be done in the limit of a sufficiently small value of α < 1.
Apart from two significant modifications, our equations for the perturbed quantities follow from the corresponding equations given in DI provided that terms involving quadratic and higher order corrections in the assumed small parameter α are discarded and that the gravitomagnetic force in the ζ component as well as relativistic post-Newtonian corrections to the r and φ-components of the Navier-Stokes equations are neglected.
The most important modification stems from the fact that in the present study an accurate description of vertical motions and displacements is needed and therefore, we take into account all contributions to the vertical velocity, v ζ . These consist of the contribution determined by time-dependent evolution of the basis vectorē ζ and the contribution determined by time dependent vertical displacements of gas elements in the disc. Accordingly, we have where U =β sin ψ −βγ cos ψ =Ψ 1 sin φ−Ψ 2 cos φ and the dot stands for the time derivative.
Only the first term in (21), which arises from the time dependence of the changing coordinate system, is normally taken into account in studies of twisted discs that use the formalism of the twisted coordinate system, see Hatchett at al 1981 and DI. However, the role of the second term, which describes the disc motion relative to the moving coordinates, is essential for our purposes.
We also take into account the time derivative of the perturbed density ρ 1 in the continuity equation and the time derivative of v ζ in the ζ-component of the Navier-Stokes equation.
These terms are proportional to a small parameter 1/(Ωt ev ) where t ev is a characteristic evolution time of the twisted disc. That is for any perturbation quantity Q, we can write |Q/Q| = O(1/(Ωt ev )). However, even though they are small, in a formal sense, when compared with leading terms in a series expansion, when retained, they allow us to regularise some otherwise formally singular expressions, that arise when dealing with the particular problem of calculating the vertical displacements in the disc, see e.g. equation (50) below.
The perturbed Navier-Stokes equations can be divided into 'horizontal' part incorporating the (r) and (φ)-components and 'vertical' part corresponding to the (ζ) component. The horizontal part follows from the equations provided in DI when the relativistic corrections are neglected and the limit of small viscosity is adopted. The (r)-component can be written where is the (r, ζ)-component of the viscosity tensor. In the expression for W, a prime denotes differentiation with respect to r.
Differentiating the (φ)-component of perturbed Navier-Stokes equations with respect to (φ) we obtain where The set of equations (21) and (24) can be used in order to express the perturbed velocity components in terms of the quantity W. As seen from equations (21) and (24), to the leading order in small parameters α ( or equivalently η) and 1/(Ωt ev ), the perturbed velocity components enter in both equations in the same combination Therefore, these equations are degenerate to leading order and accordingly the next order corrections have to be retained (Papaloizou & Pringle 1983, hereafter PP). This leads to an inverse dependence of the 'horizontal' part of the perturbed velocity on the small parameters such that both v r and v φ are ∝ min(α −1 , 1/(Ωt ev )) (see PP for details).
The (ζ)-component of the Navier-Stokes equation is written in the form where ρ 1 and P 1 are the perturbations of density and pressure, respectively. These are related through P 1 = c 2 s ρ 1 with the sound speed given by equation (10). The continuity equation is formally equivalent to perturbed continuity equation written for the case of a flat disc in cylindrical coordinate system We recall that the relation between the velocity component v ζ and the corresponding component of the Lagrangian displacement vectorξ is given by equation (22).
Apart from the terms proportional to W and U equations (23-29) are formally equivalent to the corresponding equations that would be written down directly in a fixed cylindrical coordinate system (r, φ, ζ) (PP).
The term proportional to W in equation (23) is essentially a projection of the pressure gradient. Due to the presence of perturbations with odd symmetry with respect to ζ in twisted discs, surfaces of constant pressure do not in general coincide with the orbital planes of gas elements. Therefore, there is is a projection of the pressure gradient onto such an orbital plane which is described by this term, see Fig The term proportional to 2 ∂ ∂φ U in equation (28) is made up from two contributions. The first is related to the non-holonomic character of our basis vectors. As we have mentioned above in a non-stationary situation, this leads to a non-zero vertical velocity even for a gas element being at rest with respect to the twisted coordinate system, see equation (21) above. This fact was noted for the first time by Hatchett, Begelman & Sarazin 1981. The second contribution is determined by non-inertial effects. These are taken into account by an appropriate connection coefficient in the formalism developed by Petterson 1978. Both contributions have exactly the same form. This accounts for the factor 2 in the expression 2 ∂ ∂φ U in equation (28). The set of equations (23-29) can be brought into a simpler form with help of the complex notation introduced through equation (19) for the density perturbation and the perturbations to the velocity components. Also, we assume hereafter that all perturbed quantities depend on φ and time only through an exponential factor. Thus where the frequency ω is, in general, a complex constant.
Subtracting equation (25) from equation (23) and introducing the complex notation it is easy to see that the result contains the components of the velocity perturbation only in the combination v + = v r + 2iv φ . The resulting equation can thus be represented as where we use equations (24) and (26) and set ρ ≡ ρ 0 from now on. The quantities v r and v φ can be separately expressed in terms of v + taking into account the fact that in the leading approximation in the small parameters we can equate expression (27) to zero 2 and, Now we can express v r and v φ in terms of W with help of equations (31) and (32). To do this we remark that the ζ dependence is dealt with by noting that η ∝ P and the solution is such that both v r and v φ are ∝ ζ. We then use the resulting expressions to eliminate these velocity components in the continuity equation (29) so obtaining where and Finally, equations (28) and (22) can be brought in the form and

DYNAMICS OF A TWISTED DISC
At first let us point out that the hydrodynamical equations as written in a twisted coordinate system cannot in principle provide a single dynamical equation for W without some additional considerations. Indeed, since the Euler angles β and γ are considered as dynamical variables in this formalism, the total number of dynamical variables is always larger by two than the number of equations actually needed to determine the dynamics of the system.
2 Setting (27) to zero physically reflects the fact that trajectories of gas elements in the disc may be approximated as Keplerian ellipses in the leading order. Next order correction determine parameters of the these ellipses and effect of slow precession of their main axes on time scale of order of tev , see DI and Ivanov & Illarionov 1997 for more details on physical picture of the horizontal motions in the disc. Formally, parameters and precession rate of these ellipses can be obtained from equation (31).
For example, in the case of a baratropic equation of state P = P (ρ), there are six dynamical variables used in the description of the problem, these being the three components of velocity, and the density together with the two Euler angles. However the number of dynamical equations available is four, these being the three components of the Navier-Stokes equations, and the continuity equation. Since the number of equations is smaller than the number of variables, there are in principal, different sets of variablesv(t, r), ρ 1 (t, r) and W(t, r) describing the same dynamical system and connected to each other by certain transformation laws. These transformation laws can be easily found from the condition that the velocity and density fields measured in some inertial coordinate system remain unchanged when the set of dynamical variables associated with the twisted coordinate system is transformed. Explicit expressions for these transformation laws are given in Appendix C.
It is interesting to note that the dynamical variables in the inertial coordinate system play the role of gauge independent quantities and the laws of transformation between different sets of variables in the twisted coordinate system may be regarded as gauge transformations.
This situation has many analogies in different physical systems, say in systems governed by General Relativity, such as e.g. , the theory of small cosmological perturbations, see e.g. Landau & Lifshitz 1975. In order to begin the process of obtaining a single equation describing the dynamics of a twisted disc referred hereafter as a 'dynamical equation' we substitute the density perturbation given by equation (33) where we use the hydrostatic balance equation (9) and approximate Ω + ω ≈ Ω in the first term on the right hand side. Equation (38) plays an important role in our analysis. A dynamical equation follows from (38) provided that appropriate boundary conditions at the disc surface are specified.

A simple approach to the derivation of a dynamical equation
Before discussing an approach to the problem which fully takes into account the vertical structure of the disc, we would like to show how to obtain a governing dynamical equation from equation (38) in a straightforward way. Although this approach is incomplete it allows us to obtain a dynamical equation which is approximately correct provided that the ratio of the density where the external pressure is applied to the central density, λ = ρ * /ρ c , is not too small.
To do this we integrate equation (38) over ζ taking into account the fact that the last term on the right hand side vanishes when ζ → ±ζ 0 and can be discarded. We thus obtain where Σ = +ζ 0 −ζ 0 dξρ 0 is the surface density and a typical disc height H is defined by condition H 2 = +ζ 0 −ζ 0 dζζ 2 ρ 0 /Σ. It can be easily shown that the surface terms on the left hand side are, in the low frequency limit, proportional to the Lagrangian pressure perturbation, ∆P = (dP 0 /dζ)ξ ζ + P 1 . Introducing the complex notation and taking into account that in the low frequency limit, we can write from (37) where we make use of equation (33) and use the hydrostatic balance condition (9).
The Lagrangian pressure perturbation must be equal to the radiation pressure at the surface corresponding to an optical thickness expected to be of order unity. As we discuss in Appendix A we can assume that the radiation pressure acts only on the upper surface of the disc ( i.e. where ζ → +ζ 0 ). We set, accordingly, the surface term corresponding to the lower surface of the disc to zero, and express the term corresponding to the upper surface in equation (39) with help of equation (40) through the radiation pressure term F + given by equation (A11) of Appendix A to obtain where we neglect the second term of the left hand side of (39) proportional to ω in order to get (41). This procedure is justified below, see Section 6.4 where we derive a similar equation in a more rigorous way.
Also, let us stress that the (ζ)-component of the displacement vector is assumed to be evaluated at the surface of the unperturbed disc where the optical thickness τ is of the order of unity, and let us recall that F 0 is the radiation momentum flux per unit area perpendicular to the radial direction, see Pringle 1996.
The density ρ * (τ ≈ 1) where the external radiation pressure is applied, is assumed to be always much smaller than the mid plane density ρ c .
We comment that in the absence of radiation pressure, equation (41) describes the linear dynamics of a free warped disc. When α = 0, this is wavelike, the waves being non dispersive and having local speed ΩH/2, see Nelson & Papaloizou 1999, Papaloizou & Lin 1995, Nelson & Papaloizou 2000. Furthermore we comment, without giving details, that equation (41) may also be derived following the approach discussed in those papers, so confirming the analysis based on the adoption of a twisted coordinate system presented here.
The radiation pressure term involves the combination ξ ζ /r + iW. This can be regarded as being the angular displacement of the surface of the disc. It is composed of two parts.
The first ξ ζ /r represents the displacement of the surface of the disc relative to the midplane.
The second term iW can be viewed as representing the displacement of the midplane and is described by the twisted coordinate system. When the disc behaves like a rigid ring, the contribution of ξ ζ is small.

The small parameter
The small parameter λ = ρ * /ρ c should be compared with other small parameters of the problem in order to specify the different possible regimes of dynamical evolution of the disc.
In particular, as we show below it is important to compare it with the local growth rate, in units of Ω, associated with any possible instability of the disc that is driven by radiation pressure,ω r .
As mentioned in the introduction and elaborated further below, we will see, when λ >ω r , we can neglect ξ ζ in the expression ξ ζ /r + iW entering in (41). Physically this corresponds a situation when the surfaces of constant density in the disc are approximately parallel to each other and remain at an approximately constant distance apart. In that case, for the purpose of calculation of the radiation pressure term, the disc may be considered as being composed of rigid rings having different inclinations and precession angles and the analysis undertaken in previous studies (e.g. Pringle 1996) applies. Neglecting the displacement ξ ζ in (41), we see that this equation gives the evolution of W without the need for additional considerations. For the case α > δ the frequency ω may be neglected in the expression ω = ω + iαΩ and equation is reduced to a form analogous to that obtained by Pringle 1996.
In this case Pringle 1996 indicated that the disc may be unstable in linear theory provided that the radiation pressure term is sufficiently large.
In the opposite limit δ > α, equation (41) describes wave-like propagation of twisted disturbances of the disc (see Nelson & Papaloizou 1999, Papaloizou & Lin 1995, Nelson & Papaloizou 2000 as indicated above, modified by the presence of the radiation pressure term.

The polytropic model and the dynamical equation in a 'standard' form
For the polytropic equation of state we can bring equation (41) to a further reduced form.
In this case one can obtain from equation (10) where Γ(x) is the gamma function.
Before derivation of the general dynamical equation let us now temporary assume that the disc parameters are such that λ > δ 2 /α and the term proportional to ξ ζ in (41) can be neglected and obtain an equation analogous to what has been obtained by Pringle (1996) for the polytropic model. In this case we can use equation (42) and rewrite equation (41) in where and introduce the dimensionless radiation pressure parameter defined as the ratio of the radiation momentum flux F 0 to a characteristic gas pressure in the mid plane of the disc. In general, the conditionF 0 ≪ 1 should be fulfilled for a simple description of the evolution of a twisted disc to be valid. Let us note that when γ = 5/3 the numerical factor in the second term on the left hand side of (43) is equal to 8 9π and the numerical factor in the last term is equal to 1 24 . Equation (43) describes the 'standard' dynamics of a twisted disc immersed in a radiation field.

THE SELF-CONSISTENT APPROACH
Now let us take into account the contribution of all the terms in equation (41). In order to do that we need to relate the term in this equation involving the (ζ) component of the Lagrangian displacement vector to W. Also, the term 2Ωω +ζ −ζ ρv ζ dζ was discarded without justification in order to go from equation (39) to equation (41). We provide a formal justification for this below. In order to find the vertical component of velocity, and accordingly, the vertical displacement we solve equation (38) (41) to leading order in our small parameters α and 1/(Ωt ev ) and we also specify the form of the vertical displacement at the surface where the external pressure applies.

Solution of the vertical problem for the polytropic model
In order to find a solution of (38) it is convenient to separate the vertical velocity v ζ into two parts such that v ζ = v ζ 0 + v ζ 1 . Here the first part of the decomposition, v ζ 0 , is a particular solution of (38) with boundary conditions at the disc surfaces which apply when they are free and which are regular in the limit of vanishing density. They are such that in that limit, the Lagrangian perturbation of the pressure given by equation (40) should approach zero at the disc surfaces ζ = ±ζ 0 .
The second part of the decomposition, v ζ 1 , is a solution of the homogeneous part of equation (38): This solution will be used in order to satisfy the boundary conditions where the Lagrangian perturbation of pressure is given by equation (40). We comment that expressions such as ζ → ±ζ 0 are used below to apply to the surface where ρ = ρ * , the latter being the small density where the external pressure is applied, rather than at strictly zero density.
Firstly we derive an expression for the quantity v ζ 0 . To do this we transform the right hand side of equation (38) with help of equations (10) and (11) to the form where x = ζ/ζ 0 and N = r −3/2 (dW/dr)/ω. Note that we do not assume that the disc opening angle is constant when deriving the relations (47) and (48) and these are valid for any dependence of ζ 0 on r.
The calculation of v ζ 0 is greatly facilitated by noting that the source term S/ρ has a quadratic dependence on x. Is is also easy to see by direct substitution that a solution for the vertical velocity v ζ 0 with the same quadratic dependence on x can be readily found in the form: where expressions for the quantities b 1 and b 2 directly follow from equation (38) and the form of the source term (47) b and where we introduce the dimensionless frequencyω = ω/Ω and assume that it is small: |ω| ≪ 1. In this limit it is easy to see that the quantity b 2 is ∼ω times smaller than b 1 . It is, therefore, neglected later on, and we assume that in our future analysis.

The solution of the homogeneous problem
Now let us find the homogeneous component v ζ 1 . It turns out that the solutions of equation (45) take on a very simple form when γ = 5/3. In that case the solutions can be expressed in terms of elementary functions. Therefore, we shall specialise to this case in the main text of the paper from now on. However, the final description of the system we obtain has wider applicability (see below). The form of the solutions of the homogeneous problem corresponding to general values of γ, which can be used to construct the full solution in that case is relegated to Appendix B.
When γ = 5/3 it is convenient to introduce new variables θ and Y according to the It is obvious from equation (53) that θ = 0 and θ = π/2 correspond to the disc midplane ζ = 0 and to the free surface ζ = ζ 0 , respectively. It follows from equation (10) that the density is proportional to cos 3 θ and accordingly the variable Y is proportional to ζcomponent of the gas momentum per unit volume.
In terms of the variables defined through (53), equation (45) takes the form It can be verified by direct substitution that the two independent solutions of (54) are where κ = √ 4 + 6ω. Note that the solutions (55) are valid for any value ofω. For our purposes, however, only the case |ω| ≪ 1 is significant and we accordingly set κ = 2 + 3 2ω and take into account only zeroth and first order terms when summing series in ascending powers of ω 3 .
Thus the general solution of equation (38) can be expressed as a linear combination of these independent solutions in the form where v ζ 0 is given by equation (49) with γ = 5/3 and the arbitrary constants C 1 and C 2 can be chosen to satisfy the boundary conditions (46).

A specific gauge freedom associated with the twisted coordinate system
Following the simple approach of section 5.1, a dynamical equation for the quantity W can be obtained from equation (38) by integration over ζ provided that certain terms are discarded. However, in a more exact approach, it is easy to see that this equation cannot be used for determining the evolution of W without invoking some additional considerations.
As we have mentioned above this situation comes about because there is some degree of 3 Contrary to the inhomogeneous part v ζ 0 the next order terms inω in the homogeneous part diverge near the disc surface. They can be of order of or larger than the leading terms and must, therefore, be retained.
arbitrariness associated with the specification of the twisted coordinate system used in our study.
The resulting gauge freedom may be used in order to put constraints on some dynamical variables by choosing the most appropriate gauge from physical point of view. In our case it seems reasonable to look for a gauge where the vertical component of velocity v ζ has a small absolute value. This could be specified by imposing the requirement that v ζ (ζ = 0) = 0.
As we will see below the condition (57) together with equation (38) fully specify all dynamical variables and allow us to obtain a single equation determining the dynamical evolution of W. 4

Derivation of the dynamical equation
In order to find the dynamical equation we calculate the coefficients C 1 and C 2 entering (56) with help of equations (40), (46) and (A11). Then application of the gauge condition (57) will give us the dynamical equation.
We note that only the homogeneous part v ζ 1 of the vertical velocity can lead to a non vanishing Lagrangian pressure perturbation in the limit of zero density when ζ → ±ζ 0 . Thus this component of the solution dominates near the surface. Using equations (40), (46) and (A11), in the limit of small surface density, ρ * , we get where we assume that all quantities are evaluated at ρ = ρ * .
Differentiating expression (56), using (55) remembering that ζ = ζ 0 sin θ and considering the limit ζ → ±ζ 0 and accordingly θ → ±π/2 we approximately have The second condition in (58) tells that C 2 = − 3 4 πωC 1 . Using this fact we obtain from equations (52), (55), (56) and (59) in the same limit 4 Note that the condition (57) is not unique. Another reasonable condition would be eg. the requirement that We expect, however, that any reasonable condition used to fix the gauge would give the same dynamical equation to leading order in the small parameter 1/(Ωtev ). and where we recall that λ ≡ ρ * /ρ c = cos(θ) such that the coordinate θ corresponds to the density level where ρ = ρ * .
Using equations (10) and (11) with γ = 5/3 we can represent the factor in the front of the radial derivative on the right hand side of equation (58) in the form Finally, substituting equations (60) and (61) in equation (58) and taking into account (62) we obtain an equation for the quantity C 1 ǫδ 2 Ωr 2 d dr where we introduce an important parameter This parameter determines strength of dynamical effects caused by the radiation pressure acting on the disc. When ǫ > 1 these effects may be significant. Physically, the inequality ǫ > 1 means that the ratio of the momentum flux carried by radiation, F 0 , to the characteristic mid plane gas pressure in the disc is larger than the disc opening angle. Now let us consider the gauge condition (57). From this condition and the expression (49) it follows that where we use the variable Y 1 defined in equations (53) and (55). Taking into account that Y 1 (θ = 0) = 4 and, accordingly b 1 + 4C 1 = 0 in equation (63) and using the explicit form for the quantity b 1 given in equation (50) with γ = 5/3 we obtain Equations (63) and (66) form a complete system for describing the disc dynamics. Substituting C 1 given by equation (66) into equation (63) we can obtain a single third order equation for the quantity W. However, we find it more convenient to retain the pair of equations obtained directly from (63) and (66) to describe the disc dynamics.

Description of the dynamical system
In order to bring this pair of equations into a simpler form, we introduce a new dimensionless variable x = 4ω Ωr C 1 . Using x to replace C 1 in equations (63) and (66) we obtain and where we have used the explicit expression for the quantity N defined through equation (48) and recall thatω = ω + iαΩ = Ω(ω + iα).
It is easy to see that equation (41) follows from form equations (67) and (68) provided that the (ζ)-component of the Lagrangian displacement vector in equation (41) is replaced by In order to be able to neglect this term in equation (67), we require that |x| ≪ |λW|.
When this is satisfied, x can be determined after neglecting the first term in this equation, and assuming that W varies on a length scale comparable to r, with the result that x ∼ ǫδ 2 .
Thus a description assuming that the vertical response of the disc is rigid can be recovered only when the disc parameters are such that λ ≫ ǫδ 2 . This is precisely the condition we obtained from a simple argument given in the introduction.
In the opposite limit the contribution proportional to d dr (λ −1 x) should not be neglected in equation (67) when the effects of the radiation pressure are significant. Then as λ → 0, this equation gives Using the above to eliminate x in equation (68) we obtain an equation describing stable radially propagating disturbances which attain a frequency given byω = 4 3π λ in the limit of zero radial wave number (see also the discussion below).
We comment here that although we used a polytropic model to derive equations (67) and (68) and continue with the immediate discussion, the indicated equivalent description of the dynamical system, which uses Σ and H with the numerical factor 3π/4 replaced by 2Σ/(ρ c ζ 0 ), contains no explicit dependence on the local details of the model, other than through the product of the height at which the external pressure is applied and the value of the density there, Σ and H. This shows that our formalism has a more general applicability and we have verified this to be the case by obtaining the equivalent representation, quite generally, using an alternative analysis that does not use twisted coordinate systems, see e.g. Nelson & Papaloizou 1999, Papaloizou & Lin 1995Nelson & Papaloizou 2000.

Local analysis
In general, equations (67) and (68) should be solved numerically, for a specified disc model. However, this is beyond the scope of the present paper. Here we follow Pringle 1996 and adopt the usual local approximation scheme in which it is assumed that the radial dependence of all variables is ∝ e ikr and that the radial wave-number k is such that |kr| ≫ 1. Note however, that although details must depend on global issues and boundary conditions so that one has to proceed with caution, we may also expect that this analysis gives reasonable order of magnitude estimates even when |kr| ∼ 1.
In addition, for simplicity we shall also continue to consider the polytropic model as the essential physics is contained therein, while bearing mind that results, when expressed only in terms of the local disc model parameters, Σ, H, and ρ * ζ 0 , can be straightforwardly applied to other models on the basis of the discussion given in section 6.5 above.
Setting x, W ∝ e ikr in equations (67) and (68) we obtain 1 + 2i 3 and respectively, wherek = kr is the dimensionless wave number. One can use equations (71) and (72) to obtain dispersion relationω =ω(k). However, as we have to consider several limiting cases for which there is qualitatively different dynamics, the general form has limited usefulness.

Dynamics of a twisted disc without external radiation pressure: Viscous and wave regimes
At first let us recall the basic dynamical properties of the standard twisted disc where radiation pressure effects are absent. In this case the dispersion relation can be easily obtained from equation (72) together with the condition x = 0. We then obtaiñ The behaviour of the disc is qualitatively different in the two limiting cases of sufficiently large and sufficiently small α parameter (PP, Papaloizou & Lin 1995). For the case of a large α the 'slow' mode corresponding to the the positive imaginary part in (73) determines the disc dynamics and the dispersions relation can be approximately written as Remembering that according to our convention all quantities are proportional to e −iωt ≈ e −ων t we conclude that in this case, twisted perturbations of the disc decay. In the limit of small α it is easy to see that twisted perturbations propagate over the disc with a speed of the order of a typical sound speed with little dissipation and we havẽ ω ≈ ±ω s − iα/2,ω s = δ|k| 2 √ 6 .
From the conditionω ν <ω s one can see that the 'viscous' regime of evolution of the disc is realised when For the case of a 'global' perturbation with |k| ∼ 1 the above condition is approximately reduced to requirement that the viscosity parameter α is larger than the disc opening angle δ.

The effects of radiation pressure for relatively large α
Now let us consider effects induced by external radiation pressure. At first let us assume that the condition (76) is fulfilled. In this case we have from equation (72) that approximately Substituting equation (77) in (71) we obtaiñ where we now define the time scale characterising influence of the radiation pressure,ω −1 r , such that The dimensionless time scaleω −1 r is that associated with the radiation pressure driven instability considered by Pringle 1996.
It is clear from equation (78) that the character of the dynamical evolution in fact depends on a comparison ofω r to λ = ρ * /ρ c . Asω r = ǫδ 2k , this leads directly to the condition ǫδ 2 ≪ λ being required for the assumption that the vertical response of the disc is that of a rigid body to be valid for disturbances with radial scale comparable to r. That condition was also derived on general grounds in section 6.5 and from simple arguments in the introduction and is explored further below.
First let us assume the caseω r ≪ λ. In this limit the factor in the front of the first brackets in (78) is equal to one and because by necessity λ < 1, the first term in the brackets can be neglected, and thereforeω is purely imaginary such that This gives either growth or decay of the perturbation depending on whether the sign of the expression in the brackets is positive or negative. Sinceω r is proportional tok it can be, in principal, positive whenk > 0. Thus one can infer instability of the disc as has been done by Pringle 1996. The quantity λ drops out and the disc behaves as a collection of rigid rings.
The corresponding stability criterion can be obtained from the condition that the absolute value of the first term in the brackets is larger than the second term and therefore Note that the numerical factor, ∼ 0.1, in the expression for ǫ crit depends on the vertical structure of the disc. It seems reasonable to assume, however, that it is always smaller than one for any possible vertical distribution of pressure and density.
Let let us now consider the limitω r ≫ λ. In this limit the character of the disc evolution changes drastically. This is essentially because the vertical motion can no longer be treated as uniform as would occur for a series of rigid rings. Instead of potential growth of disc perturbations, we find mainly decaying oscillatory dynamics. The corresponding dispersion relation reads As seen from (82) the real part of the frequency does not depend on the magnitude of the radiation pressure. The imaginary part can be either negative or positive. The values of the terms in the brackets, being always negative in the limit of very small λ or very largeω r and thus very large F 0 , lead to damping in those limits. In any event we might also expect that the oscillations in the disc inclination could result in strong self-shielding of the disc which can inhibit the growth of the inclination angle.

The effects of radiation pressure for small α
Let us now consider the opposite case of a small viscosity parameter α such that the inequality (76) is not satisfied. In this case, provided that the value of radiation flux is smaller than the typical mid plane pressure of the gas, a simple analysis shows that the effect of radiation pressure is to give a correction to the dispersion relation (75).
We begin by pointing out that equation (72) in this limit can be brought to the form where the last term in the bracket is a small correction. After substitution of (83) in equation (71) we can look for solution for the dispersion relation in the form;ω = ±ω s − iα +ω 1 . It is easy to see that the correctionω 1 has the form It is clear from (84) than we always have |ω 1 | < 4 9πω r . In addition we note that |ω r /ω s | ∼ ǫδ.
As follows from the definition of the parameter ǫ (see equation (64)), the right hand side of (85) is smaller than unity provided that the radiation flux is smaller than the disc midplane pressure. Assuming that this is valid the disc dynamics, in the case of the small α, is mainly determined by propagation of waves with a speed of the order of the sound speed as in the case without radiation pressure. As in the case of very small λ considered above, we again get damping in that limit. It also seems reasonable to suppose that in this situation, fast oscillations of the inclination angle would lead to a strong self-shielding of the disc thus inhibiting possible instabilities.
7.2 Physical character of the disc dynamics in the regime where λ ≪ω r As discussed above, when λ ≪ω r the character of the disc dynamics differs drastically from that of a set of rigid rings. The disc then tends to oscillate at the frequencyω ≈ 4 3π λ, see equation (82) and the discussion in section 6.5. Since this regime has not been discussed elsewhere we briefly describe its main features.
Assuming the disc to be in this regime, as indicated in section 6.5, we can neglect unity in the brackets on the right hand side of equation (71) and obtain a simple relation between x and W such that Now we use equation (86) and equation (69) together with equation (A11) of Appendix A to see that the radiation pressure term F + → 0 when λ/ω r → 0.
Therefore, in this limit the disc tends to evolve in such a way that its surface becomes almost parallel to the initial or radial direction thus strongly decreasing the amount of any increase in the amount of intercepted radiation as a result of the perturbation. The disc gas does not cross this surface and therefore, in this limit, the radiation pressure effectively provides an impenetrable wall at the density ρ * thus prohibiting the gas motion in vertical direction. However, when there is no influence of the radiation pressure on the lower disc boundary, as mentioned already in section 6.5, stable radially propagating warp like disturbances can still still exist.

Maximal potential growth/decay rate
An interesting consequence of the effects described above is that, when equation (78) is considered, there is an extremum of the imaginary part ofω,ω I as a function ofω r . The extremum is attained whenω r =ω e r = ± 3 2 λ and the corresponding value ofω I is equal tõ Thus for any value of the external radiation pressure and corresponding positive/negative ω r corresponding to the instability/decay of the disc perturbations, the growth/decay rate cannot exceed the value given by equation (87) with the signs (+) and (−) corresponding to instability and decay, respectively.
When the effect of the external radiation pressure is significant, the second term on the right hand side of (87) can be neglected and we have a very simple expression for the maximal absolute value of the potential growth/decay rate

Conditions on disc model parameters for potential instability
As indicated above, the disc dynamics in the presence of the radiation field is more complex in the case of a relatively large α for which inequality (76)  According a the linear theory of twisted disc perturbations that makes use of a local WKB analysis, instability associated with the effect of radiation pressure appears to be unaffected by the disc density ρ * at the location it is applied when |ω r | ≪ 3 2 λ. When this condition is not satisfied, the disc dynamics is strongly affected by effects associated with the non-rigid body character of the response of the disc vertical structure and the potential instability may be either inhibited or absent.
The above inequality leads to an upper limit for the value ǫ and, accordingly, on the value of the radiation momentum flux F 0 On the other hand the quantity ǫ should be larger than ǫ crit given by equation (81). Therefore, instability seems to be possible only for a certain range of values of the radiation momentum flux. From the condition ǫ crit < ǫ * we find a condition on the disc parameters for the instability to be possible in the form λα δ 2 > π 32 |k| ≈ 0.1|k|.
Assuming that this estimate is approximately valid even when |k| ∼ 1 we obtain a lower bound for the value of the density ratio λ from (90) where δ (−2) = δ/10 −2 and α (−1) = α/10 −1 . When the density ratio is very small and such that λ < λ crit , the simple theory indicates that the disc is likely to be stable for the whole range of possible values of the radiation flux F 0 .

Some conditions for non-linearity
We can obtain a value F 0 = F * of the radiation momentum flux corresponding to the characteristic value ǫ * separating two possible regimes of evolution of the disc in the presence of radiation pressure. By making use of equation (64) we find where P c = ρ c Ω 2 ξ 2 0 estimates the gas pressure in the disc mid plane. As mentioned above, in a situation where flaring of the disc is present, the magnitude of the radiation momentum flux should be smaller than P c in order for our simple theory to be approximately correct.
This sets an upper bound for the value of λ, such that λ < δ. Assuming that the radiation momentum flux, F 0 , is near to F * , the expression (92) leads to a stringent constraint on the value of the inclination angle β. In order for the linear theory to be valid, a characteristic value of the radiation pressure in the disc F ∼ F 0 r(dβ/dr) should be smaller than the value of the gas pressure at the surface where the radiation pressure is applied, P * ≡ P (ρ * ) ≈ P c λ 5/3 for the polytropic model with γ = 5/3. From the condition F < P * , assuming that β changes on scale comparable to the radius, we obtain a condition for linear theory to be valid in the form β < β * = λ 2/3 δ ≪ 1.
Note, however, that this constraint has a significant dependence on the vertical structure of the disc. Thus, for an isothermal disc where P ∝ ρ the small parameter λ 2/3 in (93) would be absent.
It is also interesting to note that for a polytropic disc, similar considerations lead to a strong constraint on the value of β when the radiation flux F 0 is of the order of F crit corresponding to onset of a radiation pressure induced instability in the disc. In this case we can repeat the above argument but using the value ǫ crit given by equation (81)  Proceeding this we find that linear theory would be valid only if β < β crit = 10λ 5/3 α/δ.
Substituting typical values δ = 10 −2 and α = 10 −1 and assuming that λ ∼ 10 −3 we get β crit = 10 −3 . Thus, for unperturbed disc models with strong density and temperature drops toward the surface boundary of the disc, the linear theory of perturbations is, strictly speaking, valid only for quite small values of the inclination angle.

The one dimensional vertically stratified slab analogue
As we discuss in more detail in the next section there is an analogy between the motion in this regime and that occurring in the simple hydrodynamic problem of free one-dimensional oscillations of a vertically stratified column of gas, with equation of state (8), gravitational acceleration g = −Ω 2 ζ, lower free boundary, and an upper rigid boundary at some density ρ * ≪ ρ 0 . The unperturbed distributions of pressure and density are given by equation (10).
Assuming that perturbed quantities depend on time through a factor ∝ e −i(Ω+ω)t the free oscillations of the gas column are described by equation (45) and therefore expressions (55) and (56) with the response v ζ 0 = 0 provide the solutions of the free problem for a polytrope with γ = 5/3. As above we shall focus on that case below.
If we assume that such a column is rotating around a central body with angular frequency Ω and m = 1, the frequency of oscillations observed in a non rotating frame will be equal to ω. As we will see later the smallest eigen value ω, for the normal mode oscillations of this problem, is precisely equal to 4λ 3π Ω -the frequency of the disc oscillations, when the radial wavelength is very long, and it is in the regime λ ≪ω r .

LINEAR AND NON LINEAR DYNAMICS OF THE VERTICALLY STRATIFIED SLAB
When the conditionω r ≫ λ is fulfilled, internal non uniform vertical motions determine the disc dynamics. As we have mentioned in the previous section, the dynamics of the disc in the long radial wavelength limit is quite similar to the dynamics of a perturbed vertically stratified one-dimensional gas column with equation of state given by (8) and gravitational acceleration g = −Ω 2 ζ.
Contrary to the disc case, the one dimensional problem can be easily studied, both analytically and numerically. The results can be used to check the validity of some of our assumptions leading to conclusions about the disc dynamics. They can also provide some understanding of the disc dynamics in the non-linear regime when the external pressure perturbation is of the order of disc pressure at the density level where the external forcing is applied.
It turns out that the analogy between the disc problem and the one dimensional column problem is practically exact in linear perturbation theory, provided that we consider perturbations of the column subject to the upper boundary condition where P 0 is the equilibrium distribution of pressure given by equation (10), v is the gas velocity in the vertical direction, P e is the externally applied pressure and a is a dimensionless parameter. But note that if unphysical behaviour is to be avoided with this boundary condition, we require a < 0, see below. All upper boundary quantities are assumed to be evaluated at some coordinate ζ * corresponding to the density ρ * ≡ ρ(ζ = ζ * ) ≪ ρ c . Note that in the limit a → ∞, we must have v(ζ = ζ * ) = 0 in order to have a finite pressure at the boundary. Then the condition (95) reduces to a rigid wall condition. The lower boundary remains free.
We go on to consider the dynamics of the vertically stratified column in detail, both analytically and numerically. As in the disc case we set γ = 5/3.

The eigenvalues for the linear modes of the one dimensional slab
Let us assume that the perturbed quantities depend on time through a factor ∝ e −iσt . In this case the free oscillations of the gas column are described by equation (45) for any value of σ provided that we substitute the factor σ 2 − Ω 2 for the factor 2Ωω in the second term on the left hand side. Clearly, the solutions to the problem can be found from equations (53) and (55) with the parameter κ expressed in terms of σ as κ = 1 + 3 σ 2 Ω 2 .

The case of two free boundaries
If we assume that velocities are regular at both boundaries, taken to be where the density vanishes, this regularity condition determines a set of eigenvalues with frequencies σ n Ω.
The set of σ n as well as the associated eigen functions can be easily determined from (55).
It is important for our discussion that the mode corresponding to the smallest eigenvalue, σ = Ω, has a velocity that is independent of ζ. 5 . Therefore this mode describes a uniform translation of the column as a whole.

The case of an upper rigid boundary
When we adopt the boundary condition (95) together with the regular boundary condition at ζ → −ζ 0 , the eigen frequency of the shift mode acquires a small correction Ω → Ω + ω.
The resulting mode may be called as a 'modified shift' mode. As we mentioned above one can consider the column to be rotating around some centre with angular frequency Ω. In the non-rotating frame, the shift mode of the free system describes a stationary displacement of the column as a whole, while for the case with a rigid boundary, the modified shift mode corresponds an oscillation of the column with characteristic time scale ∼ 1/ω.
The change to the eigenfrequency, ω, can be readily be determined by applying the boundary condition (95) to the solution given by equation (55) under the assumption that |a| and λ = ρ * /ρ c ≪ 1.
This is found to be approximately given bỹ where we recall thatω = ω/Ω. It follows from (96) that when a > 0 the perturbations of the column decay with time. The condition a < 0 leads to instability. Making the identificatioñ ω r ≡ − 3 10 a we see that equation (96) is precisely equivalent to equation (78) withω ν = 0. Therefore, the linear dynamics of the gas column with boundary condition (95) should capture the form of the vertical motions relevant to the disc dynamics.
We see from equation (96) that the real and imaginary parts ofω, namely ω R and ω I are functions of only one variable y = a/λ. Thus we have Thus, all models with different a and λ but the same ratio of these quantities have the same value of ω expressed in units of λ . As seen from equation (97) the absolute value of ω I has a maximum when y = 5 with |ω I | = 2|a| 15π . This corresponds to the condition (88).

Non linear numerical solutions for the vertically stratified slab
In our numerical work we use the simple Lagrangian staggered leap frog scheme described in Richtmyer & Morton 1967 for an ideal gas with γ = 5/3. We use 400 grid points in the vertical direction uniformly distributed with respect to the vertical coordinate in the static initial state. Artificial viscosity is used in order to handle possible shocks and it is checked that our results are essentially independent on the value of the artificial viscosity coefficient.
In order to check the scheme we compared the eigenvalues for the free problem σ n and the Rankine-Hugoniot relations for shocks with what is obtained in numerical computations. In both cases we obtained good agreement between analytical and numerical approaches.
In our computations we use the dimensionless time τ = Ωt. The density ρ, pressure p, velocity v and position of a gas element are assumed to be expressed in units of ρ c , P c , Ωζ 0 and ζ 0 , respectively. These normalisations are used to determine the normalisations of all other quantities that make them dimensionless. For example the thermal energy per unit area of the unperturbed disc, E th = 2 γ−1 ζ 0 0 dζp. Following our procedure, this quantity is expressed in units of P c ζ 0 . Thus we have E th = 5π 48 ≈ 0.33 for γ = 5/3 in these units.

Numerical check of the eigenvalues and structure of the linear eigenmodes
In order to check the validity of the linear theory we consider the case a > 0 corresponding to decay of warp like perturbations. The case a < 0 leads to a model with unphysical behaviour because it is such that, not only the modified shift mode, but also the higher order acoustic modes are unstable. The growth rates of the latter are always larger than that of the modified shift modes leading them to dominate the evolution.
For a > 0, we start our computations with initial state having the density and pressure distributions corresponding to hydrostatic equilibrium and a very small uniform vertical velocity v = 10 −4 . We evolve the system for a very long time τ ∼ 10 3 −10 4 . Then we determine power spectra for the state variables and use them to locate the peak corresponding to the modified shift mode. It is confirmed that the associated spectral feature or line has a Lorentz profile. The location of the maximum determines the real part of the eigen frequencyω R , and the width of the line determinesω I .
We consider four different values of λ = 0.032, 0.056, 0.084, 0.11 and seven different values of a for a given value of λ, such that a = 80λ/2 k , where the integer k = 1, 2..7. When a is smaller than 5λ -the value corresponding to the maximal value of ω I we expect the column dynamics to be similar to the disc dynamics in the regime described by Pringle 1996. In this case the values of velocity and the Lagrangian displacement should not depend significantly on ζ and the column response to the external pressure is similar to that of a rigid body.
However, when a > 5λ, we expect the behaviour to be similar to the disc dynamics in the The results of our computations are plotted in Figures 1-4. In Figure 1 we show the numerically obtained dependencies of ω R /λ (the dashed curves ) and ω I /λ (the dotted curves) on a/λ for the set of parameters described above. Specific symbols are associated with particular values of λ. Circles, squares, diamonds and triangles represent the specific values of λ we used in increasing order. Solid curves correspond to the analytically determined dependencies given by equation (97). One can see that the agreement between the analytic and numerical results is good for the smallest value of λ = 0.032. However, the curves deviate more significantly as λ is increased. This departure from a dependence of ω/λ only on a/λ may be explained as an effect due to the increasing importance of of higher order terms in the power series representation of ω/λ as a function of λ, for fixed a/λ, that were not taken into account in equations (60)    tically the same and such that ω I ≈ 1.65 · 10 −3 . As expected the temporal behaviour of the displacement vector is almost indistinguishable for these two cases. The dashed curve indicates the exponential dependence on time, calculated in this case, using the decay rate ω I = 1.65 · 10 −3 . It is evident that agreement between the analytical and numerical results is good. Figure 3 is similar to Figure 2 but the dependence of the displacement vector on time is shown for a = 0.16.. Because a/λ = 5 this value of a corresponds to the maximal decay rate when λ = 0.032. In this case it is given by ω I ≈ 6.7 · 10 −3 . Note again the very good agreement between the theoretical and numerical results.
In Figure 4 we show the dependence of the displacement vector on ζ calculated at the same time for a = 1.26 -solid curve, a = 0.16 -dot-dashed curve and a = 0.02 -dashed curve. One can see from this Figure that the displacement corresponding to the small value of a is practically independent of ζ. But for the large value of a it has a prominent increase toward the upper boundary of the column. Also, note that in the latter case the value of displacement vector near the upper boundary is close to zero. Thus, although the time  Table 1. The initial amplitude of velocity profiles, C, the maximal value of the relative pressure difference (P − P 0 )/P 0 ≡ ∆P/P 0 calculated for the first pulsational period and the difference between the final and initial values of the thermal energy per unit area, ∆E T h , for the different models used in the computations. Note that quantities are expressed as multiples of the dimensional units defined in the text. 1 2 3 4 5 6 7 8 C 10 −3 5 · 10 −3 10 −2 2 · 10 −2 4 · 10 −2 8 · 10 −2 1.6 · 10 −1 3.2 · 10 −1 ∆P/P 0 4 · 10 −2 0.2 0.4 1 2.5 7 20 80 ∆E T h 9 · 10 −8 1.4 · 10 −6 5.4 · 10 −6 2 · 10 −5 8.6 · 10 −5 4.4 · 10 −4 2 · 10 −3 9 · 10 −3 dependence of the displacement vectors corresponding to the cases of large and small a is practically the same, the spacial dependence of these quantities is quite different.

Non linear calculations
In order to determine how our results depend on the initial velocity amplitude, we consider the case for which the upper boundary of the column is a rigid wall, resulting in the condition that v(ζ = ζ * ) = 0. From our previous discussion it follows that this occurs when the parameter a is sufficiently large. The uniform initial velocity profile used in our previous numerical calculations is not compatible with this boundary condition. Therefore, we use a different initial form for the velocity Since the non-linear effects are expected to operate on a short dynamical time scale of order of a few periods P = 2πΩ −1 we follow the evolution for relatively small times τ end = 200.
When the rigid wall boundary condition is applied and the linear theory is valid, the column dynamics is determined by the linear superposition of strictly periodic oscillations.
Departures from this type of motion are determined by non-linear effects. One of the most important for our purposes is the possibility of shock formation near the upper boundary.
When such shocks propagate through the column they irreversibly heat up the gas thus causing dissipation of the motion induced in the column. 6 The strength of non-linear effects is obviously proportional to the induced amplitude C.
In the first row of Table 1 we show the eight values of C used in our computations and the 6 Of course, shock formation in a realistic system would be quite different from what we consider in the text. Since the upper boundary of the disc is close to the surface of unit optical depth, in a realistic situation, the shocks would be radiative. However, it is plausible that our study indicates the range of disc parameters for which shocks may be expected. corresponding models, enumerated by numbers from one to eight. In all cases we assume that the rigid boundary is situated at the value of ζ * corresponding to λ = 0.032.
A useful measure of non-linearity of the system is the ratio of the difference between the pressure at the upper boundary and its equilibrium value ∆P = P − P 0 , to P 0 . For definiteness we estimate the maximal values of this quantity during the first non-linear pulsation and display them in the second row of Table 1. As seen from this Table when ∆P/P 0 < 1, this quantity has an approximate linear dependence on C, while for larger values this dependence is somewhat steeper. The third row represents the difference between the thermal energy (per unit area) at the end of the computations at τ end = 200 and its initial value E th ≈ 0.33, in units of P c ζ 0 . As expected this quantity is approximately proportional to the square of C. Figure 6. The dependence of the entropy per unit area, S, for the six models (3 − 8). As the entropy increases while moving from the lowermost to uppermost curve, the models range from 3 − 8 with increasing values of the initial amplitude C.

Shock development and the formation of a rarefied hot atmosphere
When ∆P/P 0 < 1 the dynamics of the model is well described by the linear theory. In the opposite case, non-linear effects are substantial. The system exhibits a complicated pulsational behaviour and strong shocks are observed. This regime is illustrated in Figure   5 where the results of calculations for model 7, with large initial amplitude C = 0.16, are presented. We show the density profiles at four different times τ = 2.5 ( solid curve ), τ = 5 (dashed curve ), τ = 7.5 ( dotted curve ) and τ = 10 ( dot dashed curve ). The density distribution at τ = 2.5 is quite similar to the equilibrium density distribution. When τ = 5 and τ = 7.5 shocks are observed propagating from the upper boundary in the direction of negative ζ. These shocks heat the gas up, thus dissipating kinetic energy. Also note the density drop at the upper boundary when τ = 5 and τ = 10. We have checked that this density drop has a persistent character and the density at the boundary at τ = τ end is significantly smaller than the initial density while the pressure does not exhibit such a significant change.
Thus, the presence of shocks in the system leads to formation of a hot rarefied region near the rigid boundary. It is not clear however, how this effect would operate in twisted discs. In the latter case the boundary is defined by the condition that the optical thickness is one. The position of this boundary may move when any density redistribution becomes significant. In our case with fixed upper boundary location, the presence of this hot low density region leads to the suppression of shock formation after a time corresponding to a few orbital periods, an effect that may not occur if the upper boundary moved in such a way as to maintain the relative pressure variations.
The effect in our case is illustrated in Figure 6 where we show the dependence of entropy per unit area, S = dζρ log(p/ρ γ ) as a function of time for the models with indices larger than 2, in the natural units ρ c ζ 0 . Curves attaining larger values correspond to larger amplitudes. One can see from this Figure that when the initial amplitude is sufficiently large, i.e.
C 0.04 and accordingly ∆P/P 0 2.5, half of the entropy change is roughly determined for τ < 10 − 15. (i.e. approximately during the first two orbital periods). Then shock production becomes less efficient due to the mechanism discussed above.

Estimate of dissipation time scale
As mentioned above, it is not clear to what extent the dynamics of the column in the strongly non linear regime is similar to the non-linear dynamics of twisted discs. One can assume, e.g. that in twisted discs, the boundary where external radiation pressure applies follows a certain density level. In this situation it seems reasonable to assume that the shock formation might persist at a rate similar to what is observed at initial times τ < 10. In this case the disc perturbations could be damped on a time scale t nl which can be roughly estimated from the results presented in Figure 1.
For model 7, assuming the thermal energy production is of the order of ∆E th = 0.002, we can estimate the non linear decay time as t nl ∼ 10(E th /∆E th )Ω −1 ∼ 1.5 · 10 3 Ω −1 where we recall that E th ≈ 0.3 in our units. Thus, under these assumptions the perturbations may decay in approximately 200 orbital periods. This rate is characteristically the internal energy content within one surface scale height per orbital period and becomes noticeable when the relative pressure variation ∆P/P 0 near the upper boundary exceeds ∼ 2.5.
Let us stress, however, that although indicative, because of an obvious sensitive dependence on the density structure near the boundary, such an estimate cannot be rigorously justified within our model. Clearly more comprehensive future investigations are needed.

DISCUSSION AND CONCLUSIONS
In this paper we have presented a self-consistent calculation of the response of a twisted disc to the action of a radiation pressure force acting in upper layers of the disc. This is assumed to be due to the interception and subsequent reradiation of radiation from a central source. The radiation pressure force is assumed to be applied at a single density level ρ * corresponding to optical thickness unity. The degree of twist and warping is assumed to be small enough that linear theory can be used to calculate the response.
The analysis of Pringle 1996 modelled the disc as a set of concentric rigid rings in a state of Keplerian rotation. These were assumed to communicate with each other through the exchange angular momentum because of the action of viscous forces. Up to now there has been no consideration of effects arising from the response of the disc vertical structure to the externally applied pressure.
In this paper we have considered the warping and twisting of an accretion disc taking account the response of the disc vertical structure to an external radiation field assuming that the effect of self-shielding of radiation by the global twisted disc is not significant. We developed a description of the evolution of the disc in terms of a pair of equations governing the one dimensional evolution of the inclination as a function of radius and time (see section 6 and equations (41), (67) and (68)). This description extends earlier ones, so that the influence of radiation pressure on discs for which warps occur in the low viscosity bending wave regime, as well as the higher viscosity diffusive regime, may be considered.
We found several qualitatively different dynamical regimes that may be realised in astrophysical discs. These are related to whether the character of the response of the disc vertical structure causes significant departures from what would be obtained for a rigid body.

Conditions for a quasi-rigid response
We found that to avoid such departures, the momentum flux due to radiation from the central source F 0 , should be smaller than a critical value, F * , given by F * = (λ/δ)P c , where δ = ζ 0 /r is the disc aspect ratio, P c is the disc mid plane pressure and it has been assumed that the the ratio (λ/δ) is small (see sections 6.5 and 7.1.2). Then, when F 0 > F crit ≈ 0.1(δ/α)P c , warping instability of the disc becomes a possibility ( eg. Pringle, 1996).
Thus for radiation driven warping of a disc that behaves like an ensemble of rigid rings, we find two requirements, namely that F * > F crit , together with F 0 < F * . These conditions together imply that the ration of the surface to mid plane density should be sufficiently large and such that λ > λ crit ≈ 0.1δ 2 /α.

9.2
Large external radiation momentum flux In the opposite limit of large F 0 > max(F * , F crit ), consideration of the disc structure in the vertical direction is essential as the vertical displacement ceases to be uniform. The perturbed gas motion in the disc is mainly determined by the vertical component of velocity and the perturbed quantities tend to oscillate at a characteristic frequency ω = (ρ c ζ 0 /2Σ)λΩ, with Ω being the local Keplerian angular velocity. In this situation, vertical motion tends to be suppressed by the external pressure field on the irradiated side of the disc while warping motions persist in the bulk of the disc and on the shielded side of the disc (see sections 6.5 and 7.2).
The disc surface intercepting the radiation tends to align parallel to the radiation rays, thus decreasing the amount of intercepted radiation, while the inclination angles associated with the disc mid plane and the opposite free boundary of the disc oscillate, being approximately equal to each other. In this limit the upper surface of optical thickness one plays the role of an impenetrable wall. Thus the development of instability of the inclination angle due to radiation pressure effects (e.g. Pringle 1996) would be inhibited in this limit.
In principle, the presence of warping motions in this regime would be associated with variability on a large characteristic time scale T ch ∼ λ −1 T K , where T K and λ are some 'typical' values of the orbital period and the density ratio.

Limits on the WKB growth rate
Following Pringle 1996 we have considered possible instabilities using a WKB approach.
As this neglects potentially important global effects and boundary conditions results are not definitive, nonetheless they should give a good indication of when instability could be possible.
We find that when F 0 < F * the growth rate increases with F 0 while in the formal limit F 0 → ∞ it approaches zero. Thus in the linear theory there is an upper limit for the growth rate of the instability of the inclination angle which is always (ρ c ζ 0 /Σ)λΩ in the absence of viscosity (see section 7.3). When present, viscosity acts to damp any instability at a rate 0.1δ 2 Ω/α, leading again to the condition λ > λ crit ≈ 0.1δ 2 /α for radiation warping instability to be feasible.

Conditions for non-linearity
Our results described above rely upon the applicability of linear perturbation theory. As discussed in section 7.5, the breakdown of linear theory is expected when the Lagrangian change of pressure induced in the surface layers becomes of the order of the unperturbed pressure even though the change of inclination angle could be very small. The corresponding constraints on the inclination angle are especially strong for disc models with a significant drop of temperature toward the boundaries of the disc.

Non linear calculations for the one dimensional slab analogue
A direct numerical approach to the issues discussed above is not yet feasible due to threedimensional nature of the problem and the many physical processes involved. However, when vertical motions dominate, the problem becomes analogous to the one dimensional problem of vertical motions induced in a stratified gas column orbiting around a central source with Keplerian angular velocity.
As discussed in Section 8 the one dimensional slab gives the same fundamental period of oscillation as obtained from the full disc theory when an appropriate boundary condition on the upper surface of the column is specified. The dependence of the eigenfrequency on theoretical parameters as well as the effective presence, in the limit of small λ, of an impenetrable wall at the upper surface of the column have been checked numerically. Where they can be compared, agreement between our analytical and numerical results is very good.
We also used the one dimensional model in order to investigate possible consequences of non-linear behaviour in the system. To do this we studied the motion of the slab ensuing from the imposition of a vertical velocity profile with varying amplitude. We found that when the ratio of the Lagrangian pressure perturbation to the local value of the pressure at the upper boundary of the disc, ∆P/P * , becomes larger than 1 − 10, strong shocks propagating downward into the column are observed (see section 8.3). In principle, these shocks may lead to non-linear dissipation of energy stored in the vertical motion and in section 8.4 we made a very crude estimate of a possible time scale of 200 orbits for ∆P/P * ∼ 10. However, this result must be checked in framework of a more sophisticated numerical approach which takes into account physical processes neglected in this study, most notably the effects of radiation transfer in the upper layers of the column.

Issues for future consideration
In a fully self-consistent study the dynamical effects of radiation pressure must be studied together with effects determined by the radiation heating of the disc atmosphere. This can be done in the most convenient way within the framework of the one dimensional model discussed above.
In a realistic disc model, where effects due to the flaring of the disc lead to the interception and scattering of radiation in the disc photosphere, when axisymmetric and unperturbed, radiation heating can significantly influence on or even determine the value of the density ratio λ.
This parameter, being the ratio of the density at the absorbing surface to the mid plane density defines the boundary between different dynamical regimes for a twisted disc. In this connection it is interesting to note that in certain vertical models of X-ray heated accretion discs, the ratio λ can be quite small, being of order of 10 −4 − 10 −5 , e.g. This assumption is essentially equivalent to the statement that all vertical motions can be accounted for by a rigid tilt and it precludes internal vertical motions that vary with the vertical coordinate.
As discussed in the main text, this is not strictly valid, as in general such motions occur.
Accordingly, the density surfaces corresponding to an optical thicknesses of order of unity where the radiation pressure force is applied are not everywhere parallel to the surface locally corresponding to the disc mid-plane and the expression for the radiation pressure acting on the disc must be, accordingly modified.
In general, an expression describing the surfaces of constant optical thickness can be easily obtained from the transformation law from the Cartesian coordinate system to the twisted coordinate system written for a gas element situated near the surface of the disc and where the upper (lower) sign corresponds to the upper (lower) disc surface. Let us recall that we assume that the inclination angle β is small and derive all equations in the linear approximation for the angle β and the Lagrangian displacement vectorξ.
The unit vector perpendicular to the disc surface is given by the standard expression whereR is the radius vector with components (x, y, z).
We make the standard assumption that the radiation propagates radially before being intercepted by the disc and neglect the effect of self-shielding. In this case an amount of radiation intercepted by the disc and subsequently reradiated away is proportional to cosine of angle between the radius vector and the normal to the disc, cos Ξ ≡n ·R/|R|. For the upper branch of the surface a simple calculation gives cos Ξ = −r ∂ ∂r where we recall that W = β ′ sin ψ − βγ ′ cos ψ = Ψ 1 sin φ − Ψ 2 cos φ. As we discuss in the main text, in this paper, for simplicity, we neglect the effect of disc flaring and accordingly set ∂ ∂r ( ζ 0 r ) = 0 in equation (A4). The upper (lower ) branch of the disc surface is able to intercept the radiation of a central source only when cos Ξ < 0 ( cos Ξ > 0). Therefore the radiation pressure, F , applied to the disc surface has the form where Θ(x) is the step function, the signs (+) and (−) stand for the upper and lower branches of the disc surface, respectively, is the momentum flux per unit of surface carried by the radiation emitted by the central source, and L is the luminosity of the source. The factor (2/3) accounts for (assumed ) isotropic character of the disc reradiation.
In order to have the same symmetries as the surface pressure term (A5) the ζ-component of the displacement vector, ξ ζ , and, accordingly, the ζ-component of velocity perturbation, v ζ , must satisfy On the other hand, for a given value of φ, v ζ may be separated on even and odd components: v ζ = v ζ e + v ζ o and both components should be taken into account in order to satisfy the appropriate boundary conditions at ζ → ±ζ 0 . However, the latter requirement poses a problem with the symmetry law (A7). Indeed, while the even part of v ζ , v ζ e , agrees with (A7) the odd one is not and the standard separation of the velocity perturbations into even and odd parts although being respected by the perturbed equations of motion, is not, strictly speaking, compatible with the symmetries provided by the boundary conditions. This problem can be circumvented with the help of the assumption that the normal decomposition into even and odd parts in ζ of the perturbed vertical velocity is replaced by where v ζ e and v ζ o have even (odd ) properties with respect to the reflection ζ → −ζ. Since the surface radiation pressure vanishes when cos Ξ = 0 or equivalently sin(φ − φ 1 ) = 0 and, accordingly, we have v ζ o = 0 at the corresponding values of φ = φ 1 and φ = φ 1 + π the velocity anzatz (A8) is self-consistent. Also, it agrees with the symmetries determined by equation (A5) and it has local simmetries respected by the equations of motion.
Above and in the main text we always consider the angle φ to be in a segment [φ 1 ..φ 1 + π] defined by the condition sin(φ − φ 1 ) > 0. With help of this convention we can treat our dynamical variables as having the standard even and odd symmetries with respect to reflection ζ → −ζ. Continuation of these variables into the segment defined by condition sin(φ − φ 1 ) < 0 can be made with the help of the expression (A8). In the case when sin(φ − φ 1 ) > 0 it follows from equations (A4) and (A5) that we have and where we recall that the density ρ * corresponds to the surface of an unperturbed disc with optical thickness τ ≈ 1. Taking into account that we can write F + = F 1 + cos φ + F 2 + sin φ and ξ ζ = ξ ζ 1 cos φ + ξ ζ 2 sin φ, we can represent equation (A9) in a convenient complex form where F + = F 1 + + iF 2 + , ξ ζ = ξ ζ 1 + iξ ζ 2 and W = βe iγ = Ψ 1 + iΨ 2 . Obviously, we also have F − = 0.

APPENDIX C: THE GAUGE TRANSFORMATIONS
As we pointed out in the main text different twisted coordinate systems can describe the same physical situation provided that they are connected by a certain family of transformations (gauge transformations) which maintain fixed distributions of physical quantities characterising the system (density and velocity fields in the case of a baratropic gas) in a non rotating (say, cylindrical) coordinate system with origin located at the central source.
These transformation laws can be found by considering transformations from the cylindrical to the twisted coordinates. We derive them here assuming, for simplicity, that the equation of state is baratropic and the gas is inviscid. The latter assumption allows us to neglect the presence of the radial component of velocity in an unperturbed background state.
By construction the quantity ρ c 1 is not changed after a gauge transformation. Therefore, it may be classified as a 'gauge independent' quantity. The result (C2) can be represented in a compact form with help of the complex notation introduced in equation (19) ρ c 1 = ρ 1 + iDW.
In an absolutely analogous way we can find the law of transformation of velocity perturbations. We have and v rc = v r − iζẆ − ζΩW.
Here the quantities v φc , v zc and v rc represent the velocity perturbations in the cylindrical coordinate system. They are gauge independent by construction. The quantities v φ , v ζ and v r represent the time derivative of the corresponding coordinates of a comoving gas element.
As we discussed above these quantities do not coincide with the projections of velocity onto the orthonormal basis used in the formalism exploiting the twisted coordinate system. The relations between these projections and the time derivatives of the coordinates can be easily found from the results provided in the Appendix of Petterson 1978, his equations (A3). They can be written in the form where it is implied that the background parts of the (ζ) and (r)-components of velocity can be neglected. Substituting equations (C8) in equations (C5)-(C7) we find v zc =v ζ + rΩW, and v rc =v r − ζΩW.
Now let us suppose that there are two twisted coordinate systems with two different W 1 and W 2 describing the same physical system. That means that the quantities characterising the density perturbation and perturbations of the velocity field in the fixed cylindrical system must be the same when expressed in terms of Euler angles and corresponding perturbations belonging to the different twisted coordinate systems. From this requirement and equations (C4) and (C9)-(C11) it follows that the differences in the corresponding to these systems perturbations of density and velocity components must be proportional to the difference between W 1 and W 2 , δW = W 1 − W 2 , and δv r = ζΩδW.
Equations (C12)-(C15) define the gauge transformations between two twisted coordinate system for an inviscid disc. Viewed in these terms, the generalisation to the case of a viscid disc is straightforward provided that the background part of the (r)-component of velocity is taken into account.
One can also check that equations (C12)-(C15) provide an exact solution of the equa-tions of motion and the continuity equation for an inviscid disc being written in a twisted coordinate system, for an arbitrary dependence of W on t and r. This solution describes the so-called trivial modes corresponding to an unperturbed disc. In this case the density and velocity fields in the cylindrical coordinate system remain equal to their unperturbed values while the dynamics in the twisted coordinates is such that perturbations of density and velocity defined with respect to the twisted coordinates are exactly compensated by those associated with time dependent Euler angles.
It is important to stress that our dynamical equations (23) and (25) as well as the continuity equation (29) are not invariant under the transformations generated by (C12)-(C15) even when the term determined by viscosity are discarded. This is due to the fact that certain terms are neglected in these equations. These terms provide only next order corrections in our expansion series in small parameters when the gauge (57) is fixed. On the other hand the combination ξ ζ /r + iW entering in the expression (A11) determine an invariant quantity -the angle between the normal to the disc surface and the radial direction. This combination must be the same in all twisted coordinate systems connected by transformations (C12)-(C15). Using equations (37) and (C8) it is easy to see that this combination is proportional tô v ζ + rΩW.
From equation (C14) it follows that the expression (C16) is, indeed, a gauge invariant quantity.