A&A 449, 861-868 (2006)
DOI: 10.1051/0004-6361:20054098

Perturbations in the self-similar Bondi flow

J. Gaite

Instituto de Matemáticas y Física Fundamental, CSIC, Serrano 113bis, 28006 Madrid, Spain

Received 24 August 2005 / Accepted 12 November 2005

Abstract
Context. The question of the stability of steady spherical accretion has been studied for many years and, recently, the concept of spatial instability has been introduced.
Aims. Here we study the perturbations of steady spherical accretion flows (Bondi solutions), restricting ourselves to the case of a self-similar flow, as a case that is amenable to analytic treatment and that is of physical interest. We restrict ourselves further to its acoustic perturbations.
Methods. The radial perturbation equation can be solved in terms of Bessel functions. We study the formulation of adequate boundary conditions and decide to use no matter-flux-perturbation conditions (at the Bondi radius and at r=0). We also consider the problems of initial conditions and time evolution, and concentrate more particularly on radial perturbations.
Results. No spatial instability at r=0 is found. The time evolution is such that perturbations eventually become ergodic-like and show no trace of instability or of acquiring any remarkable pattern.

Key words: accretion, accretion disks - hydrodynamics

1 Introduction

The Bondi solutions for stationary spherical accretion onto a compact object (Bondi 1952) have been the basis for many analytical (or semi-analytical) studies of spherical accretion. Since there is a set of solutions, selecting the most adequate among them has been a moot point. Bondi suggested that stability criteria would be useful for this purpose. However, neither analytical studies (Garlick 1979; Petterson et al. 1980) nor numerical studies (Ruffert 1994) have detected any instability, in the sense that linear modes have amplitudes that increase with time. Garlick's argument is both very general and powerful, because it is based on energy conservation. Kovalenko & Eremin (1998) point out another possible type of instability, namely, one in space rather than in time. They state that the nature of this instability is given by a perturbation whose amplitude relative to the unperturbed quantity grows infinitely with decreasing radius, so that one expects that it reaches the nonlinear regime.

Apart from this spatial notion of stability, Garlick (1979) had proposed earlier "to take a more sophisticated approach to the stability problem'' (in the time domain): the perturbation, independent of its global behaviour, could be increasing in a small region (anywhere along the flow). However, he ruled this out as "unphysical''. Garlick (1979) also considered the initial value problem briefly.

Our purposes in this paper are to pay special attention to the problem of the time evolution of the perturbations from Bondi accretion and to go as far as we can with analytic methods. Given that the general problems of initial conditions and time evolution are intrinsically difficult, we are compelled to make a number of simplifications. First of all, we consider a particularly simple form of Bondi accretion: self-similar Bondi accretion. It occurs for the value of the polytropic index $\gamma= 5/3$, a limit case with physical interest (or it holds approximately as an intermediate asymptotic limit for ${\gamma < 5/3}$). As an additional simplification, we restrict ourselves to acoustic perturbations, neglecting entropic perturbations and vorticity. This is consistent, as generation of entropy and vorticity are related and are simply advected with the flow and eventually vanish (Garlick 1979). Moreover, to carry out a complete study of time evolution, we restrict ourselves to radial perturbations. In this case, the analytic treatment is simpler, as is the graphical representation of the results.

One possible rôle of entropy and vorticity perturbations has been explored recently (Foglizzo 2001). However, in accord with the traditional view, no instabilities can arise in polytropic flow in the absence of an external source of entropy or vorticity. We also do not consider any external sources.

