text:
# Series expansions for lattice Green functions

## 1 Introduction

Green functions, also known as two-point functions or propagators, are ubiquitous in physics. They appear in classical field theories, quantum mechanics and quantum field theories to name but a few areas. When a discrete space-time approach is used as a regulator, the continuum Green functions become lattice functions. In the problems of condensed matter physics both the continuum and lattice approach are natural. Random walks on a lattice provide another field of application for Green functions, where they are generating functions for the return probability of a random walker on a lattice . Given an isotropic random walk occurring in steps of one unit on a hypercubic lattice, one considers the probability of visiting a given site, or of returning to the starting site, after $`n`$ steps. This problem was first considered by Pólya  who showed that this probability tends to one only in one and two dimensions, for $`n`$ tending to infinity. This is the recurrence/transience transition which is an important result in the study of random walks.

The free propagator appears naturally in the perturbative expansion of field theories. A lattice provides a regulator for the infinities in such expansions, and the four-dimensional discrete Green functions appear naturally . The three-dimensional function appears in the study of effective three-dimensional theories . A lattice regularization has also been used in  to study non-relativistic quantum scattering in two and three space dimensions. For dimensions greater than one, the continuum free Green function is singular at the origin, while the lattice version is finite. This provides a regularization of the UV divergences. In yet another context, the propagator was used to find the resistance of a network of resistors . (See also  in relation to conformal invariance.)

Although the integral representation of the lattice Green function for the hypercubic lattice is well-known, the analytic calculation of this integral remains a challenge. A closed form appears to exist only in one dimension and for a special case in two dimensions. In higher dimensions there are several partial results . A rather recent numerical approach was used in  and developed in  for the massless free propagator. It was applied to the two-dimensional case in . Similar position space methods were considered earlier in , and later in .

In this paper I consider the free propagator as a function of the mass squared which is allowed to vary in the whole complex plane. The space coordinates are treated as parameters. Series expansions are first obtained for large values of the mass. In one and two dimensions these series are hypergeometric while in higher dimensions they are not. This large mass expansion follows from the multidimensional integral representation of the lattice Green function. However this representation does not allow one to find the series expansion around other values of the mass. To this end recurrence relations and corresponding differential equations in the mass squared are derived. This approach gives the singularities of the Green function and allows one to expand around any point. Non-trivial monodromies are found around the singular points.

The paper is organized as follow. Section 2 introduces the notation, and the well-known integral representation for the minimal anisotropic lattice Green function for the simple cubic lattice in $`d`$ dimensions. Some general properties are also discussed. The one-dimensional integral representation in terms of Bessel functions is given. A simple derivation for the location of the singularity in the complex mass plane is found. A general expansion of the lattice Green function in the inverse of the mass is obtained, and a recurrence relation between dimensions for the coefficients of this expansion is noted. In section 3 the one-dimensional case is studied in detail as a simple illustration of the methods used. This also serves as a guide to the features common to all dimensions. The same methods are applied to the two dimensional Green function and two novel expansions are found (section 4). The recurrence relation and the differential equation are then obtained for the three-dimensional case (section 5). Section 6 concludes with general remarks for four and higher dimensions. Appendix A contains a set of formulæ used in the text. In appendix B, the method of  is extended to arbitrary mass and anisotropies.

## 2 General results

Consider the hypercubic lattice in $`d`$ dimensions, with unit vectors $`\widehat{e}_j`$, $`j=1,\mathrm{},d`$. The discrete anisotropic Green equation is

$$Hf(\stackrel{}{x})\underset{j=1}{\overset{d}{}}\alpha _j\left[f(\stackrel{}{x}+\widehat{e}_j)+f(\stackrel{}{x}\widehat{e}_j)\right]=2\beta f(\stackrel{}{x})+\delta _{\stackrel{}{x},\stackrel{}{0}}$$
(1)
The integers $`x_j`$ label the lattice sites, and the anisotropies $`\alpha _j`$ are arbitrary but non-vanishing complex numbers. When the $`\alpha _j`$’s are omitted, they should be assumed to be all equal to 1. The action of the discrete anisotropic $`d`$-dimensional Laplace operator $`\mathrm{\Delta }`$ on $`f`$ is given by $`(H2_{j=1}^d\alpha _j)f(\stackrel{}{x})`$. Here the contribution of the latter term has been absorbed in $`\beta `$ which is, depending on the context, the mass squared, the energy eigenvalue or a formal expansion parameter for a generating function. This definition of $`\beta `$ gives a natural parity symmetry (6), and renders the sets of singularities symmetric with respect to the origin. The free massless scalar propagator of lattice gauge theories corresponds to the isotropic case, $`\alpha _j=1`$ for all $`j`$, and $`\beta =d`$. A completely isotropic solution of the isotropic equation satisfies

$$2df(\widehat{e}_j)=2\beta f(\stackrel{}{0})+1,j=1,\mathrm{},d$$
(2)
Adding to a solution of (1) any solution of the homogeneous equation,

$$\underset{j=1}{\overset{d}{}}\alpha _j\left[f(\stackrel{}{x}+\widehat{e}_j)+f(\stackrel{}{x}\widehat{e}_j)\right]=2\beta f(\stackrel{}{x}),$$
(3)
yields another Green function. A large class of such solutions can be written as

$$h(\stackrel{}{x}|\stackrel{}{\alpha },\beta )=_\pi ^{+\pi }\mathrm{}_\pi ^{+\pi }\frac{d^d\stackrel{}{q}}{(2\pi )^d}\mathrm{exp}(i\stackrel{}{q}\stackrel{}{x})\stackrel{~}{h}(\stackrel{}{q})\delta (\underset{j=1}{\overset{d}{}}\alpha _j\mathrm{cos}q_j\beta )$$
(4)
where the function $`\stackrel{~}{h}`$ is arbitrary but well-behaved.

One can define a decoupling point at $`\beta =0`$, for which equation (1) becomes an equation for two sub-lattices. A lattice where $`_{j=1}^dx_j`$ is even and one where $`_{j=1}^dx_j`$ is odd.

The minimal solution of (1) is given by

$$G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )=_\pi ^{+\pi }\mathrm{}_\pi ^{+\pi }\frac{d^d\stackrel{}{q}}{(2\pi )^d}\frac{\mathrm{exp}(i\stackrel{}{q}\stackrel{}{x})}{2_{j=1}^d\alpha _j\mathrm{cos}q_j2\beta \pm iϵ}$$
(5)
where $`ϵ0^+`$. The $`ϵ`$ prescription removes integration ambiguities at the poles.

When only $`d^{}`$ anisotropy parameters are equal to each other, the solution (5) is invariant under the $`2^dd^{}!`$ parity transformations and permutations of the $`x_j`$’s. For $`d^{}=d`$, the symmetry is that of the $`d`$-dimensional hypercubic group, and the functions $`G_\pm ^{(d)}`$ are completely symmetric in the absolute values of their arguments $`x_j`$. There is also a parity symmetry relating $`\beta `$ to $`\beta `$,

$$G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )=e^{i\pi \left(1+X\right)}G_{}^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )$$
(6)
where $`X_{j=1}^d|x_j|`$, and a complex conjugation symmetry

$$\left(G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )\right)^{}=G_{}^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha }^{},\beta ^{})$$
(7)
The latter symmetry shows that the two Green functions are complex conjugates of each other.

The exponentiation formula

$$\frac{n!}{A^{n+1}}=_0^{\mathrm{}}𝑑tt^n\mathrm{exp}(tA),\mathrm{Re}(A)>0,n=0,1,2,\mathrm{}$$
(8)
and the integral representations (74) for the Bessel functions yield a one-dimensional integral representation for (5):

$$G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )=\frac{(\pm i)^{1+X}}{2}_0^+\mathrm{}𝑑t\mathrm{exp}\left(\frac{tϵ}{2}it\beta \right)J_{|x_1|}(\alpha _1t)\mathrm{}J_{|x_d|}(\alpha _dt)$$
(9)
In this expression $`ϵ`$ can be set to zero. For $`\alpha _jR`$, the integrals converge in the domains $`\mathrm{Im}\beta 0`$, without a finite number of real points. Note that (9) satisfies (1) by virtue of properties (7881).

To find the singularities of $`G_\pm ^{(d)}`$ in $`\beta `$, consider first the initial integral representation (5). For $`|\beta |>_j|\alpha _j|`$, the integrand is a continuous function, without singularities, integrated over a compact domain. Thus the Green function has no singularities for these values of $`\beta `$. This includes the point at infinity. Now consider, for simplicity, the case without anisotropies ($`\alpha _j=1`$). Equation (5) implies that the possible singularities are real. When $`d=1`$, the integral $`(\text{9})`$ converges provided the oscillating cosine of the asymptotic expansion (76) is not “canceled” by $`\mathrm{exp}(it\beta )`$. For $`\beta =\pm 1`$, and only for these values, the integral has a diverging contribution of the form $`^{\mathrm{}}𝑑t/\sqrt{t}`$, which results in the branch points $`\beta =\pm 1`$. This is confirmed by the explicit expressions given in section 3. The same reasoning holds for $`d=2`$. The product of the two cosines yields one divergent contribution, $`𝑑t/t`$ for three values of $`\beta `$: $`0`$ and $`\pm 2`$. This is confirmed in section 4. For $`\beta =0`$, note that this approach also predicts a lack of divergence for points on the odd sub-lattice. One has $`\beta \mathrm{ln}\beta `$, which vanishes as $`\beta `$ tends to 0. For $`d3`$, there are enough powers of $`\sqrt{t}`$ to give a converging integral at all $`\beta `$. This corresponds to the recurrence/transience transition in the context of random walks. Let $`[n]`$ denote the integer part of $`n`$. Taking $`[(d1)/2]`$ $`\beta `$-derivatives of (9), and using the asymptotic expansion for the Bessel function, yields $`d+1`$ singularities: $`\beta _0=d,d+2,\mathrm{},d2,d`$. For $`d1`$ and odd, they are of the branch point type: $`(\beta \beta _0)^{p\frac{1}{2}}`$, where $`p`$ is a non-negative integer. For $`d2`$ and even, the singularities are logarithmic of the type: $`(\beta \beta _0)^p^{}\mathrm{ln}^q^{}(\beta \beta _0)`$, where $`p^{}`$ is a non-negative integer and $`q^{}`$ a positive integer. A slightly more complicated but essentially similar analysis applies for arbitrary anisotropies. The location of the singularities will depend explicitly on the $`\alpha _j`$’s. These results can be made rigorous through the use of a Tauberian theorem.

The $`ϵ`$ prescription becomes a convergence factor when one uses the series expansions (75) of the Bessel functions. This gives

$$G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )=\frac{1}{2\beta ^{1+X}}\underset{n=0}{\overset{\mathrm{}}{}}c_{X+2n}^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha })\beta ^{2n}$$
(10)
where

$$c_n^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha })_\pi ^{+\pi }\mathrm{}_\pi ^{+\pi }\frac{d^d\stackrel{}{q}}{(2\pi )^d}\mathrm{exp}(i\stackrel{}{q}\stackrel{}{x})\left(\underset{j=1}{\overset{d}{}}\alpha _j\mathrm{cos}q_j\right)^n$$
(11)
and

$`c_n^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha })`$ $`=`$ $`0,n=0,\mathrm{},X1`$ (12)
$`c_{X+2n}^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha })`$ $`=`$ $`\left(X+2n\right)!{\displaystyle \underset{j=1}{\overset{d}{}}}\left({\displaystyle \frac{\alpha _j}{2}}\right)^{|x_j|}{\displaystyle \underset{k_10}{}}\mathrm{}{\displaystyle \underset{k_d0}{}}{\displaystyle \underset{j=1}{\overset{d}{}}}\left({\displaystyle \frac{\alpha _j}{2}}\right)^{2k_j}`$
$`\times {\displaystyle \frac{\delta _{k_1+\mathrm{}+k_d,n}}{k_1!\mathrm{}k_d!(|x_1|+k_1)!\mathrm{}(|x_d|+k_d)!}},n0`$
In particular one finds

$$G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )\frac{X!_{j=1}^d\alpha _j^{|x_j|}}{|x_1|!\mathrm{}|x_d|!(2\beta )^{1+X}}\beta \mathrm{}$$
(14)
which shows the Green function to vanish faster than the simple estimate $`1/(2\beta )`$ obtained from the integral representation (5). A more compact form of the coefficients (2) can be obtained by the pairwise replacement of Bessel functions through identity (82). This is done for $`d=2,3,4`$ in the following sections.

The result (10) does not depend on the sign of the $`ϵ`$ prescription. This is easily understood by noticing that for $`|\beta |>_{j=1}^d|\alpha _j|`$ the denominator in (5) does not have poles and therefore $`ϵ`$ can be set to zero. The large $`\beta `$ expansion then yields (10) with the $`c_n`$ coefficients given by (11). One can also conclude that the expansion (102) converges at least for $`|\beta |>_{j=1}^d|\alpha _j|`$. Convergence at the generalized massless point ($`\beta =_{j=1}^d|\alpha _j|`$) depends on the dimensionality of the lattice. This is related to the recurrence/transience of the random walk. In one and two dimensions the series diverge at this point, despite the $`ϵ`$ prescription. In higher dimensions the series converge.

The random walk interpretation of the $`c_n`$’s is the following. For a random walker starting from the origin, let $`P_n(\stackrel{}{x})`$ be the probability of visiting the site $`\stackrel{}{x}`$ after $`n`$ unit steps on the $`d`$-dimensional hypercubic lattice. Take the anisotropies to be positive and such that $`_{j=1}^d\alpha _j=d`$. The probability of jumping from $`\stackrel{}{x}`$ to $`\stackrel{}{x}+\widehat{e}_j`$ or to $`\stackrel{}{x}\widehat{e}_j`$ is $`\frac{\alpha _j}{2d}`$. One then has: $`P_n(\stackrel{}{x})=\frac{1}{d^n}c_n^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha })`$. The vanishing of $`c_n`$ for $`n<X`$ is therefore natural since $`X`$ is the minimal number of steps required to reach point $`\stackrel{}{x}`$.

General relations between the coefficients $`c_n`$ for different dimensions can be simply found. Let $`d^{}`$ be any positive integer smaller than $`d`$. Expanding $`_{j=1}^d\alpha _j\mathrm{cos}q_j`$ using the binomial formula readily yields

