text:
# 1 Introduction

## 1 Introduction

Seeking a realistic model of our world, from the higher-dimensional standpoint, the brane-world approach is now vigorously envestigated in various ways. It started with the wall configuration in 5 dimensions(D) space-time. Basically the brane-world physics makes use of the soliton configuration. We can explore new physics, including the beyond-standard model, by making use of the distinguished properties of soliton, such as the localized zero modes, stable configuration, singularity-free behavior due to the extendedness, etc.. From the reason of improving some difficulties in 5D model, or of getting more general aspects of the brane-world physics, further-higher-dimensional models ( string-like object in 6D, monopole configuration in 7D, $`\mathrm{}`$) are also envestigated. In this paper, we examine a 6D model.

In the original model, 3-brane(s) is introduced, by hand, as a $`\delta `$-function distribution in the extra space. We do not take such approach in order to seek a non-singular solution. We induce the same configuration as a soliton (kink) solution. This situation is very similar to the relation between the Dirac and ’tHooft-Polyakov monopoles. The latter is the soliton (in the gauge+Higgs system) interpretation of the former system where the gauge potential has $`\delta `$-function distribution in the plane perpendicular to the Dirac string. We seek a soliton solution in the 6D gravity + Higgs system. The main purpose is to establish the pole solution . The most crucial point is to confirm the convergence of some infinite series appearing in the solution. It guarantees the boundary condition. In  they are solved up to the 2nd order approximation. Here we show them up to the 6th order. The result strongly shows that the series converge and the boundary conditions are sufficiently satisfied.

## 2 The Dimensional Reduction Model in Six Dimensions

In  a solution with the pole configuration in six dimensions is presented. The model is the 6D gravity coupled with the Higgs fields.

$`S[G_{AB},\mathrm{\Phi }]={\displaystyle d^6X\sqrt{G}(\frac{1}{2}M^4\widehat{R}G^{AB}_A\mathrm{\Phi }^{}_B\mathrm{\Phi }V(\mathrm{\Phi }^{},\mathrm{\Phi }))},`$
$`V(\mathrm{\Phi }^{},\mathrm{\Phi })={\displaystyle \frac{\lambda }{4}}(|\mathrm{\Phi }|^2v_{0}^{}{}_{}{}^{2})^2+\mathrm{\Lambda },`$ (1)
where $`(X^A)(x^\mu ,\rho ,\phi ),\mu =0,1,2,3`$. $`x^\mu `$’s are regarded as our world coordinates, whereas $`(\rho ,\phi )`$ the extra ones( $`0\rho <\mathrm{},0\phi <2\pi `$). $`G_{AB}`$ is the 6D metric field, $`\mathrm{\Phi }^{}`$ and $`\mathrm{\Phi }`$ are the complex scalars (Higgs fields). $`M(>0)`$ is the 6D Planck mass. $`\lambda (>0),v_0(>0)`$ and $`\mathrm{\Lambda }`$ in the potential $`V`$ are vacuum parameters. Expecting Randall-Sundrum type dimensional reduction, we assume the following form as the line element.

$`ds^2=\mathrm{e}^{2\sigma (\rho )}\eta _{\mu \nu }dx^\mu dx^\nu +d\rho ^2+\rho ^2\mathrm{e}^{2\omega (\rho )}d\phi ^2,`$
$`0\rho <\mathrm{},0\phi <2\pi ,`$ (2)
where $`\eta _{\mu \nu }=\text{diag}(1,1,1,1)`$. This is a natural 6D minimal-extension of the original 5D model by Randall-Sundrum . Two ”warp” factors, $`\mathrm{e}^{2\sigma (\rho )}`$ and $`\mathrm{e}^{2\omega (\rho )}`$, appear. For the fixed $`\rho `$ case ($`d\rho =0`$), the metric is the Weyl transformation of the product-space, the 4D Minkowski $`\times `$ the circle $`S^1`$. The coordinate $`\rho `$ can be regarded as the scaling parameter. We require the periodicity in $`\phi `$ for $`\mathrm{\Phi }(\rho ,\phi )`$ and assume the form: $`\mathrm{\Phi }=P(\rho )\mathrm{e}^{im\phi }(m=0,\pm 1,\pm 2,\mathrm{})`$. The 6D classical field equations of (1) are obtained as

$`3\sigma ^{\prime \prime }{\displaystyle \frac{\sigma ^{}}{\rho }}+\sigma ^{}\omega ^{}+\omega ^{\prime \prime }(\omega ^{})^2+2{\displaystyle \frac{\omega ^{}}{\rho }}=2M^4P_{}^{}{}_{}{}^{2},`$
$`16\sigma _{}^{}{}_{}{}^{2}+4{\displaystyle \frac{\sigma ^{}}{\rho }}4\sigma ^{}\omega ^{}+4\sigma ^{\prime \prime }=2M^4V,`$
$`4\sigma ^{\prime \prime }10\sigma _{}^{}{}_{}{}^{2}=M^4(P_{}^{}{}_{}{}^{2}m^2{\displaystyle \frac{\mathrm{e}^{2\omega }}{\rho ^2}}P^2+V).`$ (3)
We take a unit $`M=1`$ in the following for simplicity. 222 This is allowed by the following invariance of (1) or (3) , due to the simple dimensional counting: $`\mathrm{\Phi }M^2\mathrm{\Phi },(PM^2P),v_0M^2v_0,\lambda M^2\lambda ,\mathrm{\Lambda }M^6\mathrm{\Lambda },x^\mu M^1x^\mu ,\rho M^1\rho .`$

## 3 An Exact Solution for the No Potential Case

In order to see the structure of the solution, it is useful to consider a simple case, that is, the no potential case: $`\lambda =\mathrm{\Lambda }=0(V(\mathrm{\Phi }^{},\mathrm{\Phi })=0)`$. Even in this case, the scalars are still coupled with gravity. In this section we stress those properties which are shared with the general case of Sec.4. We consider the $`m=0`$ (no flux) solution. The following exact solution of the rational function type can be obtained.

$`\sigma ^{}(\rho )={\displaystyle \frac{1}{A\rho +B}},\omega ^{}(\rho )={\displaystyle \frac{4\rho +B}{\rho (A\rho +B)}},P(\rho )_{}^{}{}_{}{}^{2}={\displaystyle \frac{2(2A5)}{(A\rho +B)^2}},`$ (4)
where $`A(5/2)`$ and $`B`$ are integration constants. $`A`$ is a dimensionless constant, whereas $`B`$ is a dimensional one ( the dimension of length ). The lower bound for $`A`$ comes from the positivity of $`P_{}^{}{}_{}{}^{2}`$. This is a 2-parameters family solution (for $`\sigma ^{},\omega ^{}`$, and $`P^{}`$) at this stage. Some special cases are i) $`A=5/2,P^{}=0`$; ii) $`B=0,4\sigma ^{}=\omega ^{}=(4/A)\times \rho ^1,P_{}^{}{}_{}{}^{2}=(2(2A5)/A^2)\times \rho ^2`$; iii) $`A=4,\omega ^{}=\rho ^1`$; iv) $`A\mathrm{}`$, $`\sigma ^{}=0`$, $`\omega ^{}=0`$,$`P^{}=0`$. The integration constants are fixed by a required boundary condition.

Due to the condition $`A5/2`$, we do not find the solution which satisfies the RS-type configuration: $`\sigma ^{}`$ const.($`>0`$) as $`\rho \mathrm{}`$. However, if we restrict the region of $`\rho `$ as $`|B|A\rho `$, then we have $`\sigma ^{}1/B,\omega ^{}(1/\rho )(A4)/B,P_{}^{}{}_{}{}^{2}2(2A5)/B^2`$. This looks a RS-type metric although we cannot take $`\rho \mathrm{}`$.

The final solution for $`\sigma ,\omega `$ and $`P`$ are given as

$`\sigma (\rho )={\displaystyle \frac{1}{A}}\mathrm{ln}{\displaystyle \frac{|A\rho +B|}{D}},\omega (\rho )=\mathrm{ln}{\displaystyle \frac{|\rho ||A\rho +B|^{\frac{A4}{A}}}{E^{\frac{4}{A}}}},`$
$`P(\rho )=\pm {\displaystyle \frac{\sqrt{2(2A5)}}{A}}\mathrm{ln}{\displaystyle \frac{|A\rho +B|}{C}},`$ (5)
where $`C(>0),D(>0),E(>0)`$ are another integration constants with the dimension of length. They come from the symmetry of the constant translation ($`\sigma \sigma +\text{const.},\omega \omega +\text{const.},PP+\text{const.}`$) in eq.(3) with $`V=0`$. This final form (5) is a 5-parameters family of solutions.

The line element is given by

$`ds^2=\left\{{\displaystyle \frac{|A\rho +B|}{D}}\right\}^{2/A}\eta _{\mu \nu }dx^\mu dx^\nu +d\rho ^2+\left\{{\displaystyle \frac{|A\rho +B|}{E}}\right\}^{2(A4)/A}E^2d\phi ^2.`$
$`0<{\displaystyle \frac{2}{A}}{\displaystyle \frac{4}{5}},{\displaystyle \frac{6}{5}}{\displaystyle \frac{2(A4)}{A}}<2,`$ (6)
For the case of $`A\mathrm{}`$ the space-time becomes locally 6D Minkowski with the deficit angle $`\delta =2\pi (1A)=\mathrm{}`$. For the restricted region $`|B|A\rho `$, the above metric reduces to

$`ds^2=\mathrm{e}^{\frac{2}{B}\rho }\left\{{\displaystyle \frac{|B|}{D}}\right\}^{2/A}\eta _{\mu \nu }dx^\mu dx^\nu +d\rho ^2+\mathrm{e}^{\frac{2(A4)}{B}\rho }\left\{{\displaystyle \frac{|B|}{E}}\right\}^{2(A4)/A}E^2d\phi ^2.`$ (7)
As a form, this is a RS-type metric although we cannot keep it in $`\rho \mathrm{}`$ region.

The 6D Riemann scalar curvature is obtained as

$`\widehat{R}={\displaystyle \frac{4(2A5)}{(A\rho +B)^2}}0,`$ (8)
which is negative semi-definite. The absolute value goes to 0 as $`\rho \mathrm{}`$. For $`A=5/2`$ (the case i) above), $`\widehat{R}=0`$ is valid everywhere. The metric in this case

$`A={\displaystyle \frac{5}{2}},`$
$`ds^2=\left\{{\displaystyle \frac{|\frac{5}{2}\rho +B|}{D}}\right\}^{4/5}\eta _{\mu \nu }dx^\mu dx^\nu +d\rho ^2+\left\{{\displaystyle \frac{E}{|\frac{5}{2}\rho +B|}}\right\}^{6/5}E^2d\phi ^2.`$ (9)
For $`\mathrm{}>A>5/2`$, we can see three cases: a) $`B>0`$, There is no curvature singularity; b) $`B=0`$, The curvature is singular at $`\rho =0`$; c) $`B<0`$, The curvature singular at $`\rho =B/A`$. As for the horizon, from the condition $`0<\frac{2}{A}\frac{4}{5}`$, there is no horizon in the region $`0\rho <\mathrm{}`$. Hence, in b) and c), the singularities are naked ones.

The structure of the solution (6) is similar to the Kasner solution in 1+3 D. Ref. explains the cosmic string solutions in the Abelian Higgs model in 1+3 D. The two branches explained there, that is, the cosmic string branch and the Melvin branch look to have correspondence with the $`A\mathrm{}`$ case and $`A=5/2`$ case respectively.

## 4 Pole Solution in Six Dimensions

Now we consider for the case of the spontaneous breakdown (Higgs) potential: $`\lambda >0,v_0>0`$, and $`\mathrm{\Lambda }`$ is general(at this stage). The boundary condition for $`P`$ should be taken as

$`\rho \mathrm{},P(\rho )+v_0,`$ (10)
where we take the plus sign, using the freedom of $`\mathrm{\Phi }\mathrm{\Phi }`$ in (3), for the asymptotic value of $`P`$. We require (based on the analogy to 5D Randall-Sundrum model) that $`\sigma ^{}`$ and $`\omega ^{}`$ go to constants as $`\rho \mathrm{}`$. Then we can deduce, using (3), that

