Damien Fournier
email: fournier@mps.mpg.de

Chris S Hanson

Laurent Gizon

Hélène Barucq

Sensitivity kernels for time-distance helioseismology Efficient computation for spherically symmetric solar models

Keywords: helioseismology -Sun, oscillations -Sun, interior -methods, numerical

Context. The interpretation of helioseismic measurements, such as wave travel-time, is based on the computation of kernels that give the sensitivity of the measurements to localized changes in the solar interior. These kernels are computed using the ray or the Born approximation. The Born approximation is preferable as it takes finite-wavelength effects into account, although it can be computationally expensive. Aims. We propose a fast algorithm to compute travel-time sensitivity kernels under the assumption that the background solar medium is spherically symmetric. Methods. Kernels are typically expressed as products of Green's functions that depend upon depth, latitude, and longitude. Here, we compute the spherical harmonic decomposition of the kernels and show that the integrals in latitude and longitude can be performed analytically. In particular, the integrals of the product of three associated Legendre polynomials can be computed.

Results. The computations are fast and accurate and only require the knowledge of the Green's function where the source is at the pole. The computation time is reduced by two orders of magnitude compared to other recent computational frameworks. Conclusions. This new method allows flexible and computationally efficient calculations of a large number of kernels, required in addressing key helioseismic problems. For example, the computation of all the kernels required for meridional flow inversion takes less than two hours on 100 cores.

Introduction

The aim of time-distance helioseismology [START_REF] Duvall | [END_REF]) is to infer the subsurface structure of the Sun by measuring seismic wave travel times between any two points at the solar surface. The interpretation of these measurements requires understanding how waves propagate in the solar interior, i.e., solving the forward problem. Due to its simplicity, the ray approximation was initially used to invert flow velocities and sound-speed perturbations to a reference background model (Kosovichev 1996). It is still used today, for example to recover the meridional circulation (Rajaguru & Antia 2015). However, this approach is a high-frequency approximation that cannot be used to recover perturbations with sizes on the order of the local wavelength (Birch & Kosovichev 2000). Gizon & Birch (2002) derived a general framework for sensitivity kernels under the Born approximation and for random sources of excitation. Birch et al. (2004), Burston et al. (2015), and Böning et al. (2016) computed Born kernels using a normal-mode summation of the eigenfunctions in a solar-like stratified background. To treat axisymmetric background media (e.g., a background that includes large-scale differential rotation) and to include frequencies above the acoustic cutoff, Gizon et al. (2017) proposed solving the wave equation in frequency space by using a 2.5D finite-element solver. All these approaches are useful, but are computationally expensive. This limits their use in the interpretation of solar data as many kernels must be computed (and averaged). In some cases it is sufficient to consider perturbations to a steady spherically symmetric reference medium. The study of meridional circulation is one such application (see, e.g., Liang et al. 2017).

In this paper, we present a way to reduce the computational time of Born sensitivity kernels in a spherically symmetric background by treating the horizontal variables (the co-latitude θ and the longitude φ) analytically using the properties of the spherical harmonics. Here this approach is demonstrated using the scalar wave equation from Gizon et al. (2017), but it could be applied to the normal-mode summation method of Böning et al. (2016) or to solving the wave equation using a high-order finitedifference scheme (Mandal et al. 2017).

Born sensitivity kernels

Green's function in a spherically symmetric background

We follow the framework of Gizon et al. (2017), where the observable ψ(r, ω) at spatial location r = (r, θ, φ) and frequency ω is linked to the divergence of the displacement: ψ(r, ω) = c(r)∇ • ξ(r, ω). This scalar quantity solves

Lψ(r, ω) = s(r, ω), ( 1
)
where L is the spatial wave operator at frequency ω,

