Nonlinear variations in axisymmetric accretion
Abstract
We subject the stationary solutions of inviscid and axially symmetric rotational accretion to a timedependent radial perturbation, which includes nonlinearity to any arbitrary order. Regardless of the order of nonlinearity, the equation of the perturbation bears a form that is similar to the metric equation of an analogue acoustic black hole. We bring out the time dependence of the perturbation in the form of a Liénard system, by requiring the perturbation to be a standing wave under the second order of nonlinearity. We perform a dynamical systems analysis of the Liénard system to reveal a saddle point in real time, whose implication is that instabilities will develop in the accreting system when the perturbation is extended into the nonlinear regime. We also model the perturbation as a highfrequency travelling wave, and carry out a WentzelKramersBrillouin analysis, treating nonlinearity iteratively as a very feeble effect. Under this approach both the amplitude and the energy flux of the perturbation exhibit growth, with the acoustic horizon segregating the regions of stability and instability.
pacs:
98.62.Mw, 46.15.Ff, 47.10.Fg, 47.50.GjI Introduction
Compressible fluids, possessing angular momentum, execute rotational motion when they are drawn into the gravitational potential well of an accretor (Frank et al., 2002). This fact assumes particular significance when the accretor is a black hole, because astrophysical black holes make their presence felt only by their strong gravitation. No direct spectral information about black holes can ever be forthcoming, due to their event horizons. So evidence of black holes can only be known from the influence of their strong gravitational fields on any proximate astrophysical fluid. A system of rotational accretion onto a black hole has three principal phenomena governing the flow, namely, gravity (as implied by general relativity), transport of angular momentum, and nonlinearity. To take these up one after the other, first, when the accretor is a black hole, we concern ourselves with the physical character of curved spacetime. Here we work in the spherically symmetric Schwarzschild geometry, by having recourse to what is known as a pseudoNewtonian potential function (Paczyński and Wiita, 1980; Nowak and Wagoner, 1991; Artemova et al., 1996; Stuchlík and Kovář, 2008), something that, in a Newtonian framework, mimics the salient properties of Schwarzschild spacetime. We do not lose the essential effects of general relativity even in this Newtonianlike perspective (Carter, 1979), in which we have replaced the Schwarzschild spacetime by an equivalent potential field. Apropos of gravity, another usual assumption is that the fluid is nonselfgravitating, with its flow being axisymmetric (Frank et al., 2002).
Regarding the motion of the fluid, as it drifts into the central black hole, our next concern is the centrifugal effect in the flow, which arises because test fluid elements possess angular momentum. This then brings up the question of the transport of angular momentum. Now, in a simple Newtonian treatment, the centrifugal effect will prevent a test fluid element with nonzero angular momentum from reaching the centre of the potential well (Carter, 1979). Relativistic theory differs on this point, in that the very strong gravity of a black hole wins against centrifugal repulsion, especially in the vicinity of the event horizon (Carter, 1979). So, fluid elements with sufficient energy are able to overcome the centrifugal potential barrier and fall into the gravitating mass, whatever their angular momentum, and fluid elements with low (not merely zero) angular momentum are captured by the hole (Frank et al., 2002). In other words, if the axisymmetric flow has low angular momentum, and is driven by a very strong gravitational field (due, say, to a supermassive black hole), the radial drift will be much more conspicuous than the azimuthal drift. Now, an outward transport of angular momentum is effected through the azimuthal dynamics, because of viscous shearing between two differentially rotating adjacent layers of the fluid, something that makes for a slow Keplerian infall of accreting matter into the potential well (Frank et al., 2002). In contrast, if the angular momentum is low and the flow is highly subKeplerian, then the dynamics is dominated by the radial drift, suggesting that the time scale of the dynamic infall process is much smaller than the viscous time scale (Liang and Thomson, 1980; Abramowicz and Zurek, 1981), which pertains to the azimuthal dynamics. In that case the viscous outward transport of angular momentum is not of much consequence, and we consider the flow inviscid. This is the underlying principle of a subKeplerian inviscid disc with low angular momentum (Abramowicz and Zurek, 1981; Chakrabarti, 1989; Nakayama and Fukue, 1989; Chakrabarti, 1990; Yang and Kafatos, 1995; Pariev, 1996; Chakrabarti, 1996; Lu et al., 1997; Das, 2002; Das et al., 2003; Ray, 2003; Barai et al., 2004; Das, 2004; Abraham et al., 2006; Chaudhury et al., 2006; Ray and Bhattacharjee, 2007a; Goswami et al., 2007; Das et al., 2007; Roy and Ray, 2009; Nag et al., 2012). In our study we have adopted this model, which effectively neglects the dissipation of both energy and angular momentum, thereby treating both as constants of the motion in a perfect fluid (Liang and Thomson, 1980; Abramowicz and Zurek, 1981; Chakrabarti, 1989). The accretion process purported by a conserved subKeplerian flow can be facilitated further if the radial velocity is significant even far away from the inner region of the disc, and on large length scales of the flow, radial velocities of such order can result, if the angular momentum is low (Igumenshchev and Beloborodov, 1997; Beloborodov and Illarionov, 2001; Proga and Begelman, 2003). Highly subKeplerian flows are found in detached binary systems fed by accretion from stellar winds (Illarionov and Sunyaev, 1975; Liang and Nolan, 1984), semidetached nonmagnetic binaries (Bisikalo et al., 1998), and supermassive black holes fed by accretion from slowly rotating stellar clusters (Illarionov, 1987). Also, in geometrically thick accretion discs, turbulence may produce subKeplerian flows (Igumenshchev and Abramowicz, 1999).
Now we take up nonlinearity, which is the primary subject of our study. Mathematically, all of fluid dynamics is a nonlinear problem, and accretion flows are no exception to this rule. However, the nonlinear attribute of astrophysical accretion has not been addressed much in the literature. Here our chosen accretion model of a conserved axisymmetric flow in a pseudoNewtonian framework allows us to focus mainly on nonlinearity, by rendering the other analytical aspects of the flow simple. Nevertheless, the mathematical task remains challenging enough, involving nonlinear partial differential equations of fluid dynamics (Landau and Lifshitz, 1987). Even in the relatively simple stationary limit of these equations, we can make a case for nonlinearity. In accretion flows the usual boundary conditions are that at large distances from the accretor, the flow is subsonic, but very close to the accretor, the flow ought to become highly supersonic, a fact borne out if the accretor is a black hole (Novikov and Thorne, 1973; Shapiro and Teukolsky, 1983; Liang and Thomson, 1980). So, in the intermediate region, the bulk flow attains the speed of acoustic propagation in the fluid, and becomes transonic in character. These critical conditions are the consequences of coupled firstorder dynamical systems, crafted out of the equations governing stationary rotational accretion (MuchotrzebCzerny, 1986; Afshordi and Paczyński, 2003; Chaudhury et al., 2006; Ray and Bhattacharjee, 2007a; Goswami et al., 2007; Roy and Ray, 2009; Bhattacharjee et al., 2009; Nag et al., 2012), in a classic textbook approach to nonlinear problems Strogatz (1994); Jordan and Smith (1999). From stationary flows alone, however, we cannot fully gauge what nonlinearity implies for astrophysical accretion. So, our larger quest lies in in the dynamics. In arguing for both dynamics and nonlinearity, we just need to look at the wellknown nonlinear problem of the realizability of solutions passing through saddle points in a stationary phase plot. The very existence of these critical solutions is threatened by even an infinitesimal deviation from the precisely needed outer boundary condition to generate the solutions numerically (Ray and Bhattacharjee, 2002, 2007a). We can, however, overcome this difficulty by a nonperturbative temporal evolution of a globally subcritical solution towards the critical state. Under qualifying approximations, such an analytical study is known (Ray and Bhattacharjee, 2007a). Nevertheless, nonlinear problems never submit themselves easily to mathematical analyses, and there exist no analytical solutions of the fully nonlinear coupled field equations governing an accretion flow. So in the absence of any analytical formulation of the complete dynamics of the flow solutions, a tentative first step is usually a timedependent perturbative analysis, and one such study (Ray, 2003) concluded that perturbations on the flow do not produce any linear mode with an amplitude that grows in time, i.e., the background solutions are marginally stable under small perturbations. The stability of inviscid subKeplerian flows has been studied in various ways (Ray, 2003; Chaudhury et al., 2006; Bhattacharjee and Ray, 2007), but conclusions made about stability under linearization can scarcely be extended to conditions governed by nonlinearity. Here is our attempt to bridge the gap.
First, we implement a timedependent radial perturbation scheme on an inviscid axisymmetric accretion flow, retaining all orders of nonlinearity in the equation of the perturbation that follows. A striking feature of the equation of the perturbation is that even on accommodating nonlinearity to any arbitrary order, it conforms to the structure of the metric equation of a scalar field in Lorentzian geometry (Section III). This fluid analogue (an “acoustic black hole”), emulating many features of a general relativistic black hole, is a matter of continuing interest in fluid mechanics from diverse points of view (Unruh, 1981; Jacobson, 1991; Unruh, 1995; Visser, 1998; Bilić, 1999; Schützhold and Unruh, 2002; Singha et al., 2005; Volovik, 2005, 2006; Ray and Bhattacharjee, 2007a, b; Roy and Ray, 2007; Naskar et al., 2007; Das et al., 2007; Mach and Malec, 2008; Barceló et al., 2011; Nag et al., 2012; Robertson, 2012; Sarkar et al., 2013; Sen and Ray, 2014). Then we use the nonlinear equation of the perturbation to study the stability of globally subsonic stationary solutions under largeamplitude (nonlinear) timedependent perturbations. Regarding the nonperturbative evolution of the accreting system, it is reasonable to suggest that the initial condition of the evolution is globally subcritical, with gravity subsequently driving the solution to a critical state, sweeping through an infinitude of intermediate subcritical states. The stability of these subcritical states is essential for a smooth temporal convergence to a stable critical state. To investigate this aspect under relatively simple nonlinear conditions, we truncate all orders of nonlinearity beyond the second order in the equation of the perturbation. Following this, we integrate out the spatial dependence of the perturbation with the help of welldefined boundary conditions on globally subcritical flows. Then we extract only the timedependent part of the perturbation, and find that it has the mathematical appearance of a Liénard system (Strogatz, 1994; Jordan and Smith, 1999) (Section IV). Using standard analytical tools of dynamical systems to study the equilibrium features of this Liénard system, we show the existence of a saddle point in real time (Section V). This implies clearly that the stationary background solutions are unstable in time, if the perturbation is extended into the nonlinear regime. The destabilizing effect of nonlinearity remains qualitatively unaltered, when we go into the case of a highfrequency travellingwave perturbation, and carry out a WentzelKramersBrillouin (henceforth WKB) analysis (Section VI).
Ii The mathematical model of the flow
In models of accretion discs, imposing the condition of hydrostatic equilibrium along the vertical direction and then performing a vertical integration, result in the collapse of the vertical geometry of the flow on the equatorial plane of the disc (Frank et al., 2002). The equatorial flow is described by two coupled fields, the radial drift velocity, , and the surface density, , of which, the latter is defined by vertically integrating the volume density, , over the disc thickness, . This gives , and in terms of , the continuity equation is (Frank et al., 2002)
(1) 
The axisymmetric accretion flow is driven by the gravitational field of a centrally located nonrotating black hole, but the structure of the neighbouring geometry, through which the flow takes place, can still be set according to the Newtonian construct of space and time, using a pseudoNewtonian potential (Paczyński and Wiita, 1980; Nowak and Wagoner, 1991; Artemova et al., 1996; Stuchlík and Kovář, 2008). Our principal conclusions remain qualitatively unaffected by the choice of a particular pseudoNewtonian potential. Involving such a potential, , we can give the height of the disc, under hydrostatic equilibrium in the vertical direction, as , with the local speed of sound, , and the local Keplerian velocity, , defined, respectively, by and . The pressure, , as it has been introduced in the definition of , is expressed in terms of the volume density, , by a polytropic equation of state, , with being the polytropic exponent. In consequence of this definition of , we note that . We can also show from the first law of thermodynamics (Chandrasekhar, 1939) that varies between unity (the isothermal limit) and , which is the ratio of the two coefficients of specific heat capacity of a gas (corresponding to the adiabatic limit), i.e., . So the polytropic prescription is of a much more general scope than the simple adiabatic case, and is suited well for the study of open systems like astrophysical flows. Using the relationship between and , we write the disc height explicitly in terms of the standard fluid flow variables as
(2) 
a result that we use to recast equation (1) as
(3) 
which is one of the two mathematical conditions governing the dynamics of the coupled fields, and .
To ascertain the second condition necessary for determining the dynamics in the radial direction, we have to first look at the dynamics along the azimuthal direction. This gives the balance of the specific angular momentum of the flow as (Frank et al., 2002)
(4) 
in which is the local angular velocity of the flow. The torque, , on the righthand side of the foregoing equation is to be read as , with being the kinematic viscosity. The quantity, , on the lefthand side of equation (4) is the specific angular momentum of the flow. The inviscid rotational flow of our interest can be designed as a limiting case of equation (4), by setting . Then, making use of equation (1), the comoving derivative of , implied by the lefthand side of equation (4), can be made to vanish. As a result we get a simple integral solution, , with the constant of the motion, , being interpreted physically as the constant specific angular momentum of the flow. This interpretation allows us to set down the second mathematical condition of the flow along the radial direction. This is the condition of the radial momentum balance (the Euler equation), in which the centrifugal term, arising due to the rotational motion of the flow, is fixed as . After that, the radial momentum balance equation is written as
(5) 
On specifying the functions, and , equations (3) and (5) give a complete hydrodynamic description of the axisymmetric flow in terms of the two fields, and . By making explicit time dependence of these two fields disappear, i.e., by setting , we obtain the steady solutions of the flow. The resulting differential equations of the inviscid rotational flow, carrying only spatial derivatives, can be integrated to get the stationary global solutions. A noticeable feature of these stationary solutions is that they are invariant under the transformation , i.e., the mathematical problems of inflows () and outflows () are identical in the stationary state (Ray and Bhattacharjee, 2007a). This invariance has adverse implications for the critical flows. Critical solutions pass through saddle points in the stationary phase portrait of the flow (Ray and Bhattacharjee, 2002; Chaudhury et al., 2006; Ray and Bhattacharjee, 2007a), and we can argue that generating a stationary solution through a saddle point will be practically impossible, because it calls for an infinite precision in the required outer boundary condition (Ray and Bhattacharjee, 2002, 2007a). Nevertheless, criticality in accretion processes is not doubted (Liang and Thomson, 1980). The key to resolving this paradox lies in considering explicit time dependence in the flow, because of which, as we note from equations (3) and (5), the invariance under the transformation, , breaks down. Obviously then, we have to make a choice of inflows or outflows at the beginning (at ), and solutions generated thereafter are unaffected by the difficulties associated with a saddle point in the stationary flow.
Imposing various boundary conditions on the stationary integral solutions results in various classes of flow (Chakrabarti, 1989; Das, 2002; Das et al., 2003). Of these, the one of practical interest obeys the boundary conditions, as (the outer boundary condition) and at small values of . The inner boundary condition naturally suggests itself when it comes to accretion onto a black hole, for which, the final infall across the event horizon must be highly supersonic (Novikov and Thorne, 1973; Shapiro and Teukolsky, 1983; Liang and Thomson, 1980). So an open solution that starts with a very low velocity far away from the accretor, but has to allow a fluid element to reach the accretor with supersonic speeds, must pass through a saddle point, where the flow becomes critical. When the stationary flow attains criticality, the bulk flow velocity matches the speed with which an acoustic wave can propagate through the moving compressible fluid. In the case of a polytropic and axisymmetric rotational flow, this speed is not exactly given by the sonic condition, but differs slightly from it by a constant numerical factor, , due to the geometry of the flow (Ray, 2003; Chaudhury et al., 2006; Ray and Bhattacharjee, 2007a; Roy and Ray, 2009; Nag et al., 2012). For a conserved flow, the critical points are either saddle points, through which open solutions pass, or they are centretype points, around which the stationary solutions form closed trajectories (Chaudhury et al., 2006). The solutions that pass through the saddle points may be either homoclinic or heteroclinic (Chakrabarti, 1989; Das, 2002; Das et al., 2003). The number of critical points depends on the choice of the potential, . For the simple Newtonian potential, only two critical points result (Ray, 2003; Ray and Bhattacharjee, 2007a), one a centretype point and the other a saddle point. In studies of axisymmetric accretion onto a black hole, however, it is an expedient practice to employ a pseudoNewtonian potential to drive the flow. In that case, the number of critical points exceeds two, and by the properties of critical points (Jordan and Smith, 1999), there shall be multiple saddle points. Since open critical solutions, connecting the outer boundary of the flow to the event horizon of a black hole, have to pass through saddle points, an inflow process that is made to traverse all these saddle points is multicritical (Liang and Thomson, 1980; Abramowicz and Zurek, 1981; MuchotrzebCzerny, 1986; Chakrabarti, 1989, 1990; Yang and Kafatos, 1995; Pariev, 1996; Lu et al., 1997; Das, 2002; Barai et al., 2004; Das, 2004; Abraham et al., 2006; Das et al., 2007). In such a flow, a fluid element reaches the accretor after having travelled through more than one saddle point, and in between two successive saddle points, the flow suffers a shock (Chakrabarti, 1989; Nakayama and Fukue, 1989; Chakrabarti, 1990; Yang and Kafatos, 1995; Lu et al., 1997; Das, 2002; Das et al., 2003; Abraham et al., 2006; Das et al., 2007). It is at the discontinuity of a shock, that a solution is demoted from its supercritical state to a subcritical one, following which, the solution has to regain supercriticality by travelling through another saddle point (Chakrabarti, 1989; Das, 2002; Das et al., 2003).
Thus far, we have looked at the properties of the accretion flow from a stationary perspective, which is relatively easy to follow. It is the dynamics of the flow that poses a mathematical problem of greater complexity. In comparatively simple studies involving timedependence (Ray, 2003; Chaudhury et al., 2006), the inviscid and axisymmetric flow is found to be stable under the effect of linearized perturbations. This, however, does not say much about the nonperturbative temporal evolution of the velocity and density fields. In such a mathematical problem, we have to work with a coupled set of nonlinear partial differential equations, as implied by equations (3) and (5), but no analytical solutions exist for these coupled dynamic nonlinear equations. Now, the nonperturbative dynamic evolution of global and profiles is crucial in generating the critical flow. It is the way in which the two fields evolve visàvis each other that determines if the critical state would be achieved or not. We can envisage the dynamic process as one in which initially the velocity field, , is subcritical and nearly uniform for all values of , in the absence of any driving force. Then with the introduction of a gravitational field (at ), about whose centre, the fluid distribution may randomly possess some net angular momentum, the hydrodynamic fields, and , start evolving in time. In the regions where the temporal growth of outpaces the temporal growth of (to which is connected), and gravity (given by ) dominates the inhibitive centrifugal effects of rotation (given by ), the infall process becomes supercritical. Otherwise, it continues to remain subcritical. Under the approximation of a “pressureless” motion of a fluid in a gravitational field (Shu, 1991), qualified support for attaining criticality has come from a nonperturbative dynamic perspective (Ray and Bhattacharjee, 2007a), guided by the criterion that the total specific mechanical energy at the end of the evolution would be the same as what it was at the start of the evolution.
Iii Nonlinearity in the perturbative analysis
Equations (3) and (5) are integrated in their stationary limits, to obtain the spatial profiles, and . A standard method of perturbative analysis is to apply small (to a linear order) timedependent radial perturbations on the stationary solutions, and . By this, however, we do not gain much insight into the timedependent evolutionary aspects of the hydrodynamic flow. Our next logical act, therefore, is to incorporate nonlinearity in the perturbative method. With the inclusion of nonlinearity in progressively higher orders, the perturbative analysis incrementally approaches the actual timedependent evolution of the global solutions, after it has started with a given stationary profile at (it makes physical sense to suggest that this initial profile is spatially subcritical over the entire flow domain).
We prescribe the perturbation as and , in which the primed quantities indicate a perturbation about a stationary background. We define a new variable, , following a mathematical procedure employed previously in studies on inviscid axisymmetric accretion (Ray, 2003; Chaudhury et al., 2006; Ray and Bhattacharjee, 2007a). We can see that this variable emerges as a constant of the motion from the stationary limit of equation (3), and this constant, , can be identified closely with the matter flow rate, within a fixed geometrical factor. In terms of and , we can write this constant . On applying the perturbation scheme for and , the perturbation in , without losing anything of nonlinearity, is derived as
(6) 
in which,
(7) 
and . Equation (6) connects the perturbed quantities, , and , to one another. To get a relation between and only, we take equation (3), and apply the perturbation scheme on it. This results in
(8) 
To obtain a similar relationship solely between and , we combine the conditions given by equations (6) and (8) to get
(9) 
We stress at this point that in equations (6), (8) and (9), all orders of nonlinearity have been maintained. Now returning to equation (3), and taking its partial time derivative, an alternative form of the perturbation on the continuity condition appears as
(10) 
Equations (8) and (10) are equivalent expressions of the same principle, but its two distinct mathematical forms arise due to the axial symmetry of the flow, a feature that does not occur in spherical symmetry (Sen and Ray, 2014). Of the two forms, equation (10) is more useful than equation (8), when we take the secondorder partial time derivative of equation (5), and use the perturbation scheme on it to obtain
(11) 
We note the importance of using equation (10) in arriving at equation (11), which has the appearance of a similar equation of a nonlinear perturbation in the case of spherical symmetry (Sen and Ray, 2014). We exploit this similarity and apply all the mathematical methods of the spherically symmetric problem (Sen and Ray, 2014) to our present problem of a rotational flow. This is a crucial advantage.
Now making use of equation (9), its secondorder partial time derivative, and equation (10), we derive a fully nonlinear equation of the perturbation from equation (11), running in a symmetric form as
(12) 
in which,
(13) 
Going by the symmetry of equation (12), we can recast it in a compact form as
(14) 
with the Greek indices running from to , under the equivalence that stands for , and stands for . The derivation of equation (14) is pertinent to any kind of stationary background solution, with the only restriction being that the perturbation is radial. Further, equation (14), or equivalently, equation (12), is a nonlinear equation containing arbitrary orders of nonlinearity in the perturbative expansion. All of the nonlinearity is carried in the metric elements, . If we were to have worked with a linearized equation, then could be read from the matrix (Chaudhury et al., 2006; Ray and Bhattacharjee, 2007a),
(15) 
in which is the steadystate value of the local speed of sound.
Now, in Lorentzian geometry the d’Alembertian for a scalar field in curved space is expressed in terms of the metric, , as
(16) 
with being the inverse of the matrix implied by (Visser, 1998; Barceló et al., 2011). Comparing equations (14) and (16) with each other, we look for an equivalence between and . Clearly, equation (14) gives an expression for that is of the type given by equation (16). We extract the metrical part of equation (14) in the linear order, and the inverse of this metric implies an acoustic horizon, when . This approach is equivalent to the way in which an analogue metric can be fashioned out of a potential flow by converting its velocity field into the gradient of a scalar function, and then by perturbing the scalar function (Visser, 1998; Bilić, 1999; Das et al., 2007; Barceló et al., 2011). In contrast to this method of exploiting the potential character of the flow, our derivation of equation (14) makes use of the continuity condition. We argue that the latter method is more comprehensive because the continuity condition is based on matter conservation, which is a more forceful conservation principle than that of energy conservation (especially concerning open astrophysical flows), on which the scalarpotential approach is founded. If the flow were to have contained mechanisms for dissipation (as it is happens in models of axisymmetric accretion), then the potentialflow method would have been ineffective, but our method of making use of the matter conservation principle would still have delivered an equation of the perturbation.
A remarkable result of our analysis is that regardless of the order of nonlinearity that we may retain, the symmetric form of the Lorentzian metric equation remains unchanged, as we can see from equation (14). The preservation of this symmetry under arbitrary orders of nonlinearity is also exhibited in various other fluid systems like the hydraulic jump flow (Ray and Bhattacharjee, 2007b), the spherically symmetric outflows of nuclear matter (Sarkar et al., 2013), and the spherically symmetric inflows of astrophysical matter (Sen and Ray, 2014). While the analogue metric holds its ground in spite of nonlinearity, a serious consequence of including nonlinearity in the mathematical treatment is to compromise the notion of a static acoustic horizon in the flow. This is because a zeroorder description of , coming from equation (15), is no longer adequate. Instead, the elements, , are to be defined by equations (13), which carry timedependence in the higher nonlinear orders. This point of view agrees with a numerical study carried out by Mach and Malec (2008), who showed that sonic horizons would move about their static positions under strong perturbations, and the analogy between a sonic horizon and the event horizon of a black hole would appear limited. So the close correspondence between the physics of acoustic flows and many features of black hole physics is valid only in the linear order.
Iv Standing waves on subcritical inflows
All physically relevant stationary inflow solutions obey the outer boundary condition, as . Among these solutions, a critical inflow will reach the accretor with a high supercritical speed, but somewhere along the way, this flow will also undergo a discontinuity due to a shock. So critical inflows are highly subcritical at the outer boundary, are highly supercritical near the accretor (the inner boundary), and have a discontinuity in the interim region (Chakrabarti, 1989; Das, 2002; Das et al., 2003). There is, however, an entire class of inflow solutions that are globally subcritical, conforming to the inner boundary condition, as . For a gravitydriven evolution of an inflow solution to a critical state, the initial state of the evolution, as well as the intermediate states in the progression towards criticality, should realistically be subcritical. So the stability of globally subcritical flows has a significant bearing on how a critical solution will evolve eventually. Imposing an Eulerian perturbation on these subcritical inflows, their stability was studied under a linearized regime, and the amplitude of the perturbation, which was modelled as a standing wave, was seen to maintain a constant profile in time (Ray, 2003). In this respect we may say that the subcritical states are marginally stable. However, we have to be cautious in extending this argument when we consider nonlinearity in the perturbative effects, as it rightly ought to be done in a fluid flow problem.
Equation (12) gives a nonlinear equation of the perturbation, accommodating nonlinearity to any desired order. We design the perturbation to behave like a standing wave about a globally subcritical stationary solution, obeying the boundary condition that the spatial part of the perturbation vanishes at two radial points in the axisymmetric geometry, one at a great distance from the accretor (the outer boundary), and the other very close to it (the inner boundary). While the former boundary is a selfevident fact of the flow, there is a certain measure of difficulty in identifying the latter. The guiding principle behind the choice of the two boundaries is that the background stationary solution should be continuous in the interim region. For a globally continuous subcritical solution, that connects the outer boundary to the accretor, the inner boundary is obviously the surface of the accretor itself. If, however, even a subcritical inflow is disrupted by a shock, then the inner boundary should be the standing front of the shock itself. It is conceivable that no part of the perturbation on the background flow may percolate through the shock front, and so, the discontinuous front itself may be set as a boundary for the perturbation. Such piecewise continuity of a stability analysis, on either side of a discontinuity, is not uncommon in studies on fluid flows (Ray and Bhattacharjee, 2007b).
Our mathematical treatment involving nonlinearity is confined to the second order only (the lowest order of nonlinearity). Even simplified so, the entire procedure carries much of the complications associated with a nonlinear problem. The restriction of not going beyond the second order of nonlinearity implies that in equations (13) will contain primed quantities in their first power only. Taken together with equation (12), this will preserve of all the terms that are nonlinear in the second order. So, performing the necessary expansion of , and in equations (13), up to the first order only, and defining a new set of metric elements, , we obtain
(17) 
in which and are to be read just as in equation (14). The elements, , in equation (17) carry all the three perturbed quantities, , and . Our next task is to substitute both and in terms of , since equation (17) is over only. To make this substitution possible, first we have to use equation (6) to close in terms of and in all . While doing so, we ignore the product term of and in equation (6), because including it will raise equation (17) to the third order of nonlinearity. By the same token, we also have to take in equation (7). Once has been substituted in this manner, we write in terms of by invoking equation (8), with the reasoning that if and are both multiplicatively separable functions of space and time, with an exponential time part (all of which are standard prescriptions in any mathematical treatment on standing waves), then
(18) 
with being a function of only, a fact that simplifies much of the calculations to follow. The exact functional form of is determined from the way the spatial part of is prescribed, but on general physical grounds we reason that when , and are all real fluctuations, should likewise be real. We shall corroborate this argument independently in Section VI, but now we express the elements, , in equation (17), entirely in terms of as
(19) 
in all of which, has been introduced as a nonlinear “switch” parameter to track all the nonlinear terms. When , only linearity remains, and in this limit we converge to the familiar linear result implied by equation (15). In the opposite extreme, when , in addition to the linear effects, the lowest order of nonlinearity (the second order) becomes activated in equation (17), and the linearized stationary conditions of a sonic horizon get disturbed due to the nonlinear dependent terms. This very feature was observed numerically by Mach and Malec (2008). Equations (19) also contain the factors, , all of which read as
(20) 
Taking equations (17), (19) and (20) together, we finally obtain a nonlinear equation of the perturbation, completed up to the second order, without the loss of any relevant term.
To render equation (17) into a workable form, we first expand it in full and then divide it throughout by . In doing so, we also exploit the symmetry afforded by . The desirable form of the equation of the perturbation should be such that its leading term would be a secondorder partial time derivative of , with unity as its coefficient. To arrive at this form, an intermediate step will involve a division by , which, binomially, is the equivalent of a multiplication by , with a truncation applied thereafter. This is dictated by the simple principle that to keep only the secondorder nonlinear terms, it will suffice to retain just those terms which carry in its first power. The result of this entire exercise is
(21)  
in which, setting , what remains is the linear equation discussed in some earlier works (Ray, 2003; Chaudhury et al., 2006). We use a solution, , in equation (21), with being a real function (Polyanin and Zaitsev, 2004). Then we multiply the resulting expression throughout by and perform some algebraic simplifications by partial integrations to finally get
(22)  
in which the overdots indicate full derivatives in time. We integrate all spatial dependence out of equation (22) by using two boundary conditions, one very far from the accretor (at the outer boundary) and the other one (at the inner boundary) either very close to the accretor, or at a standing shock front where the background solution becomes discontinuous. At both of these boundary points, we constrain the perturbation to have a vanishing amplitude in time, while the background solution maintains a continuity in the interim region. The boundary conditions ensure that all the “surface” terms of the integrals in equation (22) will vanish (which is the reason for the tedious mathematical exercise to extract several “surface” terms). So after carrying out the required integration on equation (22), over the entire region trapped between the two specified boundaries, all that survives is the purely timedependent equation, having the form,
(23) 
in which the constants, , , and , are to be read as
(24) 
respectively. The form in which we have abstracted equation (23) is that of a general Liénard system (Strogatz, 1994; Jordan and Smith, 1999). All the terms in equation (23), bearing the parameter, , have arisen in consequence of nonlinearity. Setting , we do readily regain the known linear results (Ray, 2003), but to understand the role of nonlinearity in the perturbation, we have to look into the fully nonlinear import of the Liénard system in equation (23).
V Equilibrium conditions in the Liénard system
The mathematical form of a Liénard system is that of a damped nonlinear oscillator equation, going as (Strogatz, 1994; Jordan and Smith, 1999)
(25) 
in which, is a nonlinear damping coefficient (the retention of the parameter, , alongside , attests to its nonlinearity), and , is the “potential” of the system. Going by what equation (23) suggests, we realize that and , with the constant coefficients, , , and , having to be read from equations (IV). The equilibrium properties resulting from equation (25), can be known by decomposing this secondorder differential equation into a coupled firstorder dynamical system. To that end, we introduce a new variable, , and recast equation (25) as (Jordan and Smith, 1999)
(26) 
Equilibrium conditions follow when . For the coupled dynamical system in equations (V), we are led to two equilibrium points on the – phase plane. This is just how it should be because having accommodated nonlinearity up to the second order only, equations (V) will be quadratic in both and , yielding two equilibrium solutions. Labelling these equilibrium points with a superscript, we see that in one case, whereas in the other case, . So, one of the equilibrium points is located at the origin of the – phase plane, while the location of the other depends both on the sign and the magnitude of . In effect, both the equilibrium points lie on the line, , and correspond to the turning points of . Higher orders of nonlinearity will simply proliferate equilibrium points on the line, .
Having identified the position of the two equilibriums points, our next task is to examine their stability. So we subject both equilibrium points to small perturbations, and then carry out a linear stability analysis. The perturbation schemes on both and are and , respectively. Applying this scheme on equations (V), and then linearizing in and , will yield the coupled linear dynamical system,
(27) 
in which . Using solutions of the type, and , in equations (V), the eigenvalues of the Jacobian matrix of the dynamical system follow as
(28) 
with having to be evaluated at the equilibrium points. Knowing the eigenvalues, we can classify the stability of an equilibrium point by putting its coordinates in equation (28). The equilibrium point at the origin has the coordinates, . Using these coordinates in equation (28), we get the two roots of the eigenvalues as . If , then the eigenvalues will be purely imaginary quantities, and consequently, the equilibrium point at the origin of the – plane will be a centretype point (Jordan and Smith, 1999). And indeed, when the stationary inflow solution, about which the perturbation is constrained to behave like a standing wave, is subcritical over the entire region of the spatial integration, then , because in this situation, (Ray, 2003). Viewed in the – phase plane, the stationary solutions about this centretype fixed point at the origin, , look like closed elliptical trajectories, much in the manner of the phase solutions of a simple harmonic oscillator with conserved total energy. More to the point, these solutions correspond entirely to the solutions with unchanging amplitudes derived from a linear stability analysis of standing waves on subsonic flows (Ray, 2003). Thus, in a linear framework, a marginal sense of stability is insinuated by the centretype equilibrium point at the origin of the phase plane, because solutions in its neighbourhood are purely oscillatory in time, with no change in their amplitudes. While this conclusion can be made by a linearized analysis of the standing waves (Ray, 2003), it could be arrived at equally correctly by setting (the linear condition) in equation (28). An illustration of this special case is provided in Figure 1, which traces three phase solutions of the Liénard system. One of the solutions in this plot, obtained for and corresponding physically to the linear solution, is the closed elliptical trajectory about the centretype fixed point at .
Now, from dynamical systems theory, centretype points are known to be “borderline” cases (Strogatz, 1994; Jordan and Smith, 1999). In such situations, the linearized treatment will show marginally stable behaviour, but a robust stability or an instability may emerge immediately on accounting for nonlinearity (Strogatz, 1994; Jordan and Smith, 1999). This can be explained by a simple but generic argument. Close to the coordinate, , equations (V) can be approximated in the linear form as and , which, of course, gives a centretype point, just like a simple harmonic oscillator. Going further and accounting for the higherorder nonlinear terms, equations (V) can be viewed as a coupled dynamical system in the form, and . Such a system is said to be “reversible” if and , i.e., if (or ) is an odd function of , and (or ) is an even function of (Strogatz, 1994). Centretype points are robust under this reversibility requirement. A look at equations (V) immediately reveals that is not an even function of . Therefore, the centretype point obtained due to a linearized analysis of equations (V), is a fragile one. Ample evidence of this feature can be found in the behaviour of the spiralling solution in Figure 1.
The centretype point at the origin of the phase plane has confirmed the known linear results. However, all of that is the most that a simple linear stability analysis can bring forth. With nonlinearity in its lowest order, another equilibrium point is obtained, in addition to the centretype equilibrium point. This second equilibrium point is an outcome of the quadratic order of nonlinearity in the standing wave, and its coordinates in the phase plane are . Using these coordinates in equation (28), the eigenvalues become
(29) 
Noting as before, that , and that , and are all real quantities, the inescapable conclusion is that the eigenvalues, , are also real quantities, with opposite signs. In other words, the second equilibrium point is a saddle point (Jordan and Smith, 1999). The position of this equilibrium point is at the coordinate in the – phase portrait. The absolute value of the abscissa of this coordinate, , represents a critical threshold for the initial amplitude of the perturbation. If this amplitude is less than , then the perturbation will hover close to the linearized states about the centretype point, and stability shall prevail. The spiralling solution in Figure 1 gives a clear demonstration of this fact. If, however, the amplitude of the perturbation exceeds the critical value, i.e., if , then one enters the nonlinear regime, and in time the perturbation will undergo a divergence in one of its modes (for which has a positive root). This state of affairs has been depicted in the right side of the plot in Figure 1, showing a diverging phase solution. Since the eigenvalues, , have been yielded on using solutions of the type, , the folding time scale of this growing mode of the perturbation is , with having to be read from equation (29).
So, in the nonlinear regime, the simple fact that emerges is that stationary subsonic global background solutions will become unstable under the influence of the perturbation. In the vicinity of a saddle point, if the initial amplitude of the perturbation is greater than , then the solutions will continue to diverge, and higher orders of nonlinearity (starting with the third order in this case) will not smother this effect (Strogatz, 1994; Jordan and Smith, 1999). Since a saddle point cannot be eliminated by the inclusion of higher orders of nonlinearity (Jordan and Smith, 1999), all that we may hope for is that the instability may grow in time until it reaches a saturation level imposed by the higher nonlinear orders (but the instability will never be decayed down). This type of instability has a precedence in the laboratory fluid problem of the hydraulic jump (Volovik, 2006; Ray and Bhattacharjee, 2007b). While our discussion so far has dwelt on the perturbative perspective, the saddle point also has grave consequences for evolving a critical solution in the rotational flow through the dynamic process. There can be no critical solution without gravity driving the infall process. So, from a dynamic point of view, gravity starts the evolution towards the critical state from an initial subcritical state. In the absence of any analytical prescription for the fullblown nonlinear evolution, we have tried to get as close to the nonperturbative dynamics as possible, by the inclusion of progressively higher orders of nonlinearity in our perturbative treatment. This method has revealed to us a saddle point due merely to the second order of nonlinearity. We contend that if a saddle point is to be encountered in the realtime dynamics, then there should be serious obstacles in the way of reaching a stable and stationary critical end.
Vi Highfrequency travelling waves
Stability of fluids is also studied by constraining a perturbation to behave as a travelling wave. Applying the WKB approximation, some linearized studies (Ray, 2003; Chaudhury et al., 2006) established the stability of inviscid rotational accretion, with the perturbation on it being a highfrequency travelling wave. With nonlinearity lending additional effects, the stability of inviscid rotational accretion merits another look. To that end, we restructure equation (17), using the elements, , as they are given by equations (19), to get
(30) 
The foregoing equation has lost nothing in nonlinear terms. When we set in equations (19), the elements render equation (30) linear, which can then be worked upon by a solution of the form, , with and given by a converging power series (Ray, 2003; Chaudhury et al., 2006). We show presently how this linear solution is obtained, and thereafter, we view nonlinearity as a very weak effect about the linear condition, an argument by which we continue to use in equation (30) (Polyanin and Zaitsev, 2004; Bridges and LainePearson, 2004). And since nonlinearity is now feeble, we retain only its second order in equation (30). What results from it eventually, after a long series of algebraic steps, is
(31) 
in which and are operators, given by
(32) 
with operating on the linear part of equation (31), and on its nonlinear part. Obviously, equation (31) is a nonlinear ordinary differential equation in , to solve which, we first write , so that . With this form of , we expand equation (31) as
(33)  
Next, we set down as a power series in the form,
(34) 
The principle of the WKB approximation is that successive terms in equation (34) follow the condition , i.e., the power series given by converges rapidly as increases. To facilitate this outcome, we prescribe a highfrequency travelling wave, such that its wavelength, , is much smaller than any characteristic length scale in the fluid. The smallest critical radius in the multicritical accreting system is a natural choice for such a length scale.
To find a solution of equation (33) a first step is to apply the WKB approximation to its linear limit, i.e. when . In this case, we replace all in by , with the latter implying the linear solution in the series of . After that, on making use of the series given by in the linear portion of equation (33), we obtain the two successive highest order terms going as and . Collecting all the coefficients of , summing them up, and setting the sum equal to zero, give
(35) 
which is a simple quadratic, that leads to
(36) 
Working likewise with the coefficients of , gives
(37) 
from which, on making use of equation (36), we extract a solution of as
(38) 
with being a constant of integration. Both and give the two leading terms in the series of . From equations (36) and (38), respectively, these terms go asymptotically as and , given the condition that on large length scales, while approaches its constant ambient value. The next term in the series, , behaves asymptotically as . So, under the regime of a highfrequency travelling wave, all of this implies . Therefore, it suffices to consider the two leading terms only, involving and in the series expansion of .
Equations (36) and (38) give the leading linear solution of equation (33) when . To incorporate the effect of nonlinearity, i.e., when , we employ an iterative technique on equation (33), and collect the nonlinear contribution to and only. We achieve this through the following logical sequence, remembering that nonlinearity is a small effect about the linear solution.

In equation (33), the nonlinear dependent part carries the timedependent factor, . We expand it as a power series in , going as , and apply a truncation immediately after the zerothorder term. Otherwise, retention of any order higher than unity will allow the series to make a timedependent contribution to , in violation of the restriction that has to selfconsistently bear only a spatial dependence in the feebly nonlinear regime.

In the dependent part of equation (33), we replace all nonlinear by the linear , in accordance with our iterative principle, with a truncation applied immediately after .

Thereafter, in the dependent part of equation (33), the function, , is to be read as, , which we finally expand to . We have truncated the series expansion beyond the second order, because the higher orders will be irrelevant for the series, insofar as our objective is to determine the contribution of nonlinearity only to the terms carrying and .
Now, with the involvement of nonlinearity, the sum of the coefficients of in equation (33), set to zero, will read as
(39) 
in which,