format: text

text:
# Low autocorrelated multi-phase sequences

\[

## Abstract

The interplay between the ground state energy of the generalized Bernasconi model to multi-phase, and the minimal value of the maximal autocorrelation function, $`C_{max}=\mathrm{max}_K|C_K|`$, $`K=1,..N1`$, is examined analytically and the main results are: (a) The minimal value of $`\mathrm{min}_NC_{max}`$ is $`0.435\sqrt{N}`$ significantly smaller than the typical value for random sequences $`O(\sqrt{\mathrm{log}N}\sqrt{N})`$. (b) $`\mathrm{min}_NC_{max}`$ over all sequences of length $`N`$ is obtained in an energy which is about $`30\%`$ above the ground-state energy of the generalized Bernasconi model, independent of the number of phases $`m`$. (c) The maximal merit factor $`F_{max}`$ grows linearly with $`m`$. (d) For a given $`N`$, $`\mathrm{min}_NC_{max}\sqrt{N/m}`$ indicating that for $`m=N`$, $`\mathrm{min}_NC_{max}=1`$, i.e. a Barker code exits. The analytical results are confirmed by simulations.

\] In many applications of communication science , as well as a variety of other fields, it is necessary to find sequences with low autocorrelation functions. Some of these applications utilize the pulse compression feature of the low autocorrelated sequences to obtain high resolution in radars and sonars. In other applications, the shifts of such periodic sequences can be used to identify users in multi-user systems. Due to their importance, low autocorrelated sequences have evoked a wide spread interest accompanied by the development of various methods for constructing such sequences .

In order to construct a binary sequence $`S=(s_1,\mathrm{},s_N)`$ with low off peak autocorrelations, one has to define which quantity has to be minimized. The different applications divide the low autocorrelated sequences into two types in correspondence with the quantity which they minimize. The first kind of low autocorrelated sequences minimizes the Hamiltonian of the Bernasconi model  which is given by:

$$H=\underset{K=1}{\overset{N1}{}}C_K^2$$
(1)
where for the case of non-periodic boundary conditions, which is at the center of our study,

$$C_K=\underset{j=1}{\overset{NK}{}}s_js_{j+K},$$
(2)
and for the periodic case $`C_K=_{j=1}^Ns_js_{((j+K1)modN)+1}`$. Note that for random sequences the average value of $`H`$ in the non-periodic case is $`N^2/2`$, whereas for low autocorrelated sequences this energy is reduced by a merit factor $`F>1`$ to $`N^2/2F`$. However, there are applications for which the maximal off peak autocorrelation, $`C_{max}=\mathrm{max}_K|C_K|`$, has to be minimized. The second kind of low autocorrelated sequences are the solutions of this minimization problem. Note that for random sequences $`C_{max}`$ is typically $`O(\sqrt{\mathrm{log}N}\sqrt{N})`$. For a sequence of length $`N`$, the maximal possible ratio between the peak, $`|C_0|=N`$, and the maximal off-peak autocorrelations, $`|C_K|`$ with $`K=1,\mathrm{},N1`$ is $`N/1`$. The only known binary sequences with this ratio are the Barker codes of length $`2`$, $`3`$, $`4`$, $`5`$, $`7`$, $`11`$ and $`13`$. Obviously, the Barker sequences, when they exist, furnish a minimum for the two minimization problems. However, Turyn  has shown that no other binary codes such as this exist for any length less than $`144`$ or for odd length greater than $`13`$. An exact solution for these two minimization problems is known only for systems (up to $`N=59`$) which are small enough to permit an effective exhaustive search . Extrapolation of the ground state energies which were found for small systems using an exhaustive search indicates $`F_{max}=lim_N\mathrm{}\frac{N^2}{2H}=8.5`$ . Moreover, $`F_{max}`$ was conjectured by Golay  to be bounded from above by $`12.324`$.

The following questions regarding low autocorrelated sequences are still open and are at the center of our study: (a) Do the sequences that minimize $`C_{max}`$ minimize the energy of the Bernasconi model as well? (b) In case that the two minimization problems are not equivalent, how far are the energy values of the sequences which minimize $`C_{max}`$ from the values of the ground state energy? (c) Exhaustive search of sequences with length up to $`N=59`$ show that the degeneracy of the ground state is of O(1) and is bounded from above by $`8`$. Is the degeneracy of the states which minimize $`C_{max}`$ is of O(1) as well? (d) In this context, which of the two minimal quantities is easier to approach? Additional interesting questions arise when multi-phase sequences, whose terms are complex $`m`$-th roots of $`1`$ for $`m>2`$, are considered. For such sequences

$$s_l=\mathrm{exp}(i\frac{2\pi }{m}l)$$
(3)
where l=1,…,m and the correlations are defined as

$$C_K=\underset{j=1}{\overset{NK}{}}s_js_{j+K}^{}.$$
(4)
where $`s_{j+K}^{}`$ is the conjugate of $`s_{j+K}`$. The information about multi-phase low autocorrelated sequences is more limited, mainly since the configuration space grows like $`m^N`$ which makes the exhaustive search ineffective even for very small sequences. Hence, not much is known about the influence of the number of phases both on the maximal value of $`F`$ and on the minimal value of $`C_{max}`$. Is it possible to attain better solutions by increasing the number of phases?

In this study we suggest an analytical technique to calculate the quantities of interest, namely, $`\mathrm{max}_NF`$ and $`\mathrm{min}_NC_{max}`$ for both binary and multi-phase sequences. An obvious lower bound of $`\mathrm{min}_{N,F}C_{max}`$ over all sequences of length $`N`$ and a merit factor $`F`$ is obtained by assuming that all $`|C_K|`$ equal $`\sqrt{N/(2F)}`$. Nevertheless it is tempting to investigate the relations between $`F`$ and $`\mathrm{min}_{N,F}C_{max}`$ and to measure more accurately the optimal value $`\mathrm{min}_NC_{max}`$ over all sequences of length $`N`$ similarly to the aforementioned upper bound given for $`F_{max}`$ . For our analytical study we consider the non-periodic Bernasconi model. We approximate the sequences to be random and the autocorrelation functions to be independent variables with the following Gaussian distribution

$$P(C_K)=\frac{1}{\sqrt{2\pi (NK)}}\mathrm{exp}(\frac{C_K^2}{2(NK)}).$$
(5)
Using Eq. (5), the probability $`P_F(C_{max})`$ that the autocorrelations of a sequence with a merit factor $`F=N^2/2H`$ have an upper bound $`C_{max}`$ such that $`|C_K|C_{max}`$ is given by

$$P_F(C_{max})=_{C_{max}}^{C_{max}}\underset{K}{}dC_KP(C_K)\delta (\frac{1}{N}\underset{K=1}{\overset{N1}{}}C_K^2\frac{N}{2F}).$$
(6)
Since there are $`2^N`$ distinct binary sequences of length $`N`$, it is necessary that $`P_F2^N`$ in order for a sequence with the corresponding features, namely a merit factor $`F`$ and an upper bound $`C_{max}`$, to exist. Hence, by equating $`P_F(C_{max})`$ to $`2^N`$ one can find the minimal upper bound $`\mathrm{min}_{N,F}C_{max}`$ among the upper bounds of all sequences with a merit factor $`F`$. Moreover, the Gaussian distributions of $`C_K`$ imply that $`\mathrm{min}_{N,F}C_{max}`$ must be of $`O(\sqrt{N})`$ in order for $`P_F(C_{max})`$ to be equal $`2^N`$. Assigning $`\mathrm{min}_{N,F}C_{max}=B(F)\sqrt{N}`$ and inserting the integral representation of $`\delta `$ function in Eq. (6) , the saddle point method can be used to obtain the following set of equations

$`{\displaystyle \frac{\lambda }{2F}}{\displaystyle \frac{1}{4\lambda }}((12\lambda )\mathrm{ln}(12\lambda )+2\lambda )`$ (7)
$`{\displaystyle _0^1}\mathrm{ln}\text{erf}(B(F)g(y,\lambda ))𝑑y=\mathrm{ln}2,`$ (8)
$`{\displaystyle \frac{1}{2F}}{\displaystyle \frac{1}{4\lambda ^2}}(\mathrm{ln}(12\lambda )+2\lambda )+`$ (9)
$`B(F)\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle _0^1}{\displaystyle \frac{\mathrm{exp}(B(F)^2g(y,\lambda )^2)}{\text{erf}(B(F)g(y,\lambda ))\sqrt{2}g(y,\lambda )}}𝑑y=0.`$ (10)
where $`g(y,\lambda )=\sqrt{\frac{12\lambda (1y)}{2(1y)}}`$. Solving numerically this set of equations, it turns out that a solution exists only for a bounded region of $`F`$, where $`.215F12.324`$. These two limit values of $`F`$ are alternatively obtained by taking the limits of the integration in Eq. (6) to infinity. The existence of a lower bound is a consequence of the fact that as the energy becomes smaller the allowed values of the correlations decrease and the probability for such sequences (Eq. (6)) decreases as well. As was pointed out above, when this probability drops below $`2^N`$, there are no more sequences with the required energy. Similarly, for high energy values, high values of correlations are required and the corresponding probabilities become smaller again. The analysis of $`B(F)`$ inside the allowed region reveals that the minimal value of $`B(F_c)=0.435`$ is obtained for $`F_c=8.839`$ whereas the maximal merit factor $`F=12.324`$ corresponds to a higher value of $`B(F)=1`$. Thus, the merit factor of sequences which minimizes $`C_{max}`$ is around $`71\%`$ lower then $`F_{max}`$. Exhaustive search results in the region $`N[15,28]`$ indicate that the average ratio between the two merit factors, $`<<<F(\mathrm{min}C_{max})>>/F_{max}>_N`$, is $`0.6771`$ where $`<<`$ $`>>`$ denotes the average over all sequences which minimize $`C_{max}`$ whereas $`<`$ $`>_N`$ denotes the average over the different sequences length $`N`$. The deviation from the analytical result $`0.717`$ is attributed to finite size effects. Fig. 1 shows the behavior of $`B(F)`$ derived from the numerical solution of Eqs. (8) and (10). The behavior of $`P_F(C_{max})`$ divides this graph into three regimes. Starting with the large $`F`$ regime, there is only one sequence whose $`C_{max}`$ equals $`B(F)\sqrt{N}`$. This feature holds until $`F_c=8.839`$ for which the minimal $`B(F)=0.435`$ is obtained. Further decrease of $`F`$ results in an increasing number of sequences with an upper bound of $`B(F)\sqrt{N}`$. This increment terminates at $`F=0.8`$ and for smaller values of $`F`$ the number of sequences reduces to $`1`$ as $`F`$ approaches $`0.215`$. In order to compare the analytical behaviour of $`B`$ to that obtained from simulations, it is necessary to find a way to circumvent the mismatch between the analytical values of $`F`$ and those obtained from simulations. The rescaling of $`F`$ by $`F_{max}`$ yields a parameter $`F/F_{max}`$ whose range of values, $`[0,1]`$, does not depend on a specific realization. This allows a comparison between analytical results and simulations. Results of an exhaustive search over sequences of length $`N=32`$ show that the behavior of $`B`$ as a function of $`F/F_{max}`$ resembles the analytical prediction (Fig. 2).

Moreover, the existence of the three aforementioned regions of $`P_F(C_{max})`$ is confirmed by counting the number of sequences with the minimal upper bound $`B`$ for each of the micro-canonical ensembles with a merit factor $`F`$ and length $`N=32`$. The histogram is demonstrated in the inset of Fig. 2, where the merit factor is rescaled to $`F/F_{max}`$. In the small $`F`$ regime, the average number of sequences increases with $`F`$ until $`F/F_{max}0.2`$ from which the number of sequences decreases up to the plateau which starts at $`F/F_{max}0.83`$.

Now we turn to study the influence of the number of phases both on the maximal value of $`F`$ and on the minimal value of $`C_{max}`$. We assign $`D_K`$ and $`E_K`$ to be independent Gaussian variables which represent the real and the imaginary parts of $`C_K`$ respectively. In case that $`m=2^p`$ where $`p`$ is an integer number greater than $`1`$, these two parts have the same probability distribution

$$P(D_K)=\frac{1}{\sqrt{\pi (NK)}}\mathrm{exp}(\frac{D_K^2}{(NK)}).$$
(11)
Under the same assumptions that have been used in the derivation of the binary case, Eq. (6) becomes

$`P_F^m(C_{max})={\displaystyle _{C_{max}}^{C_{max}}}{\displaystyle \underset{K}{}}dD_KdE_KP(D_K)P(E_K)`$ (12)
$`\delta ({\displaystyle \frac{1}{N}}{\displaystyle \underset{K=1}{\overset{N}{}}}(D_K^2+E_K^2){\displaystyle \frac{N}{2F}}).`$ (13)
Similarly to the binary case, if $`P_F^m(C_{max})=m^N`$ then $`B(F)\sqrt{N}`$ is the minimal value of $`C_{max}`$ for which there is a sequence with a merit factor $`F`$. The same procedure that has been used to calculate $`B(F)`$ for the binary case results in the following set of equations