$$c_n^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha })=\underset{k=0}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{n}{k}\right)c_k^{(d^{})}(x_1,\mathrm{},x_d^{}|\alpha _1,\mathrm{},\alpha _d^{})c_{nk}^{(dd^{})}(x_{d^{}+1},\mathrm{},x_d|\alpha _{d^{}+1},\mathrm{},\alpha _d)$$
(15)
Note that these relations are valid for arbitrary anisotropies. One can also get other equations by expanding with the multinomial formulæ.

## 3 The one-dimensional Green function

In one dimension it is possible to obtain closed form expressions, and corresponding series expansions. The anisotropy parameter $`\alpha _1`$ is set to one as it is an irrelevant overall factor.

The closed forms are obtained by integration in the complex $`q`$-plane. The integration contour is the rectangular path $`]\pi +i\mathrm{},\pi ][\pi ,\pi ][\pi ,\pi +i\mathrm{}[]\pi +i\mathrm{},\pi +i\mathrm{}[`$, for the upper half-plane, and its reflection about the real axis for the lower half-plane. The part at infinity gives a vanishing contribution, while the two side contributions cancel each other because $`xx_1`$ is integer. One finds

$$G_\pm ^{(1)}(x|\beta )=\pm \frac{e^{\pm ik|x|}}{2i\mathrm{sin}k},\beta =\mathrm{cos}k,k]0,\pi [$$
(16)
One can compare (16) to its continuum counter-part:

$$g_\pm ^{(1)}(x|k)=_{\mathrm{}}^{\mathrm{}}\frac{dq}{2\pi }\frac{\mathrm{exp}(iqx)}{k^2q^2\pm iϵ}=\pm \frac{e^{\pm ik|x|}}{2ik}$$
(17)
One also finds

$$G_+^{(1)}(x|\beta )=G_{}^{(1)}(x|\beta )=\{\begin{array}{ccc}\hfill +\frac{e^{ik|x|}}{2i\mathrm{sin}k},k_1[0,\pi ],k_2>0& & \\ & & \\ \hfill \frac{e^{ik|x|}}{2i\mathrm{sin}k},k_1]0,\pi [,k_2<0& & \end{array}$$
(18)
where $`\beta =\mathrm{cos}(k_1+ik_2)`$ and $`k=k_1+ik_2`$. Note that the positive exponential $`\frac{e^{+k_2|x|}}{2\mathrm{sinh}k_2}`$ ($`k_2>0`$), despite satisfying equation (1), is not obtained. This was to be expected from the integral representation (5), since for $`|\beta |>1`$ the integrand has no singular point and the integral must vanish as $`|x|\mathrm{}`$.

The general expression (102) reduces to

$`G_\pm ^{(1)}(x|\beta )`$ $`=`$ $`{\displaystyle \frac{1}{(2\beta )^{|x|+1}}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(|x|+2k)!}{k!(|x|+k)!}}(2\beta )^{2k}`$
$`=`$ $`{\displaystyle \frac{1}{(2\beta )^{|x|+1}}}_2F_1({\displaystyle \frac{|x|}{2}}+{\displaystyle \frac{1}{2}},{\displaystyle \frac{|x|}{2}}+1;|x|+1;\beta ^2)`$
This series converges uniformly for $`|\beta |>1`$ and simply for $`|\beta |=1,\beta \pm 1`$; it diverges at $`\beta =\pm 1`$ and for $`|\beta |<1`$. The equality of this hypergeometric series to the closed form (18) was otherwise known . The recurrence relation for the coefficients of the Green function, with $`c_nc_n^{(1)}(x|1)`$, reads

$$(n^2x^2)c_nn(n1)c_{n2}=0$$
(20)
The function $`G_\pm ^{(1)}`$ satisfies a second-order differential equation in $`\beta `$, of the hypergeometric type:

$$(\beta ^21)y^{\prime \prime }+3\beta y^{}+(1x^2)y=0$$
(21)
The indices at the three regular singular points $`\beta =\pm 1`$ and $`\mathrm{}`$ are $`(\frac{1}{2},0)`$ and $`(1|x|,1+|x|)`$, respectively. The two branch points $`\beta =\pm 1`$ imply the existence of monodromies in the complex-mass plane. Solution (3) corresponds to $`(\beta =\mathrm{};s=1+|x|)`$. It is then natural to investigate the properties of the solution corresponding to $`(\beta =\mathrm{};s=1|x|)`$:

$`y_1(\beta )`$ $`=`$ $`2^{|x|2}\beta _2^{|x|1}F_1({\displaystyle \frac{1|x|}{2}},1{\displaystyle \frac{|x|}{2}};1|x|;\beta ^2),x0`$ (22)
$`y_{l1}(\beta )`$ $`=`$ $`{\displaystyle \frac{1}{4\beta }}{\displaystyle \frac{1}{\sqrt{1\beta ^2}}}\mathrm{ln}\left({\displaystyle \frac{1}{\beta }}\right)+{\displaystyle \frac{1}{4}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\beta ^{(2n+1)}{\displaystyle \frac{d}{ds}}\left({\displaystyle \frac{\mathrm{\Gamma }\left(n+\frac{s}{2}\right)\mathrm{\Gamma }\left(\frac{s+1}{2}\right)}{\mathrm{\Gamma }\left(n+\frac{s+1}{2}\right)\mathrm{\Gamma }\left(\frac{s}{2}\right)}}\right)_{|s=1}`$ (23)
with $`{}_{2}{}^{}F_{1}^{}(0,\frac{1}{2};0,\beta ^2)1`$. The derivative of the coefficients in (23) can be written using the function $`\psi (z)=\frac{d}{dz}\mathrm{ln}\mathrm{\Gamma }(z)`$. The appearance of the logarithm for this second solution is due the degeneracy of the indices at $`x=0`$. The hypergeometric series (22) truncates to polynomials in $`\beta ^{+1}`$. Note that this solution almost provides a Green function. The a priori arbitrary normalization of a solution was chosen to be equal to $`2^{|x|2}`$ in (22) so that the one-dimensional Green equation, $`f(x+1|\beta )+f(x1|\beta )2\beta f(x|\beta )=\delta _{x,0}`$, is satisfied by $`f(x|\beta )`$ given in (22) and $`f(0|\beta )`$ set to 0.

The expansion around $`\beta =1`$ can be carried out similarly. The two solutions are

$`y_0(\beta )`$ $`=`$ $`{}_{2}{}^{}F_{1}^{}(1|x|,1+|x|;{\displaystyle \frac{3}{2}};{\displaystyle \frac{1\beta }{2}})`$ (24)
$`=`$ $`{\displaystyle \frac{\mathrm{sin}\left(2|x|\mathrm{Arcsin}\left(\sqrt{\frac{1\beta }{2}}\right)\right)}{|x|\sqrt{1\beta ^2}}}\mathrm{for}\beta ]1,1]`$ (25)
$`y_{1/2}(\beta )`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\beta 1}}}_2F_1({\displaystyle \frac{1}{2}}|x|,{\displaystyle \frac{1}{2}}+|x|;{\displaystyle \frac{1}{2}};{\displaystyle \frac{1\beta }{2}})`$ (26)
$`=`$ $`i{\displaystyle \frac{\sqrt{2}\mathrm{cos}\left(2|x|\mathrm{Arcsin}\left(\sqrt{\frac{1\beta }{2}}\right)\right)}{\sqrt{1\beta ^2}}}\mathrm{for}\beta ]1,1[`$ (27)
These solutions are valid for all values of $`x`$. In (27) the choice $`\sqrt{1}=i`$ was made. The solution $`y_0`$ is regular at $`\beta =1`$ and is a solution to the Green equation when appropriately normalized, while $`y_{1/2}`$ is singular at $`\beta =1`$, and is a solution to the homogeneous equation (3). Linear combinations of these solutions provide the two analytic continuations to the function defined by (3):

$$\frac{1}{2\sqrt{2}}(\sqrt{2}|x|y_0\pm y_{1/2})=\pm \frac{1}{2i\mathrm{sin}k}e^{\pm ik|x|},\beta =\mathrm{cos}k,k]0,\pi [$$
(28)
Moreover one finds

$`\beta _2^{|x|1}F_1({\displaystyle \frac{1|x|}{2}},1{\displaystyle \frac{|x|}{2}};1|x|;\beta ^2)`$
$`=2^{1|x|}|x|_2F_1(1|x|,1+|x|;{\displaystyle \frac{3}{2}};{\displaystyle \frac{1\beta }{2}}),|x|1`$ (29)
This corresponds to the polynomial solutions (22) and (24), which have to match since, as polynomials, they are defined on the whole complex $`\beta `$-plane. (An amusing by-product of the foregoing analysis is the identity: $`1=_{n=1}^{\mathrm{}}\frac{(2n+1)!}{2^{3n}(2n1)(n!)^2}`$) The expansion around $`\beta =1`$ yields similar results as can be expected from parity. This symmetry is not explicit on the series representations because the expansion point is not $`\beta =0`$.

The limits $`\beta \pm 1`$ for the Green function do not exist. However the following well-defined limit

$$G^{(1)}(x)\frac{1}{2}\underset{\beta 1^{}}{lim}\left(G_+^{(1)}(x|\beta )+G_{}^{(1)}(x|\beta )\right)=\frac{1}{2}|x|$$
(30)
is also a solution to the Green equation. In fact the latter equation is, for any $`\beta `$, a one-dimensional recurrence relation which can be solved directly by the standard method. The points $`\beta =\pm 1`$ are degeneracy points and the direct solution in this case is: $`f(x)=(\pm 1)^x(f(0)\pm \frac{1}{2}|x|)`$, for $`\beta =\pm 1`$, and $`f(0)`$ arbitrary. For $`\beta =1`$, and with the choice $`f(0)=0`$, one recovers (30).

## 4 The two-dimensional Green function

For arbitrary anisotropies, the identity (82) for the product of two Bessel functions can be used to obtain the following series expansion for the Green function:

$`G_\pm ^{(2)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )`$ $`=`$ $`{\displaystyle \frac{1}{|x_2|!\mathrm{\hspace{0.17em}2}\beta ^{X+1}}}\left({\displaystyle \frac{\alpha _1}{2}}\right)^{|x_1|}\left({\displaystyle \frac{\alpha _2}{2}}\right)^{|x_2|}`$
$`\times {\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(X+2k)!}{k!(|x_1|+k)!}}_2F_1(k,|x_1|k;|x_2|+1;\left({\displaystyle \frac{\alpha _2}{\alpha _1}}\right)^2)\left({\displaystyle \frac{\alpha _1}{2\beta }}\right)^{2k}`$
When $`\alpha _1=\alpha _2`$, these anisotropies can be set to 1 and the preceding expression becomes

$`G_\pm ^{(2)}(\stackrel{}{x}|\beta )`$ $`=`$ $`{\displaystyle \frac{1}{(2\beta )^{X+1}}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(X+k+1)_k(X+2k)!}{k!(|x_1|+k)!(|x_2|+k)!}}(2\beta )^{2k}`$ (32)
$`=`$ $`{\displaystyle \frac{X!}{(2\beta )^{X+1}|x_1|!|x_2|!}}`$
$`\times _4F_3({\displaystyle \frac{X+1}{2}},{\displaystyle \frac{X+1}{2}},{\displaystyle \frac{X}{2}}+1;{\displaystyle \frac{X}{2}}+1;X+1,|x_1|+1,|x_2|+1;4\beta ^2)`$
The series (32) converges uniformly for $`|\beta |>2`$, simply for $`|\beta |=2`$ and $`\beta \pm 2`$ and diverges at $`\beta =\pm 2`$ and $`|\beta |<2`$.

The recurrence relation for the coefficients of the Green function read

$$(n^2X^2)(n^2x^2)c_n4n^2(n1)^2c_{n2}=0$$
(34)
where $`x|x_1||x_2|`$. The Green function $`G_\pm ^{(2)}`$ satisfies a $`4^{\mathrm{th}}`$ order differential equation in $`\beta `$:

$`\beta ^2(\beta ^24)y^{\prime \prime \prime \prime }+\beta (10\beta ^216)y^{\prime \prime \prime }+\left(\beta ^2(2(x_1^2+x_2^2)+25)8\right)y^{\prime \prime }`$
$`+3\beta \left(2(x_1^2+x_2^2)+5\right)y^{}+(1X^2)(1x^2)y=0`$ (35)
The four regular singular points and their indices are:

$`\beta `$ $`=`$ $`0:s=1,1,0,0`$
$`\beta `$ $`=`$ $`\pm 2:s=2,1,0,0`$
$`\beta `$ $`=`$ $`\mathrm{}:s=1+X,\mathrm{\hspace{0.33em}1}X,\mathrm{\hspace{0.33em}1}+x,\mathrm{\hspace{0.33em}1}x`$
The three finite singular points were predicted in section 2. The series (32) is the regular solution corresponding to $`(\beta =\mathrm{};s=1+X)`$. The degeneracy of the indices at the other points signals the existence of logarithmic solutions.

The radius of convergence of the series expansion around $`\beta =0`$ is equal to 2, as $`\beta =\pm 2`$ are the closest singularities. Thus this expansion and the expansion around infinity cover the whole complex $`\beta `$-plane. The recurrence relation for the coefficients of the series expansions around $`\beta =0`$ contains two terms, just like (34) for the expansion around infinity. This permits the explicit determination of the series. One also has to find the connection coefficients involved in the linear combinations of the four solutions which reproduce the Green functions at hand, or the analytic continuations of (32). In the previous section, in (28), one had $`\frac{|x|}{2}`$ and $`\pm \frac{1}{2\sqrt{2}}`$. But the determination of the appropriate functions of $`\stackrel{}{x}`$ is not an easy task and there is no general systematic method which gives a closed form result. Here it is possible to carry out this analysis completely, and I have found the following expansions and connection coefficients:

$$G_\pm ^{(2)}(\stackrel{}{x}|\beta )=[f_0(\stackrel{}{x})\pm l_0(\stackrel{}{x})]y_0(\beta )+[f_1(\stackrel{}{x})\pm l_1(\stackrel{}{x})]y_1(\beta )\pm h_0(\stackrel{}{x})y_{l0}(\beta )\pm h_1(\stackrel{}{x})y_{l1}(\beta )$$
(36)
where

$`y_0(\beta )=_4F_3({\displaystyle \frac{1+X}{2}},{\displaystyle \frac{1X}{2}},{\displaystyle \frac{1+x}{2}};{\displaystyle \frac{1x}{2}};1,{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}};{\displaystyle \frac{\beta ^2}{4}})`$ (37)
$`y_1(\beta )=\beta _4F_3({\displaystyle \frac{2+X}{2}},{\displaystyle \frac{2X}{2}},{\displaystyle \frac{2+x}{2}};{\displaystyle \frac{2x}{2}};{\displaystyle \frac{3}{2}},{\displaystyle \frac{3}{2}},1;{\displaystyle \frac{\beta ^2}{4}})`$ (38)
$`y_{l0}(\beta )=y_0(\beta )\mathrm{ln}\beta +{\displaystyle \frac{d}{ds}}_5F_4({\displaystyle \frac{s+1+X}{2}},{\displaystyle \frac{s+1X}{2}},{\displaystyle \frac{s+1+x}{2}};{\displaystyle \frac{s+1x}{2}},1;`$ (39)
$`{\displaystyle \frac{s+2}{2}},{\displaystyle \frac{s+2}{2}},{\displaystyle \frac{s+1}{2}},{\displaystyle \frac{s+1}{2}};{\displaystyle \frac{\beta ^2}{4}})_{|s=0}`$
$`y_{l1}(\beta )=y_1(\beta )\mathrm{ln}\beta +\beta {\displaystyle \frac{d}{ds}}_5F_4({\displaystyle \frac{s+1+X}{2}},{\displaystyle \frac{s+1X}{2}},{\displaystyle \frac{s+1+x}{2}};{\displaystyle \frac{s+1x}{2}},1;`$ (40)
$`{\displaystyle \frac{s+2}{2}},{\displaystyle \frac{s+2}{2}},{\displaystyle \frac{s+1}{2}},{\displaystyle \frac{s+1}{2}};{\displaystyle \frac{\beta ^2}{4}})_{|s=1}`$
are the four independent solutions of the differential equation around $`\beta =0`$. The $`s`$-derivatives can be easily expressed in terms of $`\psi (z)`$, the logarithmic derivative of the Gamma function. The logarithms are taken real for positive $`\beta `$. The six connection functions are given by:

$`f_0(\stackrel{}{x})`$ $`=`$ $`{\displaystyle \frac{1}{4}}\mathrm{cos}\left({\displaystyle \frac{\pi }{2}}X\right)\mathrm{cos}\left({\displaystyle \frac{\pi }{2}}x\right)+{\displaystyle \frac{1}{4}}\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}X\right)\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}|x|\right)`$ (41)
$`f_1(\stackrel{}{x})`$ $`=`$ $`{\displaystyle \frac{1}{2}}|x_1^2x_2^2|f_0(\stackrel{}{x})={\displaystyle \frac{1}{2}}|x|Xf_0(\stackrel{}{x})`$ (42)
$`h_0(\stackrel{}{x})`$ $`=`$ $`h_0\mathrm{cos}\left({\displaystyle \frac{\pi }{2}}X\right)\mathrm{cos}\left({\displaystyle \frac{\pi }{2}}x\right)`$ (43)
$`h_1(\stackrel{}{x})`$ $`=`$ $`{\displaystyle \frac{h_0}{2}}|x_1^2x_2^2|\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}X\right)\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}|x|\right)={\displaystyle \frac{h_0}{2}}xX\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}X\right)\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}x\right)`$ (44)
$`l_0(\stackrel{}{x})`$ $`=`$ $`\left(\psi \left({\displaystyle \frac{X+1}{2}}\right)+\psi \left({\displaystyle \frac{|x|+1}{2}}\right)\right)h_0(\stackrel{}{x})`$ (45)
$`l_1(\stackrel{}{x})`$ $`=`$ $`\left(\psi \left({\displaystyle \frac{X}{2}}+1\right)+\psi \left({\displaystyle \frac{|x|}{2}}+1\right)2\right)h_1(\stackrel{}{x})`$ (46)
$`{\displaystyle \frac{h_0}{2}}(|x|+X)\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}X\right)\mathrm{sin}\left({\displaystyle \frac{\pi }{2}}|x|\right)`$
where $`h_0=\frac{i}{2\pi }`$ These functions satisfy a number of equations of which the simplest are:

$`Hf_0(\stackrel{}{x})=\delta _{\stackrel{}{x},\stackrel{}{0}},Hh_0(\stackrel{}{x})=0,Hl_0(\stackrel{}{x})=0`$ (47)
$`Hf_1(\stackrel{}{x})=2f_0(\stackrel{}{x}),Hh_1(\stackrel{}{x})=2h_0(\stackrel{}{x}),Hl_1(\stackrel{}{x})=2l_0(\stackrel{}{x})`$ (48)
$`x_1[f_0(x_1+1,x_2)f_0(x_11,x_2)]=x_2[f_0(x_1,x_2+1)f_0(x_1,x_21)]`$ (49)
Thus $`f_0`$ is a Green function at the decoupling point while $`h_0`$ and $`l_0`$ are solutions of the homogeneous equation at the same point. The functions $`f_1`$, $`h_1`$ and $`l_1`$ appear as potentials for the sources $`f_0`$, $`h_0`$ and $`l_0`$, respectively. The last equation for $`f_0`$ is a strong form of a directional-independence relation. The $`f_0f_1`$ part in (36) is a solution, regular at $`\beta =0`$, of the Green equation. The remaining part is a singular solution of the homogeneous equation, independently from the value of $`h_0`$. Finally, $`h_0(\stackrel{}{x})y_0(\beta )+h_1(\stackrel{}{x})y_1(\beta )`$ is a regular solution of the homogeneous equation.

The parity property is satisfied with $`\mathrm{ln}(1)=\pm i\pi `$, and the complex conjugation symmetry for both $`\beta `$ and $`\beta ^{}`$ not on the branch cut.

Consider now the expansion about $`\beta =2`$, and let $`v=\beta 2`$. The index $`s`$ can take one of the three values already found: $`0`$, $`1`$ or 2. A generic solution is given by

$$y(v)=\underset{n0}{}a_n(s)v^{n+s}$$
(50)
where

$`16(n+s)^2(n+s1)(n+s2)a_n`$ (51)
$`+4(n+s1)(n+s2)\left[5n^2+(10s9)n+5+5s^29s2(x_1^2+x_2^2)\right]a_{n1}`$
$`+2(2n+2s3)(n+s2)\left[2n^2+(4s6)n+5+2s^26s2(x_1^2+x_2^2)\right]a_{n2}`$
$`+(n+s2+X)(n+s2X)(n+s2+x)(n+s2x)a_{n3}=0`$
with $`n0`$ and $`a_{<0}0`$. For $`a_0(s)1`$, the solutions

$`y_2(v)={\displaystyle \underset{n0}{}}a_n(s=2)v^{n+2}`$ (52)
$`y_1(v)={\displaystyle \underset{n0}{}}a_n(s=1)v^{n+1},a_10`$ (53)
$`y_0(v)={\displaystyle \underset{n0}{}}a_n(s=0)v^n,a_1=a_20`$ (54)
are linearly independent and regular at $`\beta =2`$. The fourth solution has the expected logarithmic singularity:

$$y_{l0}(v)=y_0(v)\mathrm{ln}v+\underset{n0}{}v^n\frac{d}{ds}a_n(s)_{|s=0}$$
(55)
The complete Green function is a priori a linear combination of the four solutions. I have not looked for the corresponding connection coefficients. However one can conclude that the logarithmic solution corresponds to a solution of the homogeneous equation (3), and can therefore be dropped without altering the Green property. The remaining piece is the natural regularization at $`\beta =2`$ of the Green function defined by (5). The situation at $`\beta =2`$ is similar.

For $`x_1=x_2`$ and arbitrary anisotropies, Montroll found an expression in terms of Legendre functions . In terms of the definitions adopted here one has:

$$G_\mathrm{M}^{(2)}((x_1,x_1)|(\alpha _1,\alpha _2),\beta )=\frac{1}{2\pi (\alpha _1\alpha _2)^{\frac{1}{2}}}Q_{|x_1|\frac{1}{2}}\left(\frac{\beta ^2(\alpha _1^2+\alpha _2^2)}{2\alpha _1\alpha _2}\right)$$
(56)
where $`Q_\nu (z)`$ is the $`\nu ^{\mathrm{th}}`$ Legendre function of the second kind. These functions have logarithmic singularities at $`z=\pm 1`$. But for $`\nu `$ a half odd-integer, only $`z=1`$ is a singularity:

$$Q_{|x_1|\frac{1}{2}}(1ϵ)\frac{1}{2}\mathrm{ln}\left(\frac{ϵ}{2}\right)\gamma \psi \left(|x_1|+\frac{1}{2}\right)+O(ϵ),ϵ0^+$$
(57)
where $`\gamma =0.577216\mathrm{}`$. Thus one finds a logarithmic divergence at $`\beta =\pm 2`$:

$`G_\mathrm{M}^{(2)}((x_1,x_1)|(1,1),\beta )={\displaystyle \frac{1}{4\pi }}\mathrm{ln}\left({\displaystyle \frac{4\beta ^2}{4}}\right)`$ (58)
$`+{\displaystyle \frac{\gamma }{2\pi }}+{\displaystyle \frac{1}{2\pi }}\psi \left(|x_1|+{\displaystyle \frac{1}{2}}\right)+O\left({\displaystyle \frac{\beta ^24}{2}}\right),\beta 2^{}\mathrm{or}\beta 2^+`$
One also has $`G_\mathrm{M}^{(2)}((x_1,x_1)|(1,1),0)=\frac{1}{4}`$

Montroll’s result should be qualified. It is general in terms of anisotropies, but partial as it applies only to the line $`x_1=x_2`$ and $`|\beta |<2`$. Also, it is only a part of the full Green function, even under these restrictions. Consider $`x_1=x_2`$ in the Green function given by (36). At $`\beta =0`$ this expression diverges, as expected from the analysis of section 2, and thus contradicts (56). However taking the half-sum of the two Green functions one finds:

$$G^{(2)}(\stackrel{}{x}|\beta )=\frac{1}{2}\left(G_+^{(2)}(\stackrel{}{x}|\beta )+G_{}^{(2)}(\stackrel{}{x}|\beta )\right)=f_0(\stackrel{}{x})y_0(\beta )+f_1(\stackrel{}{x})y_1(\beta )$$
(59)
This solution of the Green equation is regular at $`\beta =0`$, and

$$G^{(2)}((x_1,x_1)|\beta )=\frac{1}{4}\mathrm{cos}(\pi x_1)_2F_1(\frac{1}{2}+|x_1|,\frac{1}{2}|x_1|;1;\frac{b^2}{4})$$
(60)
For $`x_1=0`$, one has the following identity:

$${}_{2}{}^{}F_{1}^{}(\frac{1}{2},\frac{1}{2};1;\frac{b^2}{4})=\frac{2}{\pi }K\left(\frac{\beta ^2}{4}\right)$$
(61)
where $`K`$ is the complete elliptic function of the first kind . One also has

$$Q_{|x_1|\frac{1}{2}}\left(\frac{\beta ^22}{2}\right)=\frac{\pi }{2}\mathrm{cos}(\pi x_1)_2F_1(\frac{1}{2}+|x_1|,\frac{1}{2}|x_1|;1;\frac{b^2}{4}),|\beta |<2$$
(62)
which shows that $`G^{(2)}((x_1,x_1)|\beta )`$ and $`G_\mathrm{M}^{(2)}((x_1,x_1)|(1,1),\beta )`$ are equal, up to a factor of $`1`$. The origin of this sign in (56) is unclear. Compare now (56) to the half-sum of the two functions (32) at $`x_1=x_2`$. The former function is even in $`\beta `$ while the latter is odd. Therefore they cannot be equal, and (56) holds only for $`|\beta |<2`$.

## 5 The three-dimensional Green function

Using the identity (82) for the product of two Bessel functions one obtains the following series expansion for the three-dimensional function:

$`G_\pm ^{(3)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )`$ $`=`$ $`{\displaystyle \frac{1}{|x_2|!\mathrm{\hspace{0.33em}2}\beta ^{X+1}}}\left({\displaystyle \frac{\alpha _1}{2}}\right)^{|x_1|}\left({\displaystyle \frac{\alpha _2}{2}}\right)^{|x_2|}\left({\displaystyle \frac{\alpha _3}{2}}\right)^{|x_3|}`$
$`\times {\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(X+2k)!\left({\displaystyle \frac{\alpha _3}{2\beta }}\right)^{2k}{\displaystyle \underset{p=0}{\overset{k}{}}}{\displaystyle \frac{1}{p!(kp)!(|x_1|+p)!(|x_3|+kp)!}}`$
$`\times \left({\displaystyle \frac{\alpha _1}{\alpha _3}}\right)_2^{2p}F_1(p,|x_1|p;|x_2|+1;\left({\displaystyle \frac{\alpha _2}{\alpha _1}}\right)^2)`$
When the anisotropies $`\alpha _j`$ are equal to 1 the preceding expression simplifies to

$`G_\pm ^{(3)}(\stackrel{}{x}|\beta )`$ $`=`$ $`{\displaystyle \frac{1}{(2\beta )^{X+1}}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(X+2k)!}{(2\beta )^{2k}}}`$
$`\times {\displaystyle \underset{p=0}{\overset{k}{}}}{\displaystyle \frac{(|x_1|+|x_2|+p+1)_p}{p!(kp)!(|x_1|+p)!(|x_2|+p)!(|x_3|+kp)!}}`$
The series (5) converges for $`|\beta |3`$. For $`\stackrel{}{x}=\stackrel{}{0}`$ one can write

$`G_\pm ^{(3)}(\stackrel{}{0}|\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2k)!}{4^k(k!)^2}}u_k\beta ^{2k}`$ (65)
$`u_k`$ $`=`$ $`{\displaystyle \underset{p=0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{p}}\right)^2\left({\displaystyle \genfrac{}{}{0pt}{}{2p}{p}}\right)`$ (66)
The value of this series at $`\beta =3`$ was calculated by Watson :

$$G_\pm ^{(3)}(\stackrel{}{0}|3)=\frac{2}{\pi ^2}\left[18+12\sqrt{2}10\sqrt{3}7\sqrt{6}\right]K^2\left((2\sqrt{3})^2(\sqrt{3}\sqrt{2})^2\right)$$
(67)
The numerical value of (67) is: $`0.2527\mathrm{}`$ (The definition of  is adopted for the function $`K`$.) At the singular point $`\beta =3`$, this series converges, rather slowly, to the known value (67). The sum of the first 1001 terms gives: $`0.2502\mathrm{}`$. It is amusing to note that the $`u_k`$ ($`k1`$) appear to be divisible by 3, the number of dimensions.

It is possible to derive a recurrence relation for the coefficients of the $`3`$-dimensional Green function. Define the following even homogeneous polynomials

$`\mathrm{\Sigma }_{222}=x_1^2x_2^2x_3^2`$
$`\mathrm{\Sigma }_{422}=x_1^2x_2^2x_3^2(x_1^2+x_2^2+x_3^2)`$
$`\mathrm{\Sigma }_{2i}=x_1^{2i}+x_2^{2i}+x_3^{2i},i=1,2,3,4`$
$`\mathrm{\Sigma }_{(2i)(2i)}=x_1^{2i}x_2^{2i}+x_2^{2i}x_3^{2i}+x_1^{2i}x_3^{2i},i=1,2`$
$`\mathrm{\Sigma }_{(2i)2}=x_1^{2i}x_2^2+x_1^2x_2^{2i}+x_2^{2i}x_3^2+x_2^2x_3^{2i}+x_1^{2i}x_3^2+x_1^2x_3^{2i},i=2,3`$
For $`c_nc_n^{(3)}(\stackrel{}{x}|(1,\mathrm{},1))`$, I have found

$`(n2)(n4)[\mathrm{\Sigma }_84\mathrm{\Sigma }_{62}+6\mathrm{\Sigma }_{44}+4\mathrm{\Sigma }_{422}`$ (68)
$`4n^2(\mathrm{\Sigma }_6\mathrm{\Sigma }_{42}+10\mathrm{\Sigma }_{222})+2n^4(3\mathrm{\Sigma }_4+2\mathrm{\Sigma }_{22})4n^6\mathrm{\Sigma }_2+n^8]c_n`$
$`+4n(n1)^2(n4)[\mathrm{\Sigma }_6+\mathrm{\Sigma }_{42}+6\mathrm{\Sigma }_{222}(n^24n2)\mathrm{\Sigma }_42(3n^24n+2)\mathrm{\Sigma }_{22}`$
$`+n(n2)(5n^26n+8)\mathrm{\Sigma }_2n(n2)(n^2+2)(3n^26n+4)]c_{n2}`$
$`+2n(n1)(n2)(n3)[(n^24n+12)\mathrm{\Sigma }_4+2(5n^220n+12)\mathrm{\Sigma }_{22}`$
$`2n(n4)(7n^224n+28)\mathrm{\Sigma }_2+n(n4)(15n^496n^3+268n^2384n+248)]c_{n4}`$
$`+4n^2(n1)(n2)(n3)(n4)(n5)\left[3(n3)\mathrm{\Sigma }_2(7n^357n^2+158n162)\right]c_{n6}`$
$`+9n^2(n1)(n2)^2(n3)(n4)(n5)(n6)(n7)c_{n8}=0`$
This translates into a $`10^{\mathrm{th}}`$ order differential equation for $`y=G_\pm ^{(3)}`$:

$`\beta ^2(\beta ^21)^3(\beta ^29)y^{(10)}+\beta (\beta ^21)^2(61\beta ^4418\beta ^2+45)y^{(9)}`$ (69)
$`(\beta ^21)\left[\beta ^6(4\mathrm{\Sigma }_21433)+\beta ^4(16\mathrm{\Sigma }_2+7511)+\beta ^2(12\mathrm{\Sigma }_22673)+27\right]y^{(8)}`$
$`4\beta \left[\beta ^6(42\mathrm{\Sigma }_24167)+\beta ^4(146\mathrm{\Sigma }_2+17284)+\beta ^2(122\mathrm{\Sigma }_211683)18\mathrm{\Sigma }_2+1140\right]y^{(7)}`$
$`+[\beta ^6(6\mathrm{\Sigma }_4+4\mathrm{\Sigma }_{22}2552\mathrm{\Sigma }_2+102963)+\beta ^4(4\mathrm{\Sigma }_424\mathrm{\Sigma }_{22}+5740\mathrm{\Sigma }_2261972)`$
$`+\beta ^2(2\mathrm{\Sigma }_4+20\mathrm{\Sigma }_{22}2444\mathrm{\Sigma }_2+90750)+48\mathrm{\Sigma }_21536]y^{(6)}`$
$`\beta [\beta ^4(162\mathrm{\Sigma }_4108\mathrm{\Sigma }_{22}+17640\mathrm{\Sigma }_2337617)`$
$`+\beta ^2(76\mathrm{\Sigma }_4+392\mathrm{\Sigma }_{22}23180\mathrm{\Sigma }_2+470364)+14\mathrm{\Sigma }_4140\mathrm{\Sigma }_{22}+3812\mathrm{\Sigma }_264194]y^{(5)}`$
$`2[\beta ^4(2\mathrm{\Sigma }_62\mathrm{\Sigma }_{42}+20\mathrm{\Sigma }_{222}714\mathrm{\Sigma }_4476\mathrm{\Sigma }_{22}+28742\mathrm{\Sigma }_2278744)`$
$`+\beta ^2(2\mathrm{\Sigma }_62\mathrm{\Sigma }_{42}12\mathrm{\Sigma }_{222}+202\mathrm{\Sigma }_4+892\mathrm{\Sigma }_{22}`$
$`18870\mathrm{\Sigma }_2^{}+179166)+17\mathrm{\Sigma }_474\mathrm{\Sigma }_{22}+590\mathrm{\Sigma }_25055]y^{(4)}`$
$`4\beta [\beta ^2(16\mathrm{\Sigma }_616\mathrm{\Sigma }_{42}+160\mathrm{\Sigma }_{222}1230\mathrm{\Sigma }_4820\mathrm{\Sigma }_{22}+20776\mathrm{\Sigma }_2103135)`$
$`+8\mathrm{\Sigma }_6^{}8\mathrm{\Sigma }_{42}48\mathrm{\Sigma }_{222}+150\mathrm{\Sigma }_4+612\mathrm{\Sigma }_{22}5190\mathrm{\Sigma }_2+23022]y^{(3)}`$
$`+[\beta ^2(\mathrm{\Sigma }_84\mathrm{\Sigma }_{62}+6\mathrm{\Sigma }_{44}+4\mathrm{\Sigma }_{422}276\mathrm{\Sigma }_6+276\mathrm{\Sigma }_{42}2760\mathrm{\Sigma }_{222}`$
$`+6246\mathrm{\Sigma }_4+4164\mathrm{\Sigma }_{22}44916\mathrm{\Sigma }_2+108681)`$
$`40\mathrm{\Sigma }_6+40\mathrm{\Sigma }_{42}+240\mathrm{\Sigma }_{222}120\mathrm{\Sigma }_4^{}720\mathrm{\Sigma }_{22}+2280\mathrm{\Sigma }_24680]y^{(2)}`$
$`+9\beta [\mathrm{\Sigma }_8^{}4\mathrm{\Sigma }_{62}+6\mathrm{\Sigma }_{44}+4\mathrm{\Sigma }_{422}36\mathrm{\Sigma }_6`$
$`+36\mathrm{\Sigma }_{42}360\mathrm{\Sigma }_{222}+246\mathrm{\Sigma }_4^{}+164\mathrm{\Sigma }_{22}676\mathrm{\Sigma }_2+681)]y^{(1)}`$
$`+15\left[\mathrm{\Sigma }_8^{}4\mathrm{\Sigma }_{62}+6\mathrm{\Sigma }_{44}+4\mathrm{\Sigma }_{422}4\mathrm{\Sigma }_6+4\mathrm{\Sigma }_{42}40\mathrm{\Sigma }_{222}+6\mathrm{\Sigma }_4+4\mathrm{\Sigma }_{22}4\mathrm{\Sigma }_2+1\right]y=0`$
This equation has six regular singular points, 0, $`\pm 1`$, $`\pm 3`$ and $`\mathrm{}`$. The corresponding indices are:

$`\beta `$ $`=`$ $`0:s=7,7,6,5,5,4,3,2,1,0`$
$`\beta `$ $`=`$ $`\pm 1:s=6,5,4,3,2,1,0,{\displaystyle \frac{5}{2}}\text{,}{\displaystyle \frac{3}{2}}\text{,}{\displaystyle \frac{1}{2}}`$
$`\beta `$ $`=`$ $`\pm 3:s=8,7,6,5,4,3,2,1,0,{\displaystyle \frac{1}{2}}`$
$`\beta `$ $`=`$ $`\mathrm{}:s=5,3,\mathrm{\hspace{0.33em}1}+|x_1|+|x_2|+|x_3|,\mathrm{\hspace{0.33em}1}|x_1||x_2||x_3|,`$
$`1+|x_1|+|x_2||x_3|,\mathrm{\hspace{0.33em}1}+|x_1||x_2|+|x_3|,\mathrm{\hspace{0.33em}1}|x_1|+|x_2|+|x_3|,`$
$`1+|x_1||x_2||x_3|,\mathrm{\hspace{0.33em}1}|x_1||x_2|+|x_3|,\mathrm{\hspace{0.33em}1}|x_1|+|x_2||x_3|`$
The Green function is regular at $`\beta =\mathrm{}`$ and corresponds to the index $`1+|x_1|+|x_2|+|x_3|`$. The appearance of indices differing by integer values can result in logarithms. However the 10-term recurrence for the series expansion around $`\beta =\pm 3`$ shows that all 10 solutions do not contain logarithms. The same conclusion holds at $`\beta =\pm 1`$, with an 8-term recurrence. Therefore $`\beta =\pm 1`$ and $`\beta =\pm 3`$ are branch point singularities which are free of logarithms. At $`\beta =0`$ the indices indicate that some solutions contain logarithms. For the foregoing Green function $`\beta =0`$ is a regular point, and one should consider a linear combination of the logarithm-free solutions.

The reason for the appearance of $`\beta =0`$ as a singularity of the differential equation is unclear. Perhaps the fact that this point is a fixed point of the parity symmetry, or its status as a decoupling point may be relevant here.

At its four singular points, the Green function does not diverge as all the indices are non-negative; only the solutions which correspond to the vanishing indices give non-vanishing contributions. However, around a given singular point, one can still consider a natural regularization by dropping all the solutions in the corresponding linear combination which are associated with non-integer indices. Such solutions combine into a solution of the homogeneous equation (3). Finally note that some results were obtained by Joyce , at $`\stackrel{}{x}=\stackrel{}{0}`$ and around $`\beta =3`$.

## 6 Concluding remarks

Using equation (82) twice gives the large-$`\beta `$ series expansion of the four-dimensional Green function:

$`G_\pm ^{(4)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )`$ $`=`$ $`{\displaystyle \frac{_{j=1}^4\left(\frac{\alpha _j}{2}\right)^{|x_j|}}{|x_2|!|x_4|!}}{\displaystyle \frac{1}{2\beta ^{X+1}}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(X+2k)!\left({\displaystyle \frac{\alpha _3}{2\beta }}\right)^{2k}`$
$`\times {\displaystyle \underset{p=0}{\overset{k}{}}}{\displaystyle \frac{1}{p!(kp)!(|x_1|+p)!(|x_3|+kp)!}}\left({\displaystyle \frac{\alpha _1}{\alpha _3}}\right)^{2p}`$
$`\times _2F_1(p,|x_1|p;|x_2|+1;\left({\displaystyle \frac{\alpha _2}{\alpha _1}}\right)^2)`$
$`\times _2F_1((kp),|x_3|(kp);|x_4|+1;\left({\displaystyle \frac{\alpha _4}{\alpha _3}}\right)^2)`$
When the anisotropies are all equal to 1 the preceding expression becomes

$`G_\pm ^{(4)}(\stackrel{}{x}|\beta )`$ $`=`$ $`{\displaystyle \frac{1}{(2\beta )^{(X+1)}}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(X+2k)!}{(2\beta )^{2k}}}`$
$`\times {\displaystyle \underset{p=0}{\overset{k}{}}}{\displaystyle \frac{(|x_1|+|x_2|+p+1)_p(|x_3|+|x_4|+kp+1)_{kp}}{p!(kp)!(|x_1|+p)!(|x_2|+p)!(|x_3|+kp)!(|x_4|+kp)!}}`$
The series (6) converges for $`|\beta |4`$. For $`\stackrel{}{x}=\stackrel{}{0}`$ one can write

$`G_\pm ^{(4)}(\stackrel{}{0}|\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2k)!}{4^k(k!)^2}}v_k\beta ^{2k}`$ (72)
$`v_k`$ $`=`$ $`{\displaystyle \underset{p=0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{p}}\right)^2\left({\displaystyle \genfrac{}{}{0pt}{}{2k2p}{kp}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{2p}{p}}\right)`$ (73)
For $`\beta =4`$, this series converges to the known value of $`0.1549\mathrm{}`$ . The sum of the first 31 terms gives: $`0.1541\mathrm{}`$. Similarly to the $`u_k`$ in the preceding section, the $`v_k`$ ($`k1`$) also appear to be divisible by the number of dimensions, here equal to 4.

A derivation of the recurrence relation for the coefficients $`c_n^{(d)}(\stackrel{}{x})`$ can be done as for the lower dimensions. The order of the recurrence appears to be larger than 7. From this recurrence a differential equation can be derived. The 5 singular points are of the logarithmic type. The logarithmic solutions, at one given singular point, can be dropped leaving a regular Green function.

These features were seen to be common to the lowest dimensions. They also hold for all the higher dimensions. The coefficients $`c_n^{(d)}(\stackrel{}{x})`$ satisfy recurrence relations for all dimensions. The general form of these relations is easily inferred from the results for the lower dimensions. The coefficients appearing in the recurrence relations are polynomials in $`n`$ and the $`x_j^2`$’s. From such relations one can then derive the differential equation as was done for the lowest dimensions. Another common feature is the possibility of dropping the singular part, around one given singularity. This part is a solution of the homogeneous equation. Finally, at $`\beta =d`$, subtracting $`G_\pm ^{(d)}(\stackrel{}{0}|d)`$ from (5) provides another regularization (see  for $`d=2`$).

Determining the explicit recurrence relation and the differential equation is however not a trivial task. It is also difficult to find the $`\stackrel{}{x}`$-dependent coefficients appearing in the linear combination of the solutions around a given singularity. These techniques were applied for the lower dimensions and new results were obtained. While the low dimensions studied in this paper seem to be at the limit of tractability of these methods, the knowledge obtained about the analytic structure of the lattice Green functions in all dimensions is an important step in their study.

Acknowledgements: I would like to thank Ivan Horváth for bringing to my attention reference  and for a critical reading of the manuscript. Special thanks to Tim Newman for bringing to my attention reference , for a critical reading of the manuscript, for his lightning derivation of $`\frac{i}{2\pi }`$ and for enjoyable coffee breaks. I thank P. Arnold, P. Fendley and H. Thacker for discussions, and S. Adhikari and G. Moore for communications.

Appendix A: Bessel functions and other formulæ

The cylindrical Bessel functions $`J_n(z)`$ have both integral representations

$$J_n(z)=\frac{i^n}{\pi }_0^\pi 𝑑q\mathrm{exp}(iz\mathrm{cos}q)\mathrm{cos}(nq),nZ$$
(74)
and series expansions

$$J_n(z)=\left(\frac{z}{2}\right)^n\underset{k=0}{\overset{\mathrm{}}{}}\frac{(1)^k}{k!\mathrm{\Gamma }(n+k+1)}\left(\frac{z}{2}\right)^{2k},$$
(75)
with an infinite radius of convergence. The asymptotic behavior at infinity is given by

$$J_n(z)\sqrt{\frac{2}{\pi z}}\left(\mathrm{cos}\chi \frac{4n^21}{8z}\mathrm{sin}\chi \right),|z|\mathrm{},|\mathrm{arg}(z)|<\pi $$
(76)
where $`\chi =z\frac{\pi }{2}n\frac{\pi }{4}`$ One also has

$$_0^{\mathrm{}}J_n(x)𝑑x=1,n0$$
(77)
The Bessel functions have the following properties

$`J_n(z)=(1)^nJ_n(z)nZ`$ (78)
$`J_n(z)=(1)^nJ_n(z)nZ`$ (79)
$`J_0(0)=+1,J_n(0)=0,n0`$ (80)
$`J_{n1}(z)J_{n+1}(z)=2J_n^{}(z)nZ`$ (81)
A particular formula for the product of two Bessel functions is

$$J_m(az)J_n(bz)=\frac{\left(\frac{az}{2}\right)^m\left(\frac{bz}{2}\right)^n}{\mathrm{\Gamma }(n+1)}\underset{k=0}{\overset{\mathrm{}}{}}\frac{(1)_2^kF_1(k,mk;n+1;\frac{b^2}{a^2})}{k!\mathrm{\Gamma }(m+k+1)}\left(\frac{az}{2}\right)^{2k}$$
(82)
(A typographical error in  has been corrected.) Note that $`{}_{2}{}^{}F_{1}^{}`$ is in fact a polynomial in $`b^2/a^2`$. When $`a=b=1`$ this formula simplifies to

$$J_m(z)J_n(z)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{(1)^k(m+n+k+1)_k}{k!\mathrm{\Gamma }(m+k+1)\mathrm{\Gamma }(n+k+1)}\left(\frac{z}{2}\right)^{m+n+2k}$$
(83)
The generalized hypergeometric series $`{}_{p}{}^{}F_{q}^{}`$ are defined by :

$${}_{p}{}^{}F_{q}^{}(a_1,\mathrm{},a_p;b_1,\mathrm{},b_q;z)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{(a_1)_k\mathrm{}(a_p)_k}{(b_1)_k\mathrm{}(b_q)_k}\frac{z^k}{k!}$$
(84)
where the Pochhammer symbol $`(a)_k`$ is defined by

$$(a)_0=1,(a)_k=\frac{\mathrm{\Gamma }(a+k)}{\mathrm{\Gamma }(a)}=a(a+1)\mathrm{}(a+k1)$$
(85)
and $`pq+1`$. Define a differential operator $`\delta =z\frac{d}{dz}`$. The function (84) satisfies the differential equation

$$\left[\delta (\delta +b_11)\mathrm{}(\delta +b_q1)z(\delta +a_1)\mathrm{}(\delta +a_p)\right]y(z)=0$$
(86)
For $`p=q+1`$ this equation is Fuchsian with three regular singular points at $`0`$, $`1`$ and $`\mathrm{}`$.

Appendix B: The WLW approach for arbitrary mass

Starting from an observation of C. Vohwinkel, Lüscher and Weisz have developed a powerful algorithmic method for the numerical calculation of the massless lattice Green function in four dimensions . Their method applies immediately to any dimension. Here I show that this method generalizes to arbitrary mass and anisotropies. A conserved quantity for the two-dimensional Green function at the massless and decoupling points is also derived.

An integration by parts of the left-hand side of (87) yields the right-hand side:

$`\alpha _j\left(G_\pm ^{(d)}(\stackrel{}{x}+\widehat{e}_j|\stackrel{}{\alpha },\beta )G_\pm ^{(d)}(\stackrel{}{x}\widehat{e}_j|\stackrel{}{\alpha },\beta )\right)=x_j(\stackrel{}{x}|\stackrel{}{\alpha }),j=1,\mathrm{},d`$ (87)
$`(\stackrel{}{x}|\stackrel{}{\alpha })={\displaystyle _\pi ^\pi }\mathrm{}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^d\stackrel{}{q}}{(2\pi )^d}}\mathrm{exp}(i\stackrel{}{q}.\stackrel{}{x})\mathrm{ln}(2\beta 2{\displaystyle \underset{j=1}{\overset{d}{}}}\alpha _j\mathrm{cos}q_jiϵ)`$ (88)
Equation (1) allows one to find another expression for $``$:

$$(\stackrel{}{x}|\stackrel{}{\alpha })=\frac{2}{_{j=1}^dx_j}\left(\underset{j=1}{\overset{d}{}}\alpha _jG_\pm ^{(d)}(\stackrel{}{x}\widehat{e}_j|\stackrel{}{\alpha },\beta )\beta G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )\right)$$
(89)
provided $`_{j=1}^dx_j0`$. This gives the value of $`G_\pm ^{(d)}(\stackrel{}{x}+\widehat{e}_j|\stackrel{}{\alpha },\beta )`$ in terms of $`G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )`$ and $`G_\pm ^{(d)}(\stackrel{}{x}\widehat{e}_k|\stackrel{}{\alpha },\beta )`$. The repeated use of these recurrence relations, coupled with the $`\pm x_j`$ invariance, shows that $`G_\pm ^{(d)}(\stackrel{}{x}|\stackrel{}{\alpha },\beta )`$ is a linear combination of the $`2^d`$ values corresponding to $`x_j=0,1`$. Note that all the vertices of the unit hypercube are needed when the anisotropies are arbitrary. These $`2^d`$ values can be calculated numerically, and the particular “$`\pm `$” branch obtained depending on the given value of $`\beta `$ in the complex plane. This generalizes the approach developed in .

One can look for additional conserved quantities as was done in . However this method depends rather strongly on the dimension. For the isotropic two-dimensional case, define

$$g_0(n)=G_\pm ^{(2)}((n,0)|\beta ),g_1(n)=G_\pm ^{(2)}((n,1)|\beta ),n0$$
(90)
The Green equation gives

$$g_0(n+1)+g_0(n1)+2g_1(n)2\beta g_0(n)=0,n1$$
(91)
and an equation inferred from the above approach is

$$g_1(n+1)=\frac{2n}{n+1}\left(\beta g_1(n)g_0(n)\right)\frac{n1}{n+1}g_1(n1),n1$$
(92)
One can look for a conserved quantity in the following form

$`C(n)`$ $`=`$ $`ng_0(n)+a_1ng_1(n)+b_0(n1)g_0(n1)+b_1(n1)g_1(n1)`$ (93)
$`+c_0g_0(n)+c_1g_1(n)+d_0g_0(n1)+d_1g_1(n1),n1`$
Using (91) and (92) one finds that $`C(n)`$ is independent of $`n`$ provided

$$b_0=c_0=d_0=1,a_1=b_1=\frac{1}{\beta 1},c_1=d_1=0$$
(94)
and $`\beta =0`$ or $`\beta =2`$. This form does not allow for other values of $`\beta `$, but a conserved quantity at $`\beta =2`$ can be obtained from the one for $`\beta =2`$ through the parity symmetry. The new quantity is a priori conserved for $`n`$ odd and $`n`$ even separately. It would interesting to find out whether arbitrary values of $`\beta `$ accommodate conserved quantities.

Note that $`\beta =0,\pm 2`$ are the three singular values. Therefore $`C(n)`$ may not be well-defined. However, using the explicit expression of section 4, and taking the limit $`\beta 0`$, one finds that the infinities cancel exactly, leaving $`C_\pm C_\pm (n)=\pm \frac{i}{\pi }`$ for all $`n1`$. For the half-sum the conserved quantity is therefore $`C=0`$. The situation at $`\beta =2`$ is similar. The conserved quantities can be finite through cancellations, and the divergence corresponds to a solution of the homogeneous equation and can therefore be dropped. (See also the remark in the conclusion of  concerning this conserved quantity, and .)

Note added: Lattice Green functions also arise in the study of the statistical mechanics of the spherical model . Complex temperature singularities of this system were studied in . Lattice Green functions were also examined for the cases where factorizations in two complete elliptic integrals occur . I would like to thank P. Butera for bringing to my attention the four works cited in this note.


format: text

text:
# On the Structure of Advective Accretion Disks At High Luminosity

## 1 Introduction

Accretion discs are formed when the matter with a large angular momentum is falling into a black hole or another gravitating body. The well known objects where the accretion disks are found are protostar nebulae, binary X-ray sources, cataclysmic variables, active galactic nuclei and others. In this paper we discuss accretion disks around black holes. The standard model of geometrically thin accretion disk has been developed by Shakura (1972), Novikov & Thorne (1973) and Shakura & Sunyaev (1973), and has played a significant role in the development of accretion theory (see Pringle 1981; Frank, King, & Rain 1992; Kato, Mineshige, & Fukue 1998 for reviews). The standard model bases on the vertically averaged approach to equilibrium, and a suggestion of the local thermal balance in which the viscous heating of the gas is balanced by the local radiative cooling. Non-local effects, like the radial advection of thermal energy and the transonic nature of accretion flow, are neglected in the standard model. This simplified approach allows to reduce the general problem to a set of algebraical equations. Such a simple description becomes possible due to an approximate parameterization of the viscosity stress tensor with one non-zero component,

$$t_{r\varphi }=\alpha P,$$
$`(1)`$
suggested by Shakura (1972). The standard model gives a satisfactory appropriate solution of the problem at low accretion rates $`\dot{M}16L_{Edd}/c^2`$, where $`L_{Edd}`$ is the Eddington luminosity.

Simplified solution with inclusion of the advective terms into equations described the vertically integrated models of accretion disks was obtained by Paczyński & Bisnovatyi-Kogan (1981). This approach with some modifications have been used by many researchers to study transonic accretion flows around black holes (Muchotrzeb & Paczyński 1982; Muchotrzeb 1983; Matsumoto et al. 1984; Fukue 1987; Abramowicz et al. 1988; Chen & Taam 1993; Beloborodov 1998). The importance of the transonic nature of the accretion flows on the disk structure has been emphasized by Höshi & Shibazaki (1977), Liang & Thompson (1980) and Abramowicz & Zurek (1981), and later studied in more details by Abramowicz & Kato (1989).

Despite a significant progress in the study of optically thick accretion disks obtained during almost three decades there are a number of unsolved problems still posed in the theory. The problems are connected with a possible non-uniqueness of a solution at $`\alpha 0.01`$ and a non-standard behavior of a singular point type. It was reported by Matsumoto et al. (1984), Muchotrzeb-Czerny (1986) and Abramowicz et al. (1988) that in the case of viscosity prescription (1) the singular point changes its type from a saddle to node when one increases $`\alpha `$. The presence of the nodal-type singular point leads to creating of a possibility of multiple solutions as the authors have claimed. A similar change of the singular point type was reported by Chen & Taam (1993), who used the angular velocity gradient-dependent viscous stress,

$$t_{r\varphi }=\rho \nu r\frac{d\mathrm{\Omega }}{dr},$$
$`(2)`$
where $`\nu `$ is the kinematic viscosity coefficient defined by (14). Narayan, Kato, & Honma (1997) have compared two forms of viscosity (1) and (2) in the case of radiatively inefficient advection-dominated accretion flows. They concluded that the structures of flows corresponded to both viscosities are similar at $`\alpha <0.15`$.

In this paper we show that the mentioned problems have been created by several inconsistencies in the preceding studies. Some problems are connected with an inaccurate averaging of the equations over a disk thickness (Chen & Taam 1993), another ones appear due to an incomplete investigation of the singular points (Abramowicz et al. 1988). We have found that in the case of viscosity prescription (1) a set of equations describing the vertically averaged advective accretion disks has two singular points, independed of $`\alpha `$ and accretion rate. Note, that multiplicity of singular points in solutions for accretion flows in Paczyński-Wiita potential (3) was revealed by Fukue (1987), Chakrabarti & Molteni (1993) and Chakrabarti (1996) in a somewhat different context. We have shown, that at $`\alpha 0.01`$ the inner and outer (with respect to the black hole location) singular points are of the saddle type, and only one integral curve (“separatrix”) which crosses the inner point simultaneously crosses the outer one. This separatrix corresponds to the unique global solution which is determined by two parameters, $`\alpha `$ and $`\dot{m}=\dot{M}c^2/L_{Edd}`$, for a given black hole mass. In Figure 1a the structure of integral curves is schematically represented in the vicinity of the global solution which is shown by the thick line. At larger $`\alpha 0.1`$ the inner singular point is changed to a nodal-type one, while the outer point remains of a saddle-type. There is still one integral curve which goes continuously through both singular points providing a unique global solution, as it is shown in Figure 1b.

In the case of viscosity prescription (2) we have found that there is only one singular point which is always a saddle, and only one physical solution which passes through this point exists. Solutions which correspond to both forms of viscosity (1) and (2) are very close at low $`\alpha `$ limit, $`\alpha 0.1`$.

We have developed a numerical method to solve the set of equations describing the vertically averaged advective accretion disks. The method is based on the standard relaxation technique and explicitly uses conditions at the inner singular point and its vicinity. We have obtained these conditions by expanding a solution into power series around the singular point. Such a modification of the method allows to construct solutions which smoothly pass the singular points and satisfy the regularity conditions at these points with high computer precision in wide range of parameters $`\alpha `$ and $`\dot{m}`$.

The paper is organized as follows. In §2 we formulate a mathematical approach to the problem, write a set of equations, and formulate boundary conditions. In §3 we investigate critical points and discuss behavior of physical values in their vicinity. In §4 we describe our numerical results and discuss them in §5. Details of the numerical method and explicit expansion of physical quantities in the vicinity of the critical points are represented in Appendixes A and B, respectively.

## 2 Problem formulation

We consider a steady geometrically thin accretion disk around a non-rotating black hole. For simplicity, we use the pseudo-Newtonian approach to describe the disk structure in the vicinity of a black hole. In the approach the general relativistic effects are simulated by using Paczyński-Wiita potential (Paczyński & Wiita 1980)

$$\mathrm{\Phi }(r)=\frac{GM}{r2r_g},$$
$`(3)`$
where $`M`$ is the black hole mass and $`2r_g=2GM/c^2`$ is the gravitational radius. The disk self-gravity is neglected.

A general problem of investigation of two-dimensional structure of the accretion disks (in the radial and vertical directions) can be reduced to a one-dimension problem by averaging the disk structure in the vertical direction. In this formulation equations which are described the radial disk structure are written for the midplane density $`\rho `$, pressure $`P`$, radial velocity $`v`$ and angular velocity $`\mathrm{\Omega }`$. The mass conservation equation takes the form,

$$\dot{M}=4\pi rh\rho v,$$
$`(4)`$
where $`\dot{M}`$ is the accretion rate, $`\dot{M}>0`$, and $`h`$ is the disk half-thickness, which is expressed in terms of the isothermal sound speed $`c_s=\sqrt{P/\rho }`$ of gas,

$$h=\frac{c_s}{\mathrm{\Omega }_K}.$$
$`(5)`$
The equations of motion in the radial and azimuthal directions are

$$v\frac{dv}{dr}=\frac{1}{\rho }\frac{dP}{dr}+(\mathrm{\Omega }^2\mathrm{\Omega }_K^2)r,$$
$`(6)`$
$$\frac{\dot{M}}{4\pi }\frac{d\mathrm{}}{dr}+\frac{d}{dr}(r^2ht_{r\varphi })=0,$$
$`(7)`$
where $`\mathrm{\Omega }_K`$ is the Keplerian angular velocity, $`\mathrm{\Omega }_K^2=GM/r(r2r_g)^2`$, $`\mathrm{}=\mathrm{\Omega }r^2`$ is the specific angular momentum and $`t_{r\varphi }`$ is the ($`r`$, $`\varphi `$)-component of the viscous stress tensor. Other components of the stress tensor are assumed to be negligiblly small.

The vertically averaged energy conservation equation can be written in the form111The vertical averaging in equation (8) have been done in different way by different authors (compare e.g. Shakura & Sunyaev 1973, and Abramowicz et al. 1988). Our choice of coefficients in (9)-(11), following Chen & Taam (1993) may be not the optimal one. Aposteriory analysis had shown that using 4 in the denominator of (9) instead of 2, would be a more consistent choice, but this change has a little influence on our numerical results.,

$$F_{adv}=F^+F^{},$$
$`(8)`$
where

$$F_{adv}=\frac{\dot{M}}{2\pi r}\left[\frac{dE}{dr}+P\frac{d}{dr}\left(\frac{1}{\rho }\right)\right],$$
$`(9)`$
$$F^+=ht_{r\varphi }r\frac{d\mathrm{\Omega }}{dr},$$
$`(10)`$
$$F^{}=\frac{2aT^4c}{3\kappa \rho h},$$
$`(11)`$
are the advective energy flux, the viscous dissipation rate and the cooling rate per unit surface, respectively, $`T`$ is the midplane temperature, $`\kappa `$ is the opacity and $`a`$ is the radiation constant.

The equation of state for accretion matter consisted of a gas-radiation mixture is

$$c_s^2=T+\frac{1}{3}\frac{aT^4}{\rho },$$
$`(12)`$
where $``$ is the gas constant. The specific energy of the mixture is

$$E=\frac{3}{2}T+\frac{aT^4}{\rho }.$$
$`(13)`$
We consider two prescriptions of viscosity in our models. In one case we adopt a simple relation (1) between the viscous stress and pressure. In another case we assume the angular velocity gradient-dependent viscous stress (2), where the viscosity $`\nu `$ is taken in the form

$$\nu =\frac{2}{3}\alpha c_sh.$$
$`(14)`$
In the limit $`\mathrm{\Omega }\mathrm{\Omega }_K`$ both prescriptions (1) and (2) coincide.

Integrating equation (7) we obtain

$$r^2ht_{r\varphi }=\frac{\dot{M}}{4\pi }(\mathrm{}\mathrm{}_{in}),$$
$`(15)`$
where $`\mathrm{}_{in}`$ is an integration constant and has a meaning of the specific angular momentum of accreting matter at the black hole horizon. Depending on the used viscosity prescription (1) or (2) expression (15) results in an algebraical equation or a first order differential equation, respectively. So, in the case of viscosity prescription (1) we have only two first order differential equations (6) and (8), which require to set two parameters as boundary conditions. In the case of viscosity prescription (2) we obtain additional differential equation from (15) and have to set three parameters as boundary conditions. The integration constant $`\mathrm{}_{in}`$ is chosen to obtain a global transonic solution with a subsonic part at large radii and a supersonic part close to the black hole horizon.

## 3 Investigation of singular points

### 3.1 $`\alpha P`$ viscosity prescription

We consider first the case of viscosity prescription (1). From (15) we obtain the algebraical expression for $`\mathrm{\Omega }`$,

$$\mathrm{\Omega }=\frac{\mathrm{}_{in}}{r^2}+\alpha \frac{c_s^2}{vr}.$$
$`(16)`$
Using (16) the system of equations (6) and (8) can be reduced to the following form,

$$r\frac{v^{}}{v}=\frac{N_1}{D_1},$$
$`(17)`$
$$r\frac{c_s^{}}{c_s}=(1^2)\frac{N_1}{D_1}+1r\frac{\mathrm{\Omega }_K^{}}{\mathrm{\Omega }_K}+\frac{\mathrm{\Omega }^2\mathrm{\Omega }_K^2}{c_s^2}r^2,$$
$`(18)`$
where

$$N_1=\left(1r\frac{\mathrm{\Omega }_K^{}}{\mathrm{\Omega }_K}+\frac{\mathrm{\Omega }^2\mathrm{\Omega }_K^2}{c_s^2}r^2\right)\left(7\frac{3}{2}\beta \frac{1+\beta }{43\beta }\frac{\alpha ^2}{^2}\right)+\left(1r\frac{\mathrm{\Omega }_K^{}}{\mathrm{\Omega }_K}\right)\left(1+\frac{3}{2}\beta \frac{1\beta }{43\beta }\right)+$$

$$\alpha \frac{\mathrm{}_{in}}{vr}+\frac{1}{2}\frac{\alpha ^2}{^2}\frac{1\beta }{\dot{m}}\frac{\mathrm{\Omega }_Kr^2}{c_sr_g},$$
$`(19)`$
$$D_1=(^21)\left(7\frac{3}{2}\beta \frac{1+\beta }{43\beta }\frac{\alpha ^2}{^2}\right)\left(1+\frac{3}{2}\beta \frac{1\beta }{43\beta }+\frac{1}{2}\frac{\alpha ^2}{^2}\right).$$
$`(20)`$
In equations (17)-(20) we use the following notations: $`v^{}dv/dr`$, $`c_s^{}dc_s/dr`$, $`\mathrm{\Omega }_K^{}d\mathrm{\Omega }_K/dr`$, $`\beta =T/c_s^2`$ and $`=v/c_s`$. From (4) and (16) it follows the algebraical equation for $`\beta `$,

$$\beta ^4(1\beta )\frac{3}{4\pi }\frac{\dot{M}\mathrm{\Omega }_K^4}{arvc_s^7}=0.$$
$`(21)`$
The equation $`D_1=0`$ defines singular points of the differential equations (17) and (18), and can be reduced to the following form,

$$\left(7\frac{3}{2}\beta \frac{1+\beta }{43\beta }\right)^4\left(\alpha ^2+8\frac{3\beta ^2}{43\beta }\right)^2+\frac{\alpha ^2}{2}=0.$$
$`(22)`$
Equation (22) is a quadratic equation with respect to $`^2`$ and has two positive roots, which correspond to two singular points:

$$_{1,2}^2=\frac{1}{2}\left[\alpha ^2+8\frac{3\beta _s^2}{43\beta _s}\pm \sqrt{\left(\alpha ^2+8\frac{3\beta _s^2}{43\beta _s}\right)^22\alpha ^2\left(7\frac{3}{2}\beta _s\frac{1+\beta _s}{43\beta _s}\right)}\right]\left(7\frac{3}{2}\beta _s\frac{1+\beta _s}{43\beta _s}\right)^1,$$
$`(23)`$
where $`\beta _s`$ is the value of $`\beta `$ taken at the singular point.

In $`\alpha 1`$ limit we have:

$$_1^2=\left(8\frac{3\beta _s^2}{43\beta _s}\right)\left(7\frac{3}{2}\beta _s\frac{1+\beta _s}{43\beta _s}\right)^1$$
$`(24)`$
and

$$_2^2=\frac{\alpha ^2}{2}\left(8\frac{3\beta _s^2}{43\beta _s}\right)^2.$$
$`(25)`$
The first singular point, in which $`_s=_1`$, locates close to the black hole last stable orbit at $`r=6r_g`$. The corresponding values of $`_1`$ are $`1.118`$ and $`1.069`$ for the gas pressure supported ($`\beta =1`$) and the radiation pressure supported ($`\beta =0`$) accretion flows, respectively. This point is an analogy of the singular point in a spherical flow, where the point divides the subsonic and supersonic regions of accretion flow. The second singular point, in which $`_s=_2`$, located at larger radius, is the result of simplified viscosity prescription (1). We will use the notations $`(r_s)_{in}`$ and $`(r_s)_{out}`$ for positions of the inner and outer singular points, respectively.

At the singular points the numerator $`N_1`$ and denominator $`D_1`$ must simultaneously vanish to provide a regular behavior for a global solution. The type of the singular points must be consistent with a transonic nature of solution. For example, a spiral-type singular point does not satisfy the latter requirement, but a saddle or nodal-type point does it. Detailed analysis of topology in vicinity of singular points was done for thin accretion disks under isothermal approximation by Abramowicz & Kato (1989). They showed that saddle, nodal or spiral types are formally possible, but only saddle and nodal points are physically relevant. This study had confirmed the previously obtained numerical results by Matsumoto et al. (1984). The latter authors demonstrated in a framework of the isothermal accretion disks that the type of the inner singular point is defined by value of $`\mathrm{}_{in}`$. At larger $`\mathrm{}_{in}`$ the point is a spiral, at smaller $`\mathrm{}_{in}`$ there is no inner singular point at all, and only unique choice of $`\mathrm{}_{in}`$ of moderate values corresponds to a saddle or nodal-type singular point. The choice between saddle or nodal-type singular points can be done only by constructing a global model of the disk.

### 3.2 $`\mathrm{\Omega }`$-gradient-dependent viscous stress

In the case of viscosity prescription (2) the differential equations (6), (8) and (15) can be reduced to the following form,

$$r\frac{\mathrm{\Omega }^{}}{\mathrm{\Omega }}=\frac{3}{2}\frac{\mathrm{\Omega }_Krv}{\alpha c_s^2}\left(1\frac{\mathrm{}_{in}}{\mathrm{\Omega }r^2}\right),$$
$`(26)`$
$$r\frac{v^{}}{v}=\frac{N_2}{D_2},$$
$`(27)`$
$$r\frac{c_s^{}}{c_s}=(1^2)\frac{N_2}{D_2}+1r\frac{\mathrm{\Omega }_K^{}}{\mathrm{\Omega }_K}+\frac{\mathrm{\Omega }^2\mathrm{\Omega }_K^2}{c_s^2}r^2,$$
$`(28)`$
where $`\mathrm{\Omega }^{}d\mathrm{\Omega }/dx`$, $`\beta `$ is defined by (21) and

$$N_2=\left(1r\frac{\mathrm{\Omega }_K^{}}{\mathrm{\Omega }_K}+\frac{\mathrm{\Omega }^2\mathrm{\Omega }_K^2}{c_s^2}r^2\right)\left(7\frac{3}{2}\beta \frac{1+\beta }{43\beta }\right)+\frac{3}{4}\frac{\mathrm{\Omega }^2\mathrm{\Omega }_Kr^3v}{\alpha c_s^4}\left(1\frac{\mathrm{}_{in}}{\mathrm{\Omega }r^2}\right)^2+$$

$$\left(1r\frac{\mathrm{\Omega }_K^{}}{\mathrm{\Omega }_K}\right)\left(1+\frac{3}{2}\beta \frac{1\beta }{43\beta }\right)\frac{1\beta }{\dot{m}}\frac{\mathrm{\Omega }_Kr^2}{c_sr_g},$$
$`(29)`$
$$D_2=(^21)\left(7\frac{3}{2}\beta \frac{1+\beta }{43\beta }\right)\left(1+\frac{3}{2}\beta \frac{1\beta }{43\beta }\right).$$
$`(30)`$
There is only one singular point of equations (26)-(28) defined by the equation $`D_2=0`$. The point is an analogy to the inner singular point discussed in §3.1. To be consistent with §3.1 we use notation $`(r_s)_{in}`$ for the position of the point. At $`(r_s)_{in}`$ we have

$$_s^2=\left(8\frac{3\beta _s^2}{43\beta _s}\right)\left(7\frac{3}{2}\beta _s\frac{1+\beta _s}{43\beta _s}\right)^1.$$
$`(31)`$
Note, that the expression for $`_1`$ given by (24) coincides with (31). The latter could mean that the properties of global solutions in the case of viscosity prescriptions (1) and (2) are very similar in the inner part of flow at the limit of small viscosity, $`\alpha 1`$. Our numerical results confirm this conclusion.

Abramowicz & Kato (1989) studied analytically the type of singular point in the isothermal disks in the case of viscosity prescription (2). They showed that the point is always a saddle, and there is no case of a node. This conclusion differs from one obtained in the case of viscosity prescription (1). Our numerical models confirm this dependence of the singular point type on a form of viscosity.

## 4 Numerical results

To be used in the numerical method the sets of differential equations (17)-(18) and (26)-(28) have been re-written in the dimensionless form using the following dimensionless quantities: $`\stackrel{~}{r}=r/r_g`$, $`\stackrel{~}{v}=v/c`$, $`\stackrel{~}{c}_s=c_s/c`$, $`\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }r_g/c`$, $`\stackrel{~}{\mathrm{}}=\mathrm{}/r_gc`$, $`m=M/M_{}`$. In the subsequent discussions we will use these dimensionless quantities skipping in the notations the ‘tilde’ mark. The used method is described in Appendix A. We have calculated a number of models varied by the viscosity prescriptions and parameters $`\alpha `$ and $`\dot{m}`$. The black hole mass $`m`$ contributes into the dimensionless equations in the combinations

$$D=\frac{ac^4\kappa GM_{}}{^4}\frac{m}{\dot{m}}.$$

The parameter $`D`$ was taken to be inversely proportional to $`\dot{m}`$ with $`m=10`$, $`=1.6510^8ergg^1K^1`$ and $`\kappa =0.4cm^2g^1`$. The numerical grid covers the radial range form $`r_{in}`$ located at the inner singular point, $`r_{in}=(r_s)_{in}`$, till $`r_{out}10^4r_g`$.

We discuss first the influence of the numerical outer boundary conditions on our models. We have found that the models are insensitive to the specific values of the outer boundary conditions. By fixing $`\alpha `$ and $`\dot{m}`$ the unique global transonic solution is fully determined. This solution also uniquely determines the outer boundary values: two or three values depending on the used viscosity prescription (1) or (2), respectively. In general, we do not know a priori these ‘correct’ boundary values which the global solution passes through, and consequently, our assumed numerical boundary values are quite arbitrary. But, this values should be close enough to the ‘correct’ ones due to reason of numerical stability. Also, the ‘correct’ outer values are very close, but not exactly equal, to the values obtained from the standard model for the Keplerian accretion disks. Calculations show that all numerical solutions which have the same $`\alpha `$ and $`\dot{m}`$, but different boundary values at different $`r_{out}`$, converge to the ‘common’ solution which is not affected by the outer boundary. This ‘common’ solution represents the global solution which we seek. Significant differences between some numerical solution and the ‘common’ one (with relative errors $`10^4`$) are observed only in $`23`$ grid points before the last outer point at $`r_{out}`$. Such a behaviour of the numerical solutions can be explained by special properties of difference equations (A1) at $`\epsilon 1`$. At small $`\epsilon `$ the method is unstable.

Each model is characterized by value of $`\mathrm{}_{in}`$ \[see eq.(15)\] which has a sense of a specific angular moment of matter infalling into black hole. Figure 2 \[panels (a) and (b)\] shows the dependence of $`\mathrm{}_{in}`$ on accretion rate $`\dot{m}`$ for three values of $`\alpha =0.01`$, $`0.1`$ and $`0.5`$, and two forms of viscosity prescription (1) and (2). At low $`\dot{m}1`$ the value of $`\mathrm{}_{in}`$ is independent of $`\dot{m}`$, but weakly varies with $`\alpha `$. In the low $`\alpha `$ case, $`\alpha =0.01`$ and $`0.1`$, the values of $`\mathrm{}_{in}`$ are close to the minimum value of the Keplerian angular momentum, $`(\mathrm{}_K)_{min}=3.6742`$. At high $`\dot{m}0.1`$ the values of $`\mathrm{}_{in}`$ deviate from $`(\mathrm{}_K)_{min}`$ to larger or smaller values depending on $`\alpha `$. In the case $`\alpha =0.01`$ and $`0.1`$ one can see only minor differences between models with different forms of viscosity. But, for large $`\alpha =0.5`$ the difference in values of $`\mathrm{}_{in}`$ increases. Unfortunately, we had been able to calculated only a limited number of models in the case of viscosity prescription (2) due to technical reason (see Appendix A for details), and our comparison of both prescriptions is not complete in this respect.

Figure 2 \[panels (c) and (d)\] shows locations of the inner singular points $`(r_s)_{in}`$ as a function of $`\dot{m}`$ for different values of $`\alpha `$ and two different viscosity prescriptions. Similar to the case of $`\mathrm{}_{in}`$ the models at low $`\dot{m}`$ show a weak dependence of $`(r_s)_{in}`$ on $`\dot{m}`$. In the low $`\alpha `$ models (squares and circles in Figure 2) the values of $`(r_s)_{in}`$ are close to the location of the black hole last stable orbit at $`r=6`$. At high $`\dot{m}0.1`$ the values of $`(r_s)_{in}`$ are decreasing functions of $`\dot{m}`$ in the case of low $`\alpha =0.01`$ and $`0.1`$, and non-monotonically behave in the case of $`\alpha =0.5`$ (triangles in Figure 2).

Figure 3 shows locations of the outer singular points $`(r_s)_{out}`$ in the case of viscosity prescription (1) as a function of $`\dot{m}`$ for different values of $`\alpha `$. Values of $`(r_s)_{out}`$ are an increasing function of $`\dot{m}`$ and show a power-law behaviour at $`\dot{m}3`$. It is interesting to note that values of $`(r_s)_{out}`$ are almost independent of $`\alpha `$.

Examples of the specific angular momentum distribution $`\mathrm{}(r)`$ are shown in Figure 4 for $`\dot{m}=160`$ and three values of $`\alpha =0.01`$ (short-dashed line), $`0.1`$ (solid line) and $`0.5`$ (dotted line). The distributions correspond to viscosity prescription (1). The location of the inner singular points are indicated by the correspondent points on the curves. The Keplerian angular momentum is displayed by the long-dashed line for comparison. Only the low viscosity model with $`\alpha =0.01`$ has a super-Keplerian part in $`\mathrm{}(r)`$. Models with larger viscosity are everywhere sub-Keplerian. Note, that the singular point in the low viscosity model (short-dashed line) locates in the inner sub-Keplerian part of the disk.

Figure 5 shows the dependence of $`\beta _s`$ on $`\dot{m}`$ at the inner singular points. The change of value of $`\beta `$ form 1 to 0 corresponds to the change of a state from the gas pressure to radiative pressure dominated one. The thin disks with $`\beta 1`$ are locally stable, whereas the parts of the disk in which $`\beta 0`$ are thermally and viscously unstable at $`\dot{m}100`$ (Pringle, Rees, & Pacholczyk 1973). At larger $`\dot{m}100`$ the instability can be suppressed by the advection effect (Abramowicz et al. 1988). We have found a weak dependence of $`\beta _s(\dot{m})`$ on the assumed viscosity prescriptions.

Using analysis discussed in Appendix B we have determined a type of singular points in our numerical solutions. In Figures 2-5 the saddle-type points are indicated by the solid squares, circles and triangles. The nodal-type points are represented by the corresponding empty dots in the same figures. In the case of viscosity prescription (1) the solutions have two singular points (see §3.1). The inner singular points, $`(r_s)_{in}`$, can be saddles or nodes depending on values of $`\alpha `$ and $`\dot{m}`$. Note that the change of type from a saddle to nodal one does not introduce any features in solutions. The outer singular points, $`(r_s)_{out}`$, are always of a saddle-type. In the case of viscosity prescription (2) the solutions have only inner singular points (see §3.2) which are always of a saddle-type.

In models with low $`\dot{m}16`$ and low $`\alpha 0.1`$ the values of $`r_s`$ and $`\mathrm{}_{in}`$ are quantitatively very close to the last stable orbit location ($`r_{in}=6`$) and value of $`\mathrm{}_{in}`$ \[$`\mathrm{}_{in}=(\mathrm{}_K)_{min}`$\] assumed in the standard model (Shakura & Sunyaev 1973). The radial structure in our models is also very close to the one for the standard model in the same range of $`\dot{m}`$ and $`\alpha `$. Such a good coincidence means that the advective terms in equations (6) and (8) are negligiblly small in considered models. However, the high $`\alpha `$ models show quite significant deviation from the standard model independently of $`\dot{m}`$ (see Figure 2).

At high accretion rates, $`\dot{m}16`$, the effect of advection becomes significant. We illustrate it by calculating the luminosity $`L`$ of disk in the case of viscosity prescription (1),

$$L=4\pi _{(r_s)_{in}}^{\mathrm{}}F^{}r𝑑r=2L_{Edd}_{(\stackrel{~}{r}_s)_{in}}^{\mathrm{}}\frac{(1\beta )\stackrel{~}{c}_s\stackrel{~}{r}^{1/2}}{\stackrel{~}{r}2}𝑑\stackrel{~}{r},$$
$`(32)`$
where $`F^{}`$ is given by (11). Figure 6 shows calculated dependences of $`L/L_{Edd}`$ on $`\dot{m}`$ for different values of $`\alpha =0.01`$, $`0.1`$, $`0.5`$ (short-dashed, solid and dotted lines, respectively). There is a simple linear relation $`L/L_{Edd}=\eta \dot{m}`$ in the standard model, in which the advection is neglected. The radiative efficiency $`\eta `$ is a constant and equals $`\eta =1/16`$ in the case of gravitational potential (3). We plot this relation by the straight long-dashed line in Figure 6. One can clearly see from the figure that effect of advection results in reduction of luminosities with respect to the one for the standard model at $`\dot{m}16`$. There is a weak dependence of the luminosity on value of viscosity in disks.

Finally note, that our numerical solutions corresponding to viscosity prescription (1) have some resemblance to the results of Abramowicz et al. (1988), but they show important quantitative differences, especially for large $`\alpha `$ and $`\dot{m}`$.

## 5 Discussion

We have obtained unique solutions for structure of advective accretion disk in a wide range of accretion rates and $`\alpha `$-parameters. Both viscosity prescriptions (1) and (2) have been investigated. The solutions corresponding to both prescriptions are very close for $`\alpha 0.1`$, and begin to differ at larger $`\alpha `$. This is connected, probably with larger deviation of the angular velocity $`\mathrm{\Omega }`$ from the Keplerian one, $`\mathrm{\Omega }_K`$, leading to larger difference between $`t_{r\varphi }`$ in both prescriptions. Unfortunately, our comparison of the prescriptions is not complete due to technical problems in calculation of the high viscosity models in the case of viscosity prescription (2).

The main difference of the present study from the previous ones is in using more sophisticated numerical technique which accurately treats the regularity conditions in the inner singular point of equations. We have performed an analytical expansion at the singular point to calculate the derivatives of physical quantities. These derivatives help us to find the proper integral curve passing through the singular point. The approach allow us to avoid numerical instabilities and inaccuracies, appearing when only variables at the singular point, but not its derivatives, are included into a numerical scheme.

We have found different behaviour of integral curves depending on used viscosity prescription. In the case of viscosity prescription (1) there are two singular points located at $`(r_s)_{in}`$ and $`(r_s)_{out}`$. The inner point, $`(r_s)_{in}`$, locates close to the last stable black hole orbit (see Figure 2), and is an analogy of the singular point in spherical flow, where the point divides the subsonic and supersonic regions. The location of the outer point, $`(r_s)_{out}`$, is determined by the accretion rate (see Figure 3). At low $`\alpha 0.1`$ both points are of a saddle-type. Only one integral curve (“separatrix”) simultaneously crosses two saddle-type points, as it is shown in Figure 1a, and corresponds to the global solution which smoothly connects the supersonic innermost region of the accretion disk and the subsonic outer (formally at $`r=\mathrm{}`$) parts. For larger $`\alpha 0.1`$ the inner singular point changes its type to a node. There was suggestion by Muchotrzeb-Czerny (1986) and Abramowicz et al. (1988) that there is no unique solution in this case, because of all integral curves cross the node. Existence of a unique separatrix crossed simultaneously both singular points preserves a uniqueness of the solution in this case (see Figure 1b). The conclusion of Muchotrzeb-Czerny (1986) and Abramowicz et al. (1988) is probably connected with their neglection of outer singular points inherited to the problem.

Matsumoto et al. (1984) used slightly different form of equations (6) and (8), and they found that in this case only one singular point exists and changes type from a saddle to node. The difference with respect to our results arises because of using different form of the pressure gradient force \[the first term on the right hand side of equation (6)\]. We use the pressure and density taken at the equatorial plane in this term, whereas Matsumoto et al. (1984) used the vertically averaged quantities in it. In the latter case the term has the following form,

$$\frac{1}{\mathrm{\Sigma }}\frac{𝒫}{dr},\mathrm{where}\mathrm{\Sigma }=2_0^h\rho 𝑑z\mathrm{and}𝒫=2_0^hP𝑑z.$$
$`(33)`$
The vertically integrated approach (33) introduces the difference because in this case the free terms with $`\alpha ^2`$ in (20) are absent. Formally, it corresponds to location of the outer singular point at infinity, where $`_2=0`$. Such a visible difference is not qualitatively important for the physical solution, because conditions at the outer singular point are only shifted to infinity, and the integral curve itself has little changes. Thus, similar to our results, in the approach by Matsumoto et al. (1984) the inner critical points of both types, a saddle or node, correspond to a unique solution.

In the case of viscosity prescription (2) we have found that one singular point exist. The point is always of a saddle-type and determines a unique solution. Another results had been obtaining by Chen & Taam (1993). They also found that equations have one singular point, but the point changes its type from a saddle to nodal one, depending on $`\alpha `$ and the accretion rate. It is not clear why such a result was obtained. There are two main differences in our equations (6) and (8), and those used by Chen & Taam (1993). First, they used the same vertical averaging for the equation of motion as Matsumoto et al. (1984). Second, they used the vertically averaged energy equation which corresponds to the inappropriate polytropic relation, $`𝒫\mathrm{\Sigma }^\gamma `$, when one neglects terms corresponding to the viscous heating and radiative cooling. Here $`\gamma `$ is the effective adiabatic index, and other notations are similar to those used in (33). Our equation (8) corresponds to the correct polytropic relation, $`P\rho ^\gamma `$. It could be that one of the mentioned differences results in the change of critical point type.

Acknowledgments. This work was supported in part by RFBR through grant 99-02-18180, the Royal Swedish Academy of Sciences, the Danish Natural Science Research Council through grant No 9701841, Danmarks Grundforskningsfond through its support for establishment of the Theoretical Astrophysics Center.

## Appendix A Numerical method

We use the finite-difference method to solve the systems of ordinary differential equations discussed in $`\mathrm{\S }`$3. The method has some resemblance to that used by Igumenshchev, Abramowicz & Novikov (1998). In this approach the problem is reduced to a solution of a system of non-linear algebraical equations written for a each pair of neighboring numerical grid points. The numerical grid $`\{r_i\}`$ extends over about three orders of magnitude in the radial direction. We look for a numerical solution in which the location of the inner grid point $`r_1`$ coincides with the location of the singular point $`(r_s)_{in}`$ near black hole. In our method we approximate the differential equation $`dy/dr=f(r,y)`$ by the following finite differences,

$$\frac{y_iy_{i1}}{r_ir_{i1}}=\epsilon f_{i1}(1\epsilon )f_i,i=1,2,\mathrm{},I,$$
$`(A1)`$
where the function $`y(r)`$ must be replaced by $`v(r)`$, $`c_s^2(r)`$, and additionally by $`\mathrm{\Omega }(r)`$ in the case of viscous prescription (2), $`\epsilon `$ is a parameter, $`I`$ is the number of grid points, and the lower indices indicate the corresponding grid point. The value of parameter $`\epsilon `$ is chosen to provide a stability of the numerical scheme. We use $`\epsilon =1`$ in most of the cases.

We use the Newton-Raphson iteration scheme to solve the set of equations (A1). In general formulation equations (A1) can be represented by $`K`$ functional relations, involving $`K`$ variables $`y_k`$,

$$F_k(y_1,y_2,\mathrm{},y_K)=0,k=1,2,\mathrm{},K,$$
$`(A2)`$
or in the vector notation, $`𝐅(𝐲)=0`$. The $`(n+1)`$-iteration improvement of an approximate solution $`𝐲^n`$ of (A2) has the form,

$$𝐲^{n+1}=𝐲^n+\omega ^n\delta 𝐲^n,$$
$`(A3)`$
where $`\omega ^n`$ is a parameter, $`\omega ^n1`$, and the correction $`\delta 𝐲^n`$ is the solution of the matrix equation

$$𝐉^n\delta 𝐲^n=𝐅^n.$$
$`(A4)`$
In (A4) $`𝐉`$ is the Jacobian matrix, $`J_{lm}F_l/y_m`$. The parameter $`\omega `$ should be chosen to optimize the convergency of the iteration process to a solution. We use the following form of $`\omega `$,

$$\omega =\frac{\eta }{max(\eta ,\mathrm{\Delta })},$$
$`(A5)`$
where the parameter $`\eta =0.03`$ and $`\mathrm{\Delta }`$ is the average relative correction,

$$\mathrm{\Delta }=\frac{1}{K}\underset{k=1}{\overset{K}{}}\left|\frac{\delta y_k}{y_k}\right|.$$

In the case of equations (17), (18) we have $`K=2(I1)`$ functional relations involved $`2I`$ variables, $`v_i`$ and $`(c_s^2)_i`$. Two variables, $`v_I`$ and $`(c_s)_I`$, must be fixed as boundary conditions, when the number of independent variables equals to the number of equations, and system (A2) can be solved. In the case of equations (26)-(28) we have $`K=3(I1)`$ and three boundary values, $`v_I`$, $`(c_s)_I`$ and $`\mathrm{\Omega }_I`$, which are fixed to provide a consistent solution of (A2).

The presence of the singular points $`(r_s)_{in}`$ which coincides with $`r_1`$ introduces additional complications to our method. To satisfy to two regularity conditions,

$$D=0\mathrm{and}N=0,$$
$`(A6)`$
at $`r_1`$ we modify the iteration procedure by the following way. We add the equation $`D=0`$ to the set of equations (A1) together with a new independent variable $`\mathrm{}_{in}`$. Applying the Newton-Raphson iteration scheme we obtain a solution in which $`D=0`$ at $`r_1`$, but $`N0`$ in general. To satisfy the condition $`N=0`$ at the point $`i=1`$ the appropriate choice of $`r_1`$ must be done. We find the correct value of $`r_1`$ using the bisection method in which $`r_1`$ is changed by displacing all the grid. We apply the Newton-Raphson iterations (A4) for each new grid location.

To correctly approximate the differential equations in the first grid interval ($`r_1,r_2`$), where $`r_1=(r_s)_{in}`$, we expand solution at $`r_1`$,

$$v(r)=v_1+v_s^{}(rr_1),$$
$`(A7)`$
$$c_s^2(r)=(c_s^2)_1+(c_s^2)_s^{}(rr_1).$$
$`(A8)`$
The procedure of calculation of the coefficients $`v_s^{}`$ and $`(c_s^2)_s^{}`$ are given in Appendix B. Using (A7) and (A8) we fix the values at the second grid point as follows,

$$v_2=v_1+v_s^{}(r_2r_1),$$
$`(A9)`$
$$(c_s^2)_2=(c_s^2)_1+(c_s^2)_s^{}(r_2r_1).$$
$`(A10)`$
Relations (A9) and (A10) are used instead of the correspondent difference equations in (A1). The modified system of difference equations avoids numerical instabilities connected with the presence of the inner singular point and allows us to obtain a solution crossed continuously this point.

We have found that described numerical procedure becomes unstable in the case of equations (26)-(28) at large values of $`\alpha `$ and $`\dot{m}`$. The numerical instability arises because of influence of equation (26) and results in small-scale oscillations of all quantities. The numerical instability of similar type was found by Beloborodov (1998), who suppressed the oscillations on the level of $`10^3`$ of relative amplitude by applying a smoothing procedure to $`\mathrm{\Omega }`$. Unfortunately, such a smoothing procedure can not be included into our method because of additional coupling of equations (A1) with the regularity conditions (A6).

## Appendix B Expansion at singular point

We consider first the case of equations (26)-(28). We expand the numerator $`N_2`$ and denominator $`D_2`$ at the singular point $`r_s`$, as follows

$$N_2(r)=\left(\frac{N_2}{r}+\frac{N_2}{v}v_s^{}+\frac{N_2}{c_s^2}(c_s^2)_s^{}+\frac{N_2}{\mathrm{\Omega }}\mathrm{\Omega }_s^{}\right)(rr_s),$$
$`(B1)`$
$$D_2(r)=\left(\frac{D_2}{r}+\frac{D_2}{v}v^{}+\frac{D_2}{c_s^2}(c_s^2)^{}\right)(rr_s).$$
$`(B2)`$
In (B1) and (B2) the partial derivatives are taken at $`r=r_s`$, $`v=v_s`$, $`c_s^2=(c_s^2)_s`$ and $`\mathrm{\Omega }=\mathrm{\Omega }_s`$, and we use notations $`v_s=v(r_s)`$, $`(c_s^2)_s=c_s^2(r_s)`$ and $`\mathrm{\Omega }_s=\mathrm{\Omega }(r_s)`$. The ‘prime’ mark means the radial derivative of the correspondent quantity in (B1) and (B2). In (B2) we take into account that $`D_2/\mathrm{\Omega }=0`$. We denote for convenience,

$$\xi _s=r_s\frac{v_s^{}}{v_s},\eta _s=r_s\frac{(c_s^2)_s^{}}{(c_s^2)_s},\chi _s=r_s\frac{\mathrm{\Omega }_s^{}}{\mathrm{\Omega }_s}.$$

The value of $`\chi _s`$ can be defined with help of (26). Substituting (B1) and (B2) into equations (27) and (28) we finally obtain a quadratic equation with respect to $`\xi _s`$:

$$\left(v_s\frac{D_2}{v}+2(c_s^2)_s\frac{D_2}{c_s^2}A\right)\xi _s^2+\left[r_s\frac{D_2}{r}+2(c_s^2)_s\left(\frac{D_2}{c_s^2}B\frac{N_2}{c_s^2}A\right)v_s\frac{N_2}{v}\right]\xi _s$$

$$\left(r_s\frac{N_2}{r}+\mathrm{\Omega }_s\frac{N_2}{\mathrm{\Omega }}\chi _s+2(c_s^2)_s\frac{N_2}{c_s^2}B\right)=0,$$
$`(B3)`$
where we denote

$$A=1_s^2,\mathrm{and}B=\frac{3}{2}+\frac{r_s}{r_s2}+\frac{1}{(c_s^2)_s}\left(\mathrm{\Omega }_s^2r_s^2\frac{r_s}{(r_s2)^2}\right).$$

Having $`\xi _s`$ from (B3) one obtain $`\eta _s=A\xi _s+B`$. The found values of $`\xi _s`$, $`\eta _s`$ and $`\chi _s`$ determine derivatives $`v_s^{}`$, $`(c_s^2)_s^{}`$ and $`\mathrm{\Omega }_s^{}`$ which we use in the numerical procedure discussed in Appendix A. Equation (B3) has two roots. Which root should be used is defined by convergency of the iterations. We have found that only one of the two roots corresponds to the convergent solution. The partial derivatives of $`N_2`$ and $`D_2`$ used in (B3) can be calculated analytically or using numerical differentiation from (29) and (30). The latter method is simpler, but less accurate than the analytic one. The analytic derivation requires some efforts and produces quite long expressions which we do not present here. In the numerical procedure we have used the analytical expressions for derivatives, but in addition, we have checked them using estimates from the numerical differentiation.

To judge the type of critical point we follow the procedure described by Kato et al. (1998). We introduce a new variable $`\tau `$ defined by

$$d\tau =\frac{dr}{rD_2}.$$
$`(B4)`$
From equation (26)-(28) one obtain

$$\frac{d\mathrm{\Omega }}{d\tau }=\chi \mathrm{\Omega }D_2,\frac{dv}{d\tau }=vN_2,\frac{dc_s^2}{d\tau }=2c_s^2(AN_2+BD_2).$$
$`(B5)`$
All of the variables are now expanded around those at $`r_s`$,

$$r=r_s+\mathrm{\Delta }r,v=v_s+\mathrm{\Delta }v,c_s^2=(c_s^2)_s+\mathrm{\Delta }c_s^2,\mathrm{\Omega }=\mathrm{\Omega }_s+\mathrm{\Delta }\mathrm{\Omega }.$$
$`(B6)`$
Substituting (B6) into (B4) and (B5) and retaining only linear terms, one obtain a system of linear differential equations with respect to $`\mathrm{\Delta }r`$, $`\mathrm{\Delta }v`$, $`\mathrm{\Delta }c_s^2`$ and $`\mathrm{\Delta }\mathrm{\Omega }`$. Assuming that these quantities depend on $`\tau `$ in the form $`\mathrm{exp}(\lambda \tau )`$, one obtain the following characteristic equation which determines the eigenvalues $`\lambda `$:

$$\lambda ^2\lambda \left[r_s\frac{D_2}{r}+v_s\frac{N_2}{v}+2(c_s^2)_s\left(A\frac{N_2}{c_s^2}+B\frac{D_2}{c_s^2}\right)\right]$$

$$2Ar_s(c_s^2)_s\left(\frac{N_2}{r}\frac{D_2}{c_s^2}\frac{D_2}{r}\frac{N_2}{c_s^2}+\frac{D_2}{c_s^2}\frac{\chi _s}{r_s}\mathrm{\Omega }_s\frac{N_2}{\mathrm{\Omega }}\right)$$

$$2Bv_s(c_s^2)_s\left(\frac{D_2}{v}\frac{N_2}{c_s^2}\frac{N_2}{v}\frac{D_2}{c_s^2}\right)+$$

$$r_sv_s\left(\frac{D_2}{r}\frac{N_2}{v}\frac{N_2}{r}\frac{D_2}{v}\right)v_s\frac{D_2}{v}\chi _s\mathrm{\Omega }_s\frac{N_2}{\mathrm{\Omega }}=0.$$
$`(B7)`$
Again, all the partial derivatives in (B7) are taken at $`r=r_s`$, $`v=v_s`$, $`c_s^2=(c_s^2)_s`$ and $`\mathrm{\Omega }=\mathrm{\Omega }_s`$. If two solutions of quadratic equation (B7) are real and have different signs, $`\lambda _1\lambda _2<0`$, then the singular point is a saddle. The point is a node if two real roots of (B7) have identical signs, $`\lambda _1\lambda _2>0`$. Complex conjugate roots of (B7) corresponds to the spiral singular point.

In the case of equations (17), (18) equations (B3) and (B7) must be modified by substituting $`N_1`$ and $`D_1`$ instead of $`N_2`$ and $`D_2`$, and assuming $`N_1/\mathrm{\Omega }=0`$. The notation for $`B`$ must be also changed to

$$B=\frac{3}{2}+\frac{r_s}{r_s2}+\frac{1}{(c_s^2)_s}\left[\left(\frac{\mathrm{}_{in}}{r_s}+\alpha \frac{(c_s^2)_s}{v_s}\right)^2\frac{r_s}{(r_s2)^2}\right].$$

References

Abramowicz, M. A., & Zurek, W. H. 1981, ApJ, 246, 314

Abramowicz, M. A., & Kato, S. 1989, ApJ, 336, 304

Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646

Beloborodov, A. M. 1998, MNRAS, 297, 739

Chakrabarti, S.K. 1996, ApJ, 464, 664

Chakrabarti, S.K., & Molteni, D. 1993, ApJ, 417, 671

Chen, X., & Taam, R. E. 1993, ApJ, 412, 254

Frank, J., King, A., & Raine, D. 1992, Accretion Power in Astrophysics (Cambridge Univ. Press)

Fukue, J. 1987, PASJ, 39, 309

Höshi, R., & Shibazaki, N. 1977, Prog. Theor. Phys., 58, 1759

Igumenshchev, I. V., Abramowicz, M. A., & Novikov, I. D. 1998, MNRAS, 298, 1069

Kato, S., Fukue, J., & Mineshige, S. 1998, Black-Hole Accretion Disks (Kyoto: Kyoto Univ. Press)

Liang, E. P. T., & Thompson, K. A. 1980, ApJ, 240, 271

Matsumoto, R., Kato, S., Fukue, J., & Okazaki, A. T. 1984, PASJ, 36, 71

Muchotrzeb, B., & Paczyński, B. 1982, Acta Astr., 32, 1

Muchotrzeb, B. 1983, Acta Astr., 33, 79

Muchotrzeb-Czerny, B. 1986, Acta Astr., 36, 1

Narayan, R., Kato, S., & Honma, F. 1997, ApJ, 476, 49

Novikov, I. D., & Thorne, K. S. 1973, in C. DeWitt & B. DeWitt, eds, Black Holes. Gordon & Breach, New York, 343

Paczyński, B., & Wiita, P. J. 1980, A&A, 88, 23

Paczyński, B., & Bisnovatyi-Kogan, G.S. 1981, Acta Astr., 31, 283

Pringle, J. E., Rees, M. J., & Pacholczyk, A. G. 1973, A&A, 29, 179

Pringle, J. E. 1981, ARA&A, 19, 137

Shakura, N.I. 1972, Astron. Zh., 49, 921

Shakura, N.I., & Sunyaev, R.A. 1973, A&A, 24, 337