$`m=0\text{ (no flux)},`$
$`\rho \mathrm{},\sigma ^{}\alpha ,\omega ^{}\alpha ,\alpha =+\sqrt{{\displaystyle \frac{\mathrm{\Lambda }}{10}}}M^2.`$ (11)
This result says $`\mathrm{\Lambda }`$ should be taken negative: $`\mathrm{\Lambda }<0`$. We choose the plus sign as the asymptotic value so that the Weyl scaling factors in (2) work as shrinking the proper distance as $`\rho `$ increases. This choice guarantees the stableness of the system. As for the behavior near $`\rho =0`$(ultra-violet region), we can take, in a consistent way with (3), as

$`\rho +0,\sigma ^{}s\rho ^a,\omega ^{}w\rho ^b,Px\rho ^c`$
$`a=b=1,c=0,s0,w0,x0,`$
$`s={\displaystyle \frac{\lambda }{16}}(x^2v_{0}^{}{}_{}{}^{2})^2+{\displaystyle \frac{\mathrm{\Lambda }}{4}},`$ (12)
with the additional possible case: $`x=0`$ with $`c>1`$.

We examine the following form of solutions.

$`\sigma ^{}(\rho )=\alpha {\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{c_{2n+1}}{(2n+1)!}}\{\mathrm{tanh}(k\rho )\}^{2n+1},`$
$`\omega ^{}(\rho )=\alpha {\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d_{2n+1}}{(2n+1)!}}\{\mathrm{tanh}(k\rho )\}^{2n+1},`$
$`P(\rho )=v_0{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{e_{2n}}{(2n)!}}\{\mathrm{tanh}(k\rho )\}^{2n},`$ (13)
where $`c`$’s, $`d`$’s and $`e`$’s are coefficient-constants to be determined. $`\sigma ^{}`$ and $`\omega ^{}`$ are composed of odd powers of tanh($`k\rho `$), whereas $`P`$ is of even powers. They come from the behavior at the ultra-violet region (12). $`k`$ is a new scale parameter which shows the thickness of the pole. In the following, we take $`k=1`$ for simplicity. 333 This is because eq.(3) is invariant under the change: $`\rho k\rho ,\lambda (1/k^2)\lambda ,\mathrm{\Lambda }(1/k^2)\mathrm{\Lambda },v_0v_0`$.

The key equation for obtaining the solution of the form (13) is the following expansion formula about $`1/\rho `$. ( The factor $`1/\rho `$ appears in (3) ).

$`{\displaystyle \frac{1}{\rho }}={\displaystyle \frac{2}{\mathrm{tanh}\rho }}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{s_{2n}}{(2n)!}}\{\mathrm{tanh}\rho \}^{2n},{\displaystyle \frac{d^{2n}}{dx^{2n}}}\left({\displaystyle \frac{x}{\mathrm{ln}\frac{1+x}{1x}}}\right)|_{x=0}s_{2n},`$
$`0<\rho <\mathrm{},s_0={\displaystyle \frac{1}{2}},s_2={\displaystyle \frac{1}{3}},s_4={\displaystyle \frac{12}{5}},\mathrm{},`$ (14)
All coefficients are finite as shown above. We can take the limit $`\rho \mathrm{}`$ above, and see the infinite series converges and gives 0 : $`_{n=0}^{\mathrm{}}\frac{s_{2n}}{(2n)!}=0`$. We expect the coefficient-series appearing in (13) have the similar behaviors. It is a key, in the present treatment of the infinite series, that we take the expansion using powers of $`tanh\rho `$ not those of $`\mathrm{e}^{2\rho }`$. 444 We can not reexpress the RHS of the first equation of (14) as $`1/\rho =_{n=0}^{\mathrm{}}\frac{s_{2n}^{}}{(2n)!}\mathrm{e}^{2n\rho }`$ with finite coefficients. The same notice is said about (13).

We can show the above form of solutions indeed satisfy (3) if the coefficients satisfy some recursion relation. The first few orders are given as

$`\{\begin{array}{c}c_1=\frac{1}{4\alpha }\{\frac{1}{4}\lambda v_{0}^{}{}_{}{}^{4}(1e_{0}^{}{}_{}{}^{2})^2+\mathrm{\Lambda }\},\\ d_1=\frac{2}{3}c_1,\\ e_0:\text{free parameter},\end{array}`$ (18)
$`\{\begin{array}{c}c_3=\frac{3}{32}\frac{\lambda ^2v_{0}^{}{}_{}{}^{6}}{\alpha }e_{0}^{}{}_{}{}^{2}(1e_{0}^{}{}_{}{}^{2})^2+c_1(2+5\alpha c_1),\\ d_3=\frac{4}{3}c_1(1+5\alpha c_1),\\ e_2=\frac{1}{4}\lambda v_{0}^{}{}_{}{}^{2}e_0(1e_{0}^{}{}_{}{}^{2}).\end{array}`$ (22)
(23)
We note one free parameter $`e_0`$ appears 555 This freedom corresponds to, when $`\lambda =0`$, the constant translation symmetry ($`PP+\text{const.}`$) in eq.(3) .  and all coefficients are expressed by 4 parameters ($`\lambda ,v_0,\mathrm{\Lambda },e_0`$). The parameters, however, have 3 constraints.

$`1={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{c_{2n+1}}{(2n+1)!}},1={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d_{2n+1}}{(2n+1)!}},1={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{e_{2n}}{(2n)!}},`$ (24)
which comes from the boundary conditions (10) and (11). Therefore the present solution is a 1(=4-3) parameter family of solutions.

## 5 Evaluation of Coefficients and Analytic Result

We present here the results of concrete values of c’s, d’s and e’s for an input value $`e_0(=P(0)/v_0)=0.8`$. We solve constraints (24) by taking the first 7 terms (up to n=6th order). The most crucial point of the present solutions is to confirm the convergence of the infinite series $`_{n=0}^{\mathrm{}}\frac{c_{2n+1}}{(2n+1)!},_{n=0}^{\mathrm{}}\frac{d_{2n+1}}{(2n+1)!}`$, and $`_{n=0}^{\mathrm{}}\frac{e_{2n}}{(2n)!}`$, which guarantees the present boundary condition. The 6th-order approximation calculation gives the vacuum parameters as

$`v_0=0.9625,\mathrm{\Lambda }=2.725,\lambda =18.375.`$ (25)
666 These values should not be confused with numerical results like Sec.6. They should, in principle, be definitely determined just like the energy eigenvalues of the hydrogen atom are fixed by the boundary condition. 777 The constraint criterion is $`(1_{n=0}^6\frac{c_{2n+1}}{(2n+1)!})^2+(1_{n=0}^6\frac{d_{2n+1}}{(2n+1)!})^2+(1_{n=0}^6\frac{e_{2n}}{(2n)!})^2<10^4`$.

The obtained values of the coefficients are shown in Fig.1.

In the figure we also plot the data from the geometrical series: $`1/(1x)=1+x+\mathrm{}`$ at $`x=0.8`$. Comparing them, we can recognize the convergence of the coefficient-series (up to this approximation order). Note that all three series are ’oscillating’, which is advantageous for their convergence. Using these results, the analytical results of $`P(\rho ),\sigma ^{}(\rho )`$ and $`\omega ^{}(\rho )`$ ( (13))are shown in Fig.2.

We can see comparative behaviors of two ’warp’ factors ($`\sigma ^{},\omega ^{}`$) in the ultra violet region, $`\rho 2`$. In particular we notice a very ’delicate’ dip in the behavior of the second ’warp’ factor $`\omega ^{}`$ at $`k\rho 1`$. The (6D) Riemann scalar curvature is also shown in Fig.3.

It shows the curvature is positive both inside and outside of the pole, whereas negative (small absolute value) or zero between them. It is non-singular everywhere. There is no horizon in $`0\rho <\mathrm{}`$. The deficit angle at $`\rho =\mathrm{}`$ is $`2\pi `$.

## 6 Numerical Results

As the coupled differential equations for $`P(\rho ),\sigma ^{}(\rho )`$ and $`\omega ^{}(\rho )`$, the equations (3) with $`m=0`$, have the standard form of the numerical analysis, that is, Runge-Kutta method. We can numerically solve (3) without any ansatz about the form of the solution. In this approach, the choice of the initial values with high precision is required. In the present case we cannot take $`\rho =0`$ as the initial point because the appearance of the factor $`1/\rho `$ in (3). (Note this does not say the solution is singular at $`\rho =0`$. We seek a non-singular solution.) We take $`\rho =0.1`$ as the initial point of $`\rho `$. As for the initial values we borrow the results from Sec.5 : $`P(0.1)=0.79386,\sigma ^{}(0.1)=0.10488,\omega ^{}(0.1)=0.070046`$. The result by the Runge-Kutta calculation is shown in Fig.4.

It shows the analytic solution of Sec.5 is well reproduced in the ultraviolet region ( small-$`\rho `$ region ) but not for the infrared region. We should notice that the the consistent region (with the analytic result) , $`\rho 1.3`$, becomes definitely larger than that of the 2nd order approximation  , $`\rho 0.9`$. The numerical output data stop at $`\rho 1.7`$ with producing imaginary values. This occurs because keeping the positivity, $`P^{}(\rho )^20`$, becomes so stringent in the infrared region. The quantity becomes so small in that region and vanishes at $`\rho =\mathrm{}`$. We understand that further higher precision is required for the initial values in order to extend the consistent region furthermore.

## 7 Conclusion

We have presented the both analytical and numerical solution in the 6D reduction model. The approximation oder is $`n=6`$. Compared with the previous result of $`n=2`$ order one, the results improve definitely. Much evidence about the convergence of the coefficients series is obtained. An exact solution, for the no potential case, is also obtained.

We notice rather stable behavior of the following values under the increase of the approximation order : the vacuum parameters ($`v_0,\mathrm{\Lambda },\lambda `$); the coefficients ($`c_{2n+1},d_{2n+1},e_{2n}|n=0,1,\mathrm{},6`$); the form of $`\sigma ^{},\omega ^{}`$ and $`P`$. This strongly indicates the result of Fig.2 is quite near the exact result.

Acknowledgment

The author thanks P. van Nieuwenhuizen, M. Rocek and R. Shrock for discussions or comments about some points related to this work. Comments by A. Goldhaber, S.T. Hong, Z. Kakushadze, I. Oda, S. Vandoren and V. Zhukov are appreciated. He also thanks R. Roiban and C.M. Hung for helping him in relation to this work. The financial support by the governor of Shizuoka prefecture is greatly acknowledged. Finally he expresses gratitude to the hospitality at the C.N. Yang Institute for Physics, State University of New York at Stony Brook where this work has been done.


format: text

text:
# A new approach to the problem of modes in the Mestel disk

## 1 Introduction

It was not long after Kalnajs (kala (1971), kalb (1977), and references therein) and others had laid out a programme to perform stability analyses for galactic disks that it was applied to the Mestel disk by Zang (1976PhDT……..26Z (1976); henceforth cited as Zang) in a Ph.D. Thesis supervised by Alar Toomre. This disk, characterised by an infinitesimal height and a globally constant rotation curve, lends itself as a model for disks that show spiral structure mainly for two reasons. For one, its structure is so simple that one of its distribution functions has made it into a textbook example (Binney & Tremaine bintre (1987)); secondly, and probably more importantly, flat rotation curves seem to be rather prevalent among (high surface brightness) spiral galaxies (e.g., Persic & Salucci persal (1991)), and it has long been suspected that flat rotation curves and strong, “modal-looking” spiral structure are somehow linked (e.g., Biviano et al. biviano (1991)).

One of the drawbacks of the Mestel disk is that the orbits in it cannot be expressed in elementary functions. Zang circumvented this problem essentially through numerical integration and some elegant devices exploiting the self-similarity of the Mestel disk (i.e., orbits of a given eccentricity look like scaled copies of each other regardless of the angular momentum of a particle on them) to curb down on the computational effort. While the self-similarity is welcome in simplifying the calculation, it also implies that a self-consistent Mestel disk does not possess any scales. In particular, the absence of a time scale already leads to the expectation that discrete modes will be hard to construct.