$`{\displaystyle \frac{\lambda }{2F}}{\displaystyle \frac{1}{\lambda }}((1\lambda )\mathrm{ln}(1\lambda )+\lambda )`$ (14)
$`2{\displaystyle _0^1}\mathrm{ln}\text{erf}(B(F)g(y,\lambda ))𝑑y`$ $`=`$ $`\mathrm{ln}(m).`$ (15)
$`{\displaystyle \frac{1}{2F}}{\displaystyle \frac{1}{\lambda ^2}}(\mathrm{ln}(1\lambda )+\lambda )+`$ (16)
$`{\displaystyle \frac{2B(F)}{\sqrt{\pi }}}{\displaystyle _0^1}{\displaystyle \frac{\mathrm{exp}(B(F)^2g(y,\lambda )^2)}{\text{erf}(B(F)g(y,\lambda ))g(y,\lambda )}}𝑑y=0,`$ (17)
where $`g(y,\lambda )=\sqrt{\frac{1\lambda (1y)}{(1y)}}`$. Solving numerically these two equations, the behaviors of $`\mathrm{min}_NC_{max}/\sqrt{N}`$ and $`F_{max}`$ as a function of $`m`$ are obtained. Fig. 3 shows that $`\mathrm{min}_NC_{max}/\sqrt{N}`$ is proportional to $`1/\sqrt{m}`$ with slight deviations for small $`m`$. This relation between $`B(F_c)=\mathrm{min}_NC_{max}/\sqrt{N}`$ and $`m`$ implies that such a relation holds for $`B(F_{max})`$ as well. Since the variance of $`C_K^2`$ equals $`NK`$, their typical values drop linearly with $`K`$. Hence, $`H=N^2/2F`$ is approximated by a sum of an algebric series of $`N1`$ terms, $`C_K^2`$, with an upper bound, $`C_{max}^2`$, of $`O(N/m)`$. Assuming that $`C_K^2`$ is homogeneously distributed between $`0`$ to $`C_{max}^2`$ yields a linear increment of $`F_{max}`$ as a function of $`m`$ in agreement with the numerical solution of $`F_{max}`$ which is depicted in the inset of Fig. 4. The results of $`F_{max}`$ and $`B(F_c)=\mathrm{min}_NC_{max}/\sqrt{N}`$ as a function of $`m`$ show that the solutions of the two minimization problems are improved by increasing the number of phases $`m`$.

These results raise the question whether it is possible to increase $`m`$ such that $`\mathrm{min}_NC_{max}`$ becomes 1. An asymptotic expansion of Eqs. (15) and (17) reveals that for $`m=N`$, $`\mathrm{min}_NC_{max}=1`$. Note that $`C_{N1}`$ is always 1 and therefore $`\mathrm{min}_NC_{max}=1`$ holds for the entire regime $`m>N`$.

We used the simulated annealing method for sequences of lengths $`N=32`$ with different number of phases $`m`$, to find $`\mathrm{min}_NC_{max}`$. The relatively small sequences size was chosen to enable an appropriate scan of the configuration space in reasonable computational time. The results are exhibited in Fig. 4, and support the anticipated $`\sqrt{N/m}`$ behavior of $`\mathrm{min}_NC_{max}`$. For $`m=32`$, there is a deviation of $`\mathrm{min}_NC_{max}`$ from the analytical prediction $`\mathrm{min}_NC_{max}=1`$, probably since the simulated annealing method yields only suboptimal solutions.

Finally, we would like to examine the simulated annealing method in light of the results of this study. Simulations show that for the same running times, the simulated annealing method yields $`C_{max}`$ which is closer to its minimal value than $`_{K=1}^{N1}C_K^2`$. This can be explained by the larger degeneracy of $`C_{max}`$ compared with that of $`_{K=1}^{N1}C_K^2`$. Moreover, it turns out that in order to minimize $`C_{max}`$, it is preferable to start the searching process with the minimization of the energy function, $`_{K=1}^{N1}C_K^2`$, and then replace it with $`C_{max}`$ . In this way the system avoids the plateaus which characterize the landscape of $`C_{max}`$ in the configuration space. However, future research is necessary to find out how the results of this study can be applied to further improve the searching processes of low autocorrelated sequences.

Critical comments on the manuscript by R. Metzler are acknowledged. The research is supported in part by the Israel Academy of science.


text:
# An analytical MHD wind model with latitudinal dependences obtained using separation of the variables

## 1 Introduction

Recent detailed observations of the solar wind (e.g., via Ulysses and SoHO, Feldman et al. 1996) and other astrophysical outflows (e.g., via the Hubble Space Telescope, Biretta 1996) have highlighted the need to provide a satisfactory description of cosmic MHD outflows which deals adequately and simultaneously with their: (i) spatial 2-D (or 3-D) character, (ii) nonlinear dynamics, (iii) detailed energy balance and (iv) time-dependent nature. At the moment this task is far from being completed and still remains a challenge for the future. The present treatment is not going to change this fact that we do not have available to this day a satisfactory MHD solar/stellar wind description. Hence, by necessity, in the present paper we shall constrain our efforts only to the above items (i)-(ii).

The first models of astrophysical outflows were polytropic and 1-D (Parker 1958, 1963), or, all forces perpendicular to the equatorial plane were ignored (Weber & Davis 1967). The discovery of flows through coronal holes called for the investigation of quasi-two dimensional polytropic models at first (Kopp & Holzer 1976, Habbal & Tsinganos 1983) and fully 2-D numerical models later (Pneuman & Kopp 1971, Steinolfson et al. 1982, Nerney & Suess 1975, Sakurai 1985). Recently, time-dependent studies of thermally and/or magneto-centrifugally-driven polytropic wind-type outflows provided new and promising results in relation to the question of the degree of collimation of the outflow (Ouyed & Pudritz 1997, Keppens & Goedbloed 1999, Ustyugova et al. 1999, Tsinganos & Bogovalov 2000, Krasnopolsky et al. 2000). Non-polytropic wind modelling on the other hand, with energy and momentum addition and finite thermal conductivity has been increasingly used because observations have highlighted the fact that the acceleration of the solar wind in high-speed streams (Feldman et al. 1996, Giordano et al. 2000) does imply energy and momentum addition in the solar corona (Leer & Holzer 1980, Steinolfson 1988, Suess et al, 1996, Hansteen et al. 1997, Wang et al. 1998).

A fully 2-D MHD modelling of steady hydromagnetic wind-type outflows in open magnetic fields has been introduced (Low & Tsinganos 1986, Tsinganos & Low 1989) with a variable polytropic index $`\gamma `$ wherein it was shown that the density needs to decrease with latitude for an accelerated wind. This is because a dipolar magnetic field needs to be kept open by a pressure that must decrease towards the pole. If the density does not vary with latitude, there is a smaller pressure gradient to drive the flow near the pole, exactly where the magnetic field is open to allow the wind to escape. The resulting acceleration is too low, since gravity dominates, and the flow does not reach a high enough terminal speed. The only way out is to allow the density to increase with latitude, faster than the pressure does (Hu & Low 1989). Subsequently, solutions of the MHD equations with a latitudinally-dependent density and a helicoidal geometry of the streamlines were analysed (Tsinganos & Trussoni 1991, hereafter referred to as TT91), but the latitudinal dependence of the different quantities was assumed a priori. Under the same assumptions, this work was later continued (Trussoni & Tsinganos 1993, Trussoni et al. 1997). Other analytical solutions which have used similar assumptions have concentrated on the asymptotic analysis and its links to the intitial boundary conditions and were able to derive a criterion for the collimation of winds into jets (Sauty & Tsinganos 1994, Sauty et al. 1999).

In this paper we shall deduce from the governing MHD equations via a nonlinear separation of the variables the latitudinal variation of the density and other relevant physical quantities of the wind, instead of adopting a priori for them a specific form, as it was done in TT91 and in subsequent papers of that series. This shall also further extend previous work (Lima & Priest 1993, hereafter referred as Paper I), which deduced a general class of solutions of the hydrodynamic equations relevant to stellar winds with a helicoidal geometry. Sect. 2 is devoted to the method of solution, Sect. 3 deals with its parametric study and finally Sect. 4 ends with a discussion of the results.

## 2 Method of solution

### 2.1 Basic equations and assumptions

We shall seek 2-D solutions of the set of steady MHD equations describing the dynamical interaction of an inviscid, compressible and highly conducting plasma with an axially symmetric magnetic field created by a central rotating object. These include Maxwell’s equation for the divergence of the magnetic field, the induction equation in the limit of large conductivity and the equations for conservation of mass and momentum (Priest 1982). The pressure is related to the density and temperature through the classical ideal gas law $`p=(2k_B/m)\rho T`$ while the temperature is obtained from an energy conservation equation such as the first law of thermodynamics, including a heating/cooling term.

The system of MHD equations can be closed if the heating function is functionally related to the thermodynamic variables and the flow field. One such simple example which has been widely used in the literature is when we have a polytropic law between pressure and density. The alternative approach that we follow here is to find some other functional relationship between the heating rate, pressure, density and flow speed such that the variables are separable.

The obvious choice of coordinate system appropriate to this type of problem is that of spherical polar coordinates $`(r,\theta ,\varphi )`$, with $`\theta `$ as the co-latitude. With axisymmetry, in order to make the above equations more tractable from an analytical point of view, the simplifying assumption will be made that the meridional components of both the velocity and magnetic field can be neglected (i.e. $`V_\theta =B_\theta =0`$). Under this assumption the projection of the field lines on the poloidal plane corresponds to radial lines. Note that it has been recently shown that this assumption of radiality is a good one, at least in the case of the solar wind (Wang et al 1998, Tsinganos & Bogovalov 2000).

Then, the system of MHD equations can be re-written explicitly in the following way:

$$\frac{}{r}(B_rr^2)=0,$$
(1)
$$\frac{}{r}(rV_\varphi B_rrV_rB_\varphi )=0,$$
(2)
$$\frac{}{r}(\rho r^2V_r)=0,$$
(3)
$$\rho V_r\frac{V_r}{r}\rho \frac{V_\varphi ^2}{r}=\frac{p}{r}\frac{B_\varphi ^2}{4\pi r}\frac{B_\varphi }{4\pi }\frac{B_\varphi }{r}\frac{\rho GM}{r^2},$$
(4)
$$\rho V_\varphi ^2\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }=\frac{p}{\theta }+\frac{B_r}{4\pi }\frac{B_r}{\theta }+\frac{B_\varphi ^2}{4\pi }\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }+\frac{B_\varphi }{4\pi }\frac{B_\varphi }{\theta },$$
(5)
$$\rho V_r\frac{V_\varphi }{r}+\rho \frac{V_rV_\varphi }{r}=\frac{B_rB_\varphi }{4\pi r}+\frac{B_r}{4\pi }\frac{B_\varphi }{r}.$$
(6)
### 2.2 A technique based on a separation of the variables

At this stage, we introduce a second simplifying assumption, namely that the variables are separable in $`r`$ and $`\theta `$. This will transform the above system of partial differential equations into a system of ordinary differential equations which are analytically more tractable.

Denoting by $`r_0`$ the radius at the base of the atmosphere, we can non-dimensionalise all quantities with respect to their values at this reference level. In particular we set $`R=r/r_0`$. Using the assumption of separation of variables, we can write the radial velocity, the azimuthal velocity and the azimuthal magnetic field as, respectively,

$$V_r(R,\theta )=V_0Y(R)v_r(\theta ),$$
(7)
$$V_\varphi (R,\theta )=V_1A(R)v_\varphi (\theta ),$$
(8)
$$B_\varphi (R,\theta )=B_1M(R)b_\varphi (\theta ),$$
(9)
where $`V_0`$, $`V_1`$ and $`B_1`$ correspond to their reference values. Note that, at this stage, the functions $`v_r(\theta )`$, $`v_\varphi (\theta )`$ and $`b_\varphi (\theta )`$ are completely arbitrary. ¿From Eqs.(1), (3) we must also have

$$B_r(R,\theta )=B_0\frac{b_r(\theta )}{R^2},\rho (R,\theta )=\rho _0\frac{g(\theta )}{YR^2},$$
(10)
with $`B_0`$ and $`\rho _0`$ as the values of the radial magnetic field and density at the reference level, while $`b_r(\theta )`$ and $`g(\theta )`$ are as yet to be determined. Note that the radial magnetic field has a monopole geometry modified by the presence of $`b_r(\theta )`$.

Using the same technique described in Paper I, let us eliminate the pressure term between the $`r`$\- and $`\theta `$\- components of the momentum equation, Eqs.(4), (5), by differentiating the first one with respect to $`\theta `$ and the second one with respect to $`r`$, and then adding. The resulting expression is

