text:
# Quick and dirty methods for studying black-hole resonances.

## I Introduction

During the last few decades it has become clear that resonant phenomena play an important role in black hole physics novikov . Numerical experiments, pioneered by Vishveshwara vishu  more than thirty years ago, and continued since then by many other authors, clearly show the dominance of the so-called quasinormal mode (QNM) oscillations in most dynamical processes involving black holes. The importance of this phenomenon is further emphasised by the possibility that future observations of QNM “ringing” following, say, the merger of two black holes by the network of interferometric gravitational wave detectors may provide direct identification of black holes. This exciting prospect has motivated a large number of investigations into the nature of the QNMs. An exhaustive list of papers on the subject can be found in Refs. kokkotas ; nollert\_rev .

The spectrum of QNM frequencies has been studied extensively in the literature, especially for the Schwarzschild black hole. It is well-known that the actual mode calculation is a non-trivial task. Recall that any physically acceptable solution to the relevant (homogeneous) wave equation takes the form of a purely “ingoing” wave $`\mathrm{exp}(i\omega r_{})`$ (where $`r_{}`$ is the tortoise radial coordinate) at the event horizon and a mixture of “ingoing/outgoing” waves $`A_{\mathrm{in}}\mathrm{exp}(i\omega r_{})+A_{\mathrm{out}}\mathrm{exp}(i\omega r_{})`$ far away from the black hole (such solutions exist because the effective curvature potential falls off asymptotically towards both the horizon and infinity, i.e. when $`r_{}\pm \mathrm{}`$). A QNM solution corresponds to purely outgoing waves at infinity. In other words, to search for a QNM one must first calculate the asymptotic amplitudes $`A_{\mathrm{out}}`$ and $`A_{\mathrm{in}}`$, and then try to find a frequency that satisfies the required QNM boundary condition $`A_{\mathrm{in}}=0`$. However, direct numerical integration of the wave equation leads to serious numerical problems. Since the QNMs are expected to be damped due to gravitational-wave emission the corresponding solutions must grow exponentially as one approaches infinity and the horizon (on a $`t`$-constant hypersurface). In order to identify a QNM one must be able to filter out the solution that decays exponentially at infinity. In other words, the numerical calculation requires exponential precision.

Since the first “direct” calculation of the slowest damped Schwarszchild QNMs by Chandrasekhar and Detweiler chandra1  various alternative schemes have been devised. Detweiler detweiler1  was able to determine the fundamental Kerr QNM with adequate accuracy by examining the functional form of the asymptotic amplitudes along the real frequency axis. Ferrari and Mashhoon valeria  exploited the correspondence between resonant and bound states of the inverted potential to obtain approximate black hole modes. WKB techniques were successfully used to calculate the first few least damped modes by Schutz and Will schutz  and several subsequent authors iyer1 ; iyer2 ; iyer3 ; guinn . High precision QNMs for both Schwarzschild and Kerr black holes were first calculated by Leaver leaver1 ; leaver2  by means of a continued fraction approach. His results for the Schwarzschild case were later verified (and extended) by Andersson na2  via the phase-amplitude formalism and Nollert and Schmidt nollert1  who used a Laplace transform/analytic continuation method. Moreover, the semi-analytic phase-integral formalism froman  has proved to be a valuable tool for the calculation and interpretation of QNMs pi .

More recently, Onozawa onozawa1  extended Leaver’s calculation of the QNMs of rotating black holes and studied the behaviour of the mode frequencies as the rotation parameter is varied. The “long-lived” QNMs of a near extreme Kerr black hole were discussed by the present authors kg1 , using analytic expressions for the asymptotic amplitudes derived by Press and Teukolsky press . That study complemented a calculation by Detweiler who was the first to use these amplitudes to derive a QNM-condition for near extreme rotation, and moreover, solved it analytically for the extreme case detweiler1 . Finally, Berti and Kokkotas bk03  have used the continued fraction approach to study the behaviour of the high overtone QNMs of a Kerr black hole. Their study unveils a complicated, and quite puzzling, structure.

In the last few months there has been renewed interest in the spectrum of black hole QNMs. Motivated by the suggestion that the asymptotic behaviour of the high-overtone QNM frequencies noll ; nahigh  may be related to the so-called Barbero-Immirzi parameter which is relevant for loop quantum gravity dreyer ; motl , there has been several attempts to approximate the high-frequency part of the spectrum. Some progress has been achieved by linking the high-overtone spectrum with the black hole’s quantum properties hod . Although we will not discuss this issue further in the present article, it is an interesting problem that need to be investigated in complete detail.

The notion of Regge poles has (so far) hardly at all been discussed in the context of black holes. The Regge poles are close relatives to the QNMs jensen1 , and correspond to those complex values of the angular momentum $`l`$ for which the wavefunction is purely outgoing at infinity and purely ingoing at the horizon. The Regge pole calculation poses an easier problem than that of finding QNMs since the corresponding wavefunction is well behaved at infinity (and the horizon) for complex $`l`$ (as long as the frequency is kept real). In quantum theory, the complex angular momentum (CAM) formalism (in which Regge poles play a central role) provides an alternative (and usually more transparent) description of scattering problems nussenzveig . Some years ago, one of us applied the CAM theory to the study of scattering of massless scalar waves by a Schwarzschild black hole cam1 ; cam2 . In that study, Regge poles and the corresponding residues were calculated via the phase-integral method. The CAM theory has also been used by Chandrasekhar and Ferrari chandra2  in an investigation of axisymmetric pulsations of a relativistic star.

In this work we discuss numerical methods for calculating both QNMs and Regge poles. The key motivation for our study was the desire to devise methods that were easy to implement (hence “quick and dirty”) and yet able to produce numerically accurate results. To achieve this goal we express the Teukolsky wave function in terms of suitable phase-functions (see Section III below). The integration of the resulting first-order equations readily yields the desired ratio of asymptotic amplitudes. In order to retain sufficient numerical precision to be able to identify the ingoing wave at infinity for complex frequencies, while ensuring that the phase-functions remain smooth and well behaved, the integration is generally performed along straight lines in the complex $`r_{}`$ plane. The slope of the integration contour is chosen in such a way that the relative magnitude of the out- and ingoing wave solutions is preserved both near the horizon and at infinity. The end result is a multi-purpose scheme that (with minor modifications) facilitates the calculation of QNMs and Regge poles as well as scattering phase-shifts kg2 .

The remainder of the Paper is organized as follows: In Section II we briefly discuss the formal definition and physical interpretation of black hole resonances. Our numerical approaches to the problem are described in full detail in Section III. Section IV contains numerical results for QNMs (Section IVA) and Regge poles (Section IVB). Section V provides a summary and some final comments. Geometric units $`G=c=1`$ are adopted throughout the paper.

## II Black hole resonances: definition and physical meaning

The notion of scattering resonances is probably most familiar from quantum theory taylor . The fact that it is also relevant in black hole physics (in terms of QNMs and Regge poles) is not particularly surprising given that linear perturbations of black hole spacetimes are governed by Schrödinger type equations. The corresponding effective potential barriers are a combination of the black hole’s gravitational pull and a repulsive centrifugal term chandrabook . This means that one can formulate a scattering problem for black holes and investigate phenomena analogous to those in the quantum case. For simplicity we will mainly consider the case of scalar waves in this paper. This is not a severe restriction, however, since more “realistic” cases, eg. gravitational perturbations, can be dealt with in exactly the same way.

After adopting the Boyer-Lindquist coordinate frame, and separating the angular dependence in the standard way, we arrive at the Teukolsky equation teuk , which governs the radial component $`\psi (r)`$ of a scalar perturbation in a Kerr black hole spacetime,

$$\frac{d^2\psi }{dr_{}^2}+Q_{lm}(r_{},\omega ,a)\psi =0$$
(1)
We have also assumed a harmonic time dependence $`e^{i\omega t}`$. The tortoise coordinate $`r_{}`$ is defined as

$$\frac{dr_{}}{dr}=\frac{r^2+a^2}{\mathrm{\Delta }}$$
(2)
where $`\mathrm{\Delta }=(rr_+)(rr_{})`$ and $`r_\pm =M\pm (M^2a^2)^{1/2}`$ denote the black hole’s inner (Cauchy) and outer (event) horizon. Explicitly we have

$$r_{}=r+\frac{2Mr_+}{r_+r_{}}\mathrm{ln}\left(\frac{r}{r_+}1\right)\frac{2Mr_{}}{r_+r_{}}\mathrm{ln}\left(\frac{r}{r_{}}1\right)$$
(3)
The effective potential $`Q_{lm}`$ is given by

$$Q_{lm}(r_{},\omega ,a)=\frac{K^2+(2am\omega a^2\omega ^2E_{lm})\mathrm{\Delta }}{(r^2+a^2)^2}\frac{dG}{dr_{}}G^2$$
(4)
with $`K=(r^2+a^2)\omega am`$ and $`G=r\mathrm{\Delta }/(r^2+a^2)^2`$. It has the following asymptotic behaviour