We divide the paper into three main sections. The following section is a review of Bondi accretion and it introduces the self-similar solutions. The next section is devoted both to the study of perturbations, in particular acoustic modes, and to formulating the appropriate boundary problem of differential equations. We follow previous work but try to be more systematic. Unfortunately, the boundary problem does not lead to a Sturm-Liouville problem, even though it is connected with one. The last main section focuses on radial perturbations, as do the papers of Petterson et al. (1980) and of Theuns & David (1992). In that section, we precisely formulate and then solve the corresponding eigenvalue problem of ordinary differential equations (ODE's), finding an orthogonal system of eigenfunctions, thereby solving the problem of initial conditions. We end by discussing the relevance of our results for the problem of stability in Bondi accretion. An appendix explains the asymptotic form of small wave-length perturbations (WKB approximation).

   
2 Steady spherical accretion: the Bondi solutions

The problem of steady spherical accretion is suitable for analytical treatment and is the basis for more complicated and realistic problems (Frank et al. 1992; Shore 1992). A compact object (e.g., a star) of mass M is at rest at the origin (r=0). An infinite cloud of gas, which is at rest at infinity and has density  $\rho_\infty$ and pressure $p_\infty$, is accreted by the compact object. The motion of the gas is assumed to be spherically symmetric and steady, so the fluid partial differential equations (PDE's) become ODE's in the radius r. Furthermore, the change of mass of the compact object, the radiation effects, and the viscosity are neglected. Then the equations describing the gas are

  

$$\frac{{\rm d}}{{\rm d}r}(r^2 \rho {v}) = 0,$$ (1)

$$v \frac{{\rm d} v}{{\rm d}r} = - \frac{1}{\rho}\frac{{\rm d} p}{{\rm d}r} - \frac{GM}{r^2},$$ (2)

with $\rho$ the gas density field, v the gas velocity field, p the gas pressure, and G Newton's gravitational constant. The continuity Eq. (1) can be integrated immediately. We must add an equation of state for the fluid, which we choose to be polytropic, following Bondi (1952). Then we have a set of three equations:

   

$$4 \pi r^2 \rho v = {\cal K}_1,$$ (3)

$$v \frac{{\rm d} v}{{\rm d} r} = -\frac{GM}{r^2}- \frac{1}{\rho} \frac{{\rm d} p}{{\rm d} r},$$ (4)

$$p = {\cal K}_2~ \rho^{\gamma},$$ (5)

with ${\cal K}_1$ and ${\cal K}_2$ constants (characterising the flow), and the polytropic exponent $1 \leq \gamma \leq 5/3$. With a suitable choice of $\gamma$, the polytropic equation is equivalent to assuming that there is no heat transfer (adiabatic flow). The speed of sound, to be used later, is given by $c^2={\rm d} p/{\rm d}\rho = {\cal K}_2 \gamma~\rho^{\gamma-1}$.

By integrating the Euler Eq. (4) over r, one gets Bernouilli's equation. Furthermore, the polytropic Eq. (5) allows one to integrate ${\rm d}p$, so the problem boils down to two algebraic equations, which can be expressed in non-dimensional form by using the sound velocity in the gas at infinity $c_\infty$. Let us define the non-dimensional variables

 

$$x = \frac{r}{{\cal R}},$$ (6)

$$y = \frac{\rho}{\rho_\infty},$$ (7)

$$u = \frac{v}{c_\infty},$$ (8)

$$\lambda = \frac{{\cal K}_1}{4 \pi {\cal R}^2 \rho_\infty c_\infty};$$ (9)

that is, r is measured in units of the accretion (or Bondi) radius, ${\cal R}= GM/c^2_\infty$, v in units of the speed of sound $c_\infty$, and $\rho$ in units of  $\rho_\infty$. The non-dimensional form of the equations of motion is:

$$x^2 y u = \lambda,$$ (10)

$$\frac{u^2}{2} + \frac{y^{\gamma-1}-1}{\gamma-1} = \frac{1}{x}\cdot$$ (11)

Solving for u in the first equation and substituting in the second, we obtain

$$% \frac{\lambda^2}{2}x^{-4} y^{-2} + n\left(y^{1/n}-1\right)= x^{-1},$$ (12)

where we have introduced the adiabatic index $n = 1/(\gamma-1)$. As an equation for x, it is a quartic algebraic equation, so it admits an algebraic expression for its solutions (however complex, see Theuns & David 1992).

   
2.1 Self-similar Bondi solutions

The preceding algebraic relation between y and x (12) can be simplified in some cases. Let us assume first that the radial distance is small, namely, that $x^{-1} \gg 1$, so that the term nin Eq. (12) can be neglected (in comparison with x^-1). The condition that the density is large, i.e., $y \gg 1$, is equivalent (assuming that $\gamma$ is not too close to one). Then we have

$$% \frac{\lambda^2}{2}x^{-3} y^{-2} + n~y^{1/n} x= 1.$$ (13)

For convenience, we define the variable z = y^1/n x and rewrite this equation as

$$% \frac{\lambda^2}{2}x^{-3+2n} z^{-2n} + n~ z= 1.$$ (14)

This is still a complicated relation between x and z, so we must look for further simplification.

We see that, in the particular case n=3/2, Eq. (14) is just an algebraic equation in the sole variable z, with a solution that only depends on the non-dimensional accretion rate $\lambda$. Therefore, in this case, the non-dimensional density takes the form $y = z(\lambda)^{3/2}~ x^{-3/2}$ and the non-dimensional velocity takes the form $u = \lambda ~z(\lambda)^{-3/2}~ x^{-1/2}$. We see that both density and velocity are power laws of the radius. More complicated solutions of Eq. (14) can be considered but we will focus on the case $n=3/2 \Rightarrow\gamma=5/3$, which provides simple solutions and is realistic since it corresponds to the adiabatic flow of perfect monoatomic gases.

Then we write

  

$$\rho(r) = \alpha~ r^{-3/2},$$ (15)

$$v(r) = \beta ~r^{-1/2},$$ (16)

with $\alpha$ and $\beta$ such that

$$% \alpha \beta = \frac{{\cal K}_1}{4 \pi}\cdot$$ (17)

In this case, the speed of sound is

$$% c^2 (r)= {\cal K}_2~ \frac{5}{3}~\alpha^{2/3} r^{-1}= \sigma^2 r^{-1},$$ (18)

so the Mach number is given by

$$% {\cal M} = \frac{v(r)}{c(r)}= \frac{\beta}{\sigma},$$ (19)

which is constant. Thus each solution is characterised by either $\{\alpha, \beta, \sigma\}$ or $\{{\cal K}_1,{\cal K}_2, {\cal M}\}.$

We note that self-similar Bondi solutions have no discontinuities and are either subsonic or supersonic everywhere. On account of the boundary condition at infinity, it is natural to discard the supersonic self-similar solutions.

2.1.1 Intermediate asymptotic limit when ${\gamma < 5/3}$

It is commonly observed that nonlinear equations that do not have everywhere valid similarity solutions can nevertheless develop them in an intermediate asymptotic regime. This is a regime that takes place between two very different scales in such way that both scales can be neglected (Barenblatt 1996). In the problem of steady spherical accretion by a compact object, we have two fundamental scales, namely, the Bondi radius ${\cal R}$ and the sonic radius. The latter is defined by the implicit equation

$$% {r_{\rm s}} = \frac{GM}{2 c^2(r_{\rm s})}$$ (20)

(Frank et al. 1992; Shore 1992). This equation can be written with non-dimensional variables as

$$% {z_{\rm s}(x_{\rm s})} = y_{\rm s}(x_{\rm s})^{1/n} x_{\rm s} = \frac{1}{2},$$ (21)

where $y_{\rm s}(x_{\rm s})$ is the solution of Eq. (12) corresponding to $x_{\rm s}$. Writing Eq. (12) in terms of $x_{\rm s}$ and $z_{\rm s} = 1/2$, we have

$$% {2^{2n-1}}\lambda^2~x_{\rm s}^{-3+2n} = n~x_{\rm s} + 1 -\frac{n}{2}\cdot$$ (22)

If $2 > n > 3/2 \Leftrightarrow3/2 < \gamma< 5/3$ this equation has zero, one, or two solutions, depending on the value of the non-dimensional accretion rate $\lambda$ (given n). As $\lambda$ increases, the number of solutions increases as well, being one at the critical value  $\lambda_{\rm max}$, which is only a function of n. This critical value is called $\lambda_{\rm max}$ because it is the maximal value such that the flow is continuous (Theuns & David 1992). The value

$$\begin{eqnarray*}x_{\rm s}(\lambda_{{\rm max}}) = \frac{1}{2} - \frac{3}{4n} \end{eqnarray*}$$

is only a function of n and approaches zero as $n \rightarrow3/2$.

If n > 3/2 but not much larger ( ${\gamma < 5/3}$ but not much smaller), and $\lambda$ is not much smaller than  $\lambda_{\rm max}$, the solution of Eq. (14) for z will have a dependence on x but only a logarithmic one, such that the power laws (15) and (16) will hold with corrections proportional to  $\log (r/r_{\rm s})$. In other words, in the intermediate asymptotics $r_{\rm s} \ll r \ll {\cal R}$, the power laws are approximately correct (so the Mach number ${\cal M}$ hardly grows with the flow).

For precision, we will set $\gamma= 5/3$ in the following, keeping in mind that the results have somewhat wider applicability.

   
3 Linear perturbations of the self-similar Bondi solutions

Let us assume that the steady spherically symmetric solutions are perturbed, such that the total velocity and density fields become $\rho(r)+ \delta \rho({\mbox{\boldmath {$x$ }}},t)$ and ${\mbox{\boldmath {$v$ }}}(r) + \delta {\mbox{\boldmath {$v$ }}}({\mbox{\boldmath {$x$ }}},t)$, with $\rho(r)$ and ${\mbox{\boldmath {$v$ }}}(r)=v(r) \hat {\mbox{\boldmath {$r$ }}}$ the self-similar solutions, given by Eqs. (15) and (16) in the self-similar case ( $\hat{\mbox{\boldmath {$r$ }}}$ is the unit radial vector).

The independent perturbations are $\delta\rho$, $\delta{\mbox{\boldmath {$v$ }}}$ and the "entropy perturbation''  $\delta{\rm s}$, where the "entropy'' is defined as

$$\begin{eqnarray*}s = \log \frac{c^2}{\gamma~\rho^{\gamma-1}} \end{eqnarray*}$$

(Kovalenko & Eremin 1998). We shall assume that the perturbations are such that $\delta{\rm s} = 0$, so that they keep the polytropic equation of state. In particular, they may be adiabatic perturbations; we shall henceforth denote thermodynamic functions as if this is the case. Then the equations for linear perturbations are

  

$$\frac{\partial \delta \rho}{\partial t} + \nabla \cdot (\delta \rho \; {\mbox{\boldmath {$v$ }}} + \rho \; {\delta {\mbox{\boldmath {$v$ }}}}) = 0,$$ (23)

$$\frac{\partial {\delta {\mbox{\boldmath {$v$ }}}}}{\partial t} + ... ...ldmath {$v$ }}}\cdot\nabla)\; {\mbox{\boldmath {$v$ }}} = - \nabla \delta \chi,$$ (24)