$`{\displaystyle \frac{}{\theta }}\left(\rho {\displaystyle \frac{V_\varphi ^2}{r}}\right){\displaystyle \frac{}{\theta }}\left(\rho V_r{\displaystyle \frac{V_r}{r}}\right){\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{}{\theta }}\left({\displaystyle \frac{B_\varphi ^2}{r}}\right)`$
$`{\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{}{\theta }}\left(B_\varphi {\displaystyle \frac{B_\varphi }{r}}\right){\displaystyle \frac{}{\theta }}\left({\displaystyle \frac{\rho GM}{r^2}}\right)`$
$`{\displaystyle \frac{}{r}}\left(\rho V_\varphi ^2{\displaystyle \frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }}\right)+{\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{}{r}}\left(B_r{\displaystyle \frac{B_r}{\theta }}\right)`$
$`+{\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{}{r}}\left(B_\varphi ^2{\displaystyle \frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }}\right)+{\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{}{r}}\left(B_\varphi {\displaystyle \frac{B_\varphi }{\theta }}\right)`$ $`=0.`$
Under the assumption of separation of variables, the above equation will be transformed into an ordinary differential equation involving functions of $`R`$ alone. For that purpose, the functions of $`\theta `$ in each term will be set proportional to each other. We should keep in mind, however, that this last procedure may not give us the most general separable solution in the sense that some terms could be also set proportional to the sum of the others. Thus, from the last two terms, we may write

$$b_\varphi \frac{\mathrm{d}b_\varphi }{\mathrm{d}\theta }=ϵb_\varphi ^2\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta },$$
(12)
which implies that

$$\frac{\mathrm{d}}{\mathrm{d}\theta }\left(\frac{b_\varphi ^2}{\mathrm{sin}^{2ϵ}\theta }\right)=0,$$
(13)
and so

$$b_\varphi (\theta )=\mathrm{sin}^ϵ\theta ,$$
(14)
where $`ϵ`$ is an arbitrary constant. Note that the constant of integration has been set to unity, without loss of generality. Any other value for this constant can be incorporated into the radial dependence of $`B_\varphi (R,\theta )`$. In what follows, the choice of the constants of proportionality between different terms as well as the constants of integration will be made so as to obtain the simplest possible solutions, without loss of generality. Comparing the $`7^{th}`$ and $`8^{th}`$ terms, we can put

$$b_r\frac{\mathrm{d}b_r}{\mathrm{d}\theta }=\mu ϵb_\varphi ^2\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta },$$
(15)
which implies

$$b_r(\theta )=\sqrt{1+\mu \mathrm{sin}^{2ϵ}\theta },$$
(16)
where $`\mu `$ is a second arbitrary constant.

From the $`5^{th}`$ and $`8^{th}`$ terms in Eq.(LABEL:18) we may also obtain

$$\frac{\mathrm{d}g}{\mathrm{d}\theta }=2\delta ϵb_\varphi ^2\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta },$$
(17)
thus resulting in

$$g(\theta )=1+\delta \mathrm{sin}^{2ϵ}\theta ,$$
(18)
in which $`\delta `$ is the third arbitrary constant. As for the $`6^{th}`$ and $`8^{th}`$ terms we can write

$$gv_\varphi ^2=b_\varphi ^2,$$
(19)
or, equivalently,

$$v_\varphi (\theta )=\frac{\mathrm{sin}^ϵ\theta }{\sqrt{1+\delta \mathrm{sin}^{2ϵ}\theta }}.$$
(20)
Finally, from the $`\varphi `$-component of the momentum equation, Eq.(6), for the separation in $`R`$ and $`\theta `$ to work we must have

$$gv_rv_\varphi =b_rb_\varphi ,$$
(21)
giving

$$v_r(\theta )=\sqrt{\frac{1+\mu \mathrm{sin}^{2ϵ}\theta }{1+\delta \mathrm{sin}^{2ϵ}\theta }}.$$
(22)
Note that from the four remaining terms in Eq.(LABEL:18) — $`1^{st}`$, $`2^{nd}`$ $`3^{rd}`$ and $`4^{th}`$ — we obtain identities, provided the proportionality constants are chosen accordingly.

We have thus deduced general $`\theta `$-dependences for the density and hydromagnetic field, under the assumption of separation of variables, Eqs.(14), (16), (18), (20) and (22), in terms of three arbitrary constants: ($`\delta ,ϵ,\mu `$). These parameters effectively control the anisotropy in the outflow, in the magnetic field and in the density distribution.

First, the parameter $`\delta `$ is related to the ratio of the density at the equator ($`\theta =\pi /2`$) to that at the pole ($`\theta =0`$) and so the higher it is, the more the density distribution deviates from the spherically symmetric case ($`\delta =0`$).

Second, the parameter $`\mu `$ is related to the ratio of the radial kinetic energy density at the equator to that at the pole. Finally, the parameter $`ϵ`$ controls the width of the profile of the speed and density for some fixed variation between pole and equator.

Note that in the TT91 model, the simplest possible forms of $`b_r(\theta )`$, $`g(\theta )`$ and $`v_r(\theta )`$ were chosen a priori such that they were able to simulate existing observations of the solar wind. Their expressions constitute a special case of the forms deduced in this work, corresponding to $`ϵ=1`$ and $`\mu =1`$.

### 2.3 Angular momentum and angular velocity

Returning to the $`\varphi `$-component of the momentum equation, Eq.(6), it can be manipulated to give

$$\frac{}{r}\left(r\mathrm{sin}\theta V_\varphi r\mathrm{sin}\theta \frac{B_rB_\varphi }{4\pi \rho V_r}\right)=0.$$
(23)
The solution of this equation introduces a function of $`\theta `$ of the form

$$L(\theta )=r\mathrm{sin}\theta V_\varphi r\mathrm{sin}\theta \frac{B_rB_\varphi }{4\pi \rho V_r},$$
(24)
which is the total angular momentum per unit mass loss carried away by the wind, along each flow line $`\theta =\mathrm{const}`$ (see Weber & Davis 1967). The first term on the right-hand side of Eq.(24) is the angular momentum carried by the advection of the flow while the second term represents the torque associated with the magnetic stresses. Analogously, the solution to the induction equation, Eq.(2), introduces another function of $`\theta `$ of the form

$$\mathrm{\Omega }(\theta )=\frac{1}{r\mathrm{sin}\theta }\left(V_\varphi B_\varphi \frac{V_r}{B_r}\right),$$
(25)
which corresponds to the angular velocity of the roots of the field lines on the surface of the central object.

The azimuthal components of the velocity and magnetic fields can now be obtained via these two functions $`L(\theta )`$ and $`\mathrm{\Omega }(\theta )`$. After some straightforward manipulation we arrive at:

$$V_\varphi =\frac{r_0}{Y_{}}\left(\frac{Y_{}Y}{1M_A^2}\right)R\mathrm{sin}\theta \mathrm{\Omega }(\theta ),$$
(26)
$$B_\varphi =\frac{B_0r_0}{V_0Y_{}}\left(\frac{R^2/R_{}^21}{1M_A^2}\right)\frac{\mathrm{\Omega }(\theta )\mathrm{sin}\theta \sqrt{1+\delta \mathrm{sin}^{2ϵ}\theta }}{R},$$
(27)
where $`M_A^2=(V_r/V_A)^2`$ is the radial Alfvén Mach number which corresponds to the ratio of the radial velocity, $`V_r`$, to the radial Alfvénic velocity, $`V_A=B_r/\sqrt{4\pi \rho }`$, and $`(R_{},Y_{})`$ is the point for which $`M_A=1`$. At this point, $`Y_{}R_{}^2=(M_A^0)^2`$, where

$$M_A^0=\frac{V_0}{V_0^A}$$
(28)
is the ratio of the radial velocity to the Alfvénic radial velocity at the base of the wind. A physical interpretation of the radial Alfvén Mach number $`M_A`$ can be found if we write $`M_A=(\rho V_r^2/2)/(B_r^2/8\pi )`$, which represents the ratio of kinetic to magnetic energy density. If $`M_A1`$ magnetic energy dominates, whereas if $`M_A1`$ the dominant energy is kinetic. The radius $`R_{}`$ represents the distance at which magnetic energy ceases to dominate over the kinetic energy. The regularity condition at $`R_{}`$, $`M_A(R_{})=1`$, constrains the two functions $`L(\theta )`$ and $`\mathrm{\Omega }(\theta )`$: $`L(\theta )=\mathrm{\Omega }(\theta )r_{}^2\mathrm{sin}^2\theta `$.

If we now substitute Eqs.(20), (14) into Eqs.(8), (9), respectively, and compare the resulting expressions with Eqs.(26) and (27), we can deduce that $`V_0B_1=V_1B_0`$, together with the general form for $`\mathrm{\Omega }(\theta )`$,

$$\mathrm{\Omega }(\theta )=\frac{\lambda V_0Y_{}}{r_0}\frac{\mathrm{sin}^{ϵ1}\theta }{\sqrt{1+\delta \mathrm{sin}^{2ϵ}\theta }},$$
(29)
where $`\lambda =V_1/V_0=B_1/B_0`$ measures the ratio of azimuthal to radial velocities at the base of the wind, while

$$L(\theta )=\frac{\lambda V_0r_0}{(M_A^0)^2}\frac{\mathrm{sin}^{ϵ+1}\theta }{\sqrt{1+\delta \mathrm{sin}^{2ϵ}\theta }}.$$
(30)
The angular dependence of $`\mathrm{\Omega }(\theta )`$ as the parameters $`ϵ`$ and $`\delta `$ vary is worth some attention and is plotted in Fig. 2. This behaviour of $`\mathrm{\Omega }(\theta ;ϵ,\delta )`$ which emerges naturally from the separation of the variables in the governing equations can be compared, for example, with the rotation law of sunspots and solar photospheric magnetic fields. Except for the regions around the poles (since from Eq. (29), $`\mathrm{\Omega }(\theta )`$ goes to zero for $`\theta =0`$ and $`180^{}`$), the graphs of Figs. 2(a,b,c) for $`ϵ>1`$ show similarities with the angular dependence of $`\mathrm{\Omega }(\theta )`$ arising from the yearly averaged rotation profiles of photospheric magnetic fields (Snodgrass 1983).

### 2.4 Mass and angular momentum efflux

A quantity of interest is the mass efflux ($`\rho V_rr^2`$) or, equivalently, the mass loss rate per infinitesimal solid angle $`d\mathrm{\Sigma }`$ at the angle $`\theta `$ (see TT91). From the equation of conservation of mass Eq.(3) this has to be a function of $`\theta `$ alone, which we shall denote by $`\dot{m}(\theta )`$

$$\dot{m}(\theta )=\rho _0V_0r_0^2\sqrt{(1+\mu \mathrm{sin}^{2ϵ}\theta )(1+\delta \mathrm{sin}^{2ϵ}\theta )}.$$
(31)
Figure 3 shows the variation of the mass efflux with latitude. It vanishes at the equator only for $`\mu =1`$.

For a particular $`\delta `$, if $`|\mu |<\delta /(2\delta +1)`$ the maximum of $`\dot{m}`$ occurs for $`\theta =90^{}`$, while if $`|\mu |>\delta /(2\delta +1)`$ it occurs for $`0<\theta <90^{}`$. This is shown in Fig. 3(b) for $`\delta =4`$. We simply note in passing that the first of these angular dependences of the mass loss is reminiscent of some observed intense mass loss rates that are thought to occur through equatorial stellar winds (e.g. from Be stars observed equator-on). Also, measurements of the Lyman alpha emission by the satellites Mariner 10, Prognoz, Ulysses and SoHO (Bertaux et al. 1997) have implied that there is higher mass efflux at the equator than at the pole.

We may next introduce the angular momentum efflux $`\dot{l}(\theta )\rho V_rr^2L(\theta )`$ which is the angular momentum loss rate per infinitesimal solid angle $`d\mathrm{\Sigma }`$ at the angle $`\theta `$,

$$\dot{l}(\theta )=\frac{\lambda \rho _0V_0^2r_0^3}{(M_A^0)^2}\mathrm{sin}^{ϵ+1}\theta \sqrt{1+\mu \mathrm{sin}^{2ϵ}\theta }.$$
(32)
The variation of $`\dot{l}(\theta )`$ with latitude is shown in Fig. 4. As with the mass efflux, the angular momentum efflux vanishes at the equator for $`\mu =1`$.

In TT91 the angular momentum efflux was assumed a priori, from which the simplest possible forms for the azimuthal hydromagnetic field were obtained.

### 2.5 Hydromagnetic field and density

The azimuthal components can now be obtained

$$V_\varphi (R,\theta )=\lambda V_0\frac{R\mathrm{sin}^ϵ\theta }{\sqrt{1+\delta \mathrm{sin}^{2ϵ}\theta }}\left(\frac{Y_{}Y}{1M_A^2}\right),$$
(33)
$$B_\varphi (R,\theta )=\lambda B_0\frac{\mathrm{sin}^ϵ\theta }{R}\left(\frac{R^2/R_{}^21}{1M_A^2}\right).$$
(34)
Note that $`V_\varphi `$ is maximum at the equator. We recall that the radial hydromagnetic field and density are of the form