$$Q_{lm}\{\begin{array}{cc}k^2\text{as }rr_+(r_{}\mathrm{}),\hfill & \\ \omega ^2\text{as }r+\mathrm{}(r_{}\mathrm{}),\hfill & \end{array}$$
(5)
where $`k=\omega m\omega _+`$ with $`\omega _+=a/2Mr_+`$ denoting the angular frequency of the event horizon. The eigenvalue of the relevant angular equation is denoted by $`E_{lm}`$, cf. teuk .

As a result of (5), the solution of (1) will (to leading order) be a linear combination of simple exponentials as $`r_{}\pm \mathrm{}`$. A physically acceptable solution describes a wave propagating towards the horizon at $`r_{}\mathrm{}`$, i.e. a purely “ingoing” solution. We can normalize this solution according to

$$\psi ^{in}\{\begin{array}{cc}e^{ikr_{}}\text{as }r_{}\mathrm{},\hfill & \\ A_{\mathrm{in}}e^{i\omega r_{}}+A_{\mathrm{out}}e^{i\omega r_{}}\text{as }r_{}+\mathrm{}.\hfill & \end{array}$$
(6)
Another useful solution is purely “outgoing” at infinity,

$$\psi ^{up}e^{i\omega r_{}}$$
(7)
The Wronskian of any two linearly independent solutions of (1) is constant. Using the two solutions above we easily get

$$W(\psi ^{in},\psi ^{up})=2i\omega A_{\mathrm{in}}$$
(8)
Finally, we define the $`𝒮`$-matrix, a quantity familiar from scattering theory taylor , as the ratio

$$𝒮(\omega ,l)(1)^{\mathrm{}+1}\frac{A^{\mathrm{out}}}{A^{\mathrm{in}}}=\mathrm{exp}(2i\delta _l)$$
(9)
where $`\delta _l`$ defines the scattering phase-shifts.

One can define two classes of “resonances” that correspond to purely ingoing waves at the horizon and purely outgoing waves at infinity. Loosely speaking, such solutions are “characteristic” of the black hole spacetime since they do not depend on waves coming in from infinity. The first class of resonances are the QNMs. They correspond to complex frequency poles $`\omega _n`$ of $`𝒮`$ (for integer $`l`$). The QNM frequencies (obviously) correspond to solutions such that the Wronskian (8) vanishes. The second possible class of resonances are the Regge poles. To identify these one allows the angular momentum ($`l`$) to assume complex values (while keeping $`\omega `$ real) and searches for poles $`l_n`$ of $`𝒮`$. The close connection between these classes of solutions has recently been discussed by Décanini, Folacci and Jensen jensen1 .

It is well-known that the numerical determination of QNMs requires considerable care. Since we are discussing (stable) damped modes, the QNM frequencies will have negative imaginary parts. It follows then, that each QNM eigenfunction $`e^{i\omega _nt}\psi ^{\mathrm{in}}(\omega _n,r_{})`$ will grow exponentially both towards infinity and at the horizon. Since one must be able to filter out the ingoing-wave solution at infinity in order to identify a QNM, the calculation requires exponential precision. It should be pointed out that this problem is not physical: It arises because we are viewing only a part of the overall problem, as is clearly shown in discussions of the initial-value problem for perturbed black holes nollert1 ; kokkotas ; naprd . Those studies reveal the true physical meaning of QNMs: They describe part of the late time behaviour of the perturbing field as observed at some fixed distance from the black hole nollert1 . They can be viewed as arising as a part of the field is temporarily “trapped” near the scattering center (the black hole) and leaks out to infinity. The same picture arises in the study of QNMs of much simpler systems such as static objects scattering waves in flat spacetime jensen2 ; lax .

The QNMs arise in scattering situations that can be thought of as “time dependent” (strictly speaking, the evolution and scattering of initial data with compact support). Alternatively, one might consider a time-independent scattering “experiment”. Namely, the scattering of monochromatic plane waves (where the time dependence is trivially given by $`\mathrm{exp}(i\omega t)`$ for any given $`\omega `$). In this case, the most important “observable” quantity is the differential cross section, which is a measure of the scatterer’s size as seen by waves scattered into a certain solid angle taylor . The field at infinity consists of the incoming plane wave plus a scattered piece, and can be written

$$\psi (r_{}+\mathrm{})\psi _{\mathrm{plane}}+f(\theta )\frac{e^{i\omega r_{}}}{r}$$
(10)
where we have assumed, for simplicity, axisymmetric scattering (for a treatment of non-axisymmetric scattering by rotating black holes, see kg2 ). The cross section can be expressed in terms of the (complex-valued) scattering amplitude $`f(\theta )`$,

$$\frac{d\sigma }{d\mathrm{\Omega }}=|f(\theta )|^2$$
(11)
For scalar waves scattered off a Schwarzschild black hole one finds that matzner ,

$$f(\theta )=\frac{1}{2i\omega }\underset{l=0}{\overset{\mathrm{}}{}}(2l+1)(𝒮_l1)P_l(\mathrm{cos}\theta )$$
(12)
from which we see that all the relevant scattering information is encoded in the $`𝒮`$-matrix.

The angular scattering problem is usually studied in terms of the phase-shifts $`\delta _l`$. However, it is not unusual for the sum in (12) to be slowly converging (black-hole scattering being a example, although this difficulty can to some extent be avoided following the prescription in handler ). In many situations it is therefore desirable to have a different representation for the scattering amplitude. The most celebrated alternative is the CAM approach nussenzveig . By means of a Watson-Sommerfeld transformation the sum over $`l`$ in (12) can be rewritten as an integral along a contour in the complex $`l`$-plane cam1 . This integral naturally splits into two terms,

$$f(\theta )=f_P(\theta )+f_B(\theta )$$
(13)
The first term is a sum over the residues of the poles $`l_n`$ of $`𝒮`$; the Regge poles taylor . It is explicitly given by,

$$f_P(\theta )=\frac{i\pi }{\omega }\underset{n=0}{\overset{\mathrm{}}{}}\frac{r_n(l_n+1/2)}{\mathrm{cos}\pi (l_n+1/2)}P_{l_n}(\mathrm{cos}\theta )$$
(14)
where $`r_n`$ is the residue of each Regge pole, i.e. is defined by

$$𝒮\frac{r_n}{ll_n}$$
(15)
in the vicinity of the $`n^{\mathrm{th}}`$ pole, and $`P_l`$ is a Legendre function of complex degree. The second term in (13) is a “background integral”, which is known to be less relevant as far as diffraction effects are concerned. As was demonstrated in cam2 , only the first few poles are required for determining the scattering cross section for angles away from the forward direction (the direction of incoming plane wave). Despite the fact that the problem of scattering by rotating black holes has not yet been studied within the CAM framework, it is natural to define Kerr black hole Regge poles as singularities of $`𝒮`$ in the complex $`l`$ plane while keeping $`\omega `$ real and $`m`$ an integer.

The physical interpretation of Regge states is discussed in standard textbooks nussenzveig ; taylor , cf. also the discussion of scattering of plane waves by a Schwarzschild black hole in cam2 . By approximating the Legendre function for $`|l_n+1/2|1`$ and $`\theta `$ not too close to $`0`$ or $`\pi `$ one can show that

$$P_{l_n}(\mathrm{cos}\theta )\left[\frac{1}{2\pi (l_n+1/2)\mathrm{sin}\theta }\right]^{1/2}\{e^{i(l_n+1/2)(\pi \theta )i\pi /4}+e^{i(l_n+1/2)(\pi \theta )+i\pi /4}\}$$
(16)
From this formula it follows that each Regge pole can be interpreted as a combination of two “surface” waves traveling around the black hole in opposite directions. These waves decay exponentially as they propagate, with an “angular life” $`1/\text{Im}(l_n+1/2)`$. The real part of the Regge pole is associated with the distance $`R_n`$ from the black hole at which the angular decay occurs, Re $`(l_n+1/2)\omega R_n`$ (note that this “localization principle” is formally valid only for $`l1`$ taylor ). As can be verified using the numerical results presented in Section VI, the value of $`R_n`$ lies very close to the critical impact parameter for null geodesics (which is $`b_\mathrm{c}=3\sqrt{3}M`$ for a Schwarzschild black hole). This suggests that (as one would have expected) Regge states describe waves temporarily trapped in the vicinity of the unstable photon orbit, leaking to infinity and to the horizon, as they travel around the black hole.

## III Two “quick and dirty” schemes

In this Section we will discuss two closely related methods for determining black-hole resonances. Both methods are based on the idea of integrating slowly varying phase-functions rather than the, possibly rapidly, oscillating wavefunction itself. The representation of the solutions to Schrödinger type equations in terms of phase-functions was first discussed by Milne milne  and, a long time later, by Newman and Thorson newman . The second key ingredient in our schemes is the use of complex values of the radial coordinate. By analytically continuing the resonant solutions into the complex coordinate plane one can suppress the divergences associated with (in particular) the QNM boundary conditions. In essence, our methods are close relatives to the phase-amplitude method that was developed by one of us na2  to study of scattering resonances. The main difference is that the methods we describe here make use of much simpler integration contours, and are consequently much easier to implement.

### III.1 The Prüfer method

The family of the so-called Prüfer transformations has been widely used in the numerical analysis of Sturm-Liouville problems pryce . It has also been applied to the calculation of resonances and phase-shifts in standard quantum scattering problems, see for example Refs. pajunen1 ; pajunen2 .

Consider a Schrödinger-type equation, where either the frequency $`\omega `$ or the angular momentum $`l`$ are allowed to take complex values,

$$\frac{d^2\psi }{dx^2}+Q(x,l,\omega )\psi =0$$
(17)
We consider cases where the radial coordinate $`x`$ spans the entire real axis (in the black hole case $`x`$ corresponds to the tortoise coordinate $`r_{}`$), and the effective potential $`Q`$ is a single barrier with asymptotic behaviour

$$Q\{\begin{array}{cc}q^2\text{as }x\mathrm{}\hfill & \\ \omega ^2\text{as }x+\mathrm{}\hfill & \end{array}$$
(18)
where $`q=q(\omega )`$. It follows then, that the function $`\psi (x,l,\omega )`$ will be a linear combination of exponentials as $`x\pm \mathrm{}`$. Following the discussion in the previous Section, let us assume that we are interested in a solution to (17) which is purely “ingoing” as $`x\mathrm{}`$ and mixed ingoing/outgoing as $`x+\mathrm{}`$. That is, we want to determine the “in” solution, see eqn. (6), which can be written

$$\psi \{\begin{array}{cc}e^{iqx}\text{as }x\mathrm{},\hfill & \\ B\mathrm{sin}[\omega x+\zeta ]\text{as }x+\mathrm{},\hfill & \end{array}$$
(19)
where $`\zeta `$ and $`B`$ are some complex constants.

We can write the exact solution of (17) in the form

$$\psi (x,l,\omega )=\mathrm{exp}[P(x,l,\omega )𝑑x]$$
(20)
where $`P`$ is the logarithmic derivative of $`\psi `$ (a prime denotes derivative with respect to $`x`$),

$$\frac{\psi ^{}}{\psi }=P$$
(21)
This is our first phase-function. Obviously, it should obey the boundary condition $`Piq`$ for $`x\mathrm{}`$. Alternatively, we can express the function $`\psi `$ and its derivative via the transformation,

$`\psi (x,l,\omega )`$ $`=`$ $`B\mathrm{sin}[\omega x+\stackrel{~}{P}(x,l,\omega )]`$ (22)
$`\psi ^{}(x,l,\omega )`$ $`=`$ $`B\omega \mathrm{cos}[\omega x+\stackrel{~}{P}(x,l,\omega )]`$ (23)
where $`\stackrel{~}{P}`$ is the Prüfer phase function (which should have $`\zeta `$ as its limiting value as $`x+\mathrm{}`$).

The two phase functions satisfy non-linear first order differential equations;

$$\frac{dP}{dx}+P^2+Q=0$$
(24)
$$\frac{d\stackrel{~}{P}}{dx}+\left(\omega \frac{Q}{\omega }\right)\mathrm{sin}^2(\omega x+\stackrel{~}{P})=0$$
(25)
and they are directly related by

$$\stackrel{~}{P}=\omega x+\frac{1}{2i}\mathrm{log}\left[\frac{iP\omega }{iP+\omega }\right]$$
(26)
or

$$P=\omega \mathrm{cot}(\omega x+\stackrel{~}{P})$$
(27)
A major advantage of using phase-functions can be seen immediately from Eqns. (24) and (25), assuming for the moment that we are dealing with real $`\omega `$ and $`l`$. Using (18), we see that (24) and (25) lead to

$`{\displaystyle \frac{dP}{dx}}`$ $``$ $`0\text{as}x\mathrm{}`$ (28)
$`{\displaystyle \frac{d\stackrel{~}{P}}{dx}}`$ $``$ $`0\text{as}x+\mathrm{}`$ (29)
Hence, the functions $`P`$ and $`\stackrel{~}{P}`$ are slowly varying for large negative and positive values of $`x`$, respectively. This means that the oscillatory behaviour of the original wavefunction in those regimes \[which would cause the build up of unnecessary numerical errors in a direct integration of (17)\] has been effectively removed.

Finally, knowledge of the $`\stackrel{~}{P}`$ phase-function at $`x+\mathrm{}`$ provides us with the $`𝒮`$-matrix,

$$𝒮(l,\omega )=(1)^{l+1}e^{2i\zeta }$$
(30)
Obviously, the described method can be applied to the black-hole problem with $`xr_{}`$, $`qk`$ and (17) becoming the Teukolsky radial equation (1).

At the practical level, the numerical integration for the $`P`$ function shoud be initiated at a point $`r_{}^H`$ (say), close to the event horizon. Assuming that $`Q`$ is a slowly varying function for large $`|r_{}|`$, which is certainly the case in the applications under consideration, we can initiate the integration of (24) from the lowest order WKB approximation

$$P\pm iQ^{1/2}\text{ as }|r_{}|\mathrm{}$$
(31)
(alternative approximations, eg. a power series expansion, will work equally well). In this way we can readily find a solution that corresponds to the (analytically continued) ingoing wave boundary condition at the horizon. This solution is then integrated towards the matching point $`r_{}^{\mathrm{ma}}`$. As the integration reaches the vicinity of the potential peak (the local minimum of $`Q`$ in the Kerr case) we switch to the $`\stackrel{~}{P}`$ function, and then carry on the integration to “infinity”. It is mandatory to change phase-function if one wants to ensure that the integrated function remains smooth and well behaved. This is due to the fact that the phase-functions undergo a drastic change in their behaviour (they become rapidly oscillating and may even have singularities) as they are extended through the region of the potential peak. This is a manifestation of the so-called Stokes phenomenon — the “switching on” of small exponentials in the solution to an equation of form (17), well known from WKB-and phase-integral theory berry ; froman  (see also the discussion in Appendix 1 of cam3 ). This situation is demonstrated in Fig. 1 where we plot the functions $`P(r_{})`$ and $`\stackrel{~}{P}(r_{})`$ for typical parameters $`l=m=2`$, $`\omega M=0.5`$, $`a=0.5M`$.

The Prüfer approach has already been applied in the context of scattering of plane scalar waves by rotating black holes kg2 . In that study, where both $`\omega `$ and $`l`$ were assumed to be real, the use of the phase-functions described above facilitated the computation of accurate high frequency phase-shifts (which would be problematic if the original wave equation were to be integrated handler ). We can follow exactly the same procedure when $`l`$ is allowed to acquire complex values (as long as $`\omega `$ remains real). In this way we have calculated the first few Regge poles of both Schwarzschild and Kerr black holes. The numerical results will be discussed in the following Section.

Complications arise when $`\omega `$ becomes complex while $`l`$ is kept real (integer). From Eqn. (25) we see that the term which represents an ingoing (outgoing) wave becomes exponentially small (large) at large $`x`$ (for Im $`\omega <0`$). A similar situation appears close to the horizon. Therefore, we are still facing one of the problems associated with a direct integration of the wave equation. An efficient and elegant way to avoid this difficulty is to consider the problem in the complex $`r_{}`$-plane. The success of this approach has been repeatedly demonstrated with the application of the phase-integral/amplitude formalisms to the calculation of black-hole QNMs na2 ; na6 ; pi . The main idea is to integrate our equations not along the real axis, but instead use a suitably chosen contour $`C`$ such that the ingoing/outgoing wave solutions, eg. $`e^{\pm i\omega r_{}}`$, have comparable asymptotic magnitudes (of order unity) as $`r_{}\mathrm{}`$. The simplest choice of integration contour is a straight line with slope $``$Im $`\omega /`$ Re $`\omega `$, see Fig 2. Similarly, a line with slope $``$Im $`k/`$ Re $`k`$ will guarantee purely oscillatory behaviour near the event horizon. Such paths are parallel (for $`|r_{}|\mathrm{}`$) to the so-called anti-Stokes lines which play a crucial role in WKB and phase-integral theory froman ; na2 . The anti-Stokes lines are curves along which $`Q^{1/2}𝑑r_{}`$ is purely real and as a consequence the WKB solution of (1) is purely oscillatory. It should be noted that a similar integration contour has been used in a calculation of relativistic stellar pulsation modes na4 .

We use a parameterization

$$r_{}=r_{}^{\mathrm{ma}}+\rho e^{i\beta }$$
(32)
with $`\mathrm{}<\rho <+\mathrm{}`$ and $`\beta `$ given by,

$$\beta =\{\begin{array}{cc}Arg(\omega )\text{for }\rho >0,\hfill & \\ Arg(k)\text{for }\rho <0,\hfill & \end{array}$$
(33)
The matching point $`r_{}^{\mathrm{ma}}`$ is chosen to lie on the real axis ( $`\rho =0`$) and some experimentation shows that reliable results are obtained as long as the matching is performed near the peak of the curvature potential (the extremum of $`Q`$).

Rewriting Eqns. (24) and (25) in terms of $`\rho `$ we have

$$\frac{dP}{d\rho }+e^{i\beta }\left(P^2+Q\right)=0$$
(34)
$$\frac{d\stackrel{~}{P}}{d\rho }+e^{i\beta }\left(\omega \frac{Q}{\omega }\right)\mathrm{sin}^2(\stackrel{~}{P}+|\omega |\rho +\omega r_{}^{\mathrm{ma}})=0$$
(35)
The radial coordinate $`r`$ will also be assumed complex, and we determine it from

$$\frac{dr}{d\rho }=e^{i\beta }\frac{\mathrm{\Delta }(r)}{r^2+a^2}$$
(36)
With the help of (3) it is easy to see that

$$rr_+(1+e^{cr_{}})$$
(37)
where $`c=(r_+r_{})/2Mr_+`$ as $`\rho \mathrm{}`$, and that $`rr_{}`$ as $`\rho +\mathrm{}`$. Hence, the proposed integration contour corresponds to a straight line at large distances in the $`r`$-plane, and to a clockwise spiral as the horizon is approached. This strongly resembles the more complicated integration contour used in the phase-amplitude calculation of QNMs, see figures 1 and 2 in na2 .

Not surprisingly, some care must be taken in implementing complex integration contours in black-hole problems. For example, we need to understand the analytic properties of the potential $`Q`$. From Eqn. (4) we deduce that this quantity will have poles at $`r=\pm ia`$. Since the function $`r_{}(r)`$ is multivalued, there will be a series of poles in the left half of the $`r_{}`$-plane, symmetrically located with respect to the real axis. All of these poles have the same real part (which is zero for the non-rotating case), while their imaginary parts satisfy Im $`r_{}\text{Re}r_{}`$. In practice, the presence of these singularities provides an upper limit to the complex-rotation angle. This is easiest illustrated by the Schwarzschild case. According to (3) the point $`r=0`$ translates into

$$r_{}(0)=2M\mathrm{log}(1)=\pm 2Mi\pi (1+2k)\text{with}k=0,1,2,\mathrm{}$$
(38)
It is easy to understand that one should not rotate the integration contour beyond this angle.

Another important issue that deserves some discussion concerns the boundary conditions in the complex plane. When we are solving the system of phase equations (34) and (35) we are actually imposing the QNM boundary condition on the contour $`C`$ at $`\rho \pm \mathrm{}`$. However, the “true” boundary conditions are to be imposed on the real axis ($`r_{}\pm \mathrm{}`$). It is perfectly natural to wonder if we are, in some way, altering the problem by using the contour $`C`$. This question has been discussed in the context of WKB and phase-integral studies (and the complex $`r`$-plane) araujo ; froman . It was concluded that it is safe to analytically continue the boundary condition in the proposed way.

### III.2 QNM excitation coefficients

We have used the Prüfer transformation scheme described above to calculate both QNMs and Regge poles for scalar waves in the geometry of a rotating black hole. The method is easily extended to other classes of perturbing fields, the most interesting of which is likely to be that of gravitational perturbations. As will be discussed below the method is somewhat limited for extremely rapidly damped QNMs, and hence it cannot be used to investigate the asymptotic behaviour of the mode spectrum motl . However, from an astrophysical point of view the slowly damped modes are the most important and for these our scheme works very well, indeed.

In the case of QNMs one might also want to assess to what extent the modes are excited given a certain initial perturbation. As discussed in naprd  the answer can be quantified in an initial-data independent way in terms of the “excitation coefficient”

$$B_n=\frac{A_{\mathrm{out}}(\omega _n)}{2\omega _n\alpha _n}$$
(39)
where we have assumed that $`A_{\mathrm{in}}\alpha _n(\omega \omega _n)`$ in the vicinity of the $`n^{th}`$ mode frequency. While one can formulate expressions for this quantity within our previous scheme, we have found that the precision of the actual numerical results is somewhat unsatisfactory. To remedy this we have devised an alternative scheme. The purpose of the following discussion is to present this scheme, with the particular aim of achieving accurate results for the excitation coefficients.

The first half of the calculation remains identical to the one discussed above, and we thus obtain a solution $`P_L(r_{}^{\mathrm{ma}})`$ (say) which satisfies that boundary condition of ingoing waves at the event horizon. At spatial infinity we now obtain the appropriate solution in a similar way. Let us choose the phase of the square root in such a way that the solution corresponding to $`P_++iQ^{1/2}(r_{})`$ has the desired outgoing-wave character. Once this solution is extended to the matching point we have a second solution $`P_+(r_{}^{\mathrm{ma}})`$. Thus a resonant overall solution, a QNM, follows from the simple condition

$$P_+(r_{}^{\mathrm{ma}})=P_L(r_{}^{\mathrm{ma}})$$
(40)
Now we want to determine also the excitation coefficients. We can do this by considering the solution for a general frequency, that contains both in- and outgoing wave at infinity, rather than a specific QNM solution. Assume (as before) that $`P_++iQ^{1/2}(r_{})`$ corresponds to an outgoing wave solution, and that $`P_{}iQ^{1/2}(r_{})`$ represents an ingoing wave. This means that, after integrating these two functions to the matching point the general solution corresponds to the linear combination

$$\psi _R(r_{}^{\mathrm{ma}})=A\mathrm{exp}\left[_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}P_+𝑑r_{}\right]+B\mathrm{exp}\left[_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}P_{}𝑑r_{}\right]$$
(41)
where the integration is initiated at $`r_{}^{\mathrm{}}`$. We can write the logarithmic derivative of this function as

$$\frac{\psi _R^{}}{\psi _R}=\frac{P_++(B/A)P_{}\mathrm{exp}\left[_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}(P_{}P_+)𝑑r_{}\right]}{1+(B/A)\mathrm{exp}\left[_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}(P_{}P_+)𝑑r_{}\right]}$$
(42)
The condition for an overall solution at $`r_{}^{\mathrm{ma}}`$ is that the logarithmic derivative is continuous, i.e. that the above expression is equal to $`P_L(r_{}^{\mathrm{ma}})`$. This condition can be written

$$\frac{B}{A}=\frac{P_+P_L}{P_LP_{}}\mathrm{exp}\left[_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}(P_+P_{})𝑑r_{}\right]\text{ for }r_{}=r_{}^{\mathrm{ma}}$$
(43)
Now recalling that $`P_{}`$ corresponds to the ingoing wave solution at infinity we know that a QNM must correspond to $`B=0`$, and we see that the above condition is the same as the one we derived previously.