Zang found that the Mestel disk is a rather bewildering construct. His formalism leads to a singular integral equation that is quite difficult to treat numerically. However, the structure of its kernel as well as physical arguments suggest that there are no discrete modes apart from rather exotic non-rotating and non-growing ones, implying that either no growing modes exist or modes become unstable regardless of their growth rate and pattern speed. Thus, his investigations concentrated on a variant of the Mestel disk in which a certain part of the disk is declared “immobile”, the cut-out disk. This tampered system is furnished with a length scale, admits discrete modes and becomes treatable with his numeric technique.

But even the cut-out disk behaves contrary to the preconceptions of the mid-70ies. While Zang would have liked to see a mildly unstable two-armed mode, typically the only mode that is growing in an axisymmetrically stable disk is a one-armed mode, and for this, a stability limit cannot be found. Even when he made the inner cut-out so steep that $`m=2`$ modes emerged in axisymmetrically stable disks in a way that Toomre (1977ARA+A..15..437T (1977)) describes as unrealistic, these one-armed modes, again in Toomre’s words continued to “plague” him.

Despite these somewhat unexpected findings and the fact that is was only fragmentarily published in Toomre (1977ARA+A..15..437T (1977)) and Toomre (toomreswing (1981)), Zang’s thesis became a frequently-cited classic and the Mestel disk developed a life of its own. Recently, Evans & Read (1998a , 1998a ; see also Read (1997PhDT………1R (1997); henceforth cited as ERI and ERII) generalised Zang’s work to the entire family of physically plausible disks with surface densities following a power law, and Goodman & Evans (goodev (1999)) once more tackled the issue of rotating modes the untampered (self-consistent) Mestel disk.

This work set out as an attempt to avoid the numerical integrations that had to be employed by Zang and ERI by approximating the orbits and eventually deriving a closed kernel of the integral equation governing the problem. Among the merits of the availability of such a kernel is that an analytic expression for the kernel might allow a decision whether or not the self-consistent Mestel disk admits rotating modes. Also, it might facilitate an analysis of the processes at work at various resonances of the modes in the cut-out disk in the context of a full disk without having to take recourse to somewhat daring techniques to extend WKBJ approximations across resonances (e.g., Mark mark (1976)).

The organisation of this paper is as follows: In section 2, we derive the orbits and the Hamiltonian of the disk. Section 3 basically follows Kalnajs’ programme that will in the end yield an (almost) closed expression for the kernel of an integral equation governing the modes in the disk. To assess the validity of this kernel, we will compare the results reported in ERII with ours in section 4. Finally, we revisit the issue of modes in the self-consistent disk in section 5.

## 2 The Orbits

The equations of motion in the Mestel disk are

$`\ddot{r}`$ $`=`$ $`r\dot{\theta }^21/r`$ (1)
$`L_z`$ $`=`$ $`r^2\dot{\theta }=\mathrm{const}.`$ (2)
Here, $`r`$ and $`\theta `$ denote polar coordinates, and $`L_z`$ is the $`z`$-component of the angular momentum per unit mass. Units have been chosen such that the circular velocity $`v_c=2\pi \mathrm{G}\mathrm{\Sigma }_0=1`$, where $`\mathrm{G}`$ is the gravitational constant, and $`\mathrm{\Sigma }_0`$ is the surface density at a reference radius that can be set to unity because of the self-similarity of the disk.

We seek closed expressions for the solution of the equations (2). Since the exact solutions cannot be expressed in elementary functions, we need to employ some approximation. A straightforward linearisation in configuration space – the epicyclic approximation – fails for the Mestel disk, since the Hamiltonian derived from it does not contain contributions from non-circular motion components. Therefore we linearise in Fourier space, a procedure generally known as the method of harmonic balance (Bogoliubov & Mitropolski bogmip (1961)).

Inserting one of the equations (2) into the other yields

$$\ddot{r}=L_z^2/r^31/r.$$
(3)
We now substitute $`r(t)=R_0(1+\zeta \mathrm{cos}(\kappa t))`$, making $`R_0`$ the reference radius of a given orbit, $`\zeta `$ a measure for its deviation from circularity, and $`\kappa `$ an epicyclic frequency. After this, the first Fourier coefficients of (3) with respect to $`t`$ are

$`a_0`$ $`=`$ $`{\displaystyle \frac{L_z^2\left(1+\zeta ^2/2\right)R_0^2\left(1\zeta ^2\right)^2}{R_0^3\left(1\zeta ^2\right)^{5/2}}}`$ (4)
$`a_1`$ $`=`$ $`{\displaystyle \frac{2R_0^2\left(1\zeta ^2\right)^2\left(\sqrt{1\zeta ^2}1\right)+3L_z^2\zeta ^2}{\zeta R_0^3\left(1\zeta ^2\right)^{5/2}}}`$ (5)
$`b_1`$ $`=`$ $`0,`$ (6)
where $`a_i`$ and $`b_i`$ denote the cosine and sine coefficients, respectively.

Dropping all higher Fourier coefficients and substituting the result back into (3), one finds by comparing the terms independent of $`t`$

$$L_z=\frac{\sqrt{2}R_0(1\zeta ^2)}{\sqrt{\zeta ^2+2}},$$
(7)
and by comparing terms proportional to $`\mathrm{cos}(\kappa t)`$

$$\kappa =\frac{\sqrt{2}}{\zeta R_0}\sqrt{12\frac{\sqrt{1\zeta ^2}}{\zeta ^2+2}}.$$
(8)
The circular frequency of the guiding centre $`\mathrm{\Omega }`$ results from a Fourier transformation of the second equation in (2) with (7) inserted and is given by

$$\mathrm{\Omega }=\frac{L_z}{R_0^2\left(1\zeta ^2\right)^{3/2}}.$$
(9)
We can now derive the Hamiltonian for a mass point in the Mestel disk in action-angle coordinates by simply inserting (7), (8), and (9) into an Hamiltonian written in Cartesian coordinates. Assuming orbits not too far from circularity, one can linearise with respect to $`\zeta `$. After again dropping all Fourier components of order 2 or higher, one is left with

$$H(J,L_z)=\frac{1}{2}+\mathrm{ln}(L_z)+\frac{\sqrt{2}J}{L_z},$$
(10)
where $`J=\kappa \zeta ^2R_0^2/2`$ is the energy in the epicyclic motion divided by $`\kappa `$ and thus the second integral. This clearly is the simplest Hamiltonian that could have anything to do with a Mestel disk and indeed the Hamiltonian that would result from a classic epicyclic approximation if it worked. As a matter of fact, at this point one has $`\kappa =\sqrt{2}/R_0`$ and $`\mathrm{\Omega }=1/R_0`$, reproducing the result of an epicyclic analysis. This Hamiltonian has already been used by Collett et. al. (1997MNRAS.285…49C (1997)).

The set of approximated action-angle coordinates $`L_z`$, $`J`$, $`w_1:=\kappa t`$, $`w_2:=\mathrm{\Omega }t`$ is correct to first order in $`\zeta `$, and the canonical equations are correct to first order in $`\zeta `$ as well. As is demonstrated in Demleitner (diss (2000)), the Jacobian of the transformation between Cartesian coordinates and this set is $`1+O(\zeta ^2)`$, showing that the transformation is indeed canonical to the required order.

Because of the two approximations (only keeping the lowest Fourier coefficients and linearising with respect to $`\zeta `$), it is not *a priori* clear that the disk described by this transformation and the resulting Hamiltonian (10) does indeed behave like a Mestel disk. It is one of the purposes of this work to show *a posteriori*, at is were, that this is the case.

## 3 Derivation of the matrix equation

With the Hamiltonian and the transformation equations, we can now follow the programme outlined in the series of papers by Kalnajs cited above. The first step is to establish a base for the potential perturbation. We choose

$`\varphi _{\alpha ,m}`$ $`=`$ $`L_z^{i\alpha 1/2}\mathrm{e}^{imw_2i\omega t}`$ (12)
$`\mathrm{e}^{\left(i\sqrt{\sqrt{2}J/L_z}\left(\alpha \mathrm{cos}(w_1)\sqrt{2}m\mathrm{sin}(w_1)\right)\right).}`$
Here, $`m`$ is the number of circumferential maxima of the spiral pattern (the number of arms), and $`\alpha `$ parametrises how closely the spiral is wound, where a spiral with larger $`\alpha `$ is more trailing provided that $`\omega >0`$. This expression basically is the well-known logarithmic spiral $`\mathrm{e}^{i\alpha \mathrm{ln}(r)+im\theta i\omega t}`$ transformed into our action-angle coordinates with harmonics in the angle coordinates higher than one disregarded. The factor $`L_z^{1/2}`$ was added to ensure definiteness in the integration of the angular momentum. This analogy motivates following ERI in naming $`\alpha `$ the logarithmic wave number and $`m`$ the angular harmonic number. We assume that all potential perturbations can be represented as expansions of the form

$$\underset{m=0}{\overset{\mathrm{}}{}}_{\mathrm{}}^{\mathrm{}}\mathrm{\Phi }_m(\alpha )\varphi _{\alpha ,m}(w_1,w_2,J,L_z)𝑑\alpha .$$
(13)
The next step in Kalnajs’ programme is to solve the linearised Boltzmann equation

$$\frac{f_1}{t}+[f_1,H]=[f_0,\varphi _{\alpha ,m}],$$
(14)
written with Poisson brackets $`[.,.]`$ for an elementary perturbations $`\varphi _{\alpha ,m}(w_1,w_2,J,L_z)`$ of the form (12). Here, $`f_1`$ is the unknown perturbation of the distribution function, and $`f_0`$ an equilibrium distribution function; for $`f_0`$, we use the distribution function given by Binney & Tremaine (bintre (1987)), which, in our units and coordinates, and after applying Stirling’s formula, is

$`f_0(L_z,J)`$ $`=`$ $`{\displaystyle \frac{F}{L_z}}\mathrm{exp}\left({\displaystyle \frac{\sqrt{2}J}{\sigma ^2L_z}}\right)\mathrm{with}`$
$`F`$ $`=`$ $`{\displaystyle \frac{\sqrt{12\sigma ^2}^{11/\sigma ^2}}{2\sqrt{2}\mathrm{e}\pi ^2\sigma ^2\mathrm{G}.}}`$ (15)
Expanding the Poisson brackets in (14) yields

$`{\displaystyle \frac{f_1}{t}}+`$ $`{\displaystyle \frac{\sqrt{2}}{L_z}}{\displaystyle \frac{f_1}{w_1}}+{\displaystyle \frac{1}{L_z}}{\displaystyle \frac{f_1}{w_2}}=i{\displaystyle \frac{FL_z^{i\alpha }\mathrm{e}^{imw_2i\omega t}}{\sigma ^2L_z^{5/2}}}`$ (17)
$`\left(X(w_1)m\sigma ^2\right)\mathrm{exp}\left(iX(w_1){\displaystyle \frac{\sqrt{2}J}{\sigma ^2L_z}}\right)`$
with an auxiliary function $`X(w_1)=\sqrt{\sqrt{2}J/L_z}\left(\alpha \mathrm{cos}(w_1)\sqrt{2}m\mathrm{sin}(w_1)\right)`$. On the right-hand side of this equation, a term of order $`\zeta ^2`$ has been neglected.

To solve this equation, we perform a separation of variables by

$$f_1(w_1,w_2,t)=g_1(w_1)\mathrm{e}^{imw_2i\omega t}$$
(18)
and arrive at an inhomogeneous ordinary differential equation, the homogeneous solution of which is

$$\stackrel{~}{g}_1(w_1)=\mathrm{e}^{i(\omega L_zm)w_1/\sqrt{2}}.$$
(19)
The general solution can be obtained via a variation of constants, where the integration constant is fixed by the condition that the distribution function must be $`2\pi `$-periodic, $`g_1(w_1)=g_1(w_1+2\pi )`$.

Before writing down the result, we note that in (17) the terms with $`w_1`$ in the exponent and its coefficient are quite similar to each other. This invites a partial integration in the resulting expression, which after some algebra leads to

$`g_1(w_1)`$ $`=`$ $`C\sqrt{{\displaystyle \frac{L_z}{\sqrt{2}J}}}\mathrm{exp}({\displaystyle \frac{\sqrt{2}J}{\sigma ^2L_z}})(i\mathrm{exp}\left(iX(w_1)\right)`$ (22)
$`{\displaystyle \frac{1}{\mathrm{e}^{2\pi i\eta }1}}\left(\sqrt{2}\eta +m\sigma ^2\right)`$
$`{\displaystyle _0^{2\pi }}\mathrm{exp}(iX(w_1+w_1^{})i\eta w_1^{})dw_1^{}).`$
In this expression, we have abbreviated $`\eta =(\omega L_zm)/\sqrt{2}`$.

Substituting (22) into (18) yields the perturbation of the distribution function under an elementary potential perturbation $`\varphi _{\alpha ,m}`$. Since (14) is linear in $`\varphi _{\alpha ,m}`$, the response to an arbitrary perturbation $`\mathrm{\Phi }_1^{\mathrm{imposed}}`$ of the form (13) can immediately be computed by summing over all $`\alpha `$ and $`m`$. Thus, the solution of (14) can be written as

$$f_1^{\mathrm{response}}=\mathrm{\Phi }_1^{\mathrm{imposed}}$$
(23)
with a linear integral operator $``$.

Further pursuing Kalnajs’ programme, we attempt to find a linear combination of elementary perturbations such that $`f_1^{\mathrm{response}}`$ causes the $`\mathrm{\Phi }_1^{\mathrm{imposed}}`$ we started with. To do this, the Poisson equation has to be solved simultaneously to the Boltzmann equation. Since $``$ and the Laplacian occurring in the Poisson equation do not commute, their eigenfunctions are different, and a straightforward attempt to solve the two equations simultaneously will result in an integro-differential equation that will be very hard to solve indeed.

To avoid this, one converts the system of equations to a matrix equation. The idea here is to compute potential-density pairs, i.e., two sets of mutually orthogonal basis functions $`\varphi _{\alpha ,m}`$ and $`\mu _{\alpha ,m}`$ for each of the perturbations in density and potential where each pair $`\varphi _{\alpha ,m},\mu _{\alpha ,m}`$ satisfies the Poisson equation. When one expands the perturbations in the potential and the surface density in these functions, the Poisson equation simply reads $`\mathrm{\Phi }_m(\alpha )=M_m(\alpha )`$ in the coefficients $`\mathrm{\Phi }_m(\alpha )`$ and $`M_m(\alpha )`$ of these expansions.

The basis $`\varphi _{\alpha ,m}`$ is already fixed by (12), and the corresponding density base function can be found using the scheme of Clutton-Brock (clubro (1972)). For our calculation, however, the only important property of the potential-density pair the scalar product of $`\varphi _{\alpha ,m}`$ and $`\mu _{\alpha ,m}`$. This scalar product is of course invariant under canonical transformations, and to the the required order Clutton-Brock’s result

$$_0^{2\pi }𝑑\theta _0^{\mathrm{}}r𝑑r\varphi _{\alpha ,m}^{}\mu _{\alpha ^{},m^{}}=4\pi ^2K_m(\alpha )\delta _{m,m^{}}\delta (\alpha \alpha ^{})$$
(24)
for logarithmic spirals holds. Here,

$$K_m(\alpha )=\frac{\left|\mathrm{\Gamma }\left(\frac{3}{4}+\frac{1}{2}m+\frac{1}{2}i\alpha \right)\right|^2}{\pi \mathrm{G}\left|\mathrm{\Gamma }\left(\frac{1}{4}+\frac{1}{2}m+\frac{1}{2}i\alpha \right)\right|^2}$$
(25)
is a quantity closely related to ERI’s Kalnajs factor.

After expanding the potential and surface density perturbations in (23), one can project both sides on $`\varphi _{\alpha ^{},m^{}}`$. Exploiting (24) and inserting the (now trivial) Poisson equation, (23) takes the form

$`4\pi ^2`$ $`K_m^{}(\alpha ^{})\mathrm{\Phi }_m^{}^{\mathrm{response}}(\alpha ^{})={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\alpha {\displaystyle _0^{\mathrm{}}}𝑑L_z{\displaystyle _0^{\mathrm{}}}𝑑J`$ (27)
$`{\displaystyle _0^{2\pi }}𝑑w_1{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle _0^{2\pi }}𝑑w_2\mathrm{\Phi }_m^{\mathrm{imposed}}(\alpha )(\varphi _{\alpha ,m})\varphi _{\alpha ^{},m^{}}^{}.`$
Equating $`\mathrm{\Phi }^{\mathrm{response}}`$ and $`\mathrm{\Phi }^{\mathrm{imposed}}`$, an integral equation results the solutions of which, if they exist, are self-consistent modes.

The integrand of (27) can be written as a sum $`I_1+I_2`$ of a part without a further integral,

$`I_1`$ $`=`$ $`\mathrm{\Phi }_{\alpha ,m}{\displaystyle \frac{FL_z^{i(\alpha \alpha ^{})2}}{4\pi ^2\sigma ^2K_m^{}(\alpha ^{})}}\mathrm{e}^{i(mm^{})w_2\frac{\sqrt{2}J}{\sigma ^2L_z}}`$ (29)
$`\mathrm{e}^{i\sqrt{\frac{\sqrt{2}J}{L_z}}\left((\alpha \alpha ^{})\mathrm{cos}(w_1)\sqrt{2}(mm^{})\mathrm{sin}(w_1)\right)},`$
and one in which a further integration over the auxiliary angle $`w_1^{}`$ has to be performed,

$`I_2`$ $`=`$ $`{\displaystyle _0^{2\pi }}𝑑w_1^{}i\mathrm{\Phi }_{\alpha ,m}{\displaystyle \frac{\sqrt{2}FL_z^{i(\alpha \alpha ^{})2}\left(m\sigma ^2+\sqrt{2}\eta \right)}{8\pi ^2\sigma ^2K_m^{}(\alpha ^{})(\mathrm{e}^{2\pi i\eta }1)}}`$ (32)
$`\mathrm{e}^{i(mm^{})w_2i\eta w_1^{}\frac{\sqrt{2}J}{\sigma ^2L_z}+\sqrt{2}m^{}\mathrm{sin}(w_1)}`$
$`\mathrm{e}^{i\sqrt{\frac{\sqrt{2}J}{L_z}}\left(\alpha \mathrm{cos}(w_1^{}+w_1)\alpha ^{}\mathrm{cos}(w_1)\sqrt{2}m\mathrm{sin}(w_1^{}+w_1)\right)}.`$
The quadratures over $`w_2`$, $`w_1`$, and $`J`$, carried out in that order, are rather straightforward. Details can be found in Demleitner (diss (2000)). When integrating over $`w_1^{}`$ one has to evaluate an expression of the type

$$_0^{2\pi }\mathrm{e}^{i\tau w_1^{}}\mathrm{e}^{\nu \mathrm{cos}(vw_1^{})}𝑑w_1^{},$$
(33)
which can be done by expanding the second factor of the integrand into modified Bessel functions $`\mathrm{I}_k`$ using

$$\mathrm{e}^{\nu \mathrm{cos}(vw_1^{})}=\mathrm{I}_0(\nu )+2\underset{k=1}{\overset{\mathrm{}}{}}\mathrm{I}_k(\nu )\mathrm{cos}\left(k(vw_1^{})\right)$$
(34)
(Abramowitz & Stegun abrsteg (1972)). Unfortunately, the series introduced by this operation does not in general collapse and will be present in the kernel.

The integration over the angular momentum is the most involved one, resembling the calculations to compute the angular momentum function in Zang’s formalism. For the summands of the series resulting from the $`w_1^{}`$ integration, one has to treat integrals of the structure

$$_{\mathrm{}}^{\mathrm{}}\frac{\mathrm{e}^{(i\gamma +\lambda +1)u}}{2\left(2\xi ^2(1\mathrm{e}^u)^2\right)}𝑑u,$$
(35)
where

$$u=\mathrm{ln}(\omega L_z/m)$$
(36)
(note that this substitution only works for nonzero $`m`$ and $`\omega `$), $`\gamma =\alpha \alpha ^{}`$ and $`\xi =k/m`$. The integer $`k`$ is a summation index varying between $`0`$ and $`\mathrm{}`$, and the complete integral requires this expression to be evaluated for $`\lambda \{1,0,1\}`$.

We will compute this integral using the residue theorem, and thus the first step is to find the zeroes of the denominator of the integrand. These are at

$$u_j^{(1,2)}=\mathrm{ln}(1\pm \sqrt{2}\xi )+2\pi ij$$
(37)
for integer $`j`$. These zeroes are well-defined, since $`\xi `$, as the ratio of two integers, can never be equal to $`1/\sqrt{2}`$. Let for the moment $`\xi <1/\sqrt{2}`$ and $`\gamma >0`$.

There are three more or less technical problems in the application of the residue theorem to this integral. The first is that at least one pole is on the real axis and hence on the integration contour. The tempting idea to take the (existing) Cauchy principal value leads to erroneous results. As Landau (landau (1946)) noted in a plasma-physical context, the correct reasoning is that if a system supports a non-growing perturbation at present and started in equilibrium, the perturbation must have been a growing one at some time in the past. Since it has evolved continuously, the integration contours of the problem must have been continuously deformed as well and thus cannot have jumped across poles. A growing perturbation corresponds to an $`\omega `$ with a positive imaginary part, which with the substitution (36) means that the integration contour moves into the upper half-plane. Therefore, as $`\mathrm{}\omega 0`$, one has to leave the zeroes on the real axis out of the contour, as is sketched in Fig. 1.

The second problem is that as the size of the integration path $`R\mathrm{}`$ (cf. Fig. 1) one cannot continuously deform the contour since the poles are in the way, and the same poles make it hard to establish that the contribution of the line at $`\mathrm{}u=iR`$ vanishes. The standard technique to avoid this is to write down the series of residues and show that integral over the difference between this series and the integrand vanishes as $`R\mathrm{}`$. Since this differences has no poles, the contour can be readily deformed.

The last problem is that for $`\lambda =\pm 1`$ the contributions at $`\mathrm{}u=\pm R`$ do not necessarily vanish. Indeed, one finds that these “border terms”

$$=\{\begin{array}{cc}\frac{i}{4\gamma \xi ^22}lim_x\mathrm{}\frac{\mathrm{e}^{i\gamma x}}{\gamma };\hfill & \lambda =1\text{}u\mathrm{}\hfill \\ \frac{i}{2}lim_x\mathrm{}\frac{\mathrm{e}^{i\gamma x}}{\gamma };\hfill & \lambda =1\text{}u\mathrm{}\hfill \end{array}$$
(38)
Of course, the limits occurring in these expressions do not exist. However, it should be remembered that this integral is being evaluated as part of the kernel of an integral equation. Therefore, the limits can (and have to) be taken in the sense of a distribution limit, $`lim_{x\pm \mathrm{}}\mathrm{e}^{i\gamma x}/\gamma =\pm i\pi \delta (\gamma )`$.

With these tools, one can go about evaluating the residues of the integrand, which can easily be done using $`\mathrm{res}(f,z)=lim_{uz}f(u)(uz)`$ and an application of de l’Hospital’s rule. We will give the end result below.

A few words about the other cases are in order. First, if $`\xi >1/\sqrt{2}`$, the poles at $`u_j^{(1)}`$ move away from the real axis (cf. Fig. 1). One can keep the expression resulting for $`\xi <1/\sqrt{2}`$ by choosing an appropriate branch of the logarithm in the result. Basically, one has to see to it that the $`u_j^{(1)}`$ in Fig. 1 “shift down” rather than “up” by $`i\pi `$. To ensure this, one sets $`\mathrm{ln}(1)=i\pi `$, and with $`\mathrm{ln}(1)=0`$ a natural branch cut is along the positive imaginary axis, so that $`\frac{i\pi }{2}\mathrm{ln}(\mathrm{e}^{i\varphi })>3\pi i/2`$.

If $`\gamma <0`$, one has to integrate through the lower half-plane, and consequently the poles on the real axis are added into the residue sum. The resulting expression is, of course, identical to the one for $`\gamma >0`$.

After carrying out all integrations and the summation of $`m`$ (which is trivial because the integration over $`w_2`$ contributes a $`\delta _{m,m^{}}`$, thus decoupling different angular harmonics), (27) takes the form

$$\mathrm{\Phi }_m^{}(\alpha ^{})=(\alpha ^{})\mathrm{\Phi }_m^{}(\alpha ^{})+_{\mathrm{}}^{\mathrm{}}𝒦(\alpha ,\alpha ^{})\mathrm{\Phi }_m^{}(\alpha )𝑑\alpha .$$
(39)
The summand with $``$ originates from the terms with $`\delta `$-functions and computes to

$`(\alpha ^{})`$ $`=`$ $`{\displaystyle \frac{\sqrt{12\sigma ^2}^{11/\sigma ^2}}{2\sqrt{2}\sigma ^2}}{\displaystyle \frac{\left|\mathrm{\Gamma }\left(\frac{1}{4}+\frac{1}{2}m+\frac{1}{2}i\alpha ^{}\right)\right|^2}{\left|\mathrm{\Gamma }\left(\frac{3}{4}+\frac{1}{2}m+\frac{1}{2}i\alpha ^{}\right)\right|^2}}(1`$ (42)
$`(1+\sigma ^2)\mathrm{e}^{\sigma ^2(\alpha ^2+2m^2)/2}(\mathrm{I}_0\left({\displaystyle \frac{\sigma ^2}{2}}(\alpha ^2+2m^2)\right)`$
$`+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2m^2}{m^22k^2}}\mathrm{I}_k\left({\displaystyle \frac{\sigma ^2}{2}}(\alpha _{}^{}{}_{}{}^{2}+2m_{}^{}{}_{}{}^{2})\right))).`$
The kernel itself becomes

$`𝒦`$ $`(\alpha ,\alpha ^{})=i{\displaystyle \frac{\omega ^{i(\alpha \alpha ^{})}\sqrt{12\sigma ^2}^{11/\sigma ^2}}{2\sigma ^2\mathrm{e}^{\sigma ^2(\alpha ^2+\alpha ^2+4m^2)/4+1}\left(\mathrm{e}^{2\pi (\alpha \alpha ^{})}1\right)}}`$ (46)
$`{\displaystyle \frac{\left|\mathrm{\Gamma }\left(\frac{1}{4}+\frac{1}{2}m+\frac{1}{2}i\alpha ^{}\right)\right|^2}{\left|\mathrm{\Gamma }\left(\frac{3}{4}+\frac{1}{2}m+\frac{1}{2}i\alpha ^{}\right)\right|^2}}(m^{i(\alpha \alpha ^{})}\sigma ^2\mathrm{I}_0(\tau )`$
$`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}\mathrm{I}_k(\tau )((\sqrt{2}km^{}\sigma ^2)\mathrm{\Lambda }^k\mathrm{e}^{(i(\alpha \alpha ^{})1)\mathrm{ln}(m^{}\sqrt{2}k)}`$
$`(\sqrt{2}k+m^{}\sigma ^2)\mathrm{\Lambda }^k\mathrm{e}^{(i(\alpha \alpha ^{})1)\mathrm{ln}(m^{}+\sqrt{2}k)}),`$
where we have abbreviated

$`\tau `$ $`=`$ $`{\displaystyle \frac{\sigma ^2}{2}}\sqrt{(\alpha ^2+2m^2)(\alpha ^2+2m^2)}`$ (47)
$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{\alpha \alpha ^{}+\sqrt{2}m^{}(\sqrt{2}m^{}+i(\alpha \alpha ^{}))}{\sqrt{(2m^2+\alpha ^2)(2m^2+\alpha ^2)}}}.`$ (48)
Formally, (39) is an integral equation of the second kind. However, $`𝒦`$ is not compact – the infinite integration range and the singularity on the $`\alpha =\alpha ^{}`$ diagonal prevent this –, and thus the solution theory for it is significantly different from what one has in the Fredholm case. We will return to this issue later. For now, let us just follow Zang in noting that $`\omega `$ factors out of the kernel, and thus the value of $`\omega `$ is irrelevant for the solvability of this equation.

The special case of $`\omega =0`$ has not been treated so far. It is interesting for two reasons: For one, when searching for the limit of axisymmetric stability, it suffices to consider $`\omega =0`$, since axisymmetric patterns cannot rotate anyway and considering growth is not necessary for establishing the stability margin, and for two, “neutral” modes (for which $`\omega =0`$ even when $`m0`$) were the only modes the existence of which ERII and Zang report for the self-consistent disk.

The expression (46) cannot be taken to the limit $`\omega 0`$, because $`\omega ^{i(\alpha \beta )}`$ in general has an essential singularity at $`\omega =0`$, which according to Picard’s theorem implies that the image of an arbitrarily small neighbourhood of 0 is the entire complex plane with the possible exception of one point. Instead, one has to substitute $`\omega =0`$ before the integration over $`L_z`$. This simplifies the terms quite significantly, resulting in

$$𝒦_0(\alpha ^{})=2(\alpha ^{})$$
(49)
The integral equation has now collapsed into the algebraic equation

$`\mathrm{\Phi }_m^{}(\alpha ^{})`$ $`=`$ $`𝒦_0(\alpha ^{})\mathrm{\Phi }_m^{}(\alpha ^{}),\text{or}`$ (50)
$`0`$ $`=`$ $`𝒦_0(\alpha ^{})1,`$ (51)
which clearly simplifies solving the equation by at least an order of magnitude.

## 4 The cut-out disk

The Mestel disk has three rather bizarre properties: Its mass is infinite, the surface density is not bounded in the centre, and it does not possess a natural time-scale. All these shortcomings can be cured by a relatively simple measure that was already employed by Zang. The idea is to immobilise a radially varying part of the mass, declaring it as, say, halo or bulge. Thus one keeps the simple structure of the disk, but cuts the singularity in the centre and possibly the long mass tail out to infinity as far as the Boltzmann equation is concerned. Also, the cut-out will break the self-similarity and thus introduces a time scale.

As an added benefit, the ratio of active to total mass and thus (at radially constant velocity dispersion) Toomre’s stability parameter $`Q`$ (Toomre toomreq (1964)) becomes a function of the radius, opening the possibility to introduce a $`Q`$-barrier in the disk. In such a $`Q`$-barrier, incoming waves are refracted away from the centre and thereby evade their absorption at the inner Lindblad resonance. In common semi-phenomenological theories of spiral structure (Bertin et. al. bertin (1989)), this $`Q`$-barrier serves as a part of a feedback loop, to be closed by, e.g., the corotation at which swing amplification (Toomre toomreswing (1981)) or related mechanisms would not only provide reflection but even amplification of wave packets.

In the shearing sheet (Goldreich & Lynden-Bell 1965MNRAS.130..125G (1965), Julian & Toomre 1966ApJ…146..810J (1966)), it is possible to show rigorously that the introduction of something quite similar to a $`Q`$-barrier leads to the emergence of stationary modes where no such modes are present in the classic border-less analyses (Fuchs fuchspre (2001)).

Technically, the immobilisation of the mass is done by multiplying the distribution function with a cut-out function $`H(r)`$, with $`0H(r)1`$, while keeping the Hamiltonian and thus the rotation curve unchanged. Since the main objective of our investigation of the cut-out disk was an assessment of the quality of our approximation, we used the cut-out function used by ERI (who in turn basically used what Zang had suggested). In our case, the introduction of the cut-out has the effect of eliminating the border terms in the $`L_z`$ integration, but only if one uses Zang’s doubly cut-out disk (i.e., the cut-out function tends to zero for both $`r0`$ and $`r\mathrm{}`$). Zang’s family of cut-out functions is parametrised by two positive integers $`M`$ and $`N`$, the cut-out indices, describing the steepness of the cut-outs in the outskirts and centre of the disk, respectively. Since both of these integers appear as the order of poles of the integrand in the $`L_z`$ integration, and the residues at these poles get very messy indeed as their order increases, we settled for the case $`M=N=2`$, resulting in a cut-out function

$$H(L_z)=\frac{L_z^2L_c^2}{(L_z^2+L_0^2)(L_z^2+L_c^2)},$$
(52)
where $`L_0`$ and $`L_c`$ are two parameters qualifying the locations of the inner and outer cut-outs; following ERII, we usually use $`L_0=1`$ and $`L_c=10`$.

It should be noted that for further analytic exploration of the cut-out disk, this cut-out is not favourable, since it still leads to a very cumbersome kernel. A cut-out of $`H(L_z)=2L_zL_0/(L_0^2+L_z^2)`$ results in a much more compact kernel and reproduces the overall behaviour the cut-out disk (due to its lower active surface density, the resulting disk has an even higher stability than one with Zang’s cut-out, though).

Since the cut-out function (52) only depends on $`L_z`$, most of the calculations done for the self-consistent disk just carry over. The integration over $`L_z`$ becomes somewhat easier, since there are no border terms, but this advantage is more than offset by the sheer length of the terms and the occurrence of higher-order poles. The resulting integral equation has an “almost” compact kernel, with no singularity on the diagonal and a relatively sharp decline of the kernel as $`\alpha `$ and $`\alpha ^{}`$ approach infinity. Consequently, ERI’s naive approach of simply Fredholm-discretizing the integral equation and cutting off at some appropriate $`\alpha `$ still works well. The only difference to ERI’s method of solving this integral equation was that we found a spline-based quadrature to converge faster than the integration scheme proposed by them.

Skipping the details – which can be found in Demleitner (diss (2000)) –, let us briefly compare the stability of disks described with our approximation with Read’s findings. Since ERII mostly report results for singly cut-out disks which are not favourable for an anaylsis using our tools, we give the values from Read’s (1997PhDT………1R (1997)) more comprehensive tables. Our disks become unstable to axisymmetric ($`m=0`$) perturbations at $`\sigma _{\mathrm{min}}=0.274`$; Read only gives values for $`L_c=100`$, where we find $`\sigma _{\mathrm{min}}=0.349`$ versus Read’s $`\sigma _{\mathrm{min}}=0.327`$. This agreement to within roughly 5% is remarkable considering the relatively high temperature of this disk. The higher stability of the disk with lower $`L_c`$ is mostly due to its lower active surface density.

The one-armed ($`m=1`$) modes investigated in ERII behave somewhat erratic in that at vanishing growth rate and $`\mathrm{}\omega 0`$ in that there seem to be modes for almost arbitrary $`\sigma `$. That, however, is a consequence of the fact that ERII performed these computations without an outer cut-out. Introducing an outer cut-out, the disk seems to become stable at $`\sigma =0.57`$, although clearly at that point the assumption of a razor-thin disk becomes quite unrealistic and $`v_c\sigma `$ is severely violated.

In the – at least historically – most interesting case of two-armed ($`m=2`$) modes, our disks become unstable at $`\sigma _{\mathrm{min}}=0.214`$ with a pattern speed $`\mathrm{\Omega }_p=\omega /m=0.548`$, whereas Read finds $`\sigma _{\mathrm{min}}=0.205`$ and $`\mathrm{\Omega }_p=\omega /m=0.542`$. Thus, not only is the qualitative result that the $`M=N=2`$ cut-out disk is stable with respect to two-armed perturbations reproduced, the quantitative agreement again is far better than one might expect given the rather severe assumptions that go into our model.

Read (1997PhDT………1R (1997)) does not give stability limits for the $`m=3`$ and $`m=4`$ perturbations in doubly cut-out disks, but a comparison of the eigenvalue plots leads to the expectation that our results – $`\sigma _{\mathrm{min}}=0.141`$ at $`\mathrm{\Omega }_p=0.61`$ for $`m=3`$ and $`\sigma _{\mathrm{min}}=0.124`$ at $`\mathrm{\Omega }_p=0.588`$ for $`m=4`$ – would again reproduce the results of the more exact calculations quite well. The $`m=4`$ stability limit has to be taken with a grain of salt, since our formalism reproduces modes at high velocity dispersion as $`\omega 0`$ in accordance with Read’s findings. It is remarkable, though, that these modes are not quenched by an outer cut-out like their $`m=0`$ brethren.

The properties of the modes computed from the kernels are, by and large, identical to those discussed by ERII, and our model even shows more subtle features like the rapid overtaking of eigenvalues when $`m=4`$.

As to comparison to simulations, the fact that the expressions get very unwieldy for larger $`M`$ and $`N`$ is somewhat unfortunate. Sellwood & Evans (2001) point out that for the global mode to dominate over particle noise in simulations, it is highly desirable to maintain a high amplification, implying steep cut-outs and thus relatively large cut-out indices. The values chosen by Sellwood & Evans for their “standard” disk, $`N=4`$ and $`M=6`$, are beyond the reasonable range of applicability for our approach. However, they compute a doubly cut-out model that does not possess any global modes at $`Q=1`$ (remember that $`Q`$ refers to the axisymmetric stability of the self-consistent disk) and state that the $`N=1`$ and $`M=2`$ disk with $`L_0=1`$ becomes stable with respect to $`m=1`$ perturbations below $`L_c5`$. The rotation curve they use in this model is falling but only very slightly so, so that a disk with a completely flat rotation curve should behave analogously. Indeed, ignoring the “rogue” eigenvalues (Demleitner 2000) – that become relevant in this model –, we find that even the $`m=1`$ mode gets quenched below $`L_c=4.7`$.

Summarising, our approximation does well in reproducing the results from the previous calculations for the cut-out disks. For now, the value of this insight is rather limited, since one still has to take recourse to numerical integration to solve the integral equation, and the numerics are computationally about as expensive as those in Zang’s formalism, since the expression for the kernel is very large. Still, using a cut-out more suitable to our formalism results in more compact and manageable expressions that might facilitate a closer look at the interplay between the $`Q`$-barrier, the various resonances, and density waves in a full-disk analytical model.

## 5 The self-consistent disk

The axisymmetric stability of the self-consistent disk can be determined by solving the algebraic equation (51) with a kernel still somewhat simpler than (49) that is the result of setting $`m=0`$ in (29) and (32), yielding $`\sigma _{\mathrm{min}}=0.3860`$ with $`\alpha =3.51`$, within 2% of Read’s (1997PhDT………1R (1997)) values of $`\sigma _{\mathrm{min}}=0.3780`$ and $`\alpha =3.46`$.

Since the integral equation describing the self-consistent disk for non-axisymmetric modes is highly singular, Zang’s formalism could not provide solutions for the full problem. However, it did yield a family of rather exotic modes characterised by $`\omega =0`$ and consisting of a single logarithmic spiral – the neutral modes. It should be noted that these single logarithmic spirals are scale-free, a property that sets them apart from from any other disturbance of the form (13). In our formalism, they can be computed by solving the algebraic equation (51) with the kernel (49). The resulting behaviour again agrees very well with the findings of ERII and Zang and has to be called rather peculiar: One-armed perturbations cannot be stabilised at all, two-armed modes set in below $`\sigma =0.155`$ (or $`\sigma =0.149`$ according to Read (1997PhDT………1R (1997))), three-armed ones are extremely stable and do not set in down to at least $`\sigma =0.05`$, whereas four-armed neutral modes again cannot be stabilised at all.

The neutral modes come in pairs, i.e., if a mode with logarithmic wave number $`\alpha `$ exists, so does one with $`\alpha `$. This is a consequence of the anti-spiral theorem (Lynden-Bell & Ostriker antispir (1967)), since these neutral modes are not limits of growing modes. The marginal modes in the cut-out disk and possibly existing modes with nonzero $`\omega `$ in the self-consistent disk, on the other hand, are limits of growing modes by the application of Landau’s rule and thus not subject to the anti-spiral theorem.

ERII, from their perspective of analysing the entire family of physically viable power-law disks, briefly hypothesised that the stability properties might be related to the emergence of closed orbits of a given symmetry; these would, according to her reasoning, absorb power from perturbations and stabilise the disks when present. In our model, no closed orbits at all exist, since the ratio of epicyclic and orbital frequencies is fixed and irrational. As the model reproduces the stability properties found by ERII nevertheless, one has to look for a different explanation.

Turning to (49), one sees that the contribution from the terms except the one with the sum is fairly independent of $`m`$, whereas the sum does significantly change its value with $`m`$. The term to watch here is the $`m^2/(m^22k^2)`$. When $`m^{}\sqrt{2}k`$, the corresponding summand is much larger than its neighbours since the ratio $`I_k(x)`$ and $`I_{k+1}(x)`$ basically is $`k`$ (plus some constant) and everything else does not depend on $`k`$; when $`m`$ is not too large, that “resonant” summand will dominate the value of the entire series. The sign of this resonant summand is determined by the side relative to the nearest $`k`$ on which $`m^{}\sqrt{2}k`$ reaches its minimum – for the $`m=3`$ series, the near-resonance is at $`k=2`$ and thus one has a positive contribution with $`m^{}\sqrt{2}k=0.17`$ that enters into the kernel negatively (enhancing stability), whereas the resonance for the $`m=4`$ occurs for $`m^{}\sqrt{2}k=0.24`$ and thus destabilises the disk. The next azimuthal harmonic that has an $`m`$ close to a $`\sqrt{2}k`$ is 7, and it is again negative; not surprisingly, the $`m=7`$ harmonic is quite unstable as well.

Unfortunately, $`k`$ has no immediate meaning in physical terms, although (22), suggests that the $`k`$-th terms of the series in (49) should determine the response in the $`k`$-th harmonic in the epicyclic angle $`w_1`$. In the denominator of the resonance term $`m^22k^2`$, the factor $`\sqrt{2}`$ is simply the epicyclic frequency, so that we seem to see a resonance between the frequency with which a body in the disk encounters a forcing, $`m\mathrm{\Omega }`$, and a harmonic of the epicyclic frequency, $`k\kappa `$, which is a natural frequency a body wants to respond with. If the forcing is near resonant with the response times the epicyclic frequency but has a slightly higher angular symmetry, the response is damped, whereas in the reverse case it is amplified, a behaviour common in forced oscillations. This resembles the situation at a Lindblad resonance, with the difference that for one, the resonance is not exact (if it were, the series would evaluate to infinity), and, for second, it occurs on the entire disk.

Let us now turn to the general case of rotating modes. As has been hinted at above, the main discomfort with the notion of modes in the singular disk is that in its governing equation $`\omega `$ factors out, so that without further devices a disk is either stable or becomes unstable at all growth rates and pattern speeds at once. This seems weird enough to dismiss the possibility of such modes. On the other hand, this would imply that even a completely cold disk would be stable to all non-axisymmetric perturbations, which again would be quite an exotic situation.

Quite recently, Goodman & Evans (goodev (1999)) have tackled the problem employing the Jeans equations. The basic results of their work are that the Mestel disk’s weirdness starts in its centre, and once one fixes boundary conditions there, the Mestel disk becomes quite tame. Without that boundary condition, however, the problem is not well-posed. By the requirement that the centre does not absorb or emit energy, they curb down the two-dimensional continuum of solutions to a spectrum of one-dimensional continua, in which the ratio between growth rate and pattern speed is fixed and only the phase of $`\omega `$ remains unknown. Even this ambiguity can be removed by fixing the phase shift during the reflection of a wave at the disk’s centre, thus reducing the spectrum to a countable set of discrete points; however, there is nothing in the equilibrium model that would allow the definition of this phase change.

With a closed expression for the kernel of the Mestel disk, one is in a position to reconsider the stellar-dynamical equivalent of Goodman & Evans’ problem using the tools of the theory of singular integral equations: Does the integral equation (39) admit continuous solutions, and if, for what disk parameters?

The kernel (46) governing the self-consistent disk has a Cauchy-type singularity on the diagonal, implying that the operator $`𝒦(\alpha ,\alpha ^{})\mathrm{\Phi }(\alpha )𝑑\alpha `$ is not compact, so that the familiar Fredholm theory cannot be applied. For the treatment of integral equations of this type, Muskhelishvili (musk (1953)) and his coworkers have developed a theory. A nice summary is found in the more modern work of Polyanin & Manzhirov (polman (1998)). The following discussion is based on these works.

The integral equation (39) has the form

$$0=a(\alpha ^{})\mathrm{\Phi }(\alpha ^{})+\frac{1}{i\pi }_L\frac{𝐊(\alpha ,\alpha ^{})}{\alpha \alpha ^{}}\mathrm{\Phi }(\alpha )𝑑\alpha ,$$
(53)
where $`𝐊`$ is at least Hölder continuous, $`L`$ is the real axis (integration over the real axis is henceforth implied for integral signs), and, of course, the integral is to be understood as a Cauchy principal value. In the present case, $`a(\alpha ^{})=(\alpha ^{})1`$ and $`𝐊(\alpha ,\alpha ^{})=i\pi (\alpha \alpha ^{})𝒦(\alpha ,\alpha ^{})`$. The multiplication with $`(\alpha \alpha ^{})`$ ensures that the kernel is bounded on the diagonal, since the single zero in the denominator stemming from $`\mathrm{exp}\left(2\pi (\alpha \alpha ^{})\right)1`$ now cancels out. The only part of the kernel that could jeopardise the condition of Hölder continuity is the square root, but since only $`\alpha ^2`$ is present in the arguments of the square roots, $`𝐊`$ does indeed satisfy a Hölder condition.

It is convenient to rewrite (53) to the equivalent form

$`0`$ $`=`$ $`a(\alpha ^{})\mathrm{\Phi }(\alpha ^{})+{\displaystyle \frac{1}{i\pi }}b(\alpha ^{}){\displaystyle \frac{\mathrm{\Phi }(\alpha )}{\alpha \alpha ^{}}𝑑\alpha }`$ (55)
$`+{\displaystyle \frac{1}{i\pi }}{\displaystyle 𝐊_r(\alpha ,\alpha ^{})\mathrm{\Phi }(\alpha )𝑑\alpha },`$
where

$`b(\alpha )`$ $`=`$ $`𝐊(\alpha ^{},\alpha ^{})\text{and}`$ (56)
$`𝐊_r(\alpha ,\alpha ^{})`$ $`=`$ $`{\displaystyle \frac{𝐊(\alpha ,\alpha ^{})𝐊(\alpha ^{},\alpha ^{})}{\alpha \alpha ^{}}}.`$ (57)
After this manipulation, $`𝐊_r`$ is a compact operator, whereas

$$\widehat{𝐊}\mathrm{\Phi }:=a(\alpha ^{})\mathrm{\Phi }(\alpha ^{})+\frac{1}{i\pi }b(\alpha ^{})\frac{\mathrm{\Phi }(\alpha )}{\alpha \alpha ^{}}𝑑\alpha =0,$$
(58)
known as the dominant or characteristic part of the integral equation, decides whether or not the integral equation has a nontrivial solution. For our equation, $`b(\alpha ^{})=(\alpha ^{})`$.

To find a solution of (58), let us introduce an auxiliary function

$$\mathrm{\Psi }(z)=\frac{1}{2\pi i}\frac{\mathrm{\Phi }_m(\alpha )}{(\alpha z)}𝑑\alpha ,$$
(59)
defined for complex $`z`$. If $`\mathrm{\Phi }_m(\alpha )`$ is continuous – this is henceforth assumed –, this function is piecewise analytic in the upper and lower half planes but may be discontinuous across the real axis. The merit of this auxiliary function is that it allows a transformation of the problem of the solution of the integral equation (58) into a Riemann boundary value problem.

This can be done by applying the Sokhotski-Plemelj formulae, that, in their formulation for the real axis, state that when $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ are linked by (59),

$`\mathrm{\Psi }^+(\alpha ^{})+\mathrm{\Psi }^{}(\alpha ^{})`$ $`=`$ $`{\displaystyle \frac{1}{i\pi }}{\displaystyle \frac{\mathrm{\Phi }(\alpha )}{\alpha \alpha ^{}}𝑑\alpha }`$ (60)
$`\mathrm{\Psi }^+(\alpha ^{})\mathrm{\Psi }^{}(\alpha ^{})`$ $`=`$ $`\mathrm{\Phi }(\alpha ^{})`$ (61)
hold for all $`\alpha ^{}`$. Here, $`\mathrm{\Psi }^+(\alpha ^{})`$ and $`\mathrm{\Psi }^{}(\alpha ^{})`$ are the limiting values of $`\mathrm{\Psi }(z)`$ as $`z\alpha ^{}`$ from above and below the real axis, respectively.

Inserting (61) into (58) and collecting the terms yields

$$\mathrm{\Psi }^+(\alpha ^{})=\frac{a(\alpha ^{})b(\alpha ^{})}{a(\alpha ^{})+b(\alpha ^{})}\mathrm{\Psi }^{}(\alpha ^{})=:D(\alpha ^{})\mathrm{\Psi }^{}(\alpha ^{}).$$
(62)
Going back to the integral equation for the singular disk, one has

$$D(\alpha ^{})=1/(2(\alpha ^{})1).$$
(63)
This represents a Riemann boundary value problem, asking (in the special case relevant here) for two functions $`\mathrm{\Psi }^+`$ and $`\mathrm{\Psi }^{}`$ that are analytic on the upper and lower half planes, respectively, and that satisfy $`\mathrm{\Psi }^+(\alpha ^{})=D(\alpha ^{})\mathrm{\Psi }^{}(\alpha ^{})`$ for all real $`\alpha ^{}`$.

For reasons that will become clear below, we for the time being demand that $`D`$ has neither zeroes nor poles on the real axis. Comparing (51) and (63), one sees that this will be the case exactly when neutral modes are admitted by the disk in question, which strongly suggests that (39) has no solution in the regular case of a bounded and nonzero $`D`$, and at least one solution otherwise. With this concept mind, let us return to the analysis of the boundary value problem, for now keeping our attention on the case of a regular $`D`$.

The $`D(\alpha ^{})`$ treated here is the restriction of an analytic function (in particular, by assumption it has no poles on the real axis), and thus is guaranteed to be Hölder continuous. Furthermore, $`D(\alpha ^{})`$ is a positive real function for real $`\alpha ^{}`$. This ensures that its index, the increment of its argument over the integration contour,

$$\mathrm{Ind}(D)=\frac{1}{2\pi }d\mathrm{arg}D(\alpha ),$$
(64)
is zero. This in turn implies that $`\mathrm{ln}D`$ is a well-defined (single-valued) function. Then one can take the logarithm of both sides of (62), to arrive at

$$\mathrm{ln}\mathrm{\Psi }^+(\alpha ^{})\mathrm{ln}\mathrm{\Psi }^{}(\alpha ^{})=\mathrm{ln}D(\alpha ^{}).$$
(65)
Using the second of the Sokhotski-Plemelj formulae (61), it follows that with

$$G(z)=\frac{1}{2\pi i}\frac{\mathrm{ln}D(\alpha )}{\alpha z}𝑑\alpha $$
(66)
one can write down a solution to the boundary value problem (62)

$$X^+(\alpha ^{})=\mathrm{e}^{G^+(\alpha ^{})}X^{}(\alpha ^{})=\mathrm{e}^{G^{}(\alpha ^{})},$$
(67)
provided that the problem does have a solution.

Thus one has $`D=X^+/X^{}`$, and the boundary value problem (62) can be written as

$$\frac{\mathrm{\Psi }^+(\alpha ^{})}{X^+(\alpha ^{})}=\frac{\mathrm{\Psi }^{}(\alpha ^{})}{X^{}(\alpha ^{})}.$$
(68)
Now, the functions $`G^\pm `$ by construction have no pole in their respective domains, and therefore $`X^\pm `$ have no zeroes. If $`\mathrm{\Phi }^\pm `$ are defined via an integral as in (59), this implies that they are analytic in their respective half-planes. Thus, the functions on the two sides of (68) are analytic in their respective half-planes, and they are identical on the real axis. By Morera’s theorem this implies that the right-hand side and the left hand side represent the same analytic function, $`S(z)`$, say. Now, what are the properties of $`S(z)`$?

We already pointed out that $`S(z)`$ is analytic on the entire complex plane. At infinity, one has $`\mathrm{\Psi }^\pm (\mathrm{})=\pm \frac{1}{2}\mathrm{\Phi }(\mathrm{})`$. Since we require any solution of the original integral equation (39) to vanish at infinity, one has $`\mathrm{\Psi }^\pm (\mathrm{})=0`$ as well. Since $`D1`$ for large $`\alpha `$, $`G`$ approaches 0 and $`X^\pm `$ tends to 1. This implies that $`S`$ vanishes at infinity. Since $`S`$ is an entire function, it follows from Liouville’s theorem that $`S`$ must be identically zero. By the first Sokhotski-Plemelj formula one now has $`\mathrm{\Phi }(\alpha ^{})=S(\alpha )(X^++X^{})`$, and hence the homogenous integral equation has no nontrivial solution. This is a special case of the general result that is proven somewhat more rigorously in the literature cited above: A homogenous Cauchy integral equation with an index of zero is unsolvable.

Loosening the conditions on $`D`$, we now allow it to have “poles” at some points $`\alpha _i`$, which, for the sake of simplicity, are supposed to be poles of order one. The term “pole” has been put in quotes since it is not necessary that $`D`$ is analytic in $`\{\alpha _i\}`$, only that $`D(\alpha _i)`$ behaves like $`(\alpha ^{}\alpha _i)^1`$ at the exceptional points.

Writing

$$D(\alpha ^{})=\frac{D^{}(\alpha ^{})}{_{i=1}^n(\alpha ^{}\alpha _i)},$$
(69)
the boundary condition (62) takes the form

$$\mathrm{\Psi }^+(\alpha ^{})=\frac{D^{}(\alpha ^{})}{_{i=1}^n(\alpha ^{}\alpha _i)}\mathrm{\Psi }^{}(\alpha ^{}),$$
(70)
where $`D^{}`$ now satisfies the conditions we required of $`D`$ in (62). Solving the boundary value problem for this $`D^{}`$ and substituting the result back into (70), the equivalent of (68) now is

$$\frac{\mathrm{\Psi }^+(\alpha ^{})}{X^+(\alpha ^{})}=\frac{\mathrm{\Psi }^{}(\alpha ^{})}{X^{}(\alpha ^{})_{i=1}^n(\alpha ^{}\alpha _i)},$$
(71)
where $`X`$ and $`G`$ are again defined by (67) and (66), except that now $`D`$ is replaced by $`D^{}`$. Let us again call the analytic function defined in the left and right sides of (71) $`S(z)`$.

By the right side of (71), $`S`$ now may have poles of order 1 at $`\alpha _i`$ but still has to vanish at infinity. Applying the generalised Liouville theorem, it follows that any solution of the integral equation must be of the form

$$\mathrm{\Phi }(\alpha )=\frac{c_i}{\alpha \alpha _i}\left(X^+(\alpha )+X^{}(\alpha )(\alpha \alpha _i)\right).$$
(72)
Evidently, solutions of this type cannot vanish at $`\mathrm{}`$, since such solutions will at best converge to $`X^{}(\mathrm{})`$ times a nonzero constant as $`\alpha ^{}\mathrm{}`$.

Thus, provided that all singularities of $`S(z)`$ are poles<sup>1</sup><sup>1</sup>1This assumption is hard to prove, but generally true, which is why Pipkin (pipkin (1991)) calls it Gentlemen’s Theorem Number 1 – a theorem no gentleman would question., the self-consistent disk admits no growing modes at all. The obvious explanation for this bizarre absence of growing modes even in very cold disks is that growing modes would introduce scales into the scale-free disk.

Stepping back to the mathematics for a moment, the fact that $``$ is purely real on the real axis is crucial for that result, since even a minute complex contribution could provide for modes. For example, adding an $`iϵ`$ to $``$ would not change the result that the index (64) of the problem (62) vanishes when $`(\alpha ^{})<1/2`$ for all $`\alpha ^{}`$, so that there still would be no modes above the stability limit of the neutral modes. Below that stability limit, however, the index would be nonzero and the integral equation would admit nontrivial solutions.

Such an imaginary contribution would necessitate that the real part of the kernel (46) shows a singularity on the diagonal. We cannot offer an explanation how this might be linked to the scale-freeness of the disk at this point.

## 6 Conclusions

In this paper we have derived an analytic expression for the kernel of the Mestel disk. In the cut-out disk, it reproduces the behaviour described by Zang and ERII both qualitatively and quantitatively surprisingly well. It turns out that the expressions for the kernel are quite large in the cut-out case, at least when one uses the classic cut-out functions suggested by Zang, so the computational effort involved in finding modes is about as large with our formalism as it is with the older seminumerical scheme. On the other hand, it seems that with carefully selected cut-out functions, the kernel could be simple enough for further analytic work.

The kernel governing modes in the self-consistent disk is relatively compact. It allowed us to interpret the stability behaviour of neutral modes in terms of global resonances between the orbital (which coincides with the frequency of excitation in this case) and epicyclic frequencies. In the longstanding question of rotating modes we could establish that no rotating modes can exist above the stability limit of the neutral modes. It might seem surprising that rotating modes should know about the stability limit of the rotating modes, considering that the ratio of orbital and epicyclic frequency is no longer independent of the radius for rotating modes. However, due to the lack of a length scale in the Mestel disk, any mode is arbitrarily close to a non-rotating one.

Even below the stability limit of the neutral modes, no rotating or growing modes exist in the perfectly scale-free disk under reasonable assumptions like continuity of the coefficient function in the expansion of the potential perturbation (13). This bizarre property is probably best explained by an inability of the disk to break its self-similarity and thus by a completely artificial feature. In conclusion, it seems to us that the self-consistent disk’s peculiar properties make its investigation something of a pedagogical (Goodman & Evans goodev (1999)) exercise, whereas the cut-out disk might yet prove to be a valuable testbed for wave dynamics in a fairly realistic full-disk setting.

###### Acknowledgements.

Part of this work was supported by a grant of the Land Baden-Württemberg. M. D. wishes to thank Jenny Read for virtual but nonetheless useful discussions. An anonymous referee prompted us to clarify the summary of Goodman & Evans (goodev (1999)).


text:
# Penetration Depth and Anisotropy in MgB2

## Abstract

The penetration depth $`\lambda `$ of MgB2 was deduced from both the ac susceptibility $`\chi `$ and the magnetization $`M(H)`$ of sorted powders. The good agreement between the two sets of data without geometric correction for the grain orientation suggests that MgB2 is an isotropic superconductor.

Great interest has been raised recently by the discovery of MgB2 with Tc about $`40`$K. In its normal state, the compound appears to be a metal with a low dc resistivity, Hall coefficient,  and thermoelectric power dominated by hole carriers. The estimated long mean-free path implies that the electrical transport could well be isotropic in spite of its layer-like crystalline structure, in agreement with the band-structure calculations. Below a transition temperature of Tc $``$ 40 K, the compound seems to be a phonon-mediated BCS superconductor, as suggested by the large isotopic effect and the large negative pressure effect on Tc. The upper critical field, Hc2(T), has been directly measured above 20 K with a linear (or even an upward curvature) $`T`$-dependence. The lower critical field was either indirectly estimated based on the thermodynamics field Hc(T) above 34 K or calculated based on the non-linearity in $`M(H)`$ at low temperatures. The coherence length $`\xi _\text{o}`$(0 K), penetration depth $`\lambda `$(0 K), and Ginzburg-Landau parameter $`\kappa `$ were consequentially estimated to be 5.2 nm, 125-140 nm, and 26, respectively. All these calculations, however, were based on the assumption that MgB2 can be treated as isotropic. Otherwise, a geometric factor up to 2 would be needed to correct the random grain-orientation and to convert $`M_r`$ to $`H_c`$.

To estimate the anisotropy, $`\lambda `$ of MgB2 was directly measured using both $`ac`$ susceptibility and the non-linearity in the $`M(H)`$ of powder samples. Although the two procedures involve the anisotropy of $`\lambda `$ in very different ways, the $`\lambda `$(5 K)’s deduced from these procedures are in good agreement without geometric corrections for superconducting anisotropy. The results, therefore, suggest that the anisotropy of MgB2 is very small.

Ceramic MgB2 samples were prepared using the solid-state reaction method. Small Mg chips (99.8% pure) and B powder (99.7% pure) with a stoichiometry of Mg:B = 1:2 were sealed inside a Ta tube under an Ar atomosphere. The sealed Ta ampoule was then enclosed in a quartz tube. The assembly was heated slowly up to 950 $`{}_{}{}^{\text{o}}C`$ and was kept at this temperature for 2 hours, followed by furnace cooling. The structure was detemined by X-ray powder diffraction (XRD) using a Rigaku DMAX-IIIB diffractometer. Powder samples were prepared by sorting the pulverized powder using either sieves or the method of descending speed of the particles in acetone. No grain alignment was attempted. The grain morphology as well as the particle sizes of the powders were measured using a JEOL JSM 6400 scanning electron microscope (SEM). Magnetizations were measured in a Quantum-Design 5 T SQUID magnetometer with an ac attachement. The ac susceptibility was measured under a fixed frequency f = 413.1 Hz and an amplitude of 3 Oe.

The XRD pattern can be indexed as a hexagonal cell with lattice parameters a = 3.08 Å and c = 3.52 Å. A sharp superconducting transition was observed in both resistivity and $`ac`$ susceptibility with Tc $``$ 38 K.

The transition temperature in different fields, $`T_c(H)`$, was determined from both the dc magnetization and the ac suscepitibility measurements, which should be equivalent to a local resistivity measurement, with a dc bias of 0 - 5 T. The transition width in the ac suscepitibility is only slightly broadened under fields, i.e. from $``$ 1 K at the bias field of $`H_{dc}`$ = 0 T to $``$ 5 K at 5 T (Inset, Fig. 1). We attribute this to the flux movement under fields and take the onset temperature as $`T_c(H)`$ (Fig. 1). The data from the two methods are reasonably consistent with a slope $`dH_{c2}/dT`$ $``$ 0.4$`\pm 0.05`$ T/K near 38 K, in agreement with the results obtained by Finnemore et al.

The lower critical field was measured in powder samples with a particle size $`<`$$`\mu `$m to avoid the complications caused by intergrain coupling. No systematic variation between the deduced $`\lambda `$(T) and the particle sizes was observed. This demonstrates that the grain-boundary effect in the sample is small, a fact that is also supported by direct SEM observation of the powders.

It is known that the $`ac`$ susceptibility $`\chi `$ of a superconducting sphere of diameter $`d`$ is $`\mathrm{\Phi }(d,\lambda )`$ = $`3/(8\pi )[16(\lambda /d)coth(d/2\lambda )+12(\lambda /d)^2]`$, which reduces to $``$ $`(d/\lambda )^2`$ when $`d<2\lambda `$. Methods have been previously developed to deduce $`\lambda `$ from the $`\chi `$ of magnetically aligned powders by solving the equation $`\chi `$ = $`\mathrm{\Phi }(d,\lambda )`$ $``$ -0.002$`(d/\lambda )^2`$ in the large-$`\lambda `$ limit. However, the $`\lambda `$ so deduced will be sensitive to the uncertainty of $`\chi `$ (due to either the demagnetization factor or the superconducting volume fraction) if $`d`$ $`>`$$`>`$ 2$`\lambda `$ as in unsorted powders, which may include particles as large as 3 $`\mu `$m. In fact, the calculated $`\mathrm{ln}\lambda /\mathrm{ln}\chi `$ varies with $`d/\lambda `$ only moderately below $`d/\lambda `$ = 4, i.e from -0.5 to -0.7, but changes rapidly for larger $`d/\lambda `$. For example, $`\mathrm{ln}\lambda /\mathrm{ln}\chi `$ is -3 at $`d/\lambda `$ = 20, and a 30% uncertainty of $`\chi `$ will lead to a $`\lambda `$ anywhere between 0 and 10$`d`$. A $`d/\lambda `$ $`<`$5 is needed to obtain a 20% accuracy with an estimated 30% uncertainty in $`\chi `$. The technique can be improved by using sorted powders, which have a smaller $`d`$ and a narrower size-distribution. The powder obtained from pulverizing ceramic was thoughly mixed with acetone in a 10 ml beaker. The particles was then sorted according to the time needed for them to reach the bottom of the beaker. The sample discussed here was collection of particles deposited between 1-2 hr. Our SEM observation suggested that 99% or more particles having a size between 0.1 and 2 $`\mu `$m (Fig. 2). It should also be noted that the proposed method of calculating the effective grain-size d = $`\sqrt{(\text{d}_i^5)/(\text{d}_j^3)}`$ (where di is the diameter of individual grains) may also be questionable if $`d`$ $`>`$$`2\lambda `$. As will be shown below, a 30% error may be cuased by the approaximation alone. A regression procedure, therefore, was adopted. A $`\lambda _{raw}`$ corresponding to the draw = $`\sqrt{(\text{d}_i^5)/(\text{d}_j^3)}`$ was used as the initial value. The $`d`$ was then refined regressively as $`[\mathrm{\Phi }(\lambda ,d_i)`$d$`{}_{i}{}^{4}]/[\mathrm{\Phi }(\lambda ,d_i)`$d$`{}_{i}{}^{3}]`$. The convergence is very fast. It should be noted that the effective $`d`$ depends on $`\lambda `$, i.e. the correction varies with T and may change the T-dependence of $`\lambda `$. Our tests on YBa2Cu3O7-δ powders demonstrated that the uncertainty of the $`\lambda `$ so deduced is within 10-20% of the published data if $`d/\lambda `$ is 3 or smaller.

The SEM photo of a powder sample is shown in Fig. 2. Its $`d`$-distribution is shown in the inset of Fig. 3. A d $``$ 0.88 $`\mu `$m was obtained with draw = 1.23 $`\mu `$m. The superfluid density 1/$`\lambda ^2(T)`$ was calculated assuming an isotropic superconductivity, i.e. without corrections for the random grain orientation (Fig. 3). The 1/$`\lambda ^2(T)`$ observed roughly follows a T-dependency of \[1-($`T`$/39.4)2.7\]. It should be noted that a deviated from the fit is clear below 10 K. We are hesitate, however, to draw any conclusion about the deviation since weak-links can not be conclusively excluded at this stage. The extrapolated $`\lambda `$(0) $``$ 180 nm is slightly longer than that found by Finnemore et al.

The lower critical field $`H_{c1}`$ was also deduced as the field, where the linear $`MH`$ correlation begins to be violated, from the magnetizations of the same powder sample in an $`H`$-increase branch at 5 K (Inset, Fig. 4). Several technical difficulties in this method have been previously discussed: for instance, the intergrain coupling that may cause nonlinearity far below the intragrain $`H_{C1}`$; the surface pinning that can make the $`H_{C1}`$ value observed higher; sharp local edges, i.e. strong local demagnetizing fields, that can lead to a lower one; and the experimental resolution of the nonlinearity. Several precautions have been taken. The powder sample used here has a particle size far smaller than the average grain size, which should eliminate the effect of the intergrain coupling, as suggested by the smooth and flat $`\chi `$ observed below Tc. To improve the sensitivity of nonlinearity, the $`M(H`$, 5 K) below 50 Oe was fit as a linear function of $`H`$ using a standard least-square procedure. The deviation $`\mathrm{\Delta }`$M ($``$ 0.002 emu/cm3 below 50 Oe) from the linear fit is comparable with the experimental uncertainties in both $`M`$ and $`H`$, demonstrating the negligible effect of the residual granularity. The difference between the data and the extrapolated linear fit above 50 Oe was then calculated. The uncertainty associated with the linear fitting was marked as dashed lines in Fig. 4. To further avoid the interferences from the sharp edges of the particles, the deviation at large fields was empirically fit as $`a(HH_\text{o})^{1.8}`$ with both $`a`$ and $`H_\text{o}`$ as the free parameters (Fig. 4). We justify the fit by pointing out that the magnetization of a superconductor partially penetrated by external fields will vary as the square of the thickness penetrated, i.e $`HH_\text{o}`$ in the Bean model. It should be pointed out that the value of $`H_\text{o}`$ is not very sensitive to the index chosen. A linear fit below 200 Oe leads to only 15% change. Surface pinning is usually negligible in randomly shaped grains, and typically can only make the Hc1 observed smaller. Following the procedure, an Hc1 = $`H_\text{o}/(1g)`$ 130 Oe was obtained, where $`g`$ = 1/3 is the demagnetization factor of a sphere. This value consequently leads to a $`\lambda `$ of 203 nm, in good agreement with that from the ac $`\chi `$ within the uncertainty of the techniques. No corrections have been made to consider random grain orientations.

To further verify the result, the ac susceptibility of the same powder sample was measured at 5 K with a dc bias between 0 and 200 Oe (Fig. 5). A change of the slope was observed around H$`_\text{o}`$ $``$ 110 Oe, and the Hc1 was estimated as $``$ 160 Oe. Similar measurements have been done in several different samples and the results appear to be independent of the particle sizes.

The deduced $`\lambda `$ is slightly longer than the 140 nm and the $``$ 130 nm previously reported. The exact reason for the disagreement is not clear to us at this moment. However, the $`\lambda `$ measured here using three different methods on the same sample are self-consistent within the estimated experimental uncertainty of $`\pm `$20%.

It is interesting to note the agreement between the Hc1 from $`\chi `$ and $`M(H)`$. In a highly anisotropic layered superconductor, cuprates for example, the observed $`\chi `$ will only come from the supercurrents in the layers. The $`\chi `$ observed, therefore, should be assumed to be $`\mathrm{cos}^2\theta \mathrm{sin}\theta d\theta /\mathrm{sin}\theta d\theta `$ $``$ 1/3 of the $`3/(8\pi )[16(\lambda /d)coth(d/2\lambda )+12(\lambda /d)^2]`$. The deduced $`\lambda `$ will be 1.7 times longer if no geometric correction has been made. In general, the $`1`$/$`\lambda ^2`$ deduced from non-grain-aligned powder should be 1/3$`\lambda _{ab}^2`$+2/3$`\lambda _c^2`$ = (1/3+2/3$`\gamma )/\lambda _{ab}^2`$ in layered superconductors, where $`\lambda _{ab}`$, $`\lambda _c`$, and $`\gamma `$ are the penetration depths in and out of the layers, and the anisotropy, respectively. The Hc1 deduced from the nonlinearity of $`MH`$, on the other hand, can be even larger since the effective field perpendicular to the local layers is only a fraction, i.e. cos$`\theta `$, of the external field. The good agreement observed here, therefore, strongly suggests that MgB2 is an isotropic superconductor. The estimated $`\gamma `$ should be smaller than 1.5 assuming a experimental uncertainty of $`\pm `$20% in our $`\lambda `$-calculation. This $`\gamma `$ is far smaller than that of 10-1000 observed in various cuprates, and should be regarded as essentially isotropic.

In above data analysis, a spherical shape was assumed. The relative change of the demagnetization factor is $`\mathrm{\Delta }`$r/3r in a slightly-deformed ellipsoid with radii r+$`\mathrm{\Delta }`$r, r+$`\mathrm{\Delta }`$r/2, and r+$`\mathrm{\Delta }`$r/2. A simple calculation shows that the correction of the Hc1 will be -0.25($`\mathrm{\Delta }`$r/r)2 in the ac $`\chi `$ method, but $`\mathrm{\Delta }`$r/3r in the nonlinearity method. An average length ratio (r+$`\mathrm{\Delta }`$r)/(r-$`\mathrm{\Delta }`$r) between 0.5 to 2 (i.e. $`\mathrm{\Delta }`$r/3r$`=1/3`$), therefore, may not significantly changes above conclusion. The condition seem to be satisfied (Fig. 2).

This conclusion is in agreement with the band structure calculation, the extremely long mean-free path, the long coherence length, and the small grain-boundary effect on the supercurrents reported.

In summary, the penetration depth $`\lambda `$(T) of MgB2 was deduced from both the $`ac`$ susceptibility $`\chi `$ of powders and the nonlinearity of the $`MH`$ in the $`H`$-increase branch. The good agreement between the two methods suggests that MgB2 is an isotropic superconductor.

This work is supported in part by the NSF, the T. L. L. Temple Foundation, the John and Rebecca Moores Endowment and the State of Texas through TCSUH, and at LBNL by DOE.