$$V_r(R,\theta )=V_0Y(R)\sqrt{\frac{1+\mu \mathrm{sin}^{2ϵ}\theta }{1+\delta \mathrm{sin}^{2ϵ}\theta }},$$
(35)
$$B_r(R,\theta )=\frac{B_0}{R^2}\sqrt{1+\mu \mathrm{sin}^{2ϵ}\theta },$$
(36)
$$\rho (R,\theta )=\frac{\rho _0}{YR^2}\left(1+\delta \mathrm{sin}^{2ϵ}\theta \right).$$
(37)
### 2.6 Balance of forces

At this stage, we still have to deduce an equation for $`Y(R)`$, together with the form for the pressure $`p(R,\theta )`$ and hence the temperature $`T(R,\theta )`$. Going back to the $`r`$\- and $`\theta `$-components of the momentum equation – Eqs.(4) and (5) – and assuming that the variables $`R`$ and $`\theta `$ separate in these equations, we must have

$$Q(R,\theta )=Q_0(R)+Q_1(R)\mathrm{sin}^{2ϵ}\theta ,$$
(38)
where $`Q(R,\theta )`$ is the dimensionless pressure defined by

$$p(R,\theta )=\frac{\rho _0V_0^2}{2}Q(R,\theta ).$$
(39)
$`Q_0`$ represents the spherically symmetric part of the pressure, while $`Q_1`$ includes the effects of the anisotropy. Substitution of Eqs.(38), (39) into Eqs.(4) and (5) yields the following three equations for $`Q_0(R)`$, $`Q_1(R)`$ and $`Y(R)`$

$`Q_1(R)`$ $`=`$ $`{\displaystyle \frac{\mu }{(M_A^0)^2R^4}}+{\displaystyle \frac{\lambda ^2}{ϵY}}\left({\displaystyle \frac{YY_{}}{1M_A^2}}\right)^2`$ (40)
$`\left(1+{\displaystyle \frac{1}{ϵ}}\right){\displaystyle \frac{\lambda ^2}{(M_A^0)^2R^2}}\left({\displaystyle \frac{1R^2/R_{}^2}{1M_A^2}}\right)^2,`$
$`{\displaystyle \frac{\mathrm{d}Q_1}{\mathrm{d}R}}`$ $`=`$ $`{\displaystyle \frac{\delta \nu ^2}{YR^4}}{\displaystyle \frac{2\mu }{R^2}}{\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}R}}+{\displaystyle \frac{2\lambda ^2}{YR}}\left({\displaystyle \frac{YY_{}}{1M_A^2}}\right)^2`$ (41)
$`{\displaystyle \frac{\lambda ^2}{(M_A^0)^2R^2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}R}}\left({\displaystyle \frac{1R^2/R_{}^2}{1M_A^2}}\right)^2,`$
$`{\displaystyle \frac{\mathrm{d}Q_0}{\mathrm{d}R}}`$ $`=`$ $`{\displaystyle \frac{\nu ^2}{YR^4}}{\displaystyle \frac{2}{R^2}}{\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}R}},`$ (42)
in which $`\nu `$ is the ratio of the escape speed to the radial speed at the base of the outflow

$$\nu =\frac{V_{esc}}{V_0}=\frac{\sqrt{2GM/r_0}}{V_0}.$$
(43)
To understand the interplay between different forces involved in the mechanism of this type of wind, let us describe one by one the various terms in the above equations. All of them are written so that the pressure gradient term is isolated on the left-hand side. Equation (40) represents the force equilibrium across the field lines. On the right-hand side of this, the various terms represent, respectively, the magnetic pressure, the centrifugal force and the magnetic tension. The last two equations express the equilibrium of forces along the radial direction. Equation (41) shows the anisotropic terms, namely, the anisotropic part of the gravitational and inertial forces, the centrifugal force and the magnetic tension term, respectively. The isotropic terms are shown in Eq.(42). These are related to the gravitational and inertial forces, respectively. Equations (40) and (41) can be combined to give a single expression for $`Y(R)`$

$$\frac{\mathrm{d}Y}{\mathrm{d}R}=\frac{F(R)}{G(R)},$$
(44)
where

$`F(R)`$ $`={\displaystyle \frac{\delta \nu ^2}{YR^4}}+{\displaystyle \frac{4\mu }{(M_A^0)^2R^5}}+{\displaystyle \frac{2\lambda ^2}{ϵ}}{\displaystyle \frac{Y}{RM_A^2(1M_A^2)^2}}\times `$ (45)
$`\left[{\displaystyle \frac{(1+ϵ)M_A^2ϵ}{M_A^2}}{\displaystyle \frac{R^4}{R_{}^4}}\left((2+ϵ)M_A^2(1+ϵ)\right)\right]`$
$`G(R)`$ $`={\displaystyle \frac{2\mu (M_A^0)^2Y}{M_A^2}}{\displaystyle \frac{\lambda ^2}{ϵ(1M_A^2)^2}}\times `$ (46)
$`\left[{\displaystyle \frac{2M_A^21}{M_A^4}}{\displaystyle \frac{R^4}{R_{}^4}}1\right]`$
The above differential equation for $`Y(R)`$ requires a boundary condition. For convenience, let us choose $`Y(1)=1`$, which defines $`V_0=V_r(R=1,\theta =0)`$ from Eq.(35). ¿From the classical ideal gas law we can now express temperature as

$$T(R,\theta )=\frac{mV_0^2YR^2}{2k_B(1+\delta \mathrm{sin}^{2ϵ}\theta )}Q(R,\theta ).$$
(47)
## 3 Parametric study of the solution

### 3.1 Critical points

In order to determine the full solution we need to solve Eq.(44), (45), (46) for $`Y(R)`$. This is a first-order non-linear differential equation which can be integrated numerically, using a standard routine. A first inspection shows that this equation has a singular point at $`R=R_{}`$, $`Y=Y_{}`$, where $`M_A=1`$, or in other words where the radial velocity equals the radial Alfvénic velocity $`V_A=B_r/\sqrt{4\pi \rho }`$. The term singular point is used here in the sense that both the numerator and denominator must vanish there. In a more general geometry with meridional components, this point, known as the Alfvénic point, corresponds to the location where the poloidal speed equals the poloidal Alfvén speed. It is present in magnetic wind models (c.f. Weber & Davis 1967, Mestel 1968) and is a consequence of the steady-state assumption. It delimits the magnetically dominated region beyond which the torque exerted by the magnetic field ceases to dominate over the angular momentum carried by the fluid. The exact position of this point cannot be found analytically in this case, which complicates any type of numerical treatment as we shall see later. In a similar problem, TT91 find that expanding up to fourth order around $`(R_{},Y_{})`$, all slopes are allowed. We are thus in the presence of an improper node or star-type point. This means that no particular solution is filtered out by the presence of this singular point.

There is a second singular point, say at $`(R_X,Y_X)`$, found by satisfying simultaneously $`F=G=0`$ in Eq.(44). A first-order analysis around this point shows that only two slopes are allowed, one positive and one negative, giving rise to a saddle or $`X`$type point. Thus, this point is also a critical point in the sense that it chooses a critical solution. In a recent work, Tsinganos et al. (1996) have shown that, in general, at the critical points of axisymmetric and self-similar solutions, the component of the flow velocity perpendicular to the directions of axisymmetry (the $`\varphi `$-direction in this model) and self-similarity (the $`\theta `$-direction in this model), equals the characteristic slow/fast MHD wave speed in that direction,

$$V_r^4V_r^2(V_a^2+C_s^2)+C_s^2V_{a,r}^2=0,$$
(48)
with $`C_s`$ the sound speed and $`V_a`$, $`V_{a,r}`$ the total and radial components of the Alfvén speeds, respectively. The above equation is the equation for the speed of an MHD wave propagating in the r-direction.

Thus, at the critical point where Eq. (48) is satisfied, the $`r`$-component of the flow speed equals to the slow or the fast MHD mode wave speed in that direction. Now, in low magnetisation plasmas and in the case of parallel wave propagation, the slow mode wave speed coincides with the Alfven speed. On the other hand, in high magnetization plasmas and in the case of parallel wave propagation, it is the fast mode wave speed that coincides with the Alfven speed. In the present case we have a degeneracy in the sense that there are only two critical surfaces in view of the assumption that the critical surfaces are spherical (which seems to be a good assumption, at least in the case of the solar wind - Exarhos & Moussas 2000). Note that due to the non-polytropic assumption, the sound speed is ill-defined and we cannot easily determine whether we are in a low- or high-magnetisation regime. We conjecture that here we have a situation analoguous to low magnetisation plasmas and at the X-type critical point the super-Alfvenic radial flow speed equals to the radial component of the fast MHD mode wave speed, while the slow critical transition coincides everywhere with the Alfvénic transition because of self-similarity.

We should stress at this stage that both the Alfvénic point $`(R_{},Y_{})`$ and this $`X`$type point are (loosely) called equilibrium or critical points since they satisfy the regularity condition $`F(R)=G(R)=0`$. However, since the first one does not filter any solution (i.e. all slopes are allowed through it), while the second one (as for any saddle point) selects a particular solution sometimes referred to as the critical solution, only the X-type singularity corresponds to a true critical point.

We performed a numerical integration of Eqs.(44), (45), (46). The general topology is shown in Fig. 5 for a highly magnetized medium and for representative values of $`\mu `$ and $`ϵ`$. The position of the critical points was found by an iterative procedure. The integration yields a single solution crossing the $`X`$type point with positive slope and satisfying the boundary condition $`Y(1)=1`$. For this critical solution, the flow starts near the star with low speeds and connects to large distances where it attains large super-Alfvénic velocities. There is another critical solution, which crosses the X-type point with negative slope and is always decelerating. There is a limiting value of $`ϵ`$ above which we couldn’t find a wind solution satisfying the boundary condition $`Y(1)=1`$. For example, if $`\lambda =0.5`$, $`\nu =120`$ and $`\delta =4`$, there is no solution for $`ϵ\stackrel{_>}{_{}}5`$, if $`M_A^0=0.1`$ and for $`ϵ\stackrel{_>}{_{}}2`$, if $`M_A^0=0.01`$.

Once the positions of both $`R_{}`$ and $`R_X`$ are known, all other types of solutions can be easily found. They include breeze-type solutions that cross the Alfvénic point but not the X-type point and reach sub-Alfvénic asymptotic speeds for large distances.

Table 1 shows the positions of both singular points for $`M_A^0=0.1`$ and $`M_A^0=0.01`$ and different values of $`ϵ`$ and $`\mu `$. We have taken $`\lambda =0.5`$, $`\nu =120`$ and $`\delta =4`$.

### 3.2 Radial velocity

The initial acceleration decreases as $`M_A^0`$ decreases, Figs. 6(a, b), and also as $`\lambda `$ increases. In other words, the more the outflow is magnetised and rotating, the smaller is the initial acceleration and both the magnetic field and rotation are inhibiters of the initial acceleration of the outflow. This may be understood using a simple argument, as follows.

Assume without loss of generality that $`\delta `$ = 0, such that the weight of a parcel of plasma is the same at any angle $`\theta `$, at some fixed radial distance $`R`$. Also for simplicity assume that $`\mu =1`$ and $`ϵ=1`$. Then, because of the sinusoidal dependence of the radial flow speed with the latitude, at the equator we have static conditions such that the weight of the plasma is balanced there mainly by the centrifugal force and the radial component of the Lorenz force (the centrifugal force is always positive, i.e., outwards, while in most of the cases examined the radial component of the magnetic force is also positive). Then, moving at the pole and at the same radial distance $`R`$, we find the same inwards gravitational plasma weight (since $`\delta =0`$), but no centrifugal and radial magnetic forces acting there due to their $`\text{sin}^2\theta `$ dependence. The net result is a decelerating force, relatively to the cases $`M_A^0`$ = $`\lambda `$ = 0. For a mildly magnetized case ($`M_A^0=0.1`$), Fig. 6(b) shows that the initial acceleration is very high and the radial velocity rapidly attains its asymptotic value.

The asymptotic form of $`Y(R)`$ for large $`R`$ can be obtained after integrating Eq.(44)

$$Y^33\left(1+\frac{1}{ϵ}\right)\frac{\lambda ^2}{(M_A^0)^6|\mu |R_{}^4}\mathrm{ln}R.$$
(49)
For $`M_A^01`$, the asymptotic speed increases logarithmically with $`R`$, and $`Y`$ increases with a decrease of $`M_A^0`$ or a decrease of $`ϵ`$ and $`|\mu |`$. When $`M_A^0>1`$, the logarithmic derivative $`\mathrm{d}Y^3/\mathrm{d}\mathrm{ln}R`$, for large $`R`$, is very small (behaving like $`(M_A^0)^6`$), and $`Y`$ can be taken as constant. The flow is hydrodynamically dominated from $`R=1`$, while both singular points, $`R_{}`$ and $`R_X`$, almost coincide in the region $`R<1`$ (i.e., below the base of the wind). We recall that in the hydrodynamic case (Paper I), $`Y`$ approaches rapidly a constant value, independent of $`ϵ`$.

The cautious reader may have already noted that in Eq. (49) the asymptotic radial speed is the product of a logarithmically diverging isothermal wind speed which increases like $`(\mathrm{ln}R)^{1/3}`$ (see TT91) and Michel’s characteristic speed ($`\mathrm{\Omega }^2F_B^2/\dot{M})^{1/3}`$, where $`F_B`$ and $`\dot{M}`$ are the magnetic and mass fluxes (MacGregor 1996). This is also checked by the fact that the effective polytropic index $`\gamma (R\mathrm{})1`$, as seen in Fig. 7.

### 3.3 Temperature and effective polytropic index

In this study by not using a polytropic relationship between pressure and density ($`p\rho ^\gamma `$, with $`\gamma `$ the constant polytropic index) we avoided constraining the exchange of energy so as to keep $`\gamma `$ constant. In our approach, we can define an effective polytropic index given by

$$\gamma \left[\frac{lnp}{ln\rho }\right]_{A=\text{const.}}$$
(50)
for each field line $`A(\theta )=\text{const}`$. (Weber 1970, TT91). This effective polytropic index is no longer a constant, but instead a function of $`R`$. The form of $`\gamma (R)`$ is closely associated with the variation of temperature with distance from the central object. In particular, if $`\gamma <1`$ there is intense heating and the temperature increases, whereas if $`\gamma >1`$ there is a depletion of heating and the temperature decreases (see Fig. 7).

In Fig. 8 we have chosen to illustrate the two-dimensional character of the solution deduced in this paper, with parameters corresponding to a mildly magnetized object ($`M_A^0=0.1`$), by plotting the latitudinal dependence of the radial velocity and density. The velocity, is higher at the polar axis while the density is maximum at the equator. The angular dependence of these quantities is similar to the latitudinal distribution of the density and velocity in coronal holes.

## 4 Discussion

We have explored in this paper a new set of exact solutions of the steady axisymmetric MHD equations relevant to stellar wind problems.

Analogously to the case with no magnetic field (Paper I), we were able to obtain the latitudinal dependences of the different variables and these involve only the anisotropy parameters, $`\delta `$, $`ϵ`$ and $`\mu `$. The parameter $`\delta `$ is associated with the ratio of equatorial to polar density, $`ϵ`$ controls the width of the density and velocity profiles for a fixed variation between pole and equator, while $`\mu `$ is connected to the latitudinal anisotropy of the distribution of radial kinetic energy density. The other three parameters of the problem, which come about in the solution for the radial functions, are $`\lambda `$, which controls the strength of rotation, $`\nu `$, a measure of the escape speed from the central object and $`M_A^0`$ which is related the strength of the magnetic pressure.

The degree of collimation, together with the anisotropy in the density distribution, is effectively controlled by $`\delta `$, $`ϵ`$ and $`\mu `$. These solutions could be pertinent to the plasma dynamics in polar coronal holes, stellar and extragalactic jets and even star-forming regions. The solution with$`\mu =1`$ has zero radial velocity at the equator and a continuous radial magnetic field across the equator. A discontinuity and associated current sheet arises for $`\mu 1`$. Also, the angular velocity of the roots of the field lines at the surface of the central object is such that it resembles qualitatively solar observations only for $`ϵ>1`$.

The topology of the solutions is controlled by two distinct singular points. The first one is the usual Alfvénic point and corresponds to the location where the radial velocity of the outflow equals the radial Alfvénic velocity. All solutions pass through this high-order singularity and it corresponds to an improper node or star-type point.

There is a critical point present downstream from this Alfvén transition. It is an $`X`$-type critical (or saddle) point which is responsible for filtering a single solution corresponding to a vanishing pressure at infinity — the wind-type solution. In a similar problem, Tsinganos et al. (1996) have shown that at this new critical point, the $`r`$-component of the flow speed equals the slow/fast MHD mode wave speeds in that direction.

Because the field/streamlines are constrained to keep a helicoidal geometry, the pressure gradient for large distances has to be sufficiently large to balance the dominant tension force. Thus, in some cases the solution is not valid outside a cone around the polar axis. This limitation disappears for $`M_A^00.05`$.

There is a drastic change in the nature of the solution from high to lower values of $`M_A^0`$. The solutions for high $`M_A^0`$ manifest all the characteristics of a typical hydrodynamic wind, namely very large acceleration at the base of the atmosphere and temperature decaying very rapidly with distance. The magnetically dominated cases ($`M_A^01`$) show lower acceleration at the base of the atmosphere and an almost isothermal atmosphere at larger distances.

The density anisotropy ($`\delta `$) greatly favours the acceleration of the wind close to the base, while the strength of rotation ($`\lambda `$) and the gravitational field ($`\nu `$) slows down the initial speed. The solution far away from the source seems to be quite insensitive to $`ϵ`$. The influence of $`\mu `$, on the other hand, is such that a decrease of $`|\mu |`$ not only diminishes the relative importance of the magnetic effects by decreasing the size of the magnetic lever arm, $`R_{}`$, but also increases the value of the initial acceleration as well as the asymptotic radial velocity at large distances. The asymptotic plasma temperature increases for lower values of $`|\mu |`$.

This simple but self-consistent two-dimensional solution of the MHD equations has nevertheless some limitations. The assumption of separation of variables could be dropped, although we should not forget the nonlinear character of the governing coupled partial differential equations. The assumption of the neglect of the meridional components of the magnetic field and flow ($`V_\theta ,B_\theta `$) which limits their possible applications, could also be eliminated within the framework of the method of separation of the variables. Finally, as we’ve already discussed, the present investigation does not take into account the detailed energetics in the outflow. Nevertheless, the present study indicates that inherently non-spherically symmetric solutions of the MHD equations which do not have the previous limitations may exist. This ambitious undertaking to find them remains a challenge for the future.

###### Acknowledgements.

J. Lima wishes to thank C. Sauty for his help and support and also the staff at FORTH in Crete for their kind hospitality. This work was supported by Programa Praxis XXI of the Fundação para a Ciência e a Tecnologia, by funds from the UK Science and Engineering Research Council, the Erasmus Programme and the EEC RT Network HPRN-CT-2000-00153.


text:
# Can Schwarzschildean gravitational fields suppress gravitational waves?

## I  Introduction

Backscattering has been investigated for a long time for various wave equations (see, for instance, ). In general relativity, this topic has been studied since early 1960’s (, ). This paper continues the programme that started with the study of the backscatter of scalar  and electromagnetic fields ( and ). Here we investigate the propagation of even-parity gravitational waves in a (fixed) background Schwarzschild spacetime, assuming a nonstationary source. The discussion, however, is done without any reference to the source. We only deal with field quantities. It is assumed that initial data are those of an isolated pulse (burst) of a gravitational wave. The main question that is answered is what fraction of the initially outgoing radiation may undergo backscattering before reaching null infinity? The strength of the backscattering is assessed by bounding the fraction of the initial burst energy that will not reach a distant observer in the main pulse.

The even-parity waves are the only waves which are radiated during the axisymmetric collision of nonspinning black holes , and since in this case the Schwarzschild spacetime is a valid starting point for an approximation scheme, it gives us an opportunity to bound the strength of the phenomenon in a fairly realistic astrophysical context.

The following five sections of the paper give a theoretical description of the backscattering effect. The Sec. II brings notation and the Zerilli equation. In Sec. III, the initial data are bounded by the initial energy and the solution is sought in the form of superposition of an outgoing radiation (defined by initial data) and a backscattered term. The evolution of the backscattered term can be bounded by solutions of two differential inequalities. The bounds that are derived in Sec. IV and in the Appendix deal with a general situation; no assumption is made about the initial radiation. In Sec. V, we discuss initial data that are of compact support and in addition, the relative width of the support is small. Such data correspond to radiation that is dominated by short wavelengths. In this case stronger estimates are derived. They imply that in the limit of short wavelengths (relative width of the support tending to zero) the backscattering effect becomes negligible. In Sec. VI, the ”small relative width condition” of Sec. V is supplemented by the assumption that the initial burst is far away from the horizon.

## II Formalism

The space-time geometry is defined by a Schwarzschildean line element,

$$ds^2=(1\frac{2m}{R})dt^2+\frac{1}{1\frac{2m}{R}}dR^2+R^2d\mathrm{\Omega }^2,$$
(1)
where $`t`$ is a time coordinate, $`R`$ is a radial coordinate that coincides with the areal radius and $`d\mathrm{\Omega }^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2`$ is the line element on the unit sphere, $`0\varphi <2\pi `$ and $`0\theta \pi `$. Throughout this paper $`G`$, the Newtonian gravitational constant, and $`c`$, the velocity of light are put equal to 1.

As explained in the Introdction, we restrict ourselves, to the even-parity axial perturbations. Their propagation is ruled, in the linear approximation, by the Zerilli equation . Formulated in terms of the gauge-invariant amplitude $`\mathrm{\Psi }`$ defined by Moncrief , this equation reads, in the case of $`l=2`$ multipole ,

$$(_t^2+_r^{}^2)\mathrm{\Psi }=V\mathrm{\Psi },$$
(2)
where the potential $`V`$ is given by

$$V(R)=6(1\frac{2m}{R})^2\frac{1}{R^2}+(1\frac{2m}{R})\frac{63m^2(1+\frac{m}{R})}{2R^4(1+\frac{3m}{2R})^2},$$
(3)
and where

$$r^{}=R+2m\mathrm{ln}(\frac{R}{2m}1)$$
(4)
is the tortoise radial coordinate.

Consider a set of functions $`\mathrm{\Psi }_i(r^{}t)`$, $`i=0,1,2`$, that satisfy the following linear relations

$`_t\mathrm{\Psi }_1=3\mathrm{\Psi }_0,`$ (5)
$`_t\mathrm{\Psi }_2=\mathrm{\Psi }_1m_t\mathrm{\Psi }_1.`$ (6)
The combination

$$\stackrel{~}{\mathrm{\Psi }}\mathrm{\Psi }_0(r^{}t)+\frac{\mathrm{\Psi }_1(r^{}t)}{R}+\frac{\mathrm{\Psi }_2(r^{}t)}{R^2}$$
(7)
solves Eq. (2) in Minkowski space-time ($`m=0`$); it represents purely outgoing radiation.

Let the initial data of a solution $`\mathrm{\Psi }`$ of (2) coincide with $`\stackrel{~}{\mathrm{\Psi }}`$ at $`t=0`$. Then, initially, $`\mathrm{\Psi }`$ represents a purely outgoing wave. It should be noted that the assumption that initial data are (initially) purely outgoing is made in this paper only for the sake of clear presentation. In the linear approximation the propagation of the initially outgoing radiation is independent of whether or not ingoing radiation is present.

We decompose the sought solution $`\mathrm{\Psi }(r^{},t)`$ into the known part $`\stackrel{~}{\mathrm{\Psi }}`$ and an unknown function $`\delta `$

$$\mathrm{\Psi }=\stackrel{~}{\mathrm{\Psi }}+\delta .$$
(8)
Due to the choice of the initial data made above one has $`\delta =_t\delta =0`$, at $`t=0`$.

## III Energy estimates

Let us assume that the quadrupole initial data are defined by a smooth triad of the functions $`\mathrm{\Psi }_k`$ ($`k=0,1,2`$) with the initial support $`[a,b]`$ ($`b<\mathrm{}`$). That guarantees that the initial energy density multiplied by $`R^2`$, $`\rho =\left((_t\mathrm{\Psi })^2+(_r^{}\mathrm{\Psi })^2+V\mathrm{\Psi }^2\right)/\eta _R`$, is smooth and vanishes on the boundary $`a`$. Here $`\eta _R=12m/R`$ holds.

The energy content inside a part of a Cauchy hypersurface $`\mathrm{\Sigma }_t`$ that is exterior to a ball of a radius $`R`$ can be defined as $`E(R,t)_R^{\mathrm{}}𝑑r\rho (r,t)`$. Let us point out that in order to ensure a proper normalization of the energy flux at infinity, there should be a normalization constant in the definition of the energy $`E(R,t)`$. We decided to omit it, since later on we shall be interested only in the relative efficiency of the backscatter and thus the normalization factor cancels out. The total initial energy corresponding to the hitherto defined initial data is equal to $`E(a,0)`$.

Lemma 1. Defining

$`C_1{\displaystyle \frac{3}{2}}\sqrt{(2+\sqrt{2/3})E(a,0)},`$ (9)
$`C_2\sqrt{2E(a,0)},`$ (10)
$`C_3\sqrt{{\displaystyle \frac{2E(a,0)}{\eta _a(2\sqrt{6}+1)}}},`$ (11)
and introducing the two nonnegative functions

$`g_1(R)`$ $`=`$ $`\mathrm{ln}\left({\displaystyle \frac{2m+R}{a2m}}\right)+32m^5\left({\displaystyle \frac{1}{5(2m+R)^5}}+{\displaystyle \frac{1}{5(2m+a)^5}}\right)+`$ (14)
$`20m^4\left({\displaystyle \frac{1}{(2m+R)^4}}+{\displaystyle \frac{1}{(2m+a)^4}}\right)+80m^3\left({\displaystyle \frac{1}{3(2m+R)^3}}+{\displaystyle \frac{1}{3(2m+a)^3}}\right)+`$
$`20m^2\left({\displaystyle \frac{1}{(2m+R)^2}}+{\displaystyle \frac{1}{(2m+a)^2}}\right)+10m\left({\displaystyle \frac{1}{2m+R}}+{\displaystyle \frac{1}{2m+a}}\right)`$
and

$`g_2(R)`$ $`=`$ $`Ra+16m^4\left({\displaystyle \frac{1}{3(2m+R)^3}}+{\displaystyle \frac{1}{3(2m+a)^3}}\right)+`$ (17)
$`(16m^3)\left({\displaystyle \frac{1}{(2m+R)^2}}+{\displaystyle \frac{1}{(2m+a)^2}}\right)+`$
$`(24m^2)\left({\displaystyle \frac{1}{2m+R}}+{\displaystyle \frac{1}{2m+a}}\right)+8m\mathrm{ln}\left({\displaystyle \frac{2m+R}{2m+a}}\right),`$
the following inequalities hold at $`t=0`$ and for $`Ra`$:

$`{\displaystyle \frac{|\mathrm{\Psi }_1(R)|}{R^{3/2}}}C_1\eta _R^{3/2}\sqrt{g_1(R)},`$ (18)
$`{\displaystyle \frac{|\mathrm{\Psi }_2(R)|}{\sqrt{\eta _R}R^2}}C_2\sqrt{g_2(R)}+C_1{\displaystyle \frac{6m}{\sqrt{a}}}\sqrt{g_1(R)}\left(1\sqrt{{\displaystyle \frac{a}{R}}}\right)+`$ (19)
$`6C_3{\displaystyle \frac{m}{\sqrt{a}}}\sqrt{1\left({\displaystyle \frac{a}{R}}\right)^{2\sqrt{6}+1}}\left({\displaystyle \frac{1}{\sqrt{\eta _a}}}{\displaystyle \frac{1}{\sqrt{\frac{R}{a}\frac{2m}{a}}}}\right),`$ (20)
$`{\displaystyle \frac{|\stackrel{~}{\mathrm{\Psi }}(R)|}{\sqrt{R}}}C_3\sqrt{1\left({\displaystyle \frac{a}{R}}\right)^{2\sqrt{6}+1}},`$ (21)
$`|{\displaystyle \frac{\mathrm{\Psi }_0(R)|}{\sqrt{R}}}C_3\sqrt{1\left({\displaystyle \frac{a}{R}}\right)^{2\sqrt{6}+1}}+`$ (22)
$`C_1\eta _R^{3/2}\sqrt{g_1(R)}+C_2\sqrt{{\displaystyle \frac{g_2(R)}{R}}}+C_1{\displaystyle \frac{6m}{\sqrt{aR}}}\sqrt{g_1(R)}(1\sqrt{{\displaystyle \frac{a}{R}}})+`$ (23)
$`6C_3{\displaystyle \frac{m}{\sqrt{aR}}}\sqrt{1\left({\displaystyle \frac{a}{R}}\right)^{2\sqrt{6}+1}}\left({\displaystyle \frac{1}{\sqrt{\eta _a}}}{\displaystyle \frac{1}{\sqrt{\frac{R}{a}\frac{2m}{a}}}}\right).`$ (24)
Proof.

One can explicitly verify that

$`\eta _R\left({\displaystyle \frac{\mathrm{\Psi }_1}{R^2}}+2{\displaystyle \frac{\mathrm{\Psi }_2}{R^3}}\right)=_t\stackrel{~}{\mathrm{\Psi }}+_r^{}\stackrel{~}{\mathrm{\Psi }},`$ (25)
$`\eta _R\left(2{\displaystyle \frac{\mathrm{\Psi }_0}{R}}+{\displaystyle \frac{\mathrm{\Psi }_1}{R^2}}\right)=_t\stackrel{~}{\mathrm{\Psi }}+_r^{}\stackrel{~}{\mathrm{\Psi }}+{\displaystyle \frac{2\eta _R\stackrel{~}{\mathrm{\Psi }}}{R}}.`$ (26)
The equations (26), using the relations (7), result in

$`_R\left(\eta _R^{3/2}{\displaystyle \frac{\mathrm{\Psi }_1}{R^{3/2}}}\right)={\displaystyle \frac{3}{2R^{1/2}\eta _R^3}}\left({\displaystyle \frac{_t\stackrel{~}{\mathrm{\Psi }}}{\eta _R^{1/2}}}+\sqrt{\eta _R}_R\stackrel{~}{\mathrm{\Psi }}+{\displaystyle \frac{2\stackrel{~}{\mathrm{\Psi }}\sqrt{\eta _R}}{R}}\right),`$ (27)
$`_R\left(\eta _R^{1/2}{\displaystyle \frac{\mathrm{\Psi }_2}{R^2}}\right)={\displaystyle \frac{1}{\eta _R^2}}\left({\displaystyle \frac{_t\stackrel{~}{\mathrm{\Psi }}}{\eta __R^{1/2}}}+\sqrt{\eta _R}_R\stackrel{~}{\mathrm{\Psi }}\right)+{\displaystyle \frac{3m}{R\eta __R^{3/2}}}\left({\displaystyle \frac{\stackrel{~}{\mathrm{\Psi }}}{R}}{\displaystyle \frac{\mathrm{\Psi }_1}{R^2}}\right).`$ (28)
The integration from $`a`$ to $`R`$ and the use of the Schwarz inequality yields

$`{\displaystyle \frac{|\mathrm{\Psi }_1(R)|}{R^{3/2}}}{\displaystyle \frac{3\eta _R^{3/2}}{2}}\sqrt{(2+\sqrt{2/3})E(a,0)}\left({\displaystyle _a^R}𝑑r{\displaystyle \frac{1}{r\eta _r^6}}\right)^{1/2}.`$ (29)
Integrating $`_a^R𝑑r\frac{1}{r\eta _r^6}`$, one immediately arrives at the first of the postulated inequalities.

In order to show the third inequality, notice that $`|\stackrel{~}{\mathrm{\Psi }}(R)|R^\sqrt{6}=|_a^R𝑑r_r(\stackrel{~}{\mathrm{\Psi }}(r)r^\sqrt{6})|`$. The latter expression is bounded from above, using the Schwarz inequality, by

$`\sqrt{2{\displaystyle _a^R}𝑑r\left(\eta _r(_r\stackrel{~}{\mathrm{\Psi }})^2+6\eta _r\stackrel{~}{\mathrm{\Psi }}^2/r^2\right)}\sqrt{{\displaystyle _a^R}𝑑r\eta _r^1r^{2\sqrt{6}}}`$ (30)
$`\sqrt{{\displaystyle \frac{2E(a,0)}{\eta _a(2\sqrt{6}+1)}}}R^{\sqrt{6}+0.5}\left(1\left({\displaystyle \frac{a}{R}}\right)^{\sqrt{6}+1}\right)^{1/2},`$ (31)
where the inequality in Eq. (31) follows from the monotonicity of the energy as function of $`R`$. The first factor on the left hand side of this inequality is not greater than $`\sqrt{2E(a,0)}`$ since $`6\eta _r/r^2V(r)/\eta _r`$. The replacement of $`\eta _r^1`$ by $`\eta _a^1`$ and the integration of the other factor leads to the desired result.

The second of the equations of (28) can be integrated. The Schwarz inequality and direct integration as well as the $`\stackrel{~}{\mathrm{\Psi }}`$ and $`\mathrm{\Psi }_1`$ estimates should be used in order to get the second inequality of Lemma 1. The $`\mathrm{\Psi }_0`$-estimate, in turn, follows from the identity $`\mathrm{\Psi }_0=\stackrel{~}{\mathrm{\Psi }}\mathrm{\Psi }_1/R\mathrm{\Psi }_2/R^2`$ and the preceding estimates.

## IV The estimate of the diffused energy

Let us define the strength of the backscattered radiation that is directed inward by

$$h_{}(R,t)=\frac{1}{\eta _R}(_t+_r^{})\delta (R,t).$$
(32)
Let the outgoing null geodesic $`\stackrel{~}{\mathrm{\Gamma }}_{(R,t)}`$ originate at $`(R,t)`$. If a point lies on the initial hypersurface, then we will write $`\stackrel{~}{\mathrm{\Gamma }}_{(R,0)}\stackrel{~}{\mathrm{\Gamma }}_R`$. By $`\stackrel{~}{\mathrm{\Gamma }}_{(R_0,t_0),(R,t)}`$ will be understood a segment of $`\stackrel{~}{\mathrm{\Gamma }}_{(R_0,t_0)}`$ ending at $`(R,t)`$.

A straightforward calculation shows that the rate of the energy change along $`\stackrel{~}{\mathrm{\Gamma }}_a`$ is given by

$`(_t+_r^{})E(R,t)=`$ (33)
$`[\eta _R^2h_{}^2(R,t)+V\delta ^2(a)].`$ (34)
It is necessary to point out that in the case of the initial point $`R_0>a`$ the result would be more complicated; the differentiation of the energy along $`\stackrel{~}{\mathrm{\Gamma }}_{R_0}`$ would depend also on $`\mathrm{\Psi }_0,\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$. If, however, the outgoing null geodesics is $`\stackrel{~}{\mathrm{\Gamma }}_a`$, then it starts from $`a`$ where $`\mathrm{\Psi }_0,\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$ do vanish. Since these functions depend on the difference $`r^{}t`$, their values along outgoing geodesics are constant, and that allows one to conclude that they vanish at $`\stackrel{~}{\mathrm{\Gamma }}_a`$.

The energy loss, that is the amount of energy that diffused inward $`\stackrel{~}{\mathrm{\Gamma }}_a`$ is equal to a line integral along $`\stackrel{~}{\mathrm{\Gamma }}_a`$,

$`\delta E_aE(a,0)E_{\mathrm{}}=`$ (35)
$`{\displaystyle _a^{\mathrm{}}}dr[\eta _rh_{}^2+{\displaystyle \frac{V\delta ^2}{\eta _r}}].`$ (36)
Our goal is to find an estimate of $`\delta E_a`$ of a single pulse of radiation basing only on the information about the position and the energy of the initial pulse. Obviously, $`0\delta E_aE(a,0)`$ holds. We are interested in deriving in this section a frequency-independent bound, but later we obtain estimates that are frequency-sensitive.

$`\delta `$ is initially zero and its evolution is governed by the following equation

$`(_t^2+_r^{}^2)\delta =V\delta +(V6{\displaystyle \frac{\eta _R^2}{R^2}})\left(\mathrm{\Psi }_0+{\displaystyle \frac{\mathrm{\Psi }_1}{R}}+{\displaystyle \frac{\mathrm{\Psi }_2}{R^2}}\right)+`$ (37)
$`{\displaystyle \frac{2m\eta _R}{R^4}}\left[3\mathrm{\Psi }_1+2{\displaystyle \frac{\mathrm{\Psi }_2}{R}}\right].`$ (38)
One can define an ”energy” $`H(R,t)`$ of the field $`\delta `$ which is contained in the exterior of a sphere of radius $`R`$ as follows

$$H(R,t)=_R^{\mathrm{}}𝑑r\left(\frac{(_t\delta )^2}{\eta _r}+\eta _r(_r\delta )^2+\delta ^2\frac{V}{\eta _r}\right).$$
(39)
The rate of change of $`H`$ along $`\stackrel{~}{\mathrm{\Gamma }}_{(R,t)}`$ is given by

$`(_t+_r^{})H(R,t)=`$ (40)
$`\eta _R\left[\eta _R\left({\displaystyle \frac{_t\delta }{\eta _R}}+_R\delta \right)^2+{\displaystyle \frac{V}{\eta _R}}\delta ^2\right]`$ (41)
$`4m{\displaystyle _R^{\mathrm{}}}𝑑r\eta _r{\displaystyle \frac{_t\delta }{r^4}}\left[3\mathrm{\Psi }_1+{\displaystyle \frac{2\mathrm{\Psi }_2}{r}}+{\displaystyle \frac{63m(1+\frac{m}{r})}{4(1+\frac{3m}{2r})^2}}\left(\mathrm{\Psi }_0+{\displaystyle \frac{\mathrm{\Psi }_1}{r}}+{\displaystyle \frac{\mathrm{\Psi }_2}{r^2}}\right)\right]`$ (42)
$`4m{\displaystyle _R^{\mathrm{}}}𝑑r{\displaystyle \frac{_t\delta }{r^4}}\left[3\mathrm{\Psi }_1+{\displaystyle \frac{2\mathrm{\Psi }_2}{r}}+{\displaystyle \frac{63m(1+\frac{m}{r})}{4(1+\frac{3m}{2r})^2}}\left(\mathrm{\Psi }_0+{\displaystyle \frac{\mathrm{\Psi }_1}{r}}+{\displaystyle \frac{\mathrm{\Psi }_2}{r^2}}\right)\right].`$ (43)
Herein the inequality follows from the omission of the nonpositive boundary term. This allows one to estimate the maximal value $`H_M`$ of the $`\delta `$-energy $`H`$, namely

$$\sqrt{H_M}10.43\frac{m\sqrt{E(a,0)}}{a}+O(m^2).$$
(44)
The calculation is essentially simple, but the algebra is quite lengthy and some numerical integrations are required. Details are relegated to the Appendix. We would like to point out that the $`O(m^2)`$ terms become dominant only when the location of the initial radiation pulse is smaller than $`6.6m`$. At $`a=15m`$ the neglected terms contribute much less than the leading term proportional to $`m`$.

Now, the integration of the first part of Eq. (43) along $`\stackrel{~}{\mathrm{\Gamma }}_{(a,0)}`$ yields

$`H(\mathrm{})H(0)=`$ (45)
$`{\displaystyle _a^{\mathrm{}}}𝑑R\left[\eta _R\left({\displaystyle \frac{_t\delta }{\eta _R}}+_R\delta \right)^2+{\displaystyle \frac{V}{\eta _R}}\delta ^2\right]`$ (46)
$`{\displaystyle _a^{\mathrm{}}}𝑑R4m{\displaystyle _R^{\mathrm{}}}𝑑r{\displaystyle \frac{_t\delta }{r^4}}\left[3\mathrm{\Psi }_1+{\displaystyle \frac{2\mathrm{\Psi }_2}{r}}+{\displaystyle \frac{63m(1+\frac{m}{r})}{4(1+\frac{3m}{2r})^2}}\left(\mathrm{\Psi }_0+{\displaystyle \frac{\mathrm{\Psi }_1}{r}}+{\displaystyle \frac{\mathrm{\Psi }_2}{r^2}}\right)\right].`$ (47)
Initially, $`H`$ vanishes (both $`\delta `$ and $`_t\delta `$ vanish) and $`H`$ is manifestly nonnegative. The first integral on the right hand side of (47) is recognized to be just $`\delta E_a`$. The second integral in turn can be shown to be bounded - using the Schwarz inequality and then the results of the Appendix - by $`2\sqrt{H_M}(10.43\frac{m\sqrt{E(a,0)}}{a}+O(m^2))`$. Thus, (47) implies

$$\delta E_a2\left(10.43\frac{m\sqrt{E(a,0)}}{a}+O(m^2)\right)\sqrt{H_M}\left[54.5\left(\frac{2m}{a}\right)^2+O(m^3)\right]E(a,0);$$
(48)
the right hand side of the first inequality achieves a maximal value when $`H`$ is maximal and that implies the second inequality.

Thus in summary, for the fraction of the energy that could diffuse through the null cone $`C_a`$, it holds:

Theorem. $`\delta E_a/E(a,0)`$ satisfies the inequality

$`{\displaystyle \frac{\delta E_a}{E(a,0)}}54.5\times \left({\displaystyle \frac{2m}{a}}\right)^2+O(m^3/a^3).`$ (49)
We would like to point out that the above derivation is more efficient and simpler than the one used in  or  when $`\delta E_a`$ was estimated directly on the basis of the estimates of $`\delta `$ and $`h_{}`$. This alternative approach would require a laborious integration of the field equation and the final estimate would be much worse.

## V The wavelength of the initial radiation and the backscatter.

In this section we shall consider the backscatter of the radiation that is initially of compact support and, in addition, the condition $`(ab)/a<<1`$ is satisfied. The leading contribution \- only terms that are quadratic in $`m^2`$ \- will be found.

Under the above conditions one infers from Eq. (24) that on the initial hypersurface

$`|\mathrm{\Psi }_1(R)|b^{3/2}C_1b^{3/2}\sqrt{g_1(R)}C_1\sqrt{{\displaystyle \frac{ba}{a}}}b^{3/2}`$ (50)
and

$$|\mathrm{\Psi }_2(R)|C_2b^2\sqrt{ba}$$
(51)
are valid. With the same accuracy the inequality (73) of the Appendix reads

$`(_t+_r^{})\sqrt{H(R,t)}2m\left({\displaystyle _R^{R(b)}}𝑑r\eta _r{\displaystyle \frac{9\mathrm{\Psi }_1^2}{r^8}}\right)^{1/2}+2m\left({\displaystyle _R^{\mathrm{}}}𝑑r{\displaystyle \frac{4\mathrm{\Psi }_2^2}{r^{10}}}\right)^{1/2}`$ (52)
$`6mC_1\sqrt{{\displaystyle \frac{ba}{a}}}b^{3/2}\left({\displaystyle _R^{R(b)}}𝑑r{\displaystyle \frac{1}{r^8}}\right)^{1/2}+4mC_2b^2\sqrt{ba}\left({\displaystyle _R^{R(b)}}𝑑r{\displaystyle \frac{1}{r^{10}}}\right)^{1/2}.`$ (53)
Herein, the integration extends from $`R`$, where $`R\stackrel{~}{\mathrm{\Gamma }}_a`$, to $`R(b)`$, which is defined by $`(R(b),t)\stackrel{~}{\mathrm{\Gamma }}_b`$. One has $`R(b)R=ba`$ up to the term $`m^0`$. The integral $`_R^{R(b)}𝑑r\frac{1}{r^8}`$ is bounded from above by $`(ba)/R^8`$ and the integral $`_R^{R(b)}𝑑r\frac{1}{r^{10}}`$ is bounded from above by $`(ba)/R^{10}`$, again to lowest order in powers of $`m`$.

Thus one arrives at

$`(_t+_r^{})\sqrt{H(R,t)}4m(ba)b^{3/2}\left({\displaystyle \frac{1.5C_1}{\sqrt{a}R^4}}+{\displaystyle \frac{C_2\sqrt{b}}{R^5}}\right).`$ (54)
The integration of this inequality along the null geodesic $`\stackrel{~}{\mathrm{\Gamma }}_a`$ yields

$`\sqrt{H_M}m(ba){\displaystyle \frac{b^{3/2}}{a^{7/2}}}\left(2C_1+C_2\right)+O(m^2).`$ (55)
Taking into account the condition that $`ba<<a`$, one arrives at

$`H_M{\displaystyle \frac{4m^2}{a^2}}\left({\displaystyle \frac{ba}{a}}\right)^2\left({\displaystyle \frac{b}{a}}\right)^3\left(C_1+{\displaystyle \frac{C_2}{2}}\right)^2+O(m^3/b^3).`$ (56)
Since the amount of backscattered energy $`\delta E_a`$ is bounded from above by $`2H_M`$, as shown in the preceding section, one finally arrives at the following estimate

$`{\displaystyle \frac{\delta E_a}{E(a,0)}}{\displaystyle \frac{8m^2}{a^2}}\left({\displaystyle \frac{ba}{a}}\right)^2\times \left({\displaystyle \frac{b}{a}}\right)^3\left({\displaystyle \frac{3}{2}}\sqrt{2+\sqrt{{\displaystyle \frac{2}{3}}}}+{\displaystyle \frac{1}{\sqrt{2}}}\right)^2+O(m^3/a^3)`$ (57)
$`84{\displaystyle \frac{m^2}{a^2}}\left({\displaystyle \frac{b}{a}}\right)^3\left({\displaystyle \frac{ba}{a}}\right)^2+O(m^3/a^3).`$ (58)
If $`(ba)/a<0.1`$, then the above formula predicts

$`{\displaystyle \frac{\delta E_a}{E(a,0)}}0.84{\displaystyle \frac{m^2}{a^2}}.`$ (59)
It is clear that if the relative width of the initial pulse tends to zero then the effect becomes negligible. This can be translated into the dependence on the wavelength of the radiation : The compression of the support of a function leads to the decrease of its wavelength scale in its Fourier transform.

A careful analysis of the higher order terms would show that they give a contribution to (58) that also scales like $`\left(\frac{ba}{a}\right)^2`$. In the case when $`a2m`$, Eq. (58) would be of the form

$$\frac{\delta E_a}{E(a,0)}C(x)\left(\frac{ba}{2m}\right)^2,$$
(60)
where $`C(x)`$ is a large number and $`x2m/a`$. One can show that $`lim_{x1}C(x)=\mathrm{}`$, but on the other hand $`C(x)`$ is fixed, when $`2m/a`$ is fixed. Thus (60) implies that when $`ba`$, then the backscatter becomes negligible. Radiation that is dominated by infinitely short wavelengths does not backscatter.

## VI More estimates on high frequency radiation

We assume initial data of compact support $`[a,b]`$. The initial energy $`E(a,0)`$ (see the beginning of Sec. III) reads, expressed in terms of functions $`\mathrm{\Psi }_0,\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$, as follows

$`E(a,0)`$ $`=`$ $`{\displaystyle _a^b}dr\rho ={\displaystyle _a^b}dr[{\displaystyle \frac{6\left(r\left(r\mathrm{\Psi }_0(r)+\mathrm{\Psi }_1(r)\right)+\mathrm{\Psi }_2(r)\right)^2}{r^6}}+(\mathrm{\Psi }_0^{}(r)+{\displaystyle \frac{r\mathrm{\Psi }_1^{}(r)+\mathrm{\Psi }_2^{}(r)}{r^2}})^2+`$ (62)
$`{\displaystyle \frac{\left(2\mathrm{\Psi }_2(r)r\left(\mathrm{\Psi }_1(r)+r(r\mathrm{\Psi }_0^{}(r)+\mathrm{\Psi }_1^{}(r))+\mathrm{\Psi }_2^{}(r)\right)\right)^2}{r^6}}].`$
The radiation energy in the wave zone is known to be $`E(a,0)=C_a^b𝑑r(\mathrm{\Psi }_0^{})^2`$. This can be compatible with (62) (modulo a normalization constant, which is not relevant here), if the terms with $`(\mathrm{\Psi }_0^{})^2`$ give a leading contribution.

One notices, that if $`\mathrm{\Psi }_\mu (R)`$ ($`\mu =0,1,2`$) are of compact support, then $`|\mathrm{\Psi }_\mu (R)|=|_a^R𝑑r_r\mathrm{\Psi }_\mu (r)|\sqrt{(Ra)_a^R𝑑r_r\mathrm{\Psi }_\mu ^2(r)}`$. Combining this with (6) one arrives at

$`3|\mathrm{\Psi }_0(r)|=|\mathrm{\Psi }_1^{}(r)|3\sqrt{ba}\sqrt{{\displaystyle _a^b}𝑑r(\mathrm{\Psi }_0^{})^2},`$ (63)
$`|\mathrm{\Psi }_1(r)|=|\mathrm{\Psi }_2^{}(r)|2(ba)^{3/2}\sqrt{{\displaystyle _a^b}𝑑r(\mathrm{\Psi }_0^{})^2},`$ (64)
$`|\mathrm{\Psi }_2(r)|0.8(ba)^{5/2}\sqrt{{\displaystyle _a^b}𝑑r(\mathrm{\Psi }_0^{})^2}.`$ (65)
Taking into account (65), one concludes that if

$$C_1(ba)/a<<1,$$
(66)
($`C_1`$ is a constant of the order of 100) then the energy is well approximated by $`E(a,0)=2_a^b𝑑r(\mathrm{\Psi }_0^{})^2`$.

In such circumstances it is clear that our analysis can be greatly simplified. First of all the contribution coming from $`\mathrm{\Psi }_2`$ to the backscatter is much smaller than that due to $`\mathrm{\Psi }_1`$; notice an additional power of $`(ba)/a`$ in the relevant estimate of (65). Secondly, the estimate (50) of $`\mathrm{\Psi }_1`$ is now replaced by a stronger result

$`|\mathrm{\Psi }_0(r)|{\displaystyle \frac{1}{\sqrt{2}}}(ba)^{3/2}\sqrt{E(a,0)}.`$ (67)
The repetition of the calculation of Sec. V gives finally (taking into account the above conditions)

$`{\displaystyle \frac{\delta E_a}{E(a,0)}}\left({\displaystyle \frac{2m}{a}}\right)^2\left({\displaystyle \frac{ba}{a}}\right)^4.`$ (68)
## VII Conclusions

In our paper we derived upper bounds for the backscattering of gravitational quadrupole waves propagating outward from a central compact object. The calculations were restricted to situations where the initial configuration was either an arbitrarily shaped wave with support outside some radius $`a`$ or the wave was a sharp pulse, i. e. its extension was small compared to its initial location $`a`$. The obtained upper bounds show that, for a given central object, the backscattering is the weaker the more outside from the central object the waves start propagating, and that is also the weaker the more compact the pulses are, i. e. the higher the involved frequencies are. The both results do confirm previous completely different calculations by Price, Pullin and Kundu . Backscattering should thus be strongest for pulses which start propagating outward close to the horizon of a black hole. This claim, however, needs further investigation for the following reasons. First, we gave bounds from above and not from below for the amount of backscattered energy. Second, the linear approximation may not be accurate enough very close to the horizon.

Results of Flanagan and Hughes  and Buonanno and Damour  have shown that the merger part of the gravitational wave signal could be a significant part of the total energy emitted. The wave pulse during the merger phase can be inside 3m. For a very compact pulse located in this region the inequality (60) of Sec. V can still yield a nontrivial bound, but in the general case our estimates fail. The main reason why we are loosing much in the accuracy is that we are forced to use several times - for the sake of generality - the Schwarz inequality. The present bounds can be significantly improved if initial data are explicitly known, since in this case they can be numerically bounded by an exact expression involving the initial energy and the Schwarz inequality would be used only once. On the other hand, it has been discovered that the backscattering can be quite strong when a signal propagates from within the photon sphere .

In a forthcoming paper we shall discuss, and compare with the results of our present paper, several aspects of the backscattering of gravitational waves where the sources of the gravitational waves are taken into account.

Acknowledgements. This work has been suported in part by the KBN grant 2 PO3B 010 16. One of the authors (EM) gratefully acknowledges financial support from the DAAD during his visit in Jena.

## VIII Appendix

In order to show the estimate (44) one begins with the second inequality of (43). Notice that $`H(t=0)=0`$, since $`\delta (R,t=0)=_t\delta (R,t)|_{t=0}`$. Thence the integration of (43) along $`\stackrel{~}{\mathrm{\Gamma }}_{a,(R,t)}`$ yields

$$H_M_a^{\mathrm{}}\frac{dr}{\eta _r}|RHS(r)|,$$
(69)
where $`RHS(r)`$ stands for the right hand side of Eq. (43). Our task consists in estimating the line integral of $`|RHS(r)|`$.

In order to do this one uses the estimates of (24). The calculation is quite long and we will describe only the main points. In the first step one uses the Schwarz inequality in the right hand side of (43), in order to obtain an expression of the type

$`4m\left({\displaystyle _R^{\mathrm{}}}𝑑r{\displaystyle \frac{(_t\delta )^2}{\eta _r}}\right)^{1/2}\times \left({\displaystyle _R^{\mathrm{}}}𝑑r\eta _r{\displaystyle \frac{f^2(r)}{r^8}}\right)^{1/2},`$ (70)
where $`f(r)`$ denotes $`(3\mathrm{\Psi }_1)^2`$, $`(2\mathrm{\Psi }_2/r)^2`$, or the squares of the terms that are proportional to $`63m`$. The first integral can be bounded by $`\sqrt{H(R)}`$, therefore (43) and (70) yield

$`(_t+_r^{})\sqrt{H(R,t)}2m\left({\displaystyle _R^{\mathrm{}}}𝑑r\eta _r{\displaystyle \frac{9\mathrm{\Psi }_1^2}{r^8}}\right)^{1/2}+2m\left({\displaystyle _R^{\mathrm{}}}𝑑r{\displaystyle \frac{4\mathrm{\Psi }_2^2}{r^{10}}}\right)^{1/2}+`$ (71)
$`2m\left({\displaystyle _R^{\mathrm{}}}𝑑r\left({\displaystyle \frac{63m(1+\frac{m}{r})}{4(1+\frac{3m}{2r})^2}}\right)^2{\displaystyle \frac{\mathrm{\Psi }_0^2}{r^8}}\right)^{1/2}+2m\left({\displaystyle _R^{\mathrm{}}}𝑑r\left({\displaystyle \frac{63m(1+\frac{m}{r})}{4(1+\frac{3m}{2r})^2}}\right)^2{\displaystyle \frac{\mathrm{\Psi }_1^2}{r^{10}}}\right)^{1/2}+`$ (72)
$`2m\left({\displaystyle _R^{\mathrm{}}}𝑑r\left({\displaystyle \frac{63m(1+\frac{m}{r})}{4(1+\frac{3m}{2r})^2}}\right)^2{\displaystyle \frac{\mathrm{\Psi }_2^2}{r^{12}}}\right)^{1/2}.`$ (73)
The integrands of (73) are taken at a time $`t`$ and $`(R,t)\stackrel{~}{\mathrm{\Gamma }}_a`$ ; the integration extends over the part $`rR`$ of the Cauchy hypersurface $`\mathrm{\Sigma }_t`$. At this place one inserts the bounds on $`\mathrm{\Psi }_0`$, $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$. That requires some care; the estimates hold true on the initial hypersurface $`\mathrm{\Sigma }_0`$, while here one needs estimates on $`\mathrm{\Sigma }_t`$. This point is clarified later. It is useful to introduce dimensionless variables $`x=R/a`$ and $`\stackrel{~}{m}=m/a`$.

i) First we shall consider the contribution that is due to $`3\mathrm{\Psi }_1`$. Let $`r_0`$ be defined by $`(r,t)\stackrel{~}{\mathrm{\Gamma }}_{r_0}`$. The insertion of the bound given in (24) bounds $`2m\left(_R^{\mathrm{}}𝑑r\eta _r\frac{9\mathrm{\Psi }_1^2}{r^8}\right)^{1/2}`$ by $`6mC_1\left(_R^{\mathrm{}}𝑑r\eta _r^4\frac{g_1(r_0)}{r^5}\right)^{1/2}`$. Notice that $`g_1(r)`$ is an increasing function, therefore if one replaces $`g_1(r_0)`$ by $`g_1(r)`$, then the integral that appears here cannot be smaller. In this way one utilizes the initial information (the energy inequality (24)) in order to control the evolution. The integral in question can be performed explicitly, with the result

$`6mC_1\left({\displaystyle _R^{\mathrm{}}}𝑑r\eta _r^4{\displaystyle \frac{g_1(r)}{r^5}}\right)^{1/2}={\displaystyle \frac{6mC_1}{a^2}}\sqrt{G_1(x)}.`$ (74)
Here, it holds

$`G_1(x)={\displaystyle \frac{\stackrel{~}{m}^4(137770\stackrel{~}{m}+1880\stackrel{~}{m}^22160\stackrel{~}{m}^3+960\stackrel{~}{m}^4)}{30(1+\stackrel{~}{2}m)^5x^8}}`$ (75)
$`{\displaystyle \frac{\stackrel{~}{m}^3(9915110m+10840\stackrel{~}{m}^28880\stackrel{~}{m}^3720\stackrel{~}{m}^4+3360\stackrel{~}{m}^5)}{105(1+2m)^5x^7}}+`$ (76)
$`{\displaystyle \frac{\stackrel{~}{m}^2(298113010\stackrel{~}{m}+18440\stackrel{~}{m}^2+7920\stackrel{~}{m}^341520\stackrel{~}{m}^4+27360\stackrel{~}{m}^5)}{420(1+2\stackrel{~}{m})^5x^6}}`$ (77)
$`{\displaystyle \frac{\stackrel{~}{m}(449711370\stackrel{~}{m}21720\stackrel{~}{m}^2+133040\stackrel{~}{m}^3200240\stackrel{~}{m}^4+101600\stackrel{~}{m}^5)}{2100(1+2\stackrel{~}{m})^5x^5}}+`$ (78)
$`{\displaystyle \frac{375+4650\stackrel{~}{m}35400\stackrel{~}{m}^2+93200\stackrel{~}{m}^3110000\stackrel{~}{m}^4+49376\stackrel{~}{m}^5}{3360(1+2\stackrel{~}{m})^5x^4}}`$ (79)
$`{\displaystyle \frac{11}{336\stackrel{~}{m}x^3}}{\displaystyle \frac{11}{448\stackrel{~}{m}^2x^2}}{\displaystyle \frac{11}{448\stackrel{~}{m}^3x}}+{\displaystyle \frac{11\mathrm{ln}(\frac{x}{2\stackrel{~}{m}+x})}{896\stackrel{~}{m}^4}}`$ (80)
$`(280\stackrel{~}{m}^4640\stackrel{~}{m}^3x+560\stackrel{~}{m}^2x^2224\stackrel{~}{m}x^3+35x^4){\displaystyle \frac{\mathrm{ln}(\frac{2\stackrel{~}{m}+x}{12\stackrel{~}{m}}}{140x^8}}).`$ (81)
This rather long expression is quite well approximated by $`G_1=(1+4\mathrm{ln}x)/(16x^4)`$ if $`m/a<<1`$. The integration of (74) along a null cone $`C_a`$ is done as follows. The integral $`_1^{\mathrm{}}\sqrt{G_1(x)}`$ is bounded from above,

$`{\displaystyle _1^{\mathrm{}}}𝑑x\eta _x^1\sqrt{G_1(x)}\left({\displaystyle _1^{\mathrm{}}}𝑑xx^2G_1(x)\right)^{1/2}\left({\displaystyle _1^{\mathrm{}}}𝑑xx^2\eta _x^2\right)^{1/2}.`$ (82)
Numerical integration yields

$`6C_1{\displaystyle \frac{m}{a}}\sqrt{{\displaystyle _1^{\mathrm{}}}{\displaystyle \frac{(1+4\mathrm{ln}x)}{16x^2}}}+O((m/a)^2)`$ (83)
$`8.24{\displaystyle \frac{m\sqrt{E(a,0)}}{a}}+O((m/a)^2).`$ (84)
One can check that the neglected terms can give a contribution comparable to the leading term only at distances smaller than $`6.6m`$.

ii) The calculation concerning the contribution of the $`\mathrm{\Psi }_2`$ function is similar. The leading (proportional to $`m^0`$) term is $`_R^{\mathrm{}}𝑑r\frac{\mathrm{\Psi }_2^2}{r^{10}}`$. $`|\mathrm{\Psi }_2|`$ is bounded in terms of $`g_2(x)`$. $`g_2(x)`$ is an increasing function, and a reasoning similar to what was made when discussing $`g_1(x)`$, leads to the conclusion that one can again use the initial energy inequality given by (24). One finds that $`_R^{\mathrm{}}𝑑r\frac{\mathrm{\Psi }_2^2}{r^{10}}`$ is bounded from above by

$`{\displaystyle \frac{1}{a^4}}G_2(x):={\displaystyle \frac{1}{a^4}}{\displaystyle _x^{\mathrm{}}}𝑑y{\displaystyle \frac{g_2(y)}{y^6}}(12\stackrel{~}{m}/y)^2=`$ (85)
$`{\displaystyle \frac{4\stackrel{~}{m}^2(3+44\stackrel{~}{m}120\stackrel{~}{m}^2+96\stackrel{~}{m}^3)}{21(1+2\stackrel{~}{m})^3x^7}}`$ (86)
$`{\displaystyle \frac{2\stackrel{~}{m}(21+236\stackrel{~}{m}408\stackrel{~}{m}^2192m^3+576\stackrel{~}{m}^4)}{63(1+2\stackrel{~}{m})^3x^6}}+`$ (87)
$`{\displaystyle \frac{21+88\stackrel{~}{m}+480\stackrel{~}{m}^21968\stackrel{~}{m}^3+1760\stackrel{~}{m}^4}{105(1+2\stackrel{~}{m})^3\stackrel{~}{x}^5}}`$ (88)
$`{\displaystyle \frac{34}{105x^4}}{\displaystyle \frac{31}{630\stackrel{~}{m}x^3}}{\displaystyle \frac{31}{840\stackrel{~}{m}^2x^2}}`$ (89)
$`{\displaystyle \frac{31}{840\stackrel{~}{m}^3x}}+{\displaystyle \frac{31\mathrm{ln}(\frac{x}{x2\stackrel{~}{m}})}{1680\stackrel{~}{m}^4}}`$ (90)
$`{\displaystyle \frac{8\stackrel{~}{m}(60\stackrel{~}{m}^270\stackrel{~}{m}x+21x^2)\mathrm{ln}(\frac{2\stackrel{~}{m}+x}{12m})}{105x^7}}.`$ (91)
In the limit of $`m0`$ the function $`G_2(x)`$ coincides with $`(4+5x)/(20x^5)`$. Similarly as before, in order to get a term bounding $`\sqrt{H_M}`$, one should integrate $`\sqrt{G_2(x)}/(12\stackrel{~}{m}/x)`$ along a null cone $`C_a`$. That gives 0.15, up to terms $`O(m)`$, after manipulations similar to those done earlier. The $`O(m)`$ correction becomes dominant when $`2m/a>0.3`$. After a reasoning similar to that applied above in the case of $`\mathrm{\Psi }_1`$ one finds that the total contribution due to the bound on the $`\mathrm{\Psi }_2`$ function is equal to

$$4C_2\sqrt{_1^{\mathrm{}}𝑑x\frac{(4+5x)}{20x^3}}\frac{m\sqrt{E(a,0)}}{a}=2.19\frac{m\sqrt{E(a,0)}}{a}.$$
(92)
In summary, one obtains

$$\sqrt{H_M}10.43\frac{m\sqrt{E(a,0)}}{a}+O(m^2).$$
(93)
In the above analysis, in (73), we neglected the terms proportional to $`63m`$. They give corrections of the order $`O(m^2)`$ to the right hand side of (93). We checked that their contribution is small in the region $`a>6.6m`$. Our final result (49) tells us that $`a>\sqrt{218}m15m`$ is valid for a nontrivial estimate. Therefore, all higher order terms in (93) can be safely neglected.