Lψ := -(ω 2 + 2iωγ)ψ -2iωu • ∇ψ -c∇ • 1 ρ ∇(ρcψ) , (2)
ρ and c are the solar density and sound speed from standard solar model S (Christensen- [START_REF] Christensen-Dalsgaard | NIST Digital Library of Mathematical Functions[END_REF], γ is the attenuation, u is a background flow, and s is a stochastic source term. We assume that the sources are spatially uncorrelated and depend only on depth and frequency such that the source covariance matrix is given by

M(r, r , ω) := E[s * (r, ω)s(r , ω)] = A(r, ω)δ(r -r ), (3)
where A(r, ω) is the radial profile of the source power. The wave field ψ can be obtained using

ψ(r, ω) = G(r, r , ω)s(r , ω)ρ(r )dr , ( 4
)
where G is the Green's function:

LG

(r, r , ω) = 1 ρ(r) δ(r -r ). ( 5
)
When the background is spherically symmetric (i.e., no flow and no heterogeneity), G can be written as

G(r, r , ω) = max =0 α G (r, r , ω) m=- Y m * (θ , φ )Y m (θ, φ), (6)
where r = (r, θ, φ), r = (r , θ , φ ), Y m are the normalized spherical harmonics, α = √ 4π/(2 + 1), and G is the Legendre component of the Green's function

G (r, r , ω) = 2π 0 π 0 G(r, r , ω)P (cos θ) sin θdθdφ. (7)
Equation ( 6) can be simplified by using the addition theorem (DLMF 2017, Eq. (14.18.1)) and introducing the great-circle angle between r and r . However, Eq. ( 6) is the form required in the following sections.

Cross-covariance in a spherically symmetric background

In a spherically symmetric background the expectation value of the cross-covariance between an observation point r 1 = (r 0 , θ 1 , φ 1 ) and a point r = (r, θ, φ) is

C(r 1 , r, ω) = E[ψ * (r 1 , ω)ψ(r, ω)] = α 2 m=- Y m * (θ 1 , φ 1 )Y m (θ, φ)C (r 0 , r, ω), (8)
where

C (r 0 , r, ω) = R 0 G (r , r 0 , ω) * G (r , r, ω)A(r , ω) ρ(r ) 2 r 2 dr , (9)
and r 0 = (r 0 , 0, 0) and r = (r, 0, 0) are on the polar axis. The radius r 0 is the observation radius, for example ∼150 km above the photosphere for SDO/HMI. In obtaining Eq. ( 8), we used the property that the cross-covariance depends only on the greatcircle distance between the two points. To simplify the computations we place one point on the polar axis so that the Green's function is axisymmetric and only the mode m = 0 needs to be computed.

Using the convenient source of excitation introduced in Gizon et al. (2017), C (r 0 , r, ω) is directly linked to the imaginary part of the Green's function G (r, r 0 , ω), but this assumption is not mandatory in this paper. The important assumption concerns the covariance of the sources of excitation that needs to be of the form given by Eq. (3), such that the cross-covariance depends only on depths and the great-circle distance between the two points r 1 and r. This assumption is common in helioseismology and is generally used in forward modeling (e,g., Kosovichev et al. 2000;Böning et al. 2016;Mandal et al. 2017).

Born flow kernels

Recovering flows in the solar interior is a major goal for local helioseismology, and so we focus here on flow kernels. The method presented here can be applied to all other types of perturbations with respect to a spherically symmetric background. The Born sensitivity kernel K = (K r , K θ , K φ ) connects the travel-time perturbation δτ to the vector flow u = (u r , u θ , u φ ), such that

δτ(r 1 , r 2 ) = K(r, r 1 , r 2 ) • u(r) dr. ( 10
)
According to Gizon et al. (2017) we have

K(r,r 1 , r 2 ) = 2iρ(r) ∞ -∞ dω ωW * (r 1 , r 2 , ω) × G(r 2 , r, ω)∇C(r 1 , r, ω) -G * (r 1 , r, ω)∇C * (r 2 , r, ω) , ( 11
)
where W is a weighting function that relates a change in the cross-covariance to a change in travel-time (Gizon & Birch 2002) and

∇ = (∂ r , 1/r ∂ θ , 1/(r sin θ) ∂ φ )
is the gradient operator with respect to the scattering location r. We note that in a spherically symmetric background, seismic reciprocity implies G(r, r , ω) = G(r , r, ω) for any r and r . The reference crosscovariance also satisfies C(r , r, ω) = C(r, r , ω).

The expression for the kernel may differ when a different observable is chosen; however, the above integral will always involve the product of a Green's function with the cross-covariance. One approach to obtaining the flow kernels (Böning et al. 2016;Mandal et al. 2017) is to compute the 3D Green's function and the cross-covariance using its spherical harmonic decomposition using Eq. ( 6). A reference kernel is usually obtained for a fixed pair of observation points and later rotated to obtain kernels for other pairs of points. However, a fine resolution in θ and φ is required in order to perform this rotation accurately, which makes the computation expensive in both time and memory.

Spherical harmonic decomposition of Born flow kernels

In order to circumvent the disadvantages mentioned above (e.g., rotation) and improve accuracy, we propose a new approach based on the spherical harmonic decomposition of the kernel,

K(r, r 1 , r 2 ) = m=- K m (r, r 1 , r 2 )Y m (θ, φ), (12)
where

K m (r, r 1 , r 2 ) = 2π 0 π 0 K(r, r 1 , r 2 )Y m * (θ, φ) sin θdθdφ. ( 13
)
Decomposing G(r 1 , r, ω) and C(r 2 , r, ω) into spherical harmonics, we can obtain the spherical harmonic coefficients of each kernel.

For the u r kernel, we have

K m r (r, r 1 , r 2 ) = , α α m=-m =- I r × f r (r)Y m * (θ 2 , φ 2 )Y m * (θ 1 , φ 1 ) + g r (r)Y m (θ 1 , φ 1 )Y m (θ 2 , φ 2 ) , (14)
where

f r (r) = 2iρ(r) ∞ -∞ ωW * (ω)G (r, r 0 , ω)∂ r C (r 0 , r, ω)dω, (15) g r (r) = -2iρ(r) ∞ -∞ ωW * (ω)G * (r, r 0 , ω)∂ r C * (r 0 , r, ω)dω, (16)
and

I r = 2π 0 π 0 Y m (θ, φ)Y m (θ, φ)Y m * (θ, φ) sin θdθdφ. ( 17
)
The integral of three spherical harmonics over the unit sphere can be done analytically using the Gaunt formula (see, e.g., Edmonds 1960, Eq. (4.6.3))

I r = 4π α α α (-1) m 0 0 0 m m -m , ( 18
)
where we have used the Wigner-3j symbols (see, e.g., Edmonds 1960, p. 45). The Wigner-3j symbol vanishes when m m + m , which enables us to remove the sum over m in the equation for K u r .

It can be shown that the expression for K r , K θ , and K φ can be recast in the form

K m j (r) = , α α L m=-L I j f j (r)Y m * (θ 2 , φ 2 )Y m-m * (θ 1 , φ 1 ) + I * j g j (r)Y m * (θ 1 , φ 1 )Y m-m * (θ 2 , φ 2 ) , (19)
where j ∈ {r, θ, φ} and L = min( , ).

Proceeding in a similar way, for K r the kernel K m θ depends on the functions f θ and g θ given by

f θ (r) = 2iρ(r) ∞ -∞ ωW * (ω)G (r, r 0 , ω)C (r, r 0 , ω)dω, (20)
g θ (r) = -2iρ(r) ∞ -∞ ωW * (ω)G * (r, r 0 , ω)C * (r, r 0 , ω)dω. (21)
The horizontal integral is

I θ = 1 r 2π 0 π 0 Y m (θ, φ)∂ θ Y m (θ, φ)Y m * (θ, φ) sin θdθdφ. ( 22
)
This integral I θ is much more difficult to evaluate than I r because of the θ derivative. In order to keep only associated Legendre polynomials in I θ , we use

dP m (cos θ) dθ = 1 2 ( + m)( -m + 1)P m-1 (cos θ) -( + m + 1)( -m)P m+1 (cos θ) , (23)
where the P m l are the normalized associated Legendre polynomials. We use the convention that P m±1 = 0 if |m ± 1| > , so that Eq. ( 23) remains valid for m = ± . Integrating Eq. ( 22) over φ and using Eq. ( 23), I θ becomes

I θ = 1 2 √ 2π r -( + m )( -m + 1)J m,m -1,m + ( + m + 1)( -m )J m,m +1,m , (24)
where m = mm and

J mm m = π 0 P m (cos θ)P m (cos θ)P m (cos θ) sin θdθ. (25)
Fortunately, this integral can also be evaluated analytically. It involves a sum of products of Wigner-3j symbols (see Appendix A and Dong & Lemus 2002).

The derivation of u φ is similar to u θ and requires the evaluation of the horizontal integral

I φ = im r 2π 0 π 0 1 sin θ Y m (θ, φ)Y m (θ, φ)Y m * (θ, φ) sin θdθdφ. (26)
Using

P m (cos θ) sin θ = - 1 m        (2 + 1)( + m -1)( + m) 2 -1 P m-1 -1 (cos θ) + (2 + 1)( -m -1)( -m) 2 -1 P m+1 -1 (cos θ)        , (27)
for m 0, we obtain

I φ = i 2 √ 2π r        (2 + 1)( + m -1)( + m ) 2 -1 J m,m-m-1,m , -1, + (2 + 1)( -m -1)( -m ) 2 -1 J m,m-m+1,m , -1,        . (28)
Now that we have the equations for the kernels, let us summarize the algorithm used for the resolution: 1. Computation and storage of the Green's function G l (r, r 0 , ω)

with the source on the polar axis, as a function of depth and harmonic degree for all frequencies; 2. For each great-circle distance between r 1 and r 2 :

computation of the cross-covariance using Eq. ( 8). If the convenient source of excitation of Gizon et al. (2017) is used, the cross-covariance is directly obtained from the imaginary part of the Green's function; computation of the weighting function W; computation of the functions f j and g j ; 3. Evaluation of the integrals I j and computation of the kernel using Eq. ( 19). A summary of the different terms required to compute the different components of the flow kernels using Eq. ( 19) is given in Table 1.

We note that the algorithm presented above could of course be used to compute sensitivity kernels for the cross-covariance amplitude using the linear definition of Nagashima et al. (2017) and the appropriate choice of W. It is also possible to obtain kernels for the cross-covariance function at a given frequency by removing the weighting function. In this case, the functions g j are just the complex conjugates of f j .

Table 1. Terms required to compute the spherical harmonic coefficients of the flow sensitivity kernels K m j using Eq. ( 19)

j f j (r) I j r 2iρ(r) ∞ -∞ ωW * (ω)G (r, r 0 , ω)∂ r C (r 0 , r, ω)dω 4π α α α 0 0 0 m m -m -m θ 2iρ(r) ∞ -∞ ωW * (ω)G (r, r 0 , ω)C (r 0 , r, ω)dω 1 2 √ 2πr -( + m -m)( -m + m + 1) J m,m-m-1,m + ( + m -m + 1)( -m + m) J m,m-m+1,m φ f φ (r) = f θ (r) i 2 √ 2πr (2 +1)( -m+m-1)( -m+m) 2 -1 J m,m-m+1,m
, -1,

+ (2 +1)( +m-m-1)( +m-m) 2 -1 J m,m-m-1,m , -1,
Notes. The integrals J mm m depend only on Wigner-3j symbols and can be computed using Eq. (A.2) corresponding to the algorithm of [START_REF] Christensen-Dalsgaard | NIST Digital Library of Mathematical Functions[END_REF].

Numerical validation

To evaluate the flow kernels in this framework, the only ingredient to prescribe is the Green's function as a function of the spherical harmonic degree and depth for a source located at the pole. We compute it using 1D finite elements with the solver Montjoie [START_REF] Chabassier | High Order Finite Element Method for Solving Convected Helmholtz Equation in Radial and Axisymmetric Domains[END_REF][START_REF] Fournier | [END_REF]. The Green's functions are computed with a high enough frequency resolution to resolve the modes using the mode linewidths (5671 frequencies corresponding to 4 days of observations at 60 s cadence) and max = 400. Representations of the different components of the flow kernels between two points r 1 and r 2 centered at 40 • and separated by 42 • is shown in Figs. 1 and2. They exhibit the classical banana-doughnut shape with zero sensitivity along the ray path. Small-scale structures are visible close to the surface as we kept values of up to 400. Visually, there is no difference between the kernels computed with this new approach, the ones from Gizon et al. (2017), or the ones obtained by rotation so only one is shown here.

To allow a more quantitative approach, Fig. 3 compares kernels computed using our new approach to the approach presented by Gizon et al. (2017) where the background is axisymmetric and the Green's function is computed for each azimuthal degree m on a 2D grid. These kernels K r and K θ are averaged over longitudes ( m = 0) where r 1 is located at the pole and r 2 is at a co-latitude of 42 • (similar to Fig. (17) in Gizon et al. 2017). The results show good agreement, validating the method presented here. We note slight differences in the structure of K θ at a depth of 500 km, and attribute this to the numerics of the 2D FEM solver differing. Specifically, the 2D FEM has an inherent difficulty in computing the real part of the Green's function close to the Dirac source location (see for details [START_REF] Chabassier | High Order Finite Element Method for Solving Convected Helmholtz Equation in Radial and Axisymmetric Domains[END_REF]. In order to ensure that these small differences do not affect the interpretation of the data, we compute the travel times induced by the meridional flow model from Gizon et al. (2017). We decompose the flow in Legendre polynomial u akin to Eq. ( 7) and compute a travel time for each according to

δτ = R 0 K ,m=0 (r) • u (r) r 2 dr. ( 29
)
The bottom panel of Fig. 3 shows that this travel time as a function of is nearly indistinguishable from the travel times in Gizon et al. (2017); the differences are less than 0.5 ms.

As the computational burden of this new method depends on the maximum harmonic degree of the Green's function max and the number of required for the kernel, we illustrate the efficiency of the method on two problems of interest: kernels for meridional flow inversions and kernels for supergranulation inversions.

Fig. 3. Top and middle panels: comparison of the kernels for ¯ = 5 and 10 (blue and black, respectively) computed using the method here (solid lines) and the method of Gizon et al. (2017;dots). Bottom panel: travel times δτ due to the radial (red) and the latitudinal (green) components of the flow for each ¯ of the kernels presented here (solid line) and those of Gizon et al. (2017;dots).

Computation of kernels for meridional flow

As a first test, we compute the kernels that are required to interpret meridional flow measurements. As the flow varies slowly with latitude, we can limit the number of spherical harmonic coefficients of the kernels to ≤ 10 (see Fig. 3). As for the observations, a low-pass filter with max = 300 is applied to the Green's function. The method can be parallelized in so we use 11 cores to compute kernels up to = 10.

For meridional flow measurements, the separation distance between the source and the receiver is generally prescribed for different values of the mean latitude (see, e.g., Liang et al. 2017). For a given separation distance, we compute 15 kernels corresponding to 15 different latitudes. We then vary the separation distance in order to probe different depths.

Table 2 shows the computational times and memory requirements of the different steps of the algorithm. The computation of the Green's function for a source located at the pole can be done once and for all and stored as it is necessary to do it for every kernel. Therefore, the computation is very fast (3 seconds per frequency for max = 300) and embarrassingly parallel in frequency so it could also be recomputed every time. The computation of the frequency integrals ( f j and g j ) consists in loading the Green's function and the weighting function W and summing over frequencies. The reading of the Green's function files for all frequencies takes most of the time. The spatial integrals I r , I θ , and I φ can be computed once and stored for future use. However, the computational time is small compared to the full computation of the kernel, so we decide to recompute I j every time as the reading time can depend upon file system I/O load. The computations are parallelized in and hence need 11 cores for each step since 0 ≤ ≤ 10. The computation of the kernel for u r is faster since the computation of I r requires the evaluation of only two Wigner-3j symbols, unlike I θ . However, the major difference in computational time between K r and K θ comes from the sum in in Eq. ( 19). For K r , the sum covers the range from to + , since I r = 0 for other values due to the properties of the Wigner-3j. On the contrary, the sum in for the computation of K θ must be computed for the full range from 0 to .

Even though the computational burden of K θ is greater than K r , the total burden remains significantly lower than for other methods (see Table 3). All the kernels required to perform a meridional flow inversion can be computed within 2 h with 100 cores, and the memory requirements do not exceed 1 GB. The approach mentioned in Sect. 2.3, where the full 3D kernel is computed and rotated to obtain different latitudes, would take 11 days using 1000 cores with very significant memory requirements. In the axisymmetric approach of Gizon et al. (2017) the computational time would be about 40 days on 1000 cores for all the same set of kernels.

The computational times presented here are for point-topoint measurements; however, this framework can easily be extended to geometric averaging such as arc-to-arc measurements often performed for meridional flow measurements (e.g., Liang et al. 2017). It is only necessary to replace the product of the two spherical harmonics in Eq. ( 19) by a sum over all the points of the arc.

Computation of kernels for supergranulation

Resolving smaller scale flows such as supergranules requires a high spatial resolution, and thus the Green's function needs to include much higher harmonic degrees than for meridional flow Green's function. For example, Duvall & Hanasoge (2013) considered measurements up to max = 700, but even higher max values may be required. Furthermore, supergranulation flows have maximum power around = 120, thus the kernels should at least be computed up to = 300, or higher depending on the power distribution of the flow at large . The computational burden for these kernels is summarized in Table 4. The computation A156, page 5 of 6 of the m = 0 component of the Green's function now takes about 23 second per frequency, and the loading of the files to compute f j around 6 min. The computation of the Wigner symbols is computationally more challenging as increases, since the number of loops scales as 3 max due to loops in , , and m. While the computation of I r is still fast, the evaluation of I θ now takes 270 min on 100 cores. Computing a set of 200 kernels would take 3 days with 1000 cores, which is significantly longer than for the meridional flow kernels, but still one order of magnitude faster than the approach of Gizon et al. (2017) and with a smaller memory requirement.

Conclusions

We presented a technique that is faster than previous approaches used to compute travel-time kernels under the assumption that the background medium is spherically symmetric. This technique does not rely on the numerical computation of kernel rotations and thus the memory requirements are small. Instead, the spatial integrals are performed analytically, which also leads to higher accuracy. For example, for meridional circulation applications, the kernels can be computed one thousand times faster than with previous methods, using a tenth of the memory requirement.

Fig. 1 .

 1 Fig. 1. Slices along a constant meridian of the point-to-point 3D traveltime difference kernel for u r (left) and u θ (right). The 3D kernel for u φ is zero along this slice. The kernel is computed with r 1 and r 2 separated by 42 • , with mean latitude 40 • . The green line is the ray path between the two points and the dashed black line shows the image plane of Fig. 2.

Fig. 2 .

 2 Fig. 2. Slices of the point-to-point 3D travel-time difference kernel for u θ (left) and u φ (right) along the plane indicated in Fig. 1. The 3D kernel for u r is mostly zero within this plane. The kernel is computed with r 1 and r 2 separated by 42 • , with mean latitude 40 • . The green cross indicates the intersection of the ray path and the image plane. The dashed black line shows the image plane of Fig. 1.

Table 2 .

 2 Computational time of the different steps to obtain 15 flow kernels (15 latitudes) for a given separation distance with ≤ 10 and max = 300 using 11 cores.

	Computation steps	Time [11 cores] Memory
	Green's function	13 min	100 MB
	f r (or f θ )	9 min	1 GB
	I r	0.5 s	1 GB
	I θ	6 min	1 GB
	Sum of terms in K m r	3 min	1 GB
	Sum of terms in K m θ	32 min	1 GB

Table 3 .

 3 Comparison of the computational time and memory requirements to evaluate 225 kernels for the meridional flow.Notes. For the rotation of the 3D kernels as inBöning et al. (2016), andMandal et al. (2017), we used a horizontal grid sampled with N θ = 1001 and N φ = 2001 points.

	Method	Time [cpu hours] Memory [GB]
	This paper	170	1
	Rotation of 3D kernels	2.7 × 10 5	40
	Gizon et al. (2017)	10 6	8

Table 4 .

 4 Computational time of the different steps to obtain 15 flow kernels for a given separation distance with ≤ 300 and max = 700 using 100 cores.

	Computation steps	Time [100 cores] Memory
	Green's function	11 min	100 MB
	f r (or f θ )	6 min	1 GB
	I r	2 min	1 GB
	I θ	270 min	1 GB
	Sum of terms in K ,m r	15 h	1 GB
	Sum of terms in K ,m θ	39 h	1 GB

Acknowledgements. The computer infrastructure was provided by the German Data Center for SDO funded by the German Aerospace Center (DLR) and by the Ministry of Science of the State of Lower Saxony, Germany.

Appendix A: Algorithm for the integral of three associated Legendre polynomials

For the sake of completeness, we summarize the algorithm of [START_REF] Christensen-Dalsgaard | NIST Digital Library of Mathematical Functions[END_REF] adapted to this study. The integral of three associated Legendre polynomials,

can be computed analytically in terms of sums of products of Wigner-3j symbols

which can be evaluated analytically. In this paper, we only need this value for m 123 = ±1. As this integral is zero for odd values of 123 , due to the parity of the associated Legendre polynomials, we set 123 = 2p + 1. Then, for a given m 123 , the value of J( 123 , m 123 ) can be evaluated recursively using .4) and J(2n + 1, -1) = (2n -1) 2 4(n + 1) 2 J(2n -1, -1), (A.5)

where n = 1, 2, . . . p, together with the initial conditions J(1, 1) = -π 2 and J(1, -1) = π 4 . (A.6)