To determine the excitation coefficient requires a little more work. First we need to relate our amplitudes $`A`$ and $`B`$ to the asymptotic amplitudes $`A_{\mathrm{out}}`$ and $`A_{\mathrm{in}}`$. Recall that the asymptotic amplitudes are determined from

$$\psi A_{\mathrm{out}}e^{i\omega r_{}}+A_{\mathrm{in}}e^{i\omega r_{}}\text{ as }r_{}\mathrm{}$$
(44)
We want to compare this expression to our solutions. Since we have assumed that the WKB approximation is valid at the point where we initiate integration ($`r_{}^{\mathrm{}}`$) we can use

$$_r_{}^{\mathrm{}}^r_{}P_\pm 𝑑r_{}\pm i_r_{}^{\mathrm{}}^r_{}Q^{1/2}𝑑r_{}=\pm i\omega (r_{}r_{}^{\mathrm{}})\pm i_r_{}^{\mathrm{}}^r_{}\left(\frac{Q\omega ^2}{Q^{1/2}+\omega }\right)𝑑r_{}$$
(45)
beyond $`r_{}^{\mathrm{}}`$ (the last integral should be easy to evaluate). We then have

$$\psi _+A\mathrm{exp}\left[_r_{}^{\mathrm{}}^r_{}P_+𝑑r_{}\right]=A\mathrm{exp}\left[i\omega (r_{}r_{}^{\mathrm{}})+i_r_{}^{\mathrm{}}^r_{}\left(\frac{Q\omega ^2}{Q^{1/2}+\omega }\right)𝑑r_{}\right]\text{ for }|r_{}|>|r_{}^{\mathrm{}}|$$
(46)
And we can identify

$$A_{\mathrm{out}}=A\mathrm{exp}\left[i\omega r_{}^{\mathrm{}}+i_r_{}^{\mathrm{}}^{\mathrm{}}\left(\frac{Q\omega ^2}{Q^{1/2}+\omega }\right)𝑑r_{}\right]$$
(47)
Similarly, we get

$$A_{\mathrm{in}}=B\mathrm{exp}\left[i\omega r_{}^{\mathrm{}}i_r_{}^{\mathrm{}}^{\mathrm{}}\left(\frac{Q\omega ^2}{Q^{1/2}+\omega }\right)𝑑r_{}\right]$$
(48)
Collecting the above results, we have

$$\frac{A_{\mathrm{in}}}{A_{\mathrm{out}}}=\frac{P_+P_L}{P_LP_{}}\mathrm{exp}\left[_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}(P_+P_{})𝑑r_{}+2i\omega r_{}^{\mathrm{}}2i_r_{}^{\mathrm{}}^{\mathrm{}}\left(\frac{Q\omega ^2}{Q^{1/2}+\omega }\right)𝑑r_{}\right]$$
(49)
The final step in the derivation consists of taking the derivative of this expression with respect to $`\omega `$ and evaluating it at a QNM frequency. This immediately leads to

$$\frac{\alpha _n}{A_{\mathrm{out}}}=\frac{1}{P_LP_{}}\left[\frac{d}{d\omega }(P_+P_L)\right]\mathrm{exp}\left[_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}(P_+P_{})𝑑r_{}+2i\omega r_{}^{\mathrm{}}2i_r_{}^{\mathrm{}}^{\mathrm{}}\left(\frac{Q\omega ^2}{Q^{1/2}+\omega }\right)𝑑r_{}\right]$$
(50)
from which we can calculate the desired excitation coefficient.

To conclude, let us summarize the additional steps that are required in a calculation of the excitation coefficient once a QNM frequency has been determined. First we need to evaluate the two integrals $`_r_{}^{\mathrm{}}^{r_{}^{\mathrm{ma}}}P_\pm 𝑑r_{}`$. This is best done by introducing additional dependent variables, let us call them $`\mathrm{\Phi }_\pm `$, and then integrating the equations

$$\mathrm{\Phi }_\pm ^{}=P_\pm $$
(51)
together with (24) from $`r_{}^{\mathrm{}}`$ to $`r_{}^{\mathrm{ma}}`$. Secondly, we need the derivatives $`\frac{dP_+}{d\omega }`$ and $`\frac{dP_L}{d\omega }`$. These quantities are also best calculated by introducing additional variables, i.e. by integrating

$$\frac{d}{dr_{}}\left(\frac{dP_\pm }{d\omega }\right)=\frac{dQ}{d\omega }2P_\pm \frac{dP_\pm }{d\omega }$$
(52)
using appropriate initial data, i.e.

$$\frac{dP_\pm }{d\omega }\pm \frac{i}{2Q^{1/2}}\frac{dQ}{d\omega }$$
(53)
at $`r_{}^{\mathrm{}}`$. Thirdly, we need to evaluate

$$_r_{}^{\mathrm{}}^{\mathrm{}}\left(\frac{Q\omega ^2}{Q^{1/2}+\omega }\right)𝑑r_{}$$
(54)
This is readily done by first expanding the integrand in inverse powers of $`r`$ and then integrating. The result is convergent in all relevant cases.

### III.3 Numerical implementation

Before concluding the description of the methods we have devised to study black hole resonances, it is worth commenting on the numerical implementation of the schemes. Overall we have a system of two coupled nonlinear first-order equations to solve numerically (for either choice of phase-function), namely, eqns (34) or (35) together with (36). In reality, since we are dealing with intrinsically complex quantities, we split each equation in its real and imaginary parts. At the end of the day, we arrive at four equations to be solved. This is done by means of a standard variable-step Addams method.

To complete the solution of the problem, we combine the result of the integration of the phase functions with a complex root-finder routine (making use of Müller’s algorithm recipes ) to search the complex frequency/angular momentum plane for zeros of the ratio $`A_{\mathrm{in}}/A_{\mathrm{out}}`$. A prior knowledge of the approximate location of a resonance greatly speeds up the overall procedure.

## IV Numerical results

The main aim of this investigation was to devise simple and efficient methods for computing black-hole resonances and scattering phase-shifts. As should be clear from the previous Sections, the phase-function schemes that we are discussing are straightforward to implement. They also turn out to be numerically robust and provide results of excellent accuracy. We have tested the methods by calculating known results for both Schwarzschild and Kerr black holes. The obtained results serve as a validation of our numerical code(s). Having confirmed the reliability of our schemes we have extended the study to, in particular, the resonances of rotating black holes. Below, we will describe new results for both QNMs and Regge poles of Kerr black holes. These results shed interesting light on issues concerning the excitation of the long lived QNMs that are known to exist for near extreme Kerr black holes, the relation between the approximate QNMs for $`aM`$ used in (for example) kg1 , and also the effect of rotation on the black hole Regge poles.

### IV.1 Quasinormal modes

As a first application of our numerical methods, we have calculated the first QNMs of a Schwarzschild black hole. The obtained results agree perfectly with results in the literature na2 , and confirm that our method is reliable as long as the modes are not damped too rapidly. The method breaks down when the imaginary part of the mode frequency is about one order of magnitude larger than the real part. This means that one should not expect to be able to use our scheme to study the asymptotic behaviour of the QNM spectrum. However, all astrophysically relevant QNMs are readily calculated. We have also tested the obtained excitation coefficients by comparing results obtained for gravitational perturbations (using a code based on the Sasaki-Nakamura equation SN ) to the phase-integral results in naprd0  and the data obtained by Leaver leaver1 . A sample of results is given in Table 1. By comparing the listed results to those in Table II of naprd0  we conclude that the present results, which agree very well with those of Leaver apart from a few sign discrepancies already pointed out in naprd0 , are reliable. In fact, we believe that the present method provides the (so far) most efficient tool for investigating the QNM excitation problem.

As the black hole acquires rotation, each of the Schwarzschild QNMs (for a given $`l`$) splits into $`2l+1`$ distinct modes. We have tested the reliability of our scheme by reproducing the Kerr QNM frequencies obtained by Leaver, cf. Tables 2 and 3 in leaver1 . A sample of results for $`l=2`$, obtained by implementing our Prüfer phase function scheme for the Sasaki-Nakamura equation, are provided in Table 2. These results are in perfect agreement with those obtained by Leaver’s continued fraction calculation for $`m=0`$ and $`m=\pm 1`$. Furthermore, Table 2 is useful since the $`l=2`$ QNMs may contribute significantly to the gravitational-wave signal from, for example, a black hole born following a supernova core collapse or a compact binary merger (see finn  for discussion).

Kerr QNMs are known to exhibit interesting behaviour as $`a`$ varies detweiler1 . In particular, we see from Table 2 that rotation has a pronounced effect on modes with $`m=l>0`$. These can be thought of as “equatorial” and “co-rotating” with the black hole (cf. the symmetry properties of the associated angular functions, which in the case of scalar waves are spheroidal harmonics $`S_{lm}^{a\omega }(\theta )`$). In the limit $`aM`$ the co-rotating modes have been found to become long-lived (Im $`\omega _nM0`$) and accumulate at the frequency $`m\omega _+`$ detweiler1 ; leaver2 ; onozawa1 ; kg1 . Hence, rotation has a significant effect on these modes. In contrast, the “counter-rotating” modes ($`m<0`$) are much less affected by rotation and remain relatively close to the initial Schwarzschild QNMs for all values of $`a`$. In Fig. 3 we illustrate this behaviour by showing the change in the first few scalar field QNM frequencies as the rotation rate of the black hole is varied. The data in the figure corresponds to $`0a/M0.99`$.

The behaviour of the QNMs as the black hole approaches extreme rotation is particularly intriguing. Unfortunately, our method breaks down beyond $`a0.999M`$. The likely reason for this is the multivaluedness of the complex $`r_{}`$ coordinate, cf. for example (38), and the need to initiate the integration close to the horizon, see eqn (37). This becomes tricky as the event horizon and the inner Cauchy horizon come close together. A similar problem is known to arise in the phase-integral study of the QNMs of a Reissner-Nordström black hole na6 . This means that we cannot push our calculation all the way to the extreme limit. Of course, we can determine QNMs for all rotation rates that are likely to be astrophysically relevant. Still, the behaviour in the extreme limit is interesting from a conceptual point of view. Especially since one can approximate the QNM spectrum for near extreme black holes detweiler1 . In figure 4 we compare our numerical results for scalar field QNMs to approximate results for $`aM`$ and $`\omega m\omega _+`$. The figure shows that the numerical results for the first two $`l=m=2`$ QNMs, as extended from the Schwarzschild limit, approach the first two approximate modes as $`aM`$. This is very encouraging since it demonstrates a correspondence between two rather different computational schemes. It also provides support for any argument concerning the nature and relevance of these long-lived modes that is based on the approximate formulae. The comparison also indicates the range of validity of Detweiler’s approximate mode-condition detweiler1  in a useful way.

It has been argued that the existence of long-lived QNMs for rapidly rotating black holes could prove important for gravitational-wave astronomy. The basic idea is that a slowly damped mode will be easier to detect provided that it is initially excited to the same amplitude as a more rapidly damped mode. Since we now have the means to calculate the excitation coefficients for the Kerr QNMs we can analyse this problem quantitatively. As we have indicated in a previous paper kg1 , our results do not support the idea that the long-lived QNMs will be of particular importance. The reason for this can be seen from Table 3 , where we provide results for the fundamental $`m=\pm 2`$ QNMs for various values of the rotation parameter. From this data we see that while the absolute value of the excitation coefficient for the counter-rotating QNM only changes by about a factor of two as the black hole spins up towards the extreme limit, it drops by about two order of magnitude for the co-rotating mode. This indicates that the co-rotating long-lived QNMs may be significantly more difficult to excite in a dynamical process.

We can use the results in Table 3 to discuss the “detectability” of the long-lived QNMs. Following kg1  we express the effective gravitational-wave amplitude associated with a single QNM (with complex frequency $`\omega _n`$) as

$$h_{\mathrm{eff}}\sqrt{\frac{\text{Re }\omega _n}{\text{Im }\omega _n}}\frac{A^{\mathrm{out}}}{\alpha _n}=2\sqrt{\frac{\text{Re }\omega _n}{\text{Im }\omega _n}}\omega _nB_n$$
(55)
Given our numerical data, we can compare the result for co- and counter-rotating (scalar field) modes. Figure 5 (which is identical to Figure 1 in kg1 ) shows the outcome of this comparison. From the figure it is clear that, even though the co-rotating mode is much longer lived, its “detectability” decreases significantly as $`aM`$. We conclude that the present results suggest that the effective amplitude of the slowest damped QNM for a rapidly rotating black hole may be several orders of magnitude smaller than the amplitude of the corresponding mode of a Schwarzschild black hole, or indeed the counter-rotating Kerr modes. The possible constructive interference of the many long-lived QNMs that exist for near extreme Kerr black holes, and which could still lead to a relevant gravitational-wave signal, has been discussed in kg1 .

As a final application of our method we have traced the path of a high overtone mode ($`n=9,l=2,m=0`$) for varying black-hole spin. The result is shown in Fig. 6. The purpose of this exercise was to investigate whether the QNM trajectory exhibits the characteristic “looping” behaviour, first found in the studies of gravitational and electromagnetic field QNMs for Kerr and Reissner-Nordström black holes onozawa1 ; onozawa2 . Fig. 6 confirms that this behaviour is also present for scalar fields. It has been suggested onozawa1 ; onozawa2  that this peculiar behaviour is somehow associated with the change in the potential near to the horizon as the spin varies. However, since it is unclear how the high overtone modes are related to the potential a fully satisfactory explanation is still missing (see, however, liu  for an interesting point of view.

### IV.2 Regge poles

As mentioned in the Introduction, Regge pole calculations are comparatively easy to perform since as the associated wavefunction (6) remains well behaved as $`r_{}\pm \mathrm{}`$. This means that one does not have to deal with divergences, and the identification of in- and outgoing wave solutions is straightforward. For a given complex $`l`$, the set of equations (24) and (25) can be directly integrated along the real $`r_{}`$-axis without any serious problems.

To test our numerical code we have calculated the first few Regge poles for a Schwarzschild black hole for three different frequencies. The obtained results are listed in Table 4. As far as the Regge pole positions are concerned, the results are in agreement with the approximate phase-integral results given by Andersson cam2 . However, our calculation has confirmed an earlier suggestion by Jensen jensen3  that the phase-integral residues had a systematic error. As in the case of QNMs, we believe the present method is the most efficient and reliable way to calculate black-hole Regge poles.

We next consider the Regge poles of a Kerr black hole. This is a problem that, as far as we are aware, has not been discussed before. As in the case of QNMs, the loss of spherical symmetry in the background spacetime unfolds the $`2l+1`$ degeneracy of the non-rotating resonances. The pole trajectories (for the first two Regge poles) are illustrated in Fig. 7 for $`\omega M=1`$ and $`0a/M0.99`$, and a sample of numerical data is given in Table 5. Figure 7 shows obvious similarities with the QNMs trajectories in figure 3. This is not surprising given the close relationship between the two classes of scattering resonances jensen1 .

The most important conclusion we can draw from our results concerns the physical interpretation of the Kerr Regge poles. First of all, it is clear that $`m>0`$ $`(m<0)`$ Regge poles would describe co-rotating (counter-rotating) waves. Intuitively, we would expect that the “angular life” (the inverse of Im $`l_n`$, cf. cam2 ) of the former would be significantly increased due to rotational frame dragging. As in the case of QNMs, we find that the $`m=l`$ poles are the most “sensitive” to the black hole’s rotation. This strongly suggests that these poles can be associated with waves propagating around the black hole in the vicinity of the equatorial plane. More evidence to support this intuitive picture comes from the behaviour of Re $`l_n`$. The decrease (increase) of this quantity, as the black hole spins up, is qualitatively expected from the corresponding decrease (increase) of the critical impact parameter and unstable photon orbit for prograde (retrograde) motion in Kerr spacetime.

## V Concluding remarks

We have developed simple numerical methods for calculating black hole resonances (QNMs and Regge poles). The basic ingredients of our schemes are (i) the use of (Prüfer) phase functions in place of the original wavefunction (ii) the integration of the resulting equations along rotated contours in the $`r_{}`$-plane (when the frequency is complex), and (iii) the implementation of Müller’s algorithm for locating the complex roots of $`1/𝒮(\omega ,l)`$ function (the inverse scattering matrix element, whose zeros define the resonances).

We have applied our schemes to the problem of calculating gravitational and massless scalar field perturbation QNMs and Regge poles for rotating and non-rotating black holes. The obtained numerical data verify all existing results in the literature with high precision. Moreover, we have provided some new results. In particular, we have shown how the slowest damped Kerr QNMs match the “long-lived” modes of a near extreme black hole which can be derived by an analytic approximation to the Teukolsky equation. This ties up some loose ends as far as the Kerr QNMs are concerned and also provides important support for using the approximate results to discuss the dynamics of near extreme black holes kg1 . Our results for the excitation coefficients of Kerr QNMs add fuel to the debate concerning the detectability of gravitational waves associated with the long lived modes of a near extreme black hole. We have also provided the first results for Kerr Regge poles and discussed briefly how these results make sense given the established physical interpretation of their Schwarzschild counterparts.

In conclusion, we find that our methods can be regarded as reliable and useful, apart from (i) for very rapidly damped QNMs, and (ii) in the case of near extreme rotations ($`a/M>0.999`$ or so). This is, however, a small price to pay considering the profound simplicity and adaptability of the method. This is the main advantage of our approach. All other accurate schemes require a more sophisticated analysis: Leaver’s continued fraction method leaver1  is not easy to implement given the slow convergence of the relevant series, and the complex-coordinate phase-amplitude method developed by one of us na2  makes use of much more intricate integration contours than we used in the present analysis. Finally, we should note that, even though we considered mainly the case of scalar perturbations in this paper, our method can readily be used for the more realistic case of gravitational perturbations. We therefore believe that our “quick and dirty” approach to black hole resonances may prove useful in a variety of exciting contexts.

###### Acknowledgements.

The work of K.G. was financially supported by the State Scholarships Foundation of Greece. N.A. is a Philip Leverhulme Prize Fellow and also acknowledges support from the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137).


format: text

text:
# CI and CO in the nearby spiral galaxies IC 342 and Maffei 2

## 1 Introduction

Molecular gas is a major constituent of the interstellar medium in galaxies. Within the inner kiloparsec, many spiral galaxies exhibit a strong concentration of molecular gas towards their nucleus. We have conducted a programme to observe a number of nearby galaxies in various CO transitions and in the 492 GHz 3P13P0 CI transition in order to determine the physical condition of such central molecular gas concentrations. Results on NGC 253 (Israel, White $`\&`$ Baas 1995), NGC 7331 (Israel $`\&`$ Baas 1999), NGC 6946 and M 83 = NGC 5236 (Israel $`\&`$ Baas 2001) have already been published. In this paper, we present the results for the very nearby galaxies IC 342 and Maffei 2 whose basic properties are summarized in Table 1.

Both are major members of the IC 342/Maffei group, which is located at a distance of about 2 Mpc (Huchtmeier et al., 2000 and references therein). This appears to be the galaxy group closest to the Local Group but it is, unfortunately, located in the sky very close to the Galactic plane. Consequently, its members suffer high foreground extinction rendering several of them, including Maffei 2, all but invisible at optical wavelengths. Mainly for this reason, the distance of Maffei 2 is still quite uncertain. However, it is usually assumed that it exceeds that of IC 342. We have thus assumed a distance of 2.7 Mpc which is in line with the overall cluster value, but 50 per cent greater than the reasonably well established distance to IC 342. Both IC 342 and Maffei 2 have been studied well at longer wavelengths, including those of the CO line. Inasmuch as IC 342 is, together with NGC 253 and M 82, one of the strongest (sub)millimeter line emitters in the sky, it has been observed frequently. In the $`J`$=1–0 transition of $`{}_{}{}^{12}\mathrm{CO}`$, it was one of the first galaxies mapped at $`1^{}`$ resolution (Morris $`\&`$ Lo 1978; Rickard $`\&`$ Palmer 1981; Young $`\&`$ Scoville 1982). More extensive mapping of the $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ isotopes in both the $`J`$=1–0 and $`J`$=2–1 transitions was presented by Eckart et al. (1990) and Xie et al. (1994). Emission in the $`J`$=3–2 transition was detected and mapped by Ho et al. (1987), Steppe et al. (1990) and Irwin $`\&`$ Avery (1992), while the $`{}_{}{}^{13}\mathrm{CO}`$ isotope was detected in this transition by Wall $`\&`$ Jaffe (1990) and Mauersberger et al. (1999). Finally, the $`J`$=4–3 and, with a somewhat uncertain calibration, the $`J`$ = 6–5 $`{}_{}{}^{12}\mathrm{CO}`$ transitions were detected by Güsten et al. (1993) and Harris et al. (1991). The $`{}_{}{}^{12}\mathrm{CO}/{}_{}{}^{13}\mathrm{CO}`$ isotope ratio was determined to be $`30`$ by Henkel et al. (1998). High-resolution aperture synthesis maps of IC 342 have been published in $`J`$=1–0 $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ by Ishizuki et al. (1990), Turner $`\&`$ Hurt (1992), Wright et al., (1993) and Sakamoto et al. (1999), and in $`J`$=2–1$`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ by Turner et al., (1993) and Meier et al. (2000).

Maffei 2 was first detected in $`{}_{}{}^{12}\mathrm{CO}`$ by Rickard et al. (1977) and subsequently mapped at 23′′ resolution by Weliachew et al. (1988). Emission in the $`J`$=2–1 and $`J`$=3–2 transitions was measured by Wall et al. (1993) and Mauersberger et al. (1999). Maffei 2 was mapped in $`J`$=1–0 $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ by Ishiguro et al. (1989) and Hurt $`\&`$ Turner (1991), in $`J`$=2–1 $`{}_{}{}^{12}\mathrm{CO}`$ by Sargent et al. (1985) and in $`J`$=3–2 $`{}_{}{}^{12}\mathrm{CO}`$ by Hurt et al. (1993).

Although the molecular line emission from Maffei 2 is somewhat weaker than that from IC 342, both are frequently included in the same observing programs. Both have been observed in a great variety of other molecular species: CS (Mauersberger $`\&`$ Henkel, 1989; Mauersberger et al., 1989; Sage et al., 1990; Paglione et al., 1995), HCN and HCO+ (Nguyen-Q-Rieu et al., 1992; Downes et al., 1992; Jackson et al. 1995; Paglione et al. 1997), H13CO+ and N2H+ (Mauersberger $`\&`$ Henkel, 1991); HNC (Hüttemeister et al., 1995), CN and HC3N (Henkel et al. (1988), HNCO (Nguyen-Q-Rieu et al. 1991), OCS (Mauersberger et al., 1995) as well as H2CO and CH3OH (Hüttemeister et al., 1997). In particular the discovery of strong NH3 emission, first from IC 342 (Ho et al. 1990 and references therein), then from Maffei 2 (Henkel et al. 2000; Takano et al. 2000) showed that a significant fraction of the molecular gas in the centers of these galaxies is both dense and hot.

## 2 Observations

All observations described in this paper were carried out with the 15m James Clerk Maxwell Telescope (JCMT) on Mauna Kea (Hawaii) 111The James Clerk Maxwell Telescope is operated on a joint basis between the United Kingdom Particle Physics and Astrophysics Council (PPARC), the Netherlands Organisation for Scientific Research (NWO) and the National Research Council of Canada (NRC).. Details are given in Table 2. Up to 1993, we used a 2048 channel AOS backend covering a band of 500 MHz ($`650\mathrm{km}\mathrm{s}^1`$ at 230 GHz). After that year, the DAS digital autocorrelator system was used in bands of 500 and 750 MHz. Integration times given in Table 2 are typical values used in mapping; central positions were usually observed more than once and thus generally have significantly longer integration times. Values listed are on+off. When sufficient free baseline was available, we subtracted second or third order baselines from the profiles. In all other cases, linear baseline corrections were applied. All spectra were scaled to a main-beam brightness temperature, $`T_{\mathrm{mb}}`$ = $`T_\mathrm{A}^{}`$/$`\eta _{\mathrm{mb}}`$; relevant values for $`\eta _{\mathrm{mb}}`$ are given in Table 2. Spectra of the central positions in both galaxies are shown in Fig. 1 and summarized in Table 3. In Table 2, we have also listed the parameters describing the various maps obtained. All maps are close to fully sampled. In all maps, the mapping grid was rotated by the angle given in Table 2 so that the Y axis coincided with the galaxy major axis. The velocity-integrated maps shown in Figs. 2 and 3 have been rotated back, so that north is (again) at top and the coordinates are right ascension and declination. For IC 342, the map grid origin in the $`J`$=2–1 and $`J`$=3–2 CO maps is within a few arcseconds from the optical centre listed in Table 1, whereas the $`J`$=4–3 and \[CI\] maps find their origin close to the radio centre; the two sets of maps are thus offset from one another by $`\mathrm{\Delta }\alpha `$, $`\mathrm{\Delta }\delta `$ = +9′′, -3′′. For Maffei 2, the grid origin is again close to the optical centre; the radio centre occurs in the maps nominally at $`\mathrm{\Delta }\alpha `$, $`\mathrm{\Delta }\delta `$ = -0.8′′, +6′′.

## 3 Results

### 3.1 CO distribution

In both galaxies, molecular gas as traced by CO is strongly concentrated towards the center. The central molecular source is contained with $`R<`$ 250 pc in the case of IC 342, and $`R<200\times D/2.7`$ pc in the case of Maffei 2, where $`D`$ is the true distance of Maffei 2 in Mpc.

CO line aperture synthesis maps of IC 342 show, at resolutions of 2′′–4′′, considerable structure only hinted at in our maps (Fig. 2). Clearly, the central CO concentration is not toroidal in shape. Rather, the maps published by Sakamoto et al. (1999) and Meier et al. (2000) show individual CO peaks (cloud complexes) distributed along a perspectively foreshortened double spiral – although the foreshortening seems more than expected for an inclination of only 25. The molecular gas enhancement of spiral arms ends at about $`R`$ = 250 pc from the nucleus. The nucleus itself is located at a minimum in the CO distribution. Emission from the high-density tracer HCN offers a very similar appearance (Downes et al. 1992). Our $`J`$=2–1 12CO map shows an elongated source, in which a smoothed version of the high-resolution map by Meier et al. (2000) is easily recognized. Emission in all our CO maps, as well as the \[CI\] map peaks at the position coinciding with cloud complex B in the designation by Downes et al. (1992). This peak becomes more pronounced with increasing $`J`$ level, and in \[CI\]. This is partly the consequence of increasing angular resolution, but we also note that the peak coincides with the strongest concentration of thermal radio continuum emission (Turner $`\&`$ Ho 1983), thus with the region containing the highest spatial density of early type stars. The peak line intensities in Table 3 and the line ratios in Table 4 refer to this maximum. We thus expect them to exhibit properties commensurate with a starburst environment, i.e. those of a photon-dominated region (PDR).

The region corresponding to complexes C and D (Downes et al. 1992), referred to as ‘Eastern Ridge’ by Eckart et al. (1990), behaves differently. It is almost as pronounced as region B in the $`J`$=2–1 transition, but fades rapidly at higher $`J`$ levels. Its molecular gas must thus be cooler on the whole than that of region B. This is, in any case, consistent with the much less impressive appearance of this region in the radio continuum maps by Turner $`\&`$ Ho (1983). There appear to be relatively weak, secondary maxima in the neutral carbon emission at the northern and southern edges of the main CO distribution, in addition to the primary \[CI\] maximum at cloud B.

The major-axis position-velocity map of IC 342 covers a relatively small velocity range, in accordance with the mostly face-on orientation of the galaxy. With increasing $`J`$ level, hence increasing resolution and decreasing beamsmearing, the rotation steepens. In the $`J`$=4–3 CO p-V map the central material covers a range of about 80 km s-1 over about 20′′. The rotation in this region is characterized by a velocity gradient (cf. $`J`$=4–3 CO in Fig. 4 and $`J`$=1–0 CO in Fig. 4 of Sakamoto et al. 1999) d$`V`$/d$`\theta `$ = 7.8 $`\mathrm{km}\mathrm{s}^1`$/′′ , corresponding to d$`V`$/d$`R`$ 2.1 $`\mathrm{km}\mathrm{s}^1`$/pc in the plane of the galaxy.

The structure of the central CO source in Maffei 2 is more complex than Fig. 3 suggests. The galaxy is much more tilted than IC 342, and the integrated velocity images, showing only a resolved but featureless central source, are misleading. Maps of peak temperature $`T_{\mathrm{mb}}`$ rather than integrated-velocity temperature maps show a distinct double-peaked source (see also Hurt et al. 1993). This is in sharp contrast to IC 342, where peak-temperature and integrated-temperature maps show only marginally different morphologies. The deceptively simple appearance of Maffei 2 in Fig. 3 is easily explained by a glance at the major axis position-velocity maps in Fig. 5. The single peak in the integrated map results from the line-of-sight superposition of two distinct velocity components. This is, in fact, the same situation that we found to apply to NGC 253 (Israel et al. 1995), and it is one to watch out for in all high-inclination galaxies. The CO source in Fig. 3 becomes more compact with increasing $`J`$ level, i.e. increasing frequency. As can be seen in Fig. 5, this is mostly the consequence of increasing resolution and decreasing beamsmearing. In the peak-temperature maps (not shown) the separation of the two peaks is only 13′′ (corresponding to 170 pc). The CO distribution shown in Figs. 3 and 5 is characteristic for a toroidal distribution of molecular gas, with a 45 – 60 pc wide gap in the center. However, the ‘S’-shaped morphology of CO in high-resolution (6′′) aperture synthesis maps by Ishiguro et al. (1989) and Hurt $`\&`$ Turner (1991) strongly suggests greatly enhanced emission from molecular gas in inner spiral arms instead.

The effects of beamsmearing are quite noticeable in the p-V maps in Fig. 5. The average velocity gradient increases from 4.5 $`\mathrm{km}\mathrm{s}^1/^{\prime \prime }`$ in the $`J`$=2–1 CO map to about 18 $`\mathrm{km}\mathrm{s}^1/^{\prime \prime }`$ in the $`J`$=4–3 CO map. The latter value corresponds to d$`V`$/d$`R`$ = 1.4 $`\mathrm{km}\mathrm{s}^1`$/pc, i.e. not much different from the gradient in IC 342 (or NGC 6946 and M 83 for that matter – see Israel $`\&`$ Baas 2001). We note that the CO distribution in Fig. 5 just fills the steep leading edge of the major axis HI distribution in Maffei 2 (Fig. 9 in Hurt et al. 1996).

### 3.2 Line ratios

We have determined the intensity ratio of the observed transitions at the centers of both galaxies; all intensities are normalized to those of the $`J`$=2–1 $`{}_{}{}^{12}\mathrm{CO}`$ line. The columns in Table 4 marked ‘Nucleus’ refer to ratios in a 21′′ beam centered on the nucleus, with the exception of the \[CII\]/CO(2–1) ratio which was determined after convolution of the $`J`$=2–1 12CO measurements to the 55′′ beam of the \[CII\] observations, and the $`J=10`$ 12CO/13CO isotopic ratio which represents a weighted mean from a number of determinations at various resolutions. The \[CII\] intensities (Crawford et al. 1985; Stacey et al. 1991), were converted to velocity-integrated temperatures to obtain the line ratios in Table 4. The columns marked ‘Total Center’ refer to the intensities of the central concentration integrated over the source extent in the maps. The apparent decrease in source extent going from the lower $`J`$ levels (lower frequencies) to the higher frequencies is mostly an artifact of limited resolution. When corrected for finite beamwidth, source dimensions at e.g. the $`J`$=2–1 and $`J`$=4–3 transitions are in fact very similar. Nevertheless, the smaller area coverage in the \[CI\] and $`J`$=4–3 CO maps may ahve led to an underestimate of the total intensity. The entries in Table 4 suggest that this may indeed be the case. The $`J`$=1–0/$`J`$=2–1 ratios have relatively large uncertainties, because we could not use observations of our own, but instead have used estimated $`J`$=1–0 intensities reduced to a 21′′ beam from the references given in the table. These ratios are, in any case, close to unity.

The 12CO ratios are quite similar for both galaxies, but the intensity of the \[CI\] and \[CII\] lines relative to CO is significantly lower in Maffei 2 than in IC 342, suggesting stronger PDR effects in the latter. The intensity of the \[CII\] line indicates the presence of both high temperatures and high gas densities in the medium as the critical values for this transition are $`T_{\mathrm{kin}}`$ 91 K and $`n3500\mathrm{cm}^3`$. At the same time, such values must be reconciled with the much lower temperatures and (column) densities implied by the more modest isotopic ratios in the lower CO transitions. These isotopic ratios differ for the two galaxies. For each individual transition the isotopic ratios just agree within the errors, but the overall patterns are clearly different. For Maffei 2, the isotopic ratio increases with increasing $`J`$ level, but for IC 342 the isotopic ratio is minimal in the $`J`$=2–1 transition. These differences suggest significant differences in the molecular gas properties at the positions sampled in both galaxies

## 4 Analysis

### 4.1 Modelling of CO

The observed $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ line intensities and ratios have been modelled with the large-velocity gradient (LVG) radiative transfer models described by Jansen (1995) and Jansen et al. (1994). These provide model line intensities as a function of three input parameters: gas kinetic temperature $`T_\mathrm{k}`$, molecular hydrogen density $`n(H_2)`$ and CO column density per unit velocity ($`N(\mathrm{CO})`$/d$`V`$). By comparing model line ratios to the observed ratios we have identified the physical parameters best describing the actual conditions in the observed source. Beam-averaged properties are determined by comparing observed and model intensities. In principle, with seven measured line intensities, properties of a single gas component are overdetermined as only five independent observables are required. We found that no single-component fit could be made to the data of either IC 342 or Maffei 2.

However, good fits based on two gas components can be obtained. In order to reduce the number of free parameters, we assume identical CO isotopical abundances for both gas components and assign the specific value \[12CO\]/\[13CO\] = 40 (Mauersberger $`\&`$ Henkel 1993; Henkel et al. 1998). We identified acceptable fits by searching a grid of model parameter combinations (10 K $`T_\mathrm{k}`$ 250 K, $`10^2\mathrm{cm}^3n(\mathrm{H}_2)10^5\mathrm{cm}^3`$, $`6\times 10^{15}\mathrm{cm}^2N(CO)/dV3\times 10^{18}\mathrm{cm}^2`$) for model line ratios matching the observed set, with the relative contribution of the two components as the single free parameter. Largely as a consequence of the non-negligible finite errors in the line ratios, solutions are not unique, but rather delineate a range of values in a particular region of parameter space. To a certain extent, variations in input parameters may compensate for one another, leading to identical line ratios for somewhat different combinations of input parameters. We have rejected all solutions in which the denser gas component is also hotter than the more tenuous component, as we consider this physically unlikely on the linear scales observed. From the remainder of solutions, we have selected characteristic examples and listed these in Table 5.

The line ratios observed for IC 342 allow only a rather limited range of gas parameters. One of the two components has well-determined parameters: a column density $`N(\mathrm{CO})`$/d$`V`$ = 6–10 $`\times 10^{16}\mathrm{cm}^2`$, an $`\mathrm{H}_2`$ density of 3000 $`\mathrm{cm}^3`$) and a kinetic temperature of 100 – 150 K. The other gas component either is less dense and equally hot, or denser and much cooler with a temperature of the order of 10 – 20 K. For Maffei 2, the choice of models is less restricted. Possible solutions pair a low-density, hot component and rather large column densities with a much denser, low-column density component of uncertain temperature. Alternatively, both components have a reasonably high density of 3000 $`\mathrm{cm}^3`$, but represent either a low temperature/high column-density or a high temperature/low column-density gas.

A further check on the validity of the model results is provided by the C18O measurements obtained by Eckart et al. (1990). This isotope was not included in our modelling requirements. For isotopic ratios between 12CO/C18O = 250 (model 1) and 12CO/C18O = 225 (model 3), the parameters of our solutions reproduce the observed intensity ratios in the $`J`$=1–0 and $`J`$=2–1 transitions within 10$`\%`$.

### 4.2 Beam-averaged molecular gas properties

The relation between carbon monoxide, neutral and ionized carbon combines the observed C and C+ intensities with the chemical models by Van Dishoeck $`\&`$ Black (1988) which show a strong dependence of the $`N(\mathrm{C})/N(\mathrm{CO})`$ column density ratio on total carbon and molecular hydrogen column densities. In the analysis of \[CI\], we took the kinetic temperatures, $`\mathrm{H}_2`$ densities and filling factors resulting from the CO analysis, and then solved for column density $`N(\mathrm{CI})`$. In practice, the column density of one component is frequently well-determined, whereas that of the other is more or less degenerate. For this reason, we solved for identical velocity dispersions in the two gas components. This procedure we also followed for \[CII\]; in those cases where both model kinetic temperature and gas density are very different for the two components, we have modelled a single \[CII\] component with the higher kinetic temperature and the higher gas density of the two components. This applies specifically to models 5 and 6 (see below). Finally, we have related total carbon (i.e. C + CO) column densities to molecular hydrogen column densities by using an estimated \[C\]/\[H\] gas-phase abundance ratio.

IC 342 has a measured central abundance \[O\]/\[H\] = 2.0 $`\times `$ 10-3 (Vila-Costas $`\&`$ Edmunds 1992, Garnett et al. 1997). We are unaware of any abundance determination for Maffei 2, and therefore use the mean value for the spiral galaxies listed by Zaritzky et al. (1994), \[O\]/\[H\] = 1.7 $`\times `$ 10-3. Then, using the results obtained by Garnett et al. (1999), notably their Figs. 4 and 6, we have estimated carbon abundances \[C\]/\[H\] = 1.7$`\pm `$0.5 $`\times `$ 10-3 and 1.45$`\pm `$0.5 $`\times `$ 10-3 for IC 342 and Maffei 2 respectively. As a significant fraction of carbon is tied up in dust particles and thus unavailable in the gas-phase, we have adopted a fractional correction factor $`\delta _\mathrm{c}`$ = 0.27 (see for instance van Dishoeck $`\&`$ Black 1988), so that $`N_\mathrm{H}`$ = \[2$`N(\mathrm{H}_2)+N`$(HI)\] = A \[$`N`$(CO) + $`N`$(C)\], where A = 2200 and 2500 for IC 342 and Maffei 2 respectively. The numerical parameter A is uncertain by about a factor of two.

In Table 6, we present beam-averaged column densities for both CO and C (Co and C+), and $`\mathrm{H}_2`$ column densities derived under the assumptions just discussed. We also present the total masses estimated to be present in the central molecular concentration calculated from $`L_{tot}`$ listed in Table 3, and the face-on mass densities implied by hydrogen column density and the galaxy inclination. Notwithstanding the differences between the model cloud parameters, the beam-averaged results in Table 6 are rather similar: hydrogen column densities, masses and mass-densities are well within a factor of two from one another. Beam-averaged neutral carbon to carbon monoxide column density ratios range from $`N`$(Co)/$`N`$(CO) = 0.5–0.6 for IC 342 to $`N`$(Co)/$`N`$(CO) = 0.15–0.3 for Maffei 2. These values are well within the range found for other late-type galaxies such ad M 82, NGC 253, M 83, and NGC 6946 (White et al. 1994; Israel et al. 1995; Stutzki et al. 1997; Petitpas $`\&`$ Wilson 1998; Israel $`\&`$ Baas, 2001).

### 4.3 The center of IC 342

Even before the model fitting described in the previous section, the center of IC 342 was identified as the abode of hot and dense molecular gas, first by Martin et al. (1982) and Martin $`\&`$ Ho (1985) who found the presence of gas with kinetic temperatures estimated to be 70 K or higher. Somewhat later, measurements of the high-density tracer HCN led Downes et al. (1992) to suspect the presence of molecular gas at various densities in the range $`n(\mathrm{H}_2)=1503\times 10^4\mathrm{cm}^3`$ and kinetic temperatures in the range $`T_{\mathrm{kin}}`$ = 50 – 70 K. This gas consisted of a few dense cloud complexes and much less dense gas inbetween. Multitransition HCN measurements (Jackson et al. 1995; Paglione et al. 1997) and other molecules (Paglione et al. 1995; Mauersberger et al 1995; Hüttemeister et al 1997; Henkel et al 2000) suggest kinetic temperatures in the range of 50 - 200 K, densities of about 10$`{}_{}{}^{4}\mathrm{cm}_{}^{3}`$ as well as the presence of additional gas components at lower densities. Such temperatures and densities are not unexpected for the site of a significant burst of star formation (cf. Turner $`\&`$ Ho, 1983; 1994). The PDR model calculations by Kaufman et al. (1999), applied to the relative intensities of CO, \[CI\], \[CII\] and the far-infrared continuum of the center of IC 342 likewise suggest a temperature $`T_{\mathrm{kin}}`$ 150 K, a density around 10$`{}_{}{}^{4}\mathrm{cm}_{}^{3}`$ and an ambient radiation field strength log $`G_\mathrm{o}`$ = 2.5 ($`G_\mathrm{o}`$ expressed in units of the Habing Field, i.e. $`1.6\times 10^3`$ erg s-1 cm-2. All three model solutions in Table 5 are consistent with these values, although the higher densities suggested by the other molecular line observations are only provided by model 3. More support for model 3 follows from the C18O aperture synthesis observations carried out by Meier $`\&`$ Turner (2001) who find that the molecular cloud complexes have densities of the order of 1000–3000 $`\mathrm{cm}^3`$ and, more importantly, kinetic temperatures of 10–40 K (see also Meier et al. 2000). We thus accept model 3 as the most probable approximation of the physical conditions in the central molecular gas of IC 342. About a third of the total molecular mass resides in a dense, fairly cold component with an excitation temperature close to the kinetic temperature, and very high optical depths in all observed CO transitions; even the 13CO transitions are marginally optically thick with $`\tau `$ = 1–2. Most of the remaining mass is hot, with Tex = 50 K. The $`J`$=1–0 12CO transition is optically thin, but the higher transitions are optically thick although optical depths are not extremely high.

Obviously, actual conditions will be more varied and complex, but even then we expect the beam-averaged results to be not very different from the ones obtained here. We find a total central gas mass of about 107 $`\mathrm{M}_{}`$ which agrees well with estimates derived from e.g. HCN observations, but is much lower than the mass values quoted in previous CO studies. However, these were all based on assumed CO to $`\mathrm{H}_2`$ ratios, notably the $`\mathrm{H}_2`$ column density to CO intensity ratio $`X`$. The $`\mathrm{H}_2`$ column densities and masses presented in this paper are, in contrast, based on total gas-phase carbon amounts and the actual carbon abundance. The uncertainty in these results is essentially determined by the error in the product $`\delta _C\times `$ \[C\]/\[H\], which is hard to quantify, but which we estimate to be about a factor of two. A mass of $`10^7`$ $`\mathrm{M}_{}`$ is not very high, yet it is about 10$`\%`$ of the dynamical mass with $`R`$ = 130 pc; this renders any value much higher rather unlikely. Model 3 implies for IC 342 a value $`X=3\times 10^{19}`$ $`\mathrm{H}_2`$ mol cm-2/$`\mathrm{K}\mathrm{kms}^1`$, i.e. almost an order of magnitude lower than commonly used $`X`$ values for the Galactic Solar Neighbourhood. However, such a low value is not uncommon for a galactic center. We have found similar values for NGC 6946 and M 83 (Israel $`\&`$ Baas 2001), and such values are by now solidly established for the center of our own Milky Way (Sodroski et al. 1995; Oka et al 1998; Dahmen et al. 1998). As an important consequence of these low values for $`X`$ and beam-averaged $`N(\mathrm{H}_2)`$, the abundances of all molecular species that were derived using erroneously estimated much higher $`\mathrm{H}_2`$ column densities must be revised upwards by almost an order of magnitude. This is particularly important in the case of IC 342. Because it is one of the few galaxies bright in molecular line emission, it has extensively been used as the basis for abundance determinations of a relatively great variety of molecular species.

### 4.4 The center of Maffei 2

In a similar vein, rough estimates of conditions to be found in the center of Maffei 2 follow from observations of HCN, CS and other molecules: $`\mathrm{H}_2`$ densities in the range $`10^410^5\mathrm{cm}^3`$ and probably closer to the former than to the latter, as well as kinetic temperatures of about 100 K (Mauersberger et al. 1989; Hüttemeister et al. 1997; Paglione et al. 1997). The relative intensities of CO, CI, CII and the far-infrared continuum imply, within the context of the PDR model by Kaufman et al. (1999), somewhat higher densities close to 10$`{}_{}{}^{5}\mathrm{cm}_{}^{3}`$, temperatures of about 150 K and a radiation field log $`G_\mathrm{o}`$ = 2.3. These temperatures and densities suggest that models 5 and 6 are more appropriate than models 4 and 7 which lack the high densities that appear to be required. Models 5 and 6 are not fundamentally different. They consist of a rather low density ($`n(\mathrm{H}_2)=100\mathrm{cm}^3`$), hot ($`T_{\mathrm{kin}}`$ = 100–150 K) and very optically thick ($`\tau `$ = 15–20 for $`J`$=1–0) component and a more dense, cooler component of relatively low optical thickness. The actual temperature and density of the latter component is not firmly established. As can be seen in Table 5, temperatures between 20 and 60 K and densities between $`10^4`$ and 10$`{}_{}{}^{5}\mathrm{cm}_{}^{3}`$ can be obtained by varying the relative contributions of the two gas components, requiring only that the ratio $`T_{\mathrm{kin}}^2/n(\mathrm{H}_2)`$ is kept constant. The ionized carbon intensities can only be explained by postulating the presence of a hot and dense gas phase not sampled by the CO observations but rather resulting from extensively photodissociated and ionized gas (interface with HII regions and supernova remnants?). A similar situation was found to apply in M 83 (Israel $`\&`$ Baas, 2001).

In any case, the molecular medium in the center of Maffei 2 appears to be significantly different from that in the center of IC 342. In Maffei 2, the dense phase is warmer and probably denser, and its hot phase is much less dense. Maffei 2 is experiencing a strong central starburst comparable to the one in IC 342. A difference between the two, at least at radio wavelengths, is that Maffei 2 and M 83 have much larger amounts of nonthermal emission than IC 342 (Turner $`\&`$ Ho 1983, 1994). This suggests that the starbursts in Maffei 2 and M 83 are more evolved than the one in IC 342, perhaps explaining the molecular gas differences as well.

As the distance and metallicity of Maffei 2 have not been measured, we have had to assume plausible values. This is an additional source of uncertainty in the $`\mathrm{H}_2`$ column densities, the masses and the mass densities derived for Maffei 2 in Table 6. As the observed HI column densities are much less than the model-inferred total H column densities, $`\mathrm{H}_2`$ column densities, masses and mass densities scale linearly with the actual carbon abundance, whereas the $`\mathrm{H}_2`$ and gas masses scale with the actual distance squared. Because it is unlikely that the distance of Maffei 2 differs much from the adopted value of 2.7 Mpc, the error in the derived mass due to this uncertainty is not more than about a factor of two. The derived $`X`$ value is subject to uncertainties in the carbon abundance but not to those in the distance. The $`N(\mathrm{H}_2)`$ values listed in Table 6 imply for Maffei 2 an $`\mathrm{H}_2`$/CO conversion ratio $`X=23\times 10^{19}`$ $`\mathrm{H}_2`$ mol cm-2/$`\mathrm{K}\mathrm{kms}^1`$, so that the same comment made for IC 342 also applies to Maffei 2.

## 5 Conclusions

1. Maps of the central arcminute of the nearby starburst galaxies IC 342 and Maffei 2 in various transitions of $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$, and in \[CI\] confirm the compact nature of the central molecular gas emission in both galaxies. Most of this gas is within a 200 parsec from the nucleus. In both galaxies, the molecular gas seems to reside in bright inner spiral arms, rather than a disk or torus.
2. Relative 12CO line intensities observed in matched beams are virtually identical in IC 342 and Maffei 2. However, relative 13CO intensities, as well as the intensities of the \[CI\] and \[CII\] lines differ significantly, indicating different physical conditions for the molecular gas in the two galaxies. The observed line ratios require modelling with a multi-component molecular medium in both galaxies.
3. In IC 342, a dense component with $`n(\mathrm{H}_2)10^4\mathrm{cm}^3`$ and $`T_{kin}`$ 10 K is present together with a less dense $`n(\mathrm{H}_2)3\times 10^3\mathrm{cm}^3`$ and hotter $`T_{kin}`$ 100 K component. Total carbon column densities are about 1.5 times the CO column density. The modelling solution for IC 342 is reasonably well-defined.
4. In Maffei 2, the parameters of the two components minimally required are less clearly defined. The more tenuous component has $`n(\mathrm{H}_2)10^2\mathrm{cm}^3`$ and $`T_{\mathrm{kin}}`$ 100 - 150 K, whereas the denser component has $`n(\mathrm{H}_2)10^410^5\mathrm{cm}^3`$ and $`T_{\mathrm{kin}}`$ = 20 – 60 K. Maffei 2 seems to be more affected by CO dissociation as it has a C/CO ratio of two to three.
5. In both starburst centers a significant fraction of the molecular mass (about half to two thirds) is associated with the hot PDR phase.
6. With an estimated gas-phase \[C\]/\[H\] abundance of order 4 $`\times 10^4`$, the centers of NGC 6946 and M 83 contain within $`R`$ = 0.25 kpc similar (atomic and molecular) gas masses of about $`10^7\mathrm{M}_{}`$. Peak face-on gas mass densities are typically 70 $`\mathrm{M}_{}`$ pc-2 for IC 342 and 35 $`\mathrm{M}_{}`$ pc-2 for Maffei 2, but the results for the latter are subject to relatively large errors caused by uncertainties in its distance and its carbon abundance.

###### Acknowledgements.

We are indebted to Ewine van Dishoeck and David Jansen for providing us with their detailed radiative transfer models. We thank the JCMT personnel for their support and help in obtaining the observations discussed in this paper.