where

$$% \delta \chi = \frac{c^2 ~\delta \rho}{\rho}$$ (25)

is the enthalpy perturbation (the speed of sound is given by $c^2 = {\cal K}_2 \gamma\rho^{\gamma-1}$). Therefore, we have a system of four linear PDE's.

   
3.1 Solution for acoustic modes

For the acoustic modes, $\delta{\rm s} = 0$ and the perturbation in the velocity field can be written as $\delta {\mbox{\boldmath {$v$ }}}= \nabla \phi.$ In this case, and given that Bondi flows are curl free, we can write Eq. (24) as the following differential equation for the scalar potential field $\phi$

 

$$% \nabla \frac{\partial {\phi}}{\partial t} + \nabla ({\mbox{\boldmath {$v$ }}} \cdot \nabla \phi) = - \nabla \delta \chi.$$ (26)

This equation can be integrated to obtain (Kovalenko & Eremin 1998)

 

$$% \frac{\partial {\phi}}{\partial t} + {\mbox{\boldmath {$v$ }}} \cdot \nabla \phi = - \delta \chi,$$ (27)

which provides $\delta \chi$ and hence  $\delta\rho$ in terms of $\phi$. Equation (23) can then be written in terms of $\phi$ as follows

$$% \left[ \frac{\partial}{\partial t} + {\mbox{\boldmath {$v$ ... ...right) \phi = \frac{c^2}{\rho} \nabla \cdot (\rho \nabla \phi)$$ (28)

[Eq. (24) of Kovalenko & Eremin (1998)]. This is an equation for non-homogeneous wave propagation.

One tries a solution of Eq. (28) for the scalar potential $\phi$ by separating the variables in spherical coordinates, that is, a solution in the (complex) form

 

$$% \phi(r,\theta,\varphi,t)= R(r)~ Y_{lm}(\theta,\varphi)~{\rm e}^{-{\rm i}\omega t}.$$ (29)

Given this form of the scalar potential, the velocity perturbation is given by

$$\begin{eqnarray*}\delta{\mbox{\boldmath {$v$ }}} = \left(\delta v_r, \delta v_\theta, \delta v_\varphi\right), \end{eqnarray*}$$

with

$$\delta v_r (r,\theta,\varphi,t) = {\rm Re}\left[{\rm e}^{-{\rm i}\omega t} \frac{{\rm d}R(r)}{{\rm d}r} Y_{lm}(\theta,\varphi)\right],$$ (30)

$$\delta v_\theta (r,\theta,\varphi,t) = {\rm Re}\left[{\rm e}^{-{\rm i}\omega t}\frac{R(r)}{r} \partial_\theta Y_{lm}(\theta,\varphi)\right],$$ (31)

$$\delta v_\varphi (r,\theta,\varphi,t) = {\rm Re}\left[{\rm e}^{-{... ...ga t}\frac{R(r)}{r \sin \theta} \partial_\varphi Y_{lm} (\theta,\varphi)\right]$$ (32)

(after taking the real part of the complex quantities). Note that the angular components form a poloidal vector field and the form of the radial component ensures the absence of vorticity (Garlick 1979).

The differential equation satisfied by R(r) is

 

$$% \left(\sigma^2-\beta^2\right)r^2 R'' + \left(\frac{\sigma^2-\beta^2}{2} + 2 {\rm i} \beta \omega r^{3/2}\right) r R'$$

$$+ \left[-l(l+1)\sigma^2 + \omega^2 r^{3} + {\rm i} \beta \omega r^{3/2}\right] R = 0.$$ (33)

Following standard methods of the theory of ODE's, we write the second order equation as

 

R''(r)+ P(r) R'(r) +Q(r) R(r) = 0, (34)

with

$$P(r) = \frac{1}{2r} +\frac{2{\rm i}\beta\omega r^{1/2}}{\sigma^2-\beta^2}$$ (35)

$$Q(r) = \frac{\omega^2 r - l(l+1) \sigma^2 r^{-2} + {\rm i}\beta\omega r^{-1/2}}{\sigma^2-\beta^2}\cdot$$ (36)

The ODE can be simplified by the following change of variable:

$$% \log R(r) = \log h(r) - \frac{1}{2}\int {\rm d}r P(r),$$ (37)

which integrated yields

$$% R(r) = r^{-1/4} \exp{\frac{-2{\rm i}\beta \omega r^{3/2}}{3\left(\sigma^2-\beta^2\right)} }~ h(r).$$ (38)

In terms of this new dependent variable h(r), the differential Eq. (33) becomes

 

$$% h''(r) + \left(a_1 r^{-2} + a_2 r\right) h(r) = 0,$$ (39)

with

$$a_1 = -\frac{3}{16} - \frac{l(l+1)\sigma^2}{\sigma^2-\beta^2},$$ (40)

$$a_2 = \left(\frac{\sigma\omega}{\sigma^2-\beta^2}\right)^2\cdot$$ (41)

The preceding differential Eq. (39) has one regular singular point at r=0 and one irregular singular point at $r=\infty$. In fact, it can be transformed into the Bessel equation by making r^3/2 the dependent variable. Thus we obtain

 

$$% R(r) = r^{1/4} {\rm e}^{-{\rm i} \mu r^{3/2}} \left[C_1 J_\nu \left(\kappa r^{3/2}\right) + C_2 N_\nu\left(\kappa r^{3/2}\right)\right],$$ (42)

where

$$\begin{eqnarray*}&& \mu =\frac{2 \beta \omega}{3\left(\sigma^2-\beta^2\right)},\... ...\sqrt{[1+16l(l+1)]\sigma^2 -\beta^2}}{6\sqrt{\sigma^2-\beta^2}}, \end{eqnarray*}$$

and C₁,C₂ are, in general, two arbitrary complex numbers.

Note that the order $\nu$ of the Bessel function in solution (42) does not depend on $\beta$ and $\sigma$ separately but on its ratio, that is, the Mach number, ${\cal M} = \beta/\sigma$ (in addition to depending on the angular number l). If the values of ${\cal M}$ and l are generic such that $\nu$ is non-integer, we can write the solution of the equation as

 

$$% R(r) = r^{1/4} {\rm e}^{-{\rm i} \mu r^{3/2}} \left[C_1 J_\nu \left(\kappa r^{3/2}\right) + C_2 J_{-\nu}\left(\kappa r^{3/2}\right)\right].$$ (43)

In this form, the most general (complex) solution of the scalar potential is

 

$$% \phi(r,\theta,\varphi,t) = {\rm e}^{-{\rm i}\omega t} \; Y_{lm}... ...u \left(\kappa r^{3/2}\right) + C_2 J_{-\nu}\left(\kappa r^{3/2}\right)\right).$$ (44)

Since we assume ${\cal M}< 1$, the signs of $\mu$ and $\kappa$ are positive and $\nu$ is real.

Boundary conditions are needed in order to determine both C₁ and C₂ and the $\omega$-spectrum. Notice that the constants C₁ and C₂ need to be specified for each solution, that is, they are indexed by l,m, and the $\omega$-spectrum index.

   
3.2 Boundary conditions

The acoustic perturbations satisfy the wave Eq. (28), so a suitable boundary condition on a surface is needed. After the separation of variables in spherical coordinates, one needs to impose two boundary conditions on the radial differential equation, as this equation is of the second order (Morse & Feshbach 1953). Thus one assumes that the boundary conditions adapt to the spherical coordinates and, in particular, that the boundaries are two spherical surfaces of inner radius r₁ and outer radius r₂. The natural inner radius is either the accretor or the sonic radius. The natural candidate for outer radius is the Bondi (or accretion) radius ${\cal R}$: as r gets larger than ${\cal R}$, the thermal energy of the nearly homogeneous gas dominates over its gravitational energy and we assume that this homogeneous gas remains unperturbed.

The choice of boundary conditions of Petterson et al. (1980) (for radial perturbations) is that the flux perturbation $\delta f$ should vanish at both r₁ and r₂; since the flux of matter through a sphere of radius r is constant in the Bondi solutions (Eq. (3)), it seems natural to impose the condition that the perturbations do not change it on the inner and outer boundaries. Unlike for radial modes, for general acoustic modes the variation of the matter flux current $\mbox{\boldmath {$j$ }} = \rho{\mbox{\boldmath {$v$ }}}$, namely, $\delta \mbox{\boldmath {$j$ }} = \rho ~\delta {\mbox{\boldmath {$v$ }}} + \delta \rho ~{\mbox{\boldmath {$v$ }}}$, depends on the angles  $\theta,\varphi$. The perturbed flux of matter through a sphere of radius r is given by the angular integral of its radial component, proportional to  $\int {\rm d}\Omega~Y_{lm}(\theta,\varphi)$. This angular integral vanishes if $l \neq 0$, so this condition is only useful for radial modes.

Instead, Kovalenko & Eremin (1998) derive their boundary conditions from the law of the conservation of acoustic energy

$$% \frac{\partial E}{\partial t} = \nabla\cdot \mbox{\boldmath {$W$ }},$$ (45)

where the acoustic energy density is

$$% E = \frac{\rho}{2c^2} \left[\left(\frac{\partial\phi}{\par... ...\phi)^2 - (\mbox{\boldmath {$v$ }}\cdot\nabla \phi)^2 \right],$$ (46)

and the acoustic energy flux current is

$$% \mbox{\boldmath {$W$ }} = \frac{\rho}{c^2} \frac{\partial\p... ... \right] = - \partial_t\phi \; \delta \mbox{\boldmath {$j$ }}.$$ (47)

The acoustic energy flux through a spherical surface centred on the origin is given by

$$% r^2 \int {\rm d}\Omega\;W_r = \frac{r^2\rho}{c^2}\int {\rm... ...l_t \phi - {c^2} \partial_r\phi + v^2\partial_r \phi \right].$$ (48)

This quantity is quadratic in $\phi$ and does not necessarily vanish for any l upon performing the angular integral.

We see from Eq. (48) that sufficient conditions for the vanishing of the acoustic energy flux are either

$$% \partial_t \phi = 0$$ (49)

or

$$% r^2~\hat {\mbox{\boldmath {$r$ }}} \cdot \delta \mbox{\bold... ... \phi - {c^2} \partial_r\phi + v^2\partial_r \phi \right] = 0.$$ (50)

The latter condition implies, of course, the vanishing of its integral $\delta f$ over any angular domain, so it is the appropriate generalization of the condition of Petterson et al. (1980). Hence, we choose it as our boundary condition; namely, this quantity should vanish at both r₁ and r₂. It can be expressed only in terms of R(r) as

 

$$% \frac{r^2\rho}{c^2}\left[{\rm i} \omega v R + \left(c^2 - v^2\right) R'\right] = 0.$$ (51)

The discussion above is general, that is, it is not restricted to self-similar solutions. For self-similar solutions, the accretor and sonic radii vanish, so we take r₁ = 0; on the other hand, the Bondi radius  ${\cal R} \rightarrow\infty$. However, we need a finite outer radius to have a discrete spectrum, so we take $r_2 = {\cal R}$ to be finite, keeping in mind that the actual spectrum becomes continuous in the limit  ${\cal R} \rightarrow\infty$.

The radial Eq. (33) for R(r) is not self-adjoint (the subject of adjointness of ODE's and the Sturm-Liouville boundary problem is treated, for example, by Morse & Feshbach 1953; and Coddington & Levinson 1955). However, the transformed Eq. (39) for h(r) is patently self-adjoint. The boundary condition (51) in terms of h(r)reads

 

r^-3/4[-h + 4 r h'(r)] = 0 (52)

(at r₁ and r₂). Therefore, the boundary problem for h(r) is of Sturm-Liouville type, and the eigenvalues and eigenfunctions fulfill the corresponding properties; in particular, the eigenfunctions are orthogonal (Morse & Feshbach 1953). This is natural, because h(r)is a combination of Bessel functions (with a factor r^1/2). In contrast, the boundary problem for R(r) is not of Sturm-Liouville type and the eigenfunctions need not be orthogonal (note that the $\kappa$ dependence of $\mu$ in Eq. (42) spoils the orthogonality inherent to the Bessel functions). Having an orthogonal system of eigenfunctions is important for the problem of initial conditions, namely, for finding the coefficients of the expansion corresponding to the initial data.

We now apply Eq. (52) at r₁=0 and r₂. For the former, we need to take the limit $r \rightarrow0$ of Eq. (52). Writing it in terms of Bessel functions, with the dependent variable $z= \kappa r^{3/2}$, we have

$$% \lim_{z \rightarrow0} z^{-1/6}\big(C_1\left[J_\nu(z) + 6 z J'_\nu(z) \right] + C_2\left[J_{-\nu}(z) + 6 z J'_{-\nu}(z) \right]\big) = 0.$$ (53)

In this limit, it is sufficient to consider the first term of the Taylor expansion of the Bessel functions, namely,

$$% z^{-1/6} \left[J_\nu(z) + 6 z J'_\nu(z)\right] \approx \frac{2^{-\nu}}{\Gamma(\nu + 1)} (1+ 6 \nu) z^{\nu-1/6}.$$ (54)

When $\nu > 1/6$ (l>0) we are led to impose C₂ = 0 to remove negative exponents (divergent as $r \rightarrow0$). If $\nu = 1/6$ (corresponding to the radial modes), we must take C₁ = 0 instead.

The outer boundary condition at r₂ yields

 

$$% J_\nu (z_2) + 6 z_2 J'_\nu (z_2) = 0,$$ (55)

which provides the allowed eigenvalues $\kappa_n$, hence the $\omega_n$.

Finally, it is important to analyse the relative variations in the density and the velocity. Their behaviour in the limit $r \rightarrow0$ is:

   

$$\frac{\delta v_r}{v} \sim r^{1/2} R'(r) \sim r^{-1/4 + 3\nu/2},$$ (56)

$$\frac{\delta v_\theta}{v} \sim \frac{\delta v_\varphi}{v} \sim r^{1/2} \frac{R(r)}{r} \sim r^{-1/4 + 3\nu/2},$$ (57)

$$\frac{\delta\rho}{\rho} = \frac{1}{c^2}\left[i\omega~R(r) - v R'(r)\right] \sim r^{-1/4 + 3\nu/2},$$ (58)

for $\nu > 1/6$; for $\nu = 1/6$, ${\delta v_r}/{v} \sim {\delta\rho}/{\rho} \sim r$. All these relative variations vanish as $r \rightarrow0$, so there are no spatial instabilities at the inner radius and the linear perturbation theory is sound.

4 Radial perturbations

The simplest acoustic modes are the radial perturbations, studied by Petterson et al. (1980) and Theuns & David (1992), for general Bondi flow. Indeed, a more thorough analytical treatment is also possible in the self-similar case if we restrict ourselves to radial perturbations. Petterson et al. (1980) obtained general stability results: the problem of standing waves with vanishing flux boundary conditions only has real eigenvalues of $\omega$, and radially travelling waves can be shown not to grow (by an energy argument). Theuns & David (1992) derived an associated Sturm-Liouville boundary problem by decomposing the flux perturbation $\delta f$ as if it were composed of modulated waves with a position-dependent phase velocity. This decomposition leads to a self-adjoint equation for the amplitude. Their procedure is analogous to the transformation that led us from the radial Eq. (33) to the self-adjoint Eq. (39).

Although the associated Sturm-Liouville boundary problem yields the eigenfunctions of the original boundary problem, it does not provide us with an orthogonality relation for them, which may exist or not. Fortunately, a specific mathematical treatment of the first order equations for radial perturbations allows us to find a definite scalar product under which the eigenfunctions are orthogonal, as explained in Appendix A. Next, we show the solution to the first order boundary problem, with the orthogonality relation, and then proceed to solve the problem of initial conditions to study the time evolution.

4.1 Solution of the boundary problem

The general linear perturbation Eqs. (23) and (24) restricted to radial modes give

  

$$-{\rm i}\omega~\delta \rho + \frac{1}{r^2} \partial_r \left[ r^2 (\delta \rho \; v + \rho \; \delta v)\right] = 0,$$ (59)

$$-{\rm i}\omega~\delta v + \partial_r\left(\frac{c^2}{\rho}~\delta\rho+ v~{\delta v}\right) = 0.$$ (60)

We denote $y = (\delta f, \delta g)$, where $\delta f = r^2(\delta \rho ~ v + \rho ~ \delta v)$is the perturbation of the flow (but for the factor $4\pi$) and $\delta g = (c^2/\rho)~\delta\rho+ v~{\delta v} = {\rm i}\omega\phi$ according to Eq. (60). Solution (43) with $\nu = 1/6$ and C₁ = 0(because of the boundary condition at r=0) gives

  

$$\delta f(r) = \exp\left(-{\rm i}{\cal M}\kappa r^{3/2}\right) r^{5/4}~ {\rm i} ~ J_{5/6}\left(\kappa r^{3/2}\right),$$ (61)

$$\delta g(r) = \exp\left(-{\rm i}{\cal M}\kappa r^{3/2}\right) \frac{\sigma}{\alpha}~ r^{1/4} ~J_{-1/6}\left(\kappa r^{3/2}\right),$$ (62)

where we have used the identity J_{- 1/6}(z) + 6 z J'_{- 1/6}(z) = - 6 z J_5/6(z) to simplify $\delta f$ and divided both components by the constant total flux so that $\delta f$ becomes non-dimensional. Now we impose the boundary condition $\delta f = 0$ at r=1 (we normalize the radial variable by dividing it by $r_2 = {\cal R}$). It determines that the possible values of $\kappa$ are the zeros of the Bessel function J_5/6 and so yields the eigenvalues $\kappa_n$. Redefining the time unit such that $t \rightarrow\sigma t$, we have the modes

  

$$\delta f_n(r,t) = \exp\left[-{\rm i}\kappa_n\left({\cal M}r^{3/2}... ...\right)t\right)\right] r^{5/4}~ {\rm i} ~ J_{5/6}\left(\kappa_n r^{3/2}\right),$$ (63)

$$\delta g_n(r,t) = \exp\left[-{\rm i}\kappa_n\left({\cal M}r^{3/2}... ...t)\right] \frac{\sigma}{\alpha}~r^{1/4} ~J_{-1/6}\left(\kappa_n r^{3/2}\right),$$ (64)

where $\{\kappa_n\}_{n=1}^\infty$ are the zeros of the Bessel function J_5/6. Furthermore, one can reverse the sign of $\kappa_n$, obtaining the conjugate eigenfunction $y_n^* = (\delta f_n^*, \delta g_n^*)$, up to a constant. Hence we denote $\kappa_{-n} = -\kappa_n$ and y_-n = y_n^*. An additional constant mode is necessary, namely, y₀ = (0,1), arising from the limit of Eqs. (61) and (62) as $\kappa\rightarrow0$. It corresponds to a constant velocity potential, so the full set of modes is $\{y_n\}_{n=-\infty}^\infty$.

These modes are orthogonal with respect to the product of Eq. (A.10), that is,

 

$$\left\langle y_n, y_m \right\rangle = {\alpha} \int_0^1 \frac{1}{... ...a^2 - \beta^2} \biggl[\frac{\sigma^2}{\alpha} r^{-1/2} \delta f_n^*~ \delta f_m$$

$$\quad - \beta~r^{1/2} \left(\delta f_n^*~ \delta g_m + \delta g_n... ...+ \alpha~r^{3/2} \delta g_n^* ~\delta g_m \biggl] {\rm d}r \propto \delta_{nm}.$$ (65)

(To have a non-dimensional scalar product, we introduce the constant $\alpha$.) We have checked the orthogonality relation (65) numerically for particular values of ${\cal M}$ and for modes up to $n \simeq 100$.

4.2 Initial conditions

Having solved the boundary problem in terms of orthogonal modes, we are in a position to solve the problem of initial conditions. As initial conditions, we must take two functions $\delta f(r,0)$ and  $\delta g(r,0)$ that satisfy the boundary conditions, namely, $\delta f(0,0) = \delta f(1,0) = 0$. Then we expand $y = (\delta f, \delta g)$ in modes as

$$% y(r) = \sum_{n=-\infty}^{\infty} c_n y_n(r).$$ (66)

Since y_-n = y_n^*, we require c_-n = c_n^* to have real y. So we can also write

$$\begin{eqnarray*}y(r) = c_0 y_0 + 2~ {\rm Re}\left[\sum_{n=1}^{\infty} c_n y_n(r)\right]. \end{eqnarray*}$$

Using the scalar product of Eq. (65), we can obtain the coefficients in the expansion (66):

$$% c_n = \frac{\left\langle y_n, y \right\rangle}{\left\langle y_n, y_n \right\rangle}\cdot$$ (67)

Either $(\delta f(0,r),\delta g(0,r))$ or the set of coefficients  $\{c_n\}_{n=0}^\infty$ contain the information on both density and velocity perturbations and constitute alternate definitions of the initial conditions.

4.3 Time evolution

Given the set of coefficients $\{c_n\}_{n=0}^\infty$, the solution of the boundary problem in terms of the modes (63) and (64) can be written

$$% y(r,t) = c_0 y_0 + 2~ {\rm Re}\left[\sum_{n=1}^{\infty} c_n {\rm e}^{-{\rm i}\omega_n t} y_n(r)\right].$$ (68)

Therefore, y(t,r) adopts an expansion similar to the one of y(0,r), Eq. (66), but with time-dependent coefficients

$$\begin{eqnarray*}c_n(t) = c_n {\rm e}^{-{\rm i}\omega_n t}. \end{eqnarray*}$$

Note that the norm

$$% \langle y, y \rangle = \left\vert c_0\right\vert^2 \langle ... ...left\vert c_n\right\vert^2 \left\langle y_n, y_n \right\rangle$$ (69)

is invariant in time. We can express this norm in terms of the variables  $\delta\rho$ and $\delta v$ or, alternatively, in terms of $\partial_t \phi = -c^2 \delta\rho/\rho- v \delta v$ and $\delta v$:

 

$$% \langle y, y \rangle \int_0^1 \left(\frac{r^2 c^2}{\rho} (\delta\rho)^2 + 2 r^2 v \delta\rho\delta v + r^2 \rho(\delta v)^2 \right) {\rm d}r$$ =

$$\int_0^1 \frac{r^2 \rho}{c^2} \left[ (\partial_t \phi)^2 + \left(c^2 - v^2\right) (\delta v)^2 \right] {\rm d}r.$$ = (70)

Comparison with the expression for the energy density (46) shows that the norm is just the energy (but for a constant factor), so it is naturally invariant in time. Therefore, the mode decomposition of the norm (69) expresses the total energy as a weighted sum of the energies of all modes. The "energy spectrum'' $\{\left\vert c_n\right\vert^2\}_{n=0}^\infty$ (invariant in time) characterises the initial conditions in combination with the initial phases of the coefficients. Suitable initial conditions must be normalizable, that is, must have finite energy. From the asymptotic form (B.10), ${\langle y_n, y_n \rangle} \propto 1/(1+6n)$, we deduce that we must demand that $\lim_{n \rightarrow\infty} \left\vert c_n\right\vert^2 = 0$; otherwise, the asymptotic form of the series (69) diverges like the harmonic series, at the least.

If we define the initial conditions by the energy spectrum and the initial phases of the coefficients, the forms of $\delta f$ and $\delta g$crucially depend on the correlation between those phases. In general, the evolution given by Eq. (68) is such that the correlation between the phases decreases with time and is finally lost as $t \rightarrow\infty$. To be precise, this happens unless the eigenvalues $\omega_n$have a common ratio, in which case the solution is periodic. But this does not happen in our case. In other words, the type of motion that we find is such that it leads to an evolution that is ergodic in a restricted sense, that is, an evolution that conserves the energy spectrum but which otherwise explores the available phase space. A typical state is given by taking the coefficient phases to be random. It could happen that either $\delta f(r,t)$ or  $\delta g(r,t)$ became concentrated about a value of r at some t (and they could even diverge) but it is extremely unlikely.

We have studied the evolution of particular radial perturbations of self-similar Bondi flow with initial conditions that are sufficiently smooth. They must have an energy spectrum that decays rapidly as $n \rightarrow\infty$. As this condition is time-invariant, smoothness is preserved in time. Functions such that they are concentrated in a sort of bell shape have special interest. A Gaussian profile is not allowed because of the boundary conditions ($\delta f$ vanishes at r=0,1); but we can take polynomial bell-shaped functions, for example, $b_n(x) = [4~x(1-x)]^n, \:\:n = 1, 2, \dots$ (which are succesively more concentrated). The time evolution of initial configuration $\delta f(r,0) = b_4(r), \delta g(r,0) = 0$ is plotted for several times in Fig. 1. This function gives rise to an energy spectrum that decays very rapidly, so the modes with $n \leq 50$ correspond to a relative error in its norm (total energy) that is smaller than 10^-6. We observe that $\delta g$ initially evolves quickly, reaching a non-negligible value at r=0 in a very short time. We also observe that  $\delta f(r,0)$ is almost reproduced at $t \simeq 1.7$ but $\delta g(r,0) = 0$ does not seem to be reproduced. Actually, there is no periodicity. The reason why $\delta f$ seems periodic is that the power spectrum is very concentrated in a small number of modes that are distributed with sufficient regularity. At any rate, the phase correlation among them is destroyed for longer times, as the following plots in Fig. 1 prove. The last couple of plots show a "typical'' long-time state, with random phases.

Figure 1: Time evolution of initial configuration $\delta f(r,0) = b_4(r), \delta g(r,0) = 0$, with modes up to n=50 (corresponding to a relative error in the norm <10-6). We have taken ${\cal M}=1/2$. On the left-hand side are plots of  $\delta f(r,t)$ and on the right-hand side are plots of  $\delta g(r,t)$, for t = 0, 0.1, 1.7, 100. The bottom plots correspond to a realization of the same energy spectrum with random phases.

Open with DEXTER

5 Discussion

We have examined some of the most interesting questions regarding the instability of Bondi accretion. The first question is about the spatial stability of perturbations, in the limit  $r \rightarrow0$. For our self-similar solutions, the meaning of this limit is different from the meaning in Kovalenko & Eremin's (1998), because similarity implies that the accretor and sonic radii vanish. The behaviour of perturbations in the limit  $r \rightarrow0$ depends on the type of inner boundary condition. In this regard, we have shown that the natural boundary condition is the no-energy-flux perturbation condition, which is equivalent to either the vanishing of the velocity potential time derivative or the vanishing of the variation in the matter flux current. We have focused on the latter condition as the natural generalization of the boundary conditions of Petterson et al. (1980) for radial perturbations. Under that condition, the relative variations in the physical variables (velocity and density) are bounded in the limit  $r \rightarrow0$. Therefore, the spatial stability of perturbations is ensured. This conclusion disagrees with an assertion by Kovalenko & Eremin (1998) that some perturbations are spatially unstable, for example, a non-radial subcritical mode. We attribute this difference to their imposing the no-energy-flux condition only under time averaging (which seems incomplete).

The second question is about the time evolution of perturbations. We have seen that the evolution is ergodic-like, as the perturbations make a random walk in the section of the total phase space with given energy spectrum, losing memory of the initial conditions (the initial phases of the coefficients c_n). Regarding Garlick's proposal of perturbations that could be increasing in a small region, independent of their global behaviour, we can now conclude that this can actually happen but is extremely unlikely: indeed, we have shown that an initial condition concentrated in a small region is forgotten in the evolution and not encountered again in a random walk in the big phase space. Furthermore, if we consider the non-linear coupling among the modes, with the corresponding transfer of energy, we deduce that the energy spectrum is not really conserved, so the only truly conserved quantity is the total energy. Therefore, in the long run, the energy spectrum thermalizes and the evolution becomes fully ergodic. In addition, the presence of viscosity (or other dissipative processes) would also transform the energy of perturbations into thermal energy. We see that an energy concentration in a small region corresponds to a situation of entropy that is spontaneously decreasing, so we fully agree with Garlick that this possibility must be ruled out as unphysical.

We have restricted ourselves to acoustic perturbations of self-similar Bondi accretion, but many results must be applicable to non-self-similar Bondi accretion. For example, the stability and ergodic-like nature of the evolution of acoustic perturbations is ultimately due to the reality of the spectrum of eigenvalues $\omega_n$and to its genericity (in the sense of not having a uniform distribution).

Acknowledgements

I am grateful to Carmen Molina-París for conversations. My work is supported by the "Ramón y Cajal'' program and by grant BFM2002-01014 of the Ministerio de Educación y Ciencia.

References

• Barenblatt, G. I. 1996, Scaling, Self-Similarity, and Intermediate Asymptotics (Cambridge University Press) (In the text)
• Bondi, H. 1952, MNRAS, 112, 195 [NASA ADS] (In the text)
• Coddington, E. A., & Levinson, N. 1955, Theory of Ordinary Differential Equations (New York: Mac-Graw Hill) (In the text)
• Foglizzo, T. 2002, A&A, 368, 311 [EDP Sciences] [NASA ADS] [CrossRef]
• Frank, J., King, A., & Raine, D. 1992, Accretion Power in Astrophysics (Cambridge: Cambridge University Press) (In the text)
• Garlick, A. R. 1979, A&A, 73, 171 [NASA ADS] (In the text)
• Kovalenko, I. G., & Eremin, M. A. 1998, MNRAS, 298, 861 [NASA ADS] [CrossRef] (In the text)
• Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics, Vol. I (New York: MacGraw-Hill) (In the text)
• Petterson, J. A., Silk, J., & Ostriker, J. P. 1980, MNRAS, 191, 571 [NASA ADS] (In the text)
• Ruffert, M. 1994, ApJ, 427, 342 [NASA ADS] [CrossRef] (In the text)
• Shore, S. N. 1992, An introduction to Astrophysical Hydrodynamics (San Diego: Academic Press) (In the text)
• Theuns, T., & David, M. 1992, ApJ, 384, 587 [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

   
Appendix A: Eigenvalue problem for radial perturbations

The perturbation Eqs. (23) and (24) become, for radial modes,

  

$$\partial_t \delta \rho + \frac{1}{r^2} \partial_r \left[ r^2 (\delta \rho \; v + \rho \; \delta v)\right] = 0,$$ (A.1)

$$\partial_t {\delta v} + \partial_r\left(\frac{c^2}{\rho}~\delta\rho+ v~{\delta v}\right) = 0.$$ (A.2)

After performing the Fourier transform for the time variable,

$$-{\rm i}\omega~\delta \rho + \frac{1}{r^2} \partial_r \left[ r^2 \left(\delta \rho \; v + \rho \; \delta v\right)\right] = 0,$$ (A.3)

$$-{\rm i}\omega~\delta v + \partial_r\left(\frac{c^2}{\rho}~\delta\rho+ v~{\delta v}\right) = 0.$$ (A.4)

We define the vector $x = (r^2 \delta \rho, \delta v)$ and the matrix

$$\begin{eqnarray*}A = \left( \begin{array}{cc} v(r) & r^2~\rho (r) \\ \frac{{c(r)}^2}{r^2~\rho (r)} & v(r) \end{array}\right), \end{eqnarray*}$$

such that $\det A = v^2 - c^2 \neq 0.$ Then we can write the ODE's (59) and (60) as the vector equation

$$% \frac{{\rm d}(A\cdot x)}{{\rm d}r} = {\rm i} \omega~x.$$ (A.5)

Its adjoint equation is

$$% A^\dag\cdot\frac{{\rm d} {\bar y}}{{\rm d}r} = {\rm i}\omega^*~{\bar y},$$ (A.6)

where the adjoint matrix $A^\dag$ is, in this case, just the transpose of A, and we have written the unknown vector function as ${\bar y}$ for later convenience. The derivation of this adjoint equation holds as long as the boundary conditions are such that  ${\bar y}^*\cdot A x$ vanishes at the ends of the interval of integration that defines the scalar product of functions.

According to the general theory of ODE's (Coddington & Levinson 1955), the eigenfunctions of the original ODE and the eigenfunctions of its adjoint equation are orthogonal, that is,

$$% \int {{\bar y}_n}^* \cdot x_m ~{\rm d}r \propto \delta_{nm},$$ (A.7)

and their respective eigenvalues are conjugate. But the eigenvalues are, in fact, real in our case due to a property of the matrix A; namely, it is transformed into its adjoint by the interchange of the two components of the vectors. To make this property formal, let us define the matrix

$$\begin{eqnarray*}U = \left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right), \end{eqnarray*}$$

such that $U^\dag = U$ and U^-1 = U. Then $A^\dag = U\cdot A\cdot U$, so we obtain from Eq. (A.6)

$$% A\cdot\frac{{\rm d} (U\cdot{\bar y}_n)}{{\rm d}r} = {\rm i}\omega_n^*~U\cdot{\bar y}_n.$$ (A.8)

On the other hand, defining

$$\begin{eqnarray*}y_n = A\cdot x_n = \left(r^2 \left(\delta \rho \; v + \rho \; \delta v\right), \frac{c^2}{\rho}~\delta\rho+ v~{\delta v}\right)_n, \end{eqnarray*}$$

we can write Eq. (A.5) as

$$% A\cdot\frac{{\rm d} y_n}{{\rm d}r} = {\rm i}\omega_n~y_n.$$ (A.9)

Therefore, we can identify $y_n = U\cdot{\bar y}_n$; that is, we have an equivalence between the original equation and its adjoint. In consequence, we have that $\omega_n = \omega_n^*$. This constitutes a proof of the reality of eigenvalues that is different from those given by Petterson et al. (1980) or Theuns & David (1992).

Furthermore, given the preceding identification, we can write the orthogonality relation (A.7) as the orthogonality of eigenfunctions x_n or, alternately, of eigenfunctions y_n; we shall use the latter, namely,

$$% \int y_n^* \cdot U \cdot A^{-1}\cdot y_m ~{\rm d}r \propto \delta_{nm}.$$ (A.10)

Note that the matrix $U \cdot A^{-1}$ is self-adjoint and positive and plays the rôle of a metric in the space of vectors y.

As we have pointed out above, the adjoint equation is valid if ${\bar y}^*\cdot A ~x = {\bar y}^*\cdot y = {y^1}^* y^2 + {y^2}^* y^1 = 0$at the ends of the integration interval (using the components of y = (y¹,y²)). Therefore, sufficient boundary conditions are given by either component y¹ or y² vanishing at both ends. We choose the former to vanish, as explained in the main text.

   
Appendix B: Asymptotic form of acoustic modes

B.1 Connection with the WKB solution

When the wave-length is small compared with the typical dimensions of the zone in which waves propagate, the WKB solution is a good approximation. The WKB solution of the differential equations for radial modes was obtained by Petterson et al. (1980):

$$\delta f(r) = \frac{1}{\sqrt{c~v}} \left[ A_1 \exp \left({\rm i}\... ...ght) + A_2 \exp \left({\rm i}\omega\int\! \frac{{\rm d}r}{v-c} \right) \right],$$ (B.1)

$$\delta\rho(r) = \frac{1}{r^2\sqrt{c~v}} \left[\frac{A_1}{v+c} \ex... ...2}{v-c} \exp\left({\rm i}\omega\int\! \frac{{\rm d}r}{v-c} \right) \right]\cdot$$ (B.2)

In this equation, the two functions that form the linear combination with (arbitrary) coefficients A₁ and A₂ represent the outgoing and incoming waves, respectively. Substituting the self-similar expressions of v and c and performing the integrals, we obtain the corresponding expressions of the WKB solutions  $\delta f(r)$ and  $\delta\rho(r)$. From them we can further obtain

 

$$% \delta g(r) = \frac{\sigma}{\alpha}~ r^{-1/2} \left[ A_1 \exp \... ...sigma)} - A_2 \exp \frac{2{\rm i}\omega~ r^{3/2}}{3(\beta-\sigma)} \right]\cdot$$ (B.3)

This formula is valid for large $\omega$ or, equivalently, large r.

The two component waves are analogous to the functions $J_{\nu}$ and $N_{\nu}$ in the radial solution (42). To find the precise relation between them, it is convenient to express the radial solution (42) in terms of the complex combinations of the Bessel functions, namely, the Hankel functions  $H^{(1)}_{\nu}$ and  $H^{(2)}_{\nu}$ (Morse & Feshbach 1953). The asymptotic forms of the Hankel functions for large values of their arguments is

  

$$H^{(1)}_{\nu}(x) = \sqrt{\frac{2}{\pi x}} \exp\left[{\rm i}\left(x-\frac{\nu\pi}{2}- \frac{\pi}{4}\right)\right] + O\left(x^{-3/2}\right),$$ (B.4)

$$H^{(2)}_{\nu}(x) = \sqrt{\frac{2}{\pi x}} \exp\left[-{\rm i}\left(x-\frac{\nu\pi}{2}- \frac{\pi}{4}\right)\right] + O\left(x^{-3/2}\right).$$ (B.5)

Expressing Eq. (42) in terms of Hankel functions (with two new constants C'₁ and C'₂) and using the preceding asymptotic forms, we have

$$% R(r) \!\sim\! r^{1/4} {\rm e}^{-{\rm i} \mu r^{3/2}} r^{-3/4} \... ...kappa r^{3/2}\right) \!+\! C'_2 \exp\left(-{\rm i}\kappa r^{3/2}\right)\right],$$ (B.6)

still with $\mu =\frac{2 \beta \omega}{3(\sigma^2-\beta^2)}$, $\kappa= \frac{2 \sigma\omega}{3 (\sigma^2 - \beta^2)}.$ Note that this asymptotic form becomes $\nu$-independent and therefore l-independent. Given that

$$\begin{eqnarray*}-\mu \pm \kappa= - \frac{2 \beta \omega}{3\left(\sigma^2-\beta... ...\sigma^2-\beta^2\right)} = \frac{2 \omega}{3(\beta\pm \sigma)}, \end{eqnarray*}$$

we have

$$% R(r) \sim r^{-1/2} \left[C'_1 \exp \frac{2 {\rm i}\omega~ r^{3/... ...igma)} + C'_2 \exp\frac{2 {\rm i}\omega~ r^{3/2}}{3(\beta- \sigma)}\right]\cdot$$ (B.7)

This last form is equivalent to the WKB form (B.3) but is valid for non-radial modes as well.

B.2 Asymptotic form of radial eigenfunctions

For $n \gg 1$ we can substitute in Eqs. (63) and (64) the asymptotic forms of the Bessel functions for high values of their arguments (Morse & Feshbach 1953):

$$\begin{eqnarray*}J_\nu(x) = \sqrt{\frac{2}{\pi x}} \cos\left(x-\frac{\nu\pi}{2}- \frac{\pi}{4}\right) + O\left(x^{-3/2}\right) \end{eqnarray*}$$

[consistently with Eqs. (B.4) and (B.5)]. In particular,

$$\begin{eqnarray*}&& J_{5/6}\left(\kappa_n r^{3/2}\right) = \sqrt{\frac{2}{\pi \k... ...} r^{-3/4} \cos\left(\kappa_n r^{3/2} -\frac{\pi}{6}\right)\cdot \end{eqnarray*}$$

Replacing $\cos\left(\kappa_n r^{3/2} - 2\pi/3\right)$ with $\sin\left(\kappa_n r^{3/2} - \pi/6\right)$, we can make use of the boundary condition $\delta f(1,t) = 0$ to obtain the spectrum $\kappa_n \approx \left(n + 1/6\right)\pi, \; n= 1,2,\ldots$. The corresponding proper functions are

 

$$\delta f_n(r,t) \; \approx \; \exp\left[-{\rm i}\pi\left(n + \frac{1}{6}\right) \left({\cal M}r^{3/2}+\frac{3}{2} \left(1 -{\cal M}^2\right)t\right)\right]$$

$$\times r^{1/2}~ {\rm i} ~\frac{\alpha}{\sigma}~ \frac{1}{\pi} ~{\... ... \left[\left(n \!+\! \frac{1}{6}\right)\pi~ r^{3/2} \!-\! \frac{\pi}{6}\right],$$ (B.8)

 

$$% \delta g_n(r,t) \; \approx \; \exp\left[-{\rm i}\pi\left(n + \frac{1}{6}\right) \left({\cal M}r^{3/2}+\frac{3}{2} \left(1 -{\cal M}^2\right)t\right)\right]$$

$$\times r^{-1/2} ~\frac{1}{\pi} ~{\sqrt{\frac{2}{n + 1/6}}}~ \cos \left[\left(n + \frac{1}{6}\right)\pi~ r^{3/2} - \frac{\pi}{6}\right]\cdot$$ (B.9)

In Fig. B.1 we plot two low radial modes ( ${\cal M}=1/2$), that is, their real parts, to compare them with their asymptotic expressions according to Eqs. (B.8) and (B.9). We note that the difference between the actual and the asymptotic forms diminishes as n or r grows, but it is small even in the most unfavourable case, namely, n=1 and r < 0.5.

Figure B.1: Real parts of radial modes with n=1 and n=4 (solid lines) compared with their asymptotic trigonometric forms (dashed lines). Note that the difference diminishes as r grows.

The simpler form in terms of trigonometric functions allows us to obtain an analytic form of the orthogonality of eigenfunctions:

 

$$\left\langle y_n, y_m \right\rangle \approx 12 \frac{\int\limits_0^1 {\rm e}^{\frac{{\rm i} }{2}\left( m \!-\... ...left( 1 + 6~m \right) \left(1 + 6~n \right)}}~ \left(1 - {{\cal M}}^2 \right) }$$

$$\delta_{nm}~ \frac{8}{\left( 1 + 6~n \right) {\pi }^2 \left(1 - {{\cal M}}^2\right)}\cdot$$ = (B.10)

In particular, we have a useful approximation of the norm. We have found that the norm given by this formula is in fact already a good approximation for n =1, with a relative error smaller than 0.001(for ${\cal M}=1/2$), which rapidly diminishes as n increases.