A new method to calculate normal modes

interpreted with the aid of an acoustic potential which characterizes vertical propagation of sound waves. In addition to the efﬁciency in the convergence to the eigenfrequencies, numerical tests show strong numerical stability of the method. It stems from the stability in the relaxation method because of the similarity in algebraic structure. For those reasons, we propose the method as an efﬁcient way in calculating synthetic sesimograms, barograms and ionograms for recently observed phenomena relating with coupling between the solid earth, the oceans and the atmosphere.

ionospheric sounding technique. Those phenomena all indicate that we must include the atmosphere in the total system of the oscillations to describe elasto-gravitational phenomena correctly.
Once we include the atmosphere in the earth model, we face troublesome problems because the structure of the atmosphere varies daily and annually. In fact, Nishida et al. (2000) reported that amplitudes of 0 S 29 exhibit annual variation of about 40 per cent although amplitudes of the other modes located far from the acoustic branch vary by about only 10 per cent. This phenomenon may be attributed to a change in the resonant structure between the solid mode and the acoustic mode. To study such a phenomenon, we must calculate modal solutions for a variety of earth models with different atmospheric structure. The efficiency of calculation becomes a significant issue in such a situation. Another problem associated with the atmosphere is leakage of wave energy into space. The atmosphere has neither rigid lid nor free surface. Roughly speaking, waves having shorter wavelengths than the pressure scale height easily propagate into space and are dissipated there. Eigenfrequency and eigenfunctions of such oscillations must have complex numbers because of the energy leaks even if a perfectly elastic earth model is used. Searching for eigenfrequencies in the complex frequency plane is another difficulty in the mode calculations.
On the other hand, numerical calculation of normal modes of the Earth is one of the classical problems in seismology. Many methods have been developed and are explained in many textbooks (e.g. Takeuchi & Saito 1972;Woodhouse 1988;Saito 1988). These methods are based on numerical integration of the differential equations of motion (shooting methods). The methods are efficient for searching real eigenfrequencies of non-dissipative cases but have some troubles to seek complex eigenfrequencies of dissipative or leaky oscillations. In addition, to assure accuracy of calculations, we must calculate the minor of the integral matrix (Gilbert & Backus 1966). In solar seismology, the Henyey-type relaxation method is used to calculate acoustic and gravity modes of the Sun and stars (e.g. Henyey et al. 1964;Unno et al. 1989). The latter method is easily extended to complex cases and is generally more stable for eigenvalue problems of a self-gravitating gaseous sphere whose density and pressure vary exponentially near its surface. The efficiency of the method, however, depends on initial guesses of eigenfrequencies and eigenfunctions.
Due to the above difficulties in the direct methods of convergence in the complex plane, Lognonné et al. (1998) and Artru et al. (2001Artru et al. ( , 2004) used a variational method to calculate normal modes of the coupled system of the solid earth and the atmosphere. A nuisance in this method is the calculation of many basis functions which are normal mode eigenfunctions of the earth model with atmosphere but with free boundary conditions, and it needs mapping a function with free boundary condition into a function with an open boundary condition. And couplings between modes in the variational calculation is usually restricted to several dozen modes. It means that the degree of freedom expressing the synthesized eigenfunctions is usually less than the degree of freedom of the model space (i.e. the total number of grids used for expressing the earth model and basis functions). If one is interested in a few particular modes such as 0 S 29 and 0 P 29 , the calculation of a lot of basis functions is expensive in computational costs. More direct solvers for normal modes are desirable in such a problem.
Here, we propose a new method to overcome the deficiency of the above methods. This method is formulated with a similar treatment of the global matrix governing oscillations as that in the Henyey-type relaxation method, but it solves the problem more directly like the Runge-Kutta method without initial guesses of eigenfunctions. In Section 2, we will describe the mathematical basis of the method. In Section 3, some results of numerical tests are shown, and good convergence of complex eigenfrequencies and strong stability in our method are shown. In Section 4, we discuss the merit of our method in comparison with the other methods. Finally, we conclude the proposed method is effective to calculate normal modes including solid modes, acoustic modes and gravity modes of the total system of the solid earth, the oceans and the atmosphere. The method can be applied to other planets, satellites and stars having a mixed structure of solid, oceans and atmosphere such as Jupiter, Venus, Mars and extrasolar planets.

Equations of motion
The equations of motion of an elasto-gravitating spherically symmetric transversely isotropic body are described in many textbooks (e.g. Takeuchi & Saito 1972;Dahlen & Tromp 1998). Here we do not describe them in detail and we abbreviate them as d dr where r is radius, y is the vector composed of eigenfunctions {U (r ), X (r ), V (r ), Y (r ), φ(r ), ψ(r )} for spheroidal oscillations and {W (r ), Z (r )} for toroidal oscillations, A is the motion matrix given by the coefficient matrix of eqs (82) in a solid medium and (86) in a fluid medium for spheroidal oscillations and eq. (76) for toroidal oscillations of Takeuchi & Saito (1972). The definition of the symbols U, V and W for displacement and X , Y and Z for vertical traction are the same as those in Woodhouse & Dahlen (1978). φ is the perturbation of gravitational potential and is −y 5 , and ψ is −y 6 of Takeuchi & Saito (1972). Thus terms in the fifth and sixth rows of A have different signs from those of Takeuchi & Saito (1972) in our convention. In a fluid layer, shear is lost (Y = 0) and the balance of the horizontal momentum gives where ω is the angular frequency of oscillations, ρ is the density in the fluid, g is the gravitational acceleration and r is the radius. The rank of the matrix A is reduced to 4 in fluid media.
In the case of models having anelastic properties such as PREM (Dziewonski & Anderson 1981), the elastic constants used in the spherically symmetric transversely isotropic model become complex numbers and can be given by where i is the imaginary unit, the subscript "1" means values taken at a period of 1 s, T is the period of oscillations, Q κ and Q µ are quality factors of bulk compression and shear, respectively, and The terms including ln T express the physical dispersion adopted by Dziewonski & Anderson (1981) based on the argument of Liu et al. (1976).
If the body has an atmosphere, it is convenient to use different variables {r √ ρU, r p / √ ρ, r √ ρϕ, r √ ρψ} since pressure and density decrease exponentially with altitude. They are slightly different from those used in Lognonné et al. (1998) and are natural since is proportional to the kinetic wave energy density propagating upward, whose magnitude vary mildly while the eigenfunctions themselves may vary exponentially. Matrix A for those variables becomes where N B is the Brunt-Väisälä frequency, c is the sound velocity, H p and H ρ are scale heights of pressure and density, respectively, γ is the adiabatic exponent, l is the angular degree and G is the gravitational constant. p is the Eulerian perturbation of pressure and relates to the normal traction X as follows: The derivation of matrix A in eq. (4) is straight from eq. (86) of Takeuchi & Saito (1972) which describes the equations of motion in fluid media. The effect of the viscosity is important on description of waves in the ionosphere (Artru et al. 2001). In this paper, however, we neglect the viscosity since we will focus on the long period acoustic modes for which the energy leakage due to open boundary condition is a major factor of dissipation of those modes. But note that the results of complex eigenfrequencies in this paper must be slightly modified by the viscosity.
In addition, we scale the constants and the variables in the equations of motion to be non-dimensional moderate values. The length is scaled by R that is the radius of the model, mass is by M that is the total mass of the model and angular frequency is by ω 0 = √ g 0 /R where g 0 is the surface gravity given by GM/R 2 . With these scales, density is scaled by ρ 0 = 3M/4π R 3 , pressure and stress are by p 0 = ρ 0 g 0 R, gravitational potential is by ϕ 0 = c 2 0 = g 0 R where c 0 is the scale of velocity, and 4π G then becomes 3. The last relation is the direct consequence of the definitions of ρ 0 and g 0 . All constants and variables are so scaled by the above quantities in the numerical calculations described in Section 3. Note that we use the same scales in an atmosphere where the scaled density ρ/ρ 0 does not necessarily take a moderate value and terms proportional to ρ can be neglected at a very high altitude.
To solve eq. (1) numerically, we approximate it by difference equations in the second order where y n and A n are evaluated at radius r = r n (n = 1, · · · , N ) and r = r n+1 − r n is the grid interval, in the same manner as Unno et al. (1989). Here we denote the total number of grid points by N . Rearrangement of the terms gives C n y n + D n+1 y n+1 = 0 where In this method, equations in the form of (5) are used with the following boundary conditions to obtain normal mode solutions.

Boundary conditions
Here we describe somewhat detailed expressions of boundary conditions in our method, with which the equations of motion can describe normal mode solutions completely.

Inner boundary conditions
Eqs (95) for toroidal modes and (98) and (100) for spheroidal modes of Takeuchi & Saito (1972) are adopted for the inner boundary conditions. They are the elementary solutions of an isotropic homogeneous elastic-gravitational sphere. The modal solutions can be written as a superposition of the elementary solutions y (i) which are regular at the centre of the body, for spheroidal modes. The coefficients a i are unknown and can be determined by the outer boundary conditions in the shooting method. Here we treat it in a different way. They can be eliminated in the above equation. For example, an elimination gives  where U (i) , . . . , ψ (i) are the components of the ith elementary solution. Then the boundary conditions can be written as where B is a 3 × 6 matrix and is given by for the above case. When we put the inner boundary into a liquid layer, the shear is lost (Y = 0) and the differential equation for Y becomes trivial. There are only two elementary solutions and the size of B becomes a 2 × 4 matrix.
For toroidal oscillations, the condition is very simple and is where W (1) and Z (1) are the components of the elementary solutions, or when the inner boundary conditions are set at a fluid-solid interface like the core-mantle boundary. In both cases, the conditions can be expressed by the same form as eq. (6) and the boundary matrix B becomes a 1 × 2 matrix.

Outer boundary conditions
If the outer boundary conditions are free surface conditions, then the boundary conditions can be written in the same form as eq. (6) with This merely means the free conditions If the free surface is covered by an ocean, then they are If the outer boundary are situated in the atmosphere, the boundary conditions are somewhat complicated but still can be written in the same form as eq. (6). The elementary solutions in the upper atmosphere can be approximately obtained using the characteristic solutions of matrix A of eq. (4) assuming ρ → 0 and parameters appearing in the elements of A are nearly constants (a similar discussions can be found in Unno et al. 1989;Lognonné et al. 1998). Some details of the solutions are described in the appendix.

Internal boundary conditions
Realistic earth and planetary models have inner boundaries where density and elastic properties vary discontinuously. The continuation of solutions is trivial for the solid-solid interfaces and is where the subscripts -and + denote values on the lower and upper sides of the interfaces of discontinuities, respectively. I is the identity matrix. In the case of a fluid-solid boundary, the continuation conditions are  and in the case of a solid-fluid boundary, they are  Note that the above equations have no constraints on V and Y in the fluid layer. But this is irrelevant since Y vanishes and V is not independent (see eq. 2). The number of conditions to determine the eigenfunctions is neither too high nor too low in the fluid layer, of course.
In the case of an interface with an atmosphere, the continuation conditions become slightly complicated because of the different variables in the atmosphere and it is  for an ocean-atmosphere interface. Again V and Y do not matter in the ocean layer and the atmospheric layer if viscous shear can be neglected. The extension to the case of a viscous fluid is straightforward and is not discussed in this article.
If we give different grid numbers both sides of the discontinuities, say r n = r − and r n+1 = r + where r − = r + , then the above internal boundary conditions have the same form as the difference equations of motion (see eq. 5). Hence we make no distinction between the internal boundary equations and the equation of motion as to their notation in the following discussion.

Method for solutions
The difference eq. (5) and the boundary conditions (6) describe the complete system of oscillations  where submatrices B 1 and B N denote the inner boundary conditions and the outer boundary conditions, respectively, C n and D n denote the equations of motion, and they are functions of frequency ω. Blank parts in the global matrix are all filled with zero. If a given frequency ω is an eigenfrequency, the determinant of the global matrix is zero and non-trivial solutions y n can exist. If ω is not an eigenfrequency, the global matrix can be inverted and we have only the trivial solution y n = 0. As doing in the Henyey type of relaxation method (e.g. Unno et al. 1989), we express the global matrix using different constituent sub matrices as  where S 1 is composed of B 1 and the upper half of C 1 , the upper and lower halves of Q 1 are a block filled with zero and the upper half of D 2 , respectively, and the upper and lower halves of P 2 are the lower half of C 1 and filled with zero respectively. We define, in the same manner, where superscripts 'u' and 'l' mean the upper and lower halves of the given sub matrices. S N is composed of the lower half of D N and B N . Hence, sub matrices P n , S n and Q n are all square matrices of 6 × 6 in solid layers and 4 × 4 matrices in the fluid layers. Note that the matrix in eq. (8) is just the matrix in eq. (7) and the difference between them is only in expression.
The global matrix is a square matrix of 6N × 6N for spheroidal oscillations of a solid sphere. Solving the eigenvalue problem requires high numerical costs and is generally difficult to solve if the total number of grids N is large. However, it is a banded matrix whose non-zero values are confined to submatrices P n , S n and Q n which are calculated for a given ω, and we can recursively solve eq. (8) as follows: The first row of the block matrices in eq. (8) is from which we have y 1 = −S −1 1 Q 1 y 2 . By substituting this into the second row of eq. (8), we obtain a relation between y 2 and y 3 , that is −P 2 S −1 1 Q 1 + S 2 y 2 + Q 2 y 3 = 0, and we get Proceeding with the same process, we obtain recursive relations where R n = −X −1 n Q n and X n = P n R n−1 + S n with the initial term given by where which is calculated by the recursive relation of X n and R n . If X N has its inverse for the given ω, then y N = 0 and all y n are zero owing to eq. (9). On the other hand, if the determinant of X N is zero and X N has no inverse for the given ω, eq. (10) has a non trivial solution of y N and all y n are evaluated by (9) recursively. In this case, ω must be one of the eigenfrequencies of the system and y n give the eigenfunctions at r = r n . Therefore the characteristics of sub matrix X N determine the characteristics of the original global matrix, and the problem becomes one of seeking ω for which det X N = 0. Note that the numerical cost in this recursive calculation is proportional not to N 2 but N since the number of recursive steps is only N . We give some comments on the matrix R n . From eq. (5), we obtain y n = −C −1 n D n+1 y n+1 ≡ R n y n+1 . Although this equation looks like eq. (9), they are very different. As seen from the definition of C n and D n , they can be inverted. Accordingly R n has an inverse. On the other hand, R n in eq. (9) has no inverse because the determinant of Q n in R n is zero. It can propagate y n downward only while R n can propagate them both downward and upward. Eq. (9) requires y n and y n+1 to be the solution satisfying both the inner boundary conditions and the equations of motion whereas the above equation makes them to be a solution of the equations of motion alone. The one-way propagation in our method is not deficient but produces numerical stability. The stability in use of R n will be shown numerically in the next section.
Next we consider how we seek eigenfrequencies of the system. Here, we denote an eigenvalue as a sum of a trial value ω 2 0 and its correction δω 2 , i.e.
Then eq. (10) can be expanded up to first order and we have This is a generalized eigenvalue problem in which δω 2 and y N is the eigenvalue and the eigenvector of the equation, respectively. We can expect that the correction δω 2 so obtained is significant if the trial value is close to true eigenfrequency. Eqs (11) and (12) can be solved iteratively until |δω 2 /ω 2 | < ( 1) by substituting ω 2 0 + δω 2 into ω 2 0 in the next step. The partial derivative of X N with respect to ω 2 can also be calculated in a recursive way using ∂X n ∂ω 2 = ∂P n ∂ω 2 R n−1 + P n ∂R n−1 ∂ω 2 + ∂S n ∂ω 2 and ∂R n−1 ∂ω 2 = X −1 at the same time when we calculate X n and R n . We need to keep only values of the partial derivatives of X n at the current step in the computer's memory. On the other hand, we need memory to keep all R n to calculate eigenfunctions using (9). The required memory is also proportional only to N .

N U M E R I C A L E X A M P L E S
In the previous section, we developed a new method to calculate eigenfrequency and eigenfunctions of normal modes of a spherically symmetric elasto-gravitating sphere. In this section, we describe some numerical examples which show the performance of this method.
In Table 1, we compare eigenperiods T, Q values and group velocities u of several spheroidal modes of PREM (Dziewonski & Anderson 1981) calculated by our method with those values calculated by DISPER80 (Saito 1988), and displacement eigenfunctions of those modes are shown in Fig. 1. In both methods, we used the same grid model for PREM and grid intervals ( r = r n+1 − r n ) are taken to be about 1 km (i.e. the total number of grids N is 6381). In this test calculation, we neglect the atmosphere. Our T and Q values are obtained by , respectively. Group velocity u is calculated by where I 1 and I 3 are defined in eqs (180) and (188) of Takeuchi & Saito (1972), respectively. In our notation, they are and Note that the eigenfrequency and eigenfunctions of a mode are complex because of the anelasticity of PREM in our method. So we take the real parts in the calculations of T and u. We obtained good agreement between results of both methods. In DISPER80, attenuation (Q) is treated as a small perturbation in a perfectly elastic earth model. On the other hand, our method directly includes anelasticity in the formulation of the eigenvalue problem. To see the overall accuracy of the calculation, we can calculate a relative error defined by ε = I 2 ω 2 I 1 − 1 , Table 1. Periods T, quality factors Q and group velocities u of several fundamental spheroidal modes with different angular degrees of 2, 20, 40, 80 and 160, and those of several overtone modes (n = 1, 3, 7, 13 and 25) with angular degree 2 calculated by the proposed method are compared with those (T D80 , Q D80 and u D80 ) calculated by DISPER80 (Saito 1988). In the method, Q D80 is calculated by the 1st order perturbation theory (e.g. Takeuchi & Saito 1972) using real eigenfrequency and real eigenfunctions. T RP and u RP are periods and group velocities calculated by our method without imaginary parts of elastic constants and Q RP is calculated by the first-order theory using the real eigenfunctions. Note that we use PREM as an earth model without atmosphere in this calculation.
Mode  for the free surface boundary condition (see eq. 181 of Takeuchi & Saito 1972), where As for the calculation of Table 1, the relative errors of the fundamental spheroidal modes 0 S 2 , 0 S 20 , 0 S 40 , 0 S 80 and 0 S 160 are 3.7 × 10 −8 , 1.4 × 10 −7 , 1.4 × 10 −6 , 8.4 × 10 −6 and 3.6 × 10 −5 , respectively, and those of the overtone modes 1 S 2 , 3 S 2 , 7 S 2 , 13 S 2 and 25 S 2 are 3.4 × 10 −8 , 2.7 × 10 −8 , 8.7 × 10 −8 , 1.7 × 10 −8 and 1.5 × 10 −7 , respectively. The fundamental spheroidal mode with large angular degree 0 S 160 has a somewhat worse ε due to the relatively coarse grid intervals in comparison with the skin depth of the wave energy, which is confined to the shallower part of the Earth. Further calculations for 0 S l with l ≤ 10 000 show that the error is roughly proportional to l 2 . So, with increasing in l, we need finer grid intervals (∝ l −2 ) in the shallower part of the Earth. Since the skin depth is roughly proportional to l −1 , the cost of calculation of 0 S l proportionally increases roughly with l or the frequency f by omitting calculation in the deeper part of the model. On the other hand, the error of overtone modes proportionally increase with the overtone number n or f of modes with a fixed l since r /λ becomes large (λ is the vertical wave length). Therefore, the cost of calculation of overtones also roughly increases with f under the condition of a given accuracy in the calculations. In Fig. 2, we show the behaviour of the correction δω of the eigenfrequency ω for the fundamental mode 0 S 29 (see eq. 11). Each arrow in the plots is 0.3 × δω calculated for ω at its root position. If the corrections δω are exact, the head of an arrow locates at the eigenfrequency of 0 S 29 . To prevent overlaps of arrows at the eigenfrequency, we multiplied δω by 0.3 for better visualization. In Fig. 2(a), we observe quite better behaviour of δω around 0 S 29 and 1 S 29 except near the saddle point between them. Arrows in the vicinity of 0 S 29 point straight to its eigenfrequency (Fig. 2b). Dashed lines with small solid squares are two actual paths of convergence to 0 S 29 with starting points indicated by large solid squares near the saddle point. In both cases, we can reach the convergent point within 8 steps with a given precision of = In both cases, we can reach the eigenfrequency within eight steps with a precision of = |δω 2 /ω 2 | < 10 −12 . (b) An expansion view of (a) in the vicinity of 0 S 29 . Its complex eigenfrequency is shown by a solid circle. Linearity of δω near 0 S 29 is quite good. All the arrows point to the eigenfrequency straightly. (c) The same as (b) but the earth model has an atmosphere (PREM+NRLMSIS-00). The left solid circle is the eigenfrequency of the fundamental acoustic mode 0 P 29 and the right one is 0 S 29 . In this case, the two modes are nearly degenerate. The behaviour of δω is somewhat complicated. However, arrows clearly point to the two complex eigenfrequencies and we can still distinguish 0 S 29 from 0 P 29 . |δω 2 /ω 2 | < 10 −12 . In Fig. 2(c), we show the result for a model with an atmosphere. The atmospheric model is a globally averaged one over longitude, latitude and local time based on NRLMSIS-00 (Picone et al. 2002) calculated for the atmosphere in July (see Fig. 5a). In this paper, the grid intervals r n = r n+1 − r n are also 1 km in the atmosphere. We put the open boundary condition at an altitude of 800 km. With the realistic atmosphere, the eigenfrequency of the fundamental acoustic mode 0 P 29 locates beside 0 S 29 . The behaviour of δω seems to be somewhat complicated in this case. However, arrows (δω) in the plot consistently point to the eigenfrequencies of 0 P 29 and 0 S 29 , and we can distinguish one from the other. In fact, starting from the same points as those in Fig. 5(a), we can reach the eigenfrequency of 0 S 29 within nine steps with the same precision. If we start from a different point to the left, we can easily reach the eigenfrequency of 0 P 29 similarly.
Next, we check the validity of the open boundary condition explained in the appendix. We calculate 0 S 29 and 0 P 29 for models having different heights of atmosphere. In Fig. 3, we compare displacement eigenfunctions of them for a model with a 400 km atmosphere and those for a model with an 800 km atmosphere. Both atmospheres have the same structure below 400 km. Acoustic waves with the same frequency f and wave number l as 0 P 29 are relatively well trapped in the lower 100 km part of the atmosphere (see Fig. 5b). For those waves, eigenfunctions are consistently retrieved for the different heights of the atmosphere. These results show the validity of the open boundary condition we adopted. However, modes not well trapped in the atmosphere such as 2 P 29 and 3 P 29 whose Q values are about 10 are somewhat sensitive to the height of model atmosphere. For example, the real part of the eigenfrequency of 3 P 29 is 5.10 mHz for the 400 km atmosphere but it is 5.18 mHz for the 800 km atmosphere. The difference of them becomes a few per cent on this leaky mode. On the other hand, the difference of the well trapped mode 0 P 29 whose Q is about 181 is less than 0.1 per cent. In Fig. 4, we compare eigenfunctions of the tsunami mode O 29 , the fundamental acoustic mode 0 P 29 , overtones of acoustic modes 1 P 29 , 2 P 29 and 3 P 29 for the atmosphere of 800 km height. Eigenfrequencies of those are 0.124, 3.688, 4.442, 4.876 and 5.183 mHz, and Q values are 3347, 181, 22, 10 and 12, respectively. The starting grid of those calculations is the same as Fig. 1, that is the second grid at r ∼ 1 km in the inner core. Eigenfunctions of those modes in the deeper part of the Earth are negligibly small but the proposed method can calculate them very stably from the core of the Earth to the top of the atmosphere. As an extreme case, eigenfunctions of 0 S 10000 can be calculated from the point near the centre of the Earth, of which eigenfrequency is 0.3758 Hz and Q is 14365. The high Q is due to the ocean and the skin depth of about 5 km. So far, we have shown the method can be applicable to calculations of modes for spherically symmetric earth models with a realistic atmosphere. One-dimensional (1-D) structure of the real earth gradually varies in the horizontal direction and in time (see Fig. 5). In such a situation, local eigenfrequencies and eigenfunctions are necessary to describe seismic waves (e.g. chapter 16 of Dahlen & Tromp 1998). As an example of such modal calculations, we examine the effect of annual variations in the 1-D structure of the atmosphere on complex eigenfrequencies. The temperature of the atmosphere slightly changes with season. The acoustic property and complex eigenfrequencies also vary. In such a situation, the proposed method has an advantage. Since the structure gradually changes, variations in the complex eigenfrequencies are also gradual. If one know or can calculate local eigenfrequencies for a 1-D structure, they can be used as good guesses of trial values for the complex eigenfrequencies for another successive 1-D structures and they are usually expected to converge in a few iterations. In Fig. 5(b), we compare acoustic potentials of the atmospheres in different seasons for angular degree 29. See discussion in the appendix on the definition of the acoustic potential. In Fig. 5, eigenfrequencies and Q values of 0 S 29 and 0 P 29 calculated for the structure of each month are shown. Those of the solid modes are almost unaltered all the year round as expected. The wave energy of the solid mode is, of course, mainly in the solid part of the Earth that is unchanged. On the other hand, the eigenfrequency of 0 P 29 whose energy is mostly in the atmosphere varies annually. The maximum occurs in the summer when the eigenfrequencies of the two modes are close. This is because of From the left, tsunami mode O 29 , the fundamental acoustic mode 0 P 29 , the overtones of acoustic modes 1 P 29 , 2 P 29 and 3 P 29 . In the atmosphere, the eigenfunctions multiplied by √ ρ/ρ a × r/R are plotted as in Fig. 3. These modes have most energy in the ocean and/or the atmosphere. [Hz] [Hz] [Hz] Altitude [km] 0 P 29 Figure 5. Atmospheric structure used in this study. (a) temperature, pressure, sound velocity and potential walls of acoustic waves ( f a ) and gravity waves ( f g ) with l = 29 are plotted as functions of altitude (z = r − R), respectively. (b) Comparison between acoustic potentials of the January atmosphere and the July atmosphere. The potential barrier of January around 100 km is higher and thinner than that of July. The inset is an expansion of the lower 200 km of the potentials where we compare the acoustic potentials of degree 29 of January, April, July and October. The level of eigenfrequency of 0 P 29 is shown as the horizontal dotted line. The acoustic waves with the frequency are partially reflected near the mesopause because of the higher barrier around 100 km. the protrusion of the acoustic potential around 90 km in July which shorten the acoustic path length between the surface of the ocean and the turning point nearby. The behaviour of Q seems to be somewhat complicated. It shows biannual variations. In the winter, the potential barrier near 100 km is thinner but higher. In the summer, the barrier becomes somewhat lower but grows slightly fat. When the barrier becomes higher and/or fatten at a given level, reflectivity of acoustic waves due to the barrier becomes larger and Q of the trapped oscillations becomes larger. On the other hand, the barrier becomes somewhat lower and thinner in the spring and the autumn, and Q of the acoustic mode becomes lower in those seasons. In this way, the proposed method can calculate complex eigenfrequencies and eigenfunctions efficiently for models having gradual changes in 1-D structure.
These variations in f and Q of 0 P 29 might explain the observation of Nishida et al. (2000) in which they reported the continuous excitation of 0 S 29 annually varies with an amplitude of about 40 per cent. If this is valid, the excess amplitudes observed in 0 S 29 is partly due to a contribution of the acoustic mode whose responses on the ground can vary annually. The proposed method can be useful in the excitation problem of the continuous free oscillations. We will discuss this issue in a separate paper.

D I S C U S S I O N
When we approximate the structure of the Earth and planets by a finite number N of discretized spherical shells, the eigenvalue problem becomes one to seek MN eigenvalues for a given horizontal wave number (angular degree) l, where M is 6 for spheroidal oscillations 2 for toroidal oscillations. To assure accuracy in the calculation, N must be large, i.e. we need fine structure in the description of models. Then the eigenvalue problem of degree MN becomes one hard to be solved since we must solve a polynomial equation of order MN . In this method, the global matrix has a banded structure so that all non-zero elements concentrate in a 3M band centred at the diagonal elements. Taking advantage of this, we obtain a much smaller eigenvalue problem whose eigenvalue is the correction δω 2 of the assumed eigenvalue ω 2 . This small eigenvalue problem can be solved quite easily. However, δω 2 is not good if assumed value of ω is far from one of eigenfrequencies. In this case, we must iterate the calculation until = |δω 2 /ω 2 | becomes sufficiently small by substituting ω 2 + δω 2 into ω 2 in the next step. In Fig. 2, we showed the behaviour of δω 2 around eigenfrequency of 0 S 29 . They are quite gentle and the convergence to the eigenfrequency from the frequencies near the saddle point with 1 S 29 can be achieved within eight iterations. When we add an atmosphere to the model, the fundamental acoustic mode 0 P 29 appears in the vicinity of 0 S 29 . Even in this complicated case, the method can reach the eigenfrequencies of 0 P 29 and 0 S 29 within nine steps. This good behaviour of convergence to a complex eigenfrequency is one of the remarkable points of this method, while the convergence in the shooting method is quite slow (e.g. Watada 1995).
The proposed method is formulated based on a similar treatment of the global matrix governing the system of oscillations to that in the Henyey-type relaxation method (e.g. Unno et al. 1989) used in solar seismology. In the latter method, initial guesses of eigenfunctions are also needed in the calculation. The mathematical structure of the method is essentially as follows. Eq. (8) can be symbolically expressed by where H is the global matrix and Y is the vector whose components are a series of eigenfunctions. With initial guesses of eigenfrequency ω 2 0 and eigenfunctions Y 0 and their corrections δω 2 and δY, we can write it as H ω 2 0 + δω 2 [Y 0 + δY] = 0, and then, expanding up to first order, we obtain δω 2 ∂H ∂ω 2 Y 0 + HδY = −HY 0 . With an appropriate normalization condition that one of elements of Y is unity, the above equation can be solved for δω 2 and δY. It is not an eigenvalue problem but merely solving linear algebraic equations. If the normalization condition is set on the vertical eigenfunction at the last grid, the coefficient matrix of the above linear algebraic equations is almost H except for the column corresponding to the normalization which are replaced by ∂H ∂ω 2 Y 0 . The calculation of δω 2 and δY can be performed in a recursive way like eq. (9) with extra terms, where r n is a relating matrix with the residuals HY 0 and ∂H ∂ω 2 Y 0 . Details in the method are described in Unno et al. (1989) and are not repeated here. We stress the point that the algebraic structure of our method is very similar to that of the Henyey-type relaxation method. The relaxation method is relatively better than the shooting method in the case of maintaining a decaying exponential in the presence of growing exponentials (Press et al. 1992). The stability observed in our method can be attributed to the stability of the relaxation method since the solvability depends on the structure of global matrix H or the property of R n in both methods. In addition, R n in eq. (9) and a product n R n has zero determinant and has no inverse although it looks like a propagator matrices. The propagation by R n is restricted in one direction. This restriction is not deficient and is used to avoid numerical inaccuracy that emerges in the shooting method (Gilbert & Backus 1966) and produces good stability in our method. Shibahashi & Osaki (1981) used a complex number D(ω) which is essentially same as det X N (see eq. (10)) to find non-adiabatic eigenmodes whose eigenfrequencies have large imaginary parts for stars on the left-hand side of the Cepheid instability strip in the HR diagram. When we calculate D(ω) along a closed curved line around a complex eigenfrequency, the trace of D(ω) becomes a closed curve around the origin. We can find the location of the eigenfrequency of such a mode to examine whether it encloses the origin or not. In our method, we do not use det X N itself. Instead, we expand X N with respect to ω 2 and obtain an eigenvalue problem of δω 2 . Since δω 2 is the correction to a given ω 2 , we can reach the complex eigenvalue automatically with several iterations as long as the behaviour of δω 2 is smooth. Accordingly our method requires less trials than the method of Shibahashi & Osaki (1981). Lognonné et al. (1998) used a variational method to calculate normal modes of the coupled system of the solid earth and the atmosphere. It is another way to calculate complex eigenfrequencies and eigenfunctions directly. However, this method requires the calculation of many basis functions that are the normal mode functions of an earth model with an atmosphere but with free boundary conditions, and it needs mapping a function with free boundary condition into a function with open boundary condition. If one is interested in a few particular modes such as 0 S 29 and 0 P 29 , calculation of many basis functions and solving variational problem are computationally expensive. Especially if we want to discuss effects of seasonal variation of the atmosphere, we must solve basis functions for each atmospheric model. As we saw in Fig. 6, our method can obtain the complex eigenfrequency and the eigenfunctions of a target mode without any calculations of basis functions for a variety of atmospheres and remapping eigenfunctions to the open boundary conditions. Similarly, 1-D structures of the atmosphere and the solid earth varies horizontally too although horizontal variations is generally less than vertical ones. In such a situation, local eigenfrequencies and eigenfunctions are necessary for calculation of wave propagations (e.g. Dahlen & Tromp 1998;Okamoto & Tanimoto 2003). The method developed in this paper is effective in this situation since eigenfrequencies gradually vary with the local vertical structure. Local eigenfrequencies at a location are good guesses of trial values in another eigenvalue problem for 1-D structure at successive locations and fast convergence of the eigenfrequencies is generally expected.

C O N C L U S I O N
We developed a new numerical method of normal mode calculations of the Earth and planets. It is formulated based on a similar treatment of the global matrix which govern the oscillation system to that in the Henyey-type relaxation method used in solar seismology (e.g. Unno et al. 1989). Results calculated by this method are compared with those by DISPER80, and we found good agreement between them.
In general, the search for an eigenvalue in the complex plane becomes quite difficult in comparison with the search for an eigenvalue along the real axis. In the present method, we adopt Newton's method to find the next estimate of the eigenvalue. The rate of convergence depends on the linearity of the correction δω 2 around the target eigenfrequency. Although the linearity varies from mode to mode and with the starting point of calculation, we can reach the eigenvalue in a dozen of steps from a distant initial point. Even for models with an atmosphere, we can clearly obtain different eigenfrequencies of the fundamental spheroidal mode 0 S 29 and the fundamental acoustic mode 0 P 29 which are nearly degenerate as shown in Figs 2 and 5.
Since the stability of this method is based on an algebraic structure similar to the relaxation method, we can calculate modes whose energy is confined in the very shallow part of the Earth from the vicinity of the centre of the Earth to the top of the atmosphere. In the shooting method, we need the calculation of the mth minor propagators to keep numerical accuracy in computing the dispersion relation (e.g. Gilbert & Backus 1966;Takeuchi & Saito 1972). In our method, we can obtain shallower modes with large angular degree l without calculation of the minor propagators.
In conclusion, we propose the method developed in this article as a quite efficient and stable method for calculating normal modes of the Earth, planets and stars. Especially, when we include an atmosphere in the target model, the method is generally better and a more direct way than the previous methods. Since this method provides a direct and effective calculation for eigenmodes of an anelastic earth model, it is useful in calculation of synthetic seismograms using the normal mode summation technique for the anelastic case (e.g. Lognonné 1989;Lognonné 1991) in which anelastic normal modes are not normalized in the usual way.

A C K N O W L E D G M E N T S
I thank Prof Y. Fukao for reading this manuscript and giving me some helpful comments. I also thank Profs Y. Osaki, H. Shibahashi, T. Sekii and M. Takata for inviting me to their seminar to talk about this study and significant discussions on this subject. This work is supported by JSPS grant in aid (Kiban-Kenkyu (C) No. 15540406). I also thank Dr T. Okamoto for reading this manuscript carefully and polishing my english. We also thank Prof Lognonné and Prof. Tanimoto for critical comments on this manuscript. ω 2 > ω 2 a . In the same way as ω 2 a , ω 2 g can be regarded as the potential for gravity waves since they can propagate vertically only when ω 2 < ω 2 g . Therefore, we chose λ so that λ < 0 for ω 2 > ω 2 a (acoustic waves) and λ > 0 for ω 2 < ω 2 g (gravity waves) at the outer boundary. If one neglects the first term of the definition of D, ω 2 a is just c 2 k 2 and ω 2 g is just N 2 B as expected. For waves with long horizontal wavelengths, the first term is not neglected. Then the acoustic potential relates approximately to the ratio of sound velocity c and the scale height H, that is the critical frequency (e.g. Gough 1990). Note c 2 /H 2 ∝ 1/T . The simple interpretation of the critical frequency is that only waves with longer vertical wavelength than the minimum of the scale height can be trapped in the atmosphere. It can be seen from D(ω 2 ) = 0 for which we obtain if ω 2 >N 2 B and k is small. The elementary solutions given by eqs (A1) and (A2) are used in the calculation of the open boundary conditions in the same manner as the derivation of eq. (6).