Computational helioseismology in the frequency domain: acoustic waves in axisymmetric solar models with flows
^{1} MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: gizon@mps.mpg.de
^{2} Institut für Astrophysik, GeorgAugustUniversität Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
^{3} Magique3D, Inria Bordeaux SudOuest, Université de Pau et des Pays de l’Adour, 64013 Pau, France
email: helene.barucq@inria.fr
^{4} Magique3D, Inria Bordeaux SudOuest, Université de Bordeaux, 33400 Talence, France
^{5} Institut für Numerische und Angewandte Mathematik, GeorgAugustUniversität Göttingen, Lotzestraße 18, 37083 Göttingen, Germany
Received: 4 August 2016
Accepted: 2 November 2016
Context. Local helioseismology has so far relied on semianalytical methods to compute the spatial sensitivity of wave travel times to perturbations in the solar interior. These methods are cumbersome and lack flexibility.
Aims. Here we propose a convenient framework for numerically solving the forward problem of timedistance helioseismology in the frequency domain. The fundamental quantity to be computed is the crosscovariance of the seismic wavefield.
Methods. We choose sources of wave excitation that enable us to relate the crosscovariance of the oscillations to the Green’s function in a straightforward manner. We illustrate the method by considering the 3D acoustic wave equation in an axisymmetric reference solar model, ignoring the effects of gravity on the waves. The symmetry of the background model around the rotation axis implies that the Green’s function can be written as a sum of longitudinal Fourier modes, leading to a set of independent 2D problems. We use a highorder finiteelement method to solve the 2D wave equation in frequency space. The computation is embarrassingly parallel, with each frequency and each azimuthal order solved independently on a computer cluster.
Results. We compute traveltime sensitivity kernels in spherical geometry for flows, sound speed, and density perturbations under the first Born approximation. Convergence tests show that travel times can be computed with a numerical precision better than one millisecond, as required by the most precise traveltime measurements.
Conclusions. The method presented here is computationally efficient and will be used to interpret traveltime measurements in order to infer, e.g., the largescale meridional flow in the solar convection zone. It allows the implementation of (fullwaveform) iterative inversions, whereby the axisymmetric background model is updated at each iteration.
Key words: Sun: helioseismology / Sun: oscillations / Sun: interior
© ESO, 2017
1. Introduction
Timedistance helioseismology and related techniques are methods for probing the complex dynamics of the solar convection zone (Duvall et al. 1993; Gizon & Birch 2005; Gizon et al. 2010). Information is encoded at the solar surface in the twopoint crosscovariance function of the random solar oscillations. The crosscovariance function tells us about the travel time of wave packets between any two locations, in either direction. Flows break the timesymmetry of the crosscovariance function and thus leave a signature in the observations.
Two topics of current interest include the study of meridional circulation (e.g., Zhao et al. 2013; Liang & Chou 2015) and convective flows (e.g., Hanasoge et al. 2012; Langfellner et al. 2015) using, in particular, the SDO/HMI space observations. In both cases, low flow velocities are involved, resulting in traveltime perturbations that are less than a second. Hence it is not surprising that answers vary among investigators as they choose different interpretations and modeling procedures (for a review see Hanasoge et al. 2016).
Helioseismic studies consist of several steps: measuring solar oscillations, processing and averaging the observations to extract the seismic data (e.g., wave travel times), and interpreting the seismic data using forward and inverse methods to estimate solar internal properties. In this paper we mostly consider the forward problem, i.e. the computation of synthetic seismic data for a given solar model, which is a necessary step for reliable interpretations. A short discussion of the iterative inverse problem is given at the end of the paper.
A framework was proposed by Gizon & Birch (2002) to compute perturbations to the crosscovariance function caused by weak heterogeneities. This framework has proven to be useful for local helioseismology (Birch et al. 2004, 2007; Birch & Gizon 2007; Jackiewicz et al. 2008; Švanda et al. 2011; Hanasoge et al. 2011; Burston et al. 2015; Böning et al. 2016). However, the computational expense has been a limiting factor in applications. Advances in the fields of Earth seismology, exploration geophysics and ocean acoustics provide some guidance for improvements of the computational methods and the PDEconstrained formulations of the forward and inverse methods. Hanasoge et al. (2011) took a step in the direction of incorporating some of these ideas in the time domain within the framework of Gizon & Birch (2002). Of particular relevance to the present work is the proposal by NissenMeyer et al. (2014) to consider 3D wave propagation in relevant axisymmetric background media, in order to reduce the computational cost. Using a spectral decomposition in the azimuthal direction, the forward problem separates into many independent 2D problems, which can be solved to study the spatial sensitivity of seismic travel times to 3D heterogeneities (van Driel & NissenMeyer 2014; van Driel et al. 2015; Bottero et al. 2016). In addition, theoretical studies have led to a better understanding of the connection between the crosscovariance function and the Green’s function in media permeated by random sources of excitation (Snieder & Larose 2013, and references therein).
In this paper we consider a number of current challenges for computational local helioseismology and propose some solutions. Some of these challenges are listed here:

Oscillation power spectrum. The purpose of helioseismology isto infer a model for the solar interior (density, sound speed, flows,etc., as functions of position) that is consistent with the seismicdata. In global mode helioseismology, the emphasis is on thecomparison between the model and observed mode frequencies.In contrast, local helioseismology requires models for the modeamplitudes and line profiles in addition to the mode frequencies. Thus a solar model must also describe wave excitation anddamping. Most reference solar models are 1D and their qualitycan be assessed by comparing the model and observed oscillationpower spectra. In linear inversions for local helioseismology, thereference oscillation power spectrum has to be very solarlike.

Flexibility. The semianalytical methods used so far to solve the forward problem (Gizon & Birch 2002) have proved useful (e.g., Burston et al. 2015) but lack flexibility. To first order (first Born approximation), the perturbation to the crosscovariance function is given by interaction integrals over products of Green’s functions. Various approximations are made to speed up the computation of these integrals, for example using a local Cartesian geometry and neglecting some anisotropic effects. A generalization to more realistic setups, including the treatment of various instrumental (e.g., varying pointspread function) and geometrical effects (e.g., foreshortening, lineofsight projection) is very cumbersome (e.g., Jackiewicz et al. 2007). How to include these complex effects in forward models in a reliable and efficient way is an open question.

Representation of the Green’s function. Often Green’s functions are computed using a normalmode expansion (Birch et al. 2004). This approach is problematic for a heterogeneous 3D reference model, for which eigenmodes are not readily available. Even if they were known, the summation formula would be computationally expensive. An additional complication is the treatment of the continuous spectrum above the acoustic cutoff frequency (5.3 mHz).

Timedomain simulations. Alternatively, the Green’s function may be computed numerically by solving the equations of motion in the time domain (e.g., Hanasoge & Duvall 2007; Cameron et al. 2008). However, this requires a stabilization of the background model by changing the buoyancy frequency (e.g., Schunker et al. 2011; Papini et al. 2014). Unless this operation is performed first, the linearized equations allow convective modes that grow exponentially. Unfortunately the mode frequencies are seriously affected by the stabilization of the solar model and move too far from the solar observations (Papini et al. 2014).

Computational challenge and inverse problem. Modern helioseismic observations consist of large data sets (long time series of 16 megapixel images). Crosscorrelations span a huge 5D space. Extracting the relevant information from this data set requires a good strategy and an efficient forward solver. This is especially true for iterative inversions, in which the forward solver is run for each update of the model of the solar interior.
We propose to address the above issues by carrying out the following steps:

Rewriting the perturbation to the crosscovariance function interms of the Green’s function, G, and the expectation value of thecrosscovariance function, , in the reference model. In thisformalism G and are the fundamental quantities that enter anyforward calculation in local helioseismology. The problem thenbecomes deterministic. Additionally, many systematic effectscan be accounted for by treating them as numerical operations on Gand .

Computing G and in the frequencydomain. The advantage of this approach is that frequencies are independent of each other for wave propagation in a steady background (many important problems in helioseismology fall in this category). This allows trivial parallelization of the computation in a multicore environment. Furthermore, the computation can easily be restricted to the frequency range of interest, at the necessary frequency resolution.

Adopting a flexible geometrical setup using a finiteelement discretization of the problem. This will enable us to treat problems in which the required spatial resolution depends on both the vertical and horizontal coordinates. Spherical geometry is naturally possible in this setup.

Considering a simplified problem of scalar acoustics in a stratified axially symmetric reference model for the sake of simplicity. Because of the axial symmetry the problem separates into many independent 2D problems, one for each azimuthal order. The parallelization in both frequency and azimuthal order allows very efficient computation.

Demonstrating that there is a choice of source covariance such that the expectation value of the crosscovariance is directly related to the Green’s function, even in the presence of a background flow. In this setting, the Green’s function is the only fundamental quantity that needs to be computed.
We use the above approach to compute model power spectra, timedistance diagrams, and traveltime sensitivity kernels. An application including a largescale meridional flow will be presented. Even though we consider a scalar observable in this paper, the proposed approach is intended to be generalized to more realistic observables in the future.
2. Pure acoustics in the frequency domain
2.1. Solar oscillations
In an inertial frame, the displacement ξ(r,t) of small amplitude waves obeys an equation of the form (1)where ℋ is a spatial operator and u(r) is a background flow. Wave attenuation is accounted for by temporal convolution (⊗) with the function γ(r,t). The function f(r,t) represents forcing by turbulent convection and is thus a realization drawn from a random process. The spatial differential operator ℋ is Hermitian for appropriate choices of the boundary condition and inner product (see LyndenBell & Ostriker 1967). An approximation of ℋ that captures the physics of acoustic oscillations is: (2)where ρ and c are the density and sound speed in the reference model. Assuming that the medium is steady, the problem can be written in frequency space as (3)using the Fourier convention (4)If, in addition, the source of excitation is statistically stationary, then the frequency components are independent random variables and the frequency components of the wave displacement, , are also statistically independent of each other. Thus, in frequency space, the problem separates into a set of independent boundaryvalue problems, one for each frequency.
In the following sections we drop the hat on top of Fourier transforms to simplify the notation.
2.2. Scalar wave equation
We further specify the wave equation to be solved. There are two possible choices for reducing the equation for the wave displacement to a scalar equation. One possibility is to rewrite Eq. (3)in terms of the quantity ψ = c∇·ξ. Neglecting the gravity terms, the secondorder terms in γ and u, and the cross term in γu in Eq. (3), and then taking the divergence, we obtain (5)For slow variations of u, c and γ compared to the wavelength, the wave equation for ψ simplifies to (6)with (7)The random source of excitation, s(r,ω) = c(r)∇·f(r,ω), and the attenuation, γ(r,ω) > 0, depend on frequency.
Another possible choice in order to obtain a scalar equation is to use the ansatz ξ = ρ^{1}∇(ρcψ). Under the assumption that the forcing term ρf is curl free, this leads to the same scalar equation as above.
When ρ and c are solarlike, the scalar Eq. (6)captures most of the interesting physics of p modes. For example, a solarlike density gives the correct acoustic cutoff frequency of 5.3 mHz. However, this equation is not relevant for modeling f and g modes.
Here, we choose to specify the source of excitation through the covariance function (8)where ∗ takes the complex conjugate and the functions ϵ(r) and Π(ω) control the spatial and frequency dependencies of the source covariance respectively. In the above expression we assumed spatially uncorrelated sources. This formulation is standard in local helioseismology (e.g., Birch et al. 2004).
2.3. Boundary condition
The geometrical setup is shown in Fig. 1. We work in a sphericalpolar coordinate system (r,θ,φ). The solar (or stellar) model is cylindrically symmetric around the ẑaxis. The computational boundary S encloses the solar photosphere. The wave equation is supplemented by a boundary condition at S on ψ. We denote by the unit normal vector to S and by ∂_{n}ψ the normal derivative of ψ.
Fig. 1 Geometrical setup for a background medium symmetric about an axis ẑ. The thick contour delineates the boundary, S, of the computational domain V. A position vector r in V is specified by sphericalpolar coordinates (r, θ, φ). The seismic observable ψ(r) is measured near the photosphere. 

Open with DEXTER 
A freesurface boundary condition (ψ ∝ ∇·ξ = 0, Dirichlet) is often used in helioseismology (Bogdan et al. 1996; ChristensenDalsgaard 2003, etc.). However, it is more appropriate to choose a transparent boundary condition to model the continuous spectrum above 5.3 mHz and to allow highfrequency waves to escape (e.g., Kumar et al. 1990). The use of perfectly matched layers is one way to approximate a transparent boundary (e.g., Hanasoge et al. 2011). Here, we choose the simplest approximation to the Sommerfeld outgoing radiation condition, (9)where k_{n} is the local wavenumber normal to the boundary. In order for this radiation boundary condition to hold true, the sound speed and density must be constant near the computational boundary, which is not the case in the Sun. In practice, we extend the solar reference model with a layer that smoothly transitions to constant density and sound speed (see Sect. 7.1 for details).
2.4. Hermiticity
In this section and the following, we demonstrate some important properties of the scalar wave Eq. (6), which will be used in the later sections. These basic properties would remain the same in the full vector equations.
One important property of the equations that describe nonattenuating waves in the Sun is the Hermiticity of the wave operator, which implies that the eigenvalues are real. LyndenBell & Ostriker (1967) proved Hermiticity of the wave operator for a selfgravitating fluid bounded by a free surface.
Here we ask if this condition is satisfied for the timeharmonic scalar wave equation defined by Eq. (6). The scalar wave field belongs to a complex Hilbert space. For any two functions ψ and ϕ in this space we define the inner product (10)where V is the computational domain. Performing integration by parts twice implies (11)We see that the spatial operator H is Hermitian if ψ and ϕ vanish on the boundary. However, for a radiation boundary condition, the surface integral remains and H is not Hermitian. The eigenvalues of H are no longer real.
The other terms in Eq. (6)are the damping and advection terms. Using integration by parts once, we ascertain that the damping operator is antiHermitian and the advection operator is Hermitian: The last result is a consequence of mass conservation, ∇·(ρu) = 0, for a steady flow that does not cross the computational boundary ( on S).
Combining Eqs. (6), and (11)–(13)yields (14)where (15)is the operator obtained by switching the sign of γ in L or, equivalently, by switching the sign of u and taking the complex conjugate. We note that for a radiation boundary condition, L^{†} is not the adjoint of L because the surface integral remains.
2.5. Green’s function and generalized seismic reciprocity
In geophysics, the principle of seismic reciprocity states that the same signal should be recorded if the locations of a source and a receiver are exchanged. Although seismic reciprocity is not preserved in the presence of a flow, a generalization of reciprocity can be obtained.
Let us introduce the Green’s function as the solution to (16)where δ is the Dirac delta function and G satisfies the same boundary condition as the wave field ψ. Consider a source at r_{1} in the physical domain and a source at r_{2} in the domain in which the flow has opposite sign: (17)where G(·;−u) is the Green’s function associated with the operator L_{(− u)}, with the same boundary condition as G. Multiplying the first equation by G(r,r_{2},ω;−u) and the second by G(r,r_{1},ω), integrating each over position r, and subtracting the two equations, we obtain (18)Following the same logic as earlier (integration by parts), (19)for any two functions ψ and ϕ in the space of solutions. This time the surface integral vanishes for either type of boundary condition (Dirichlet or radiation). In particular, (20)This relationship shows that the righthand side of Eq. (18)vanishes and so we find (21)This is a generalization of seismic reciprocity: the Green’s function is unchanged upon exchanging source and receiver and changing the sign of the flow. The conclusion holds even though the damping operator is not Hermitian.
3. Crosscovariance function and travel times
For the sake of simplicity, let us suppose that ψ can be directly measured on the solar surface. This choice is not unreasonable as ψ is proportional to the pressure fluctuations and is thus related to intensity fluctuations, an observable quantity. This choice is not at all a limitation of the method: other observables can be derived from ψ, for example a proxy for the lineofsight velocity.
At frequency ω, consider the crosscovariance in Fourier space as the product of the wave field at two locations of measurement, (22)In terms of the Green’s function, we have (23)where the source s(r,ω) is a realization of a random process. Under the assumption of spatially uncorrelated sources (see Eq. (8)), the crosscovariance can be written as one volume integral: (24)Following Gizon & Birch (2002), we define the perturbation to the travel time δτ between points r_{1} and r_{2} as (25)where (26)C_{ref}is a reference crosscovariance function, and w(t) is a temporal window function that selects a particular section of the data. For example, we can choose w(t) = exp [−(t−t_{g})^{2}/ 2σ^{2}] where t_{g}> 0 is a group time and σ sets the width of the window.
4. New formulation of travel time sensitivity kernels
In the presence of subsurface perturbations to the reference model the travel times computed from Eq. (25)will be affected. An understanding of the spatial sensitivity of the measurements to structural perturbations and flows requires the development of travel time sensitivity kernels. In this section we develop a new formulation of the traveltime sensitivity kernels in terms of the Green’s function and the crosscovariance computed in the reference model.
4.1. Perturbations to the medium
In addressing how travel times are sensitive to changes in the solar reference model, we consider changes in c(r), ρ(r), and u(r), as well as spatial changes to the attenuation γ(r) and the source amplitude ϵ(r). The vector flow (27)is specified by the three components u_{k} on the basis { ê_{k} } of unit vectors , and for the sphericalpolar coordinate system. In short, the physical variables of interest can be combined into a set (28)We are looking for traveltime sensitivity kernels K_{α} such that the traveltime perturbation δτ can be written as (29)for infinitesimal perturbations δq_{α}. For some variables there is a surface integral on the computational boundary, because integration by parts is needed to write traveltime perturbations in the form of Eq. (29). For these variables, we assume that the perturbations to the model vanish on the computational boundary (high in the atmosphere) and we drop the surface terms.
In order to find the kernels, the traveltime perturbation (30)is written in terms of the firstorder perturbation to the crosscovariance: (31)The next step is then to express δψ in terms of the δq_{α}.
4.2. Perturbation to the wave field
To first order (first Born approximation) the operator defined in Eq. (6)acts on the perturbed wave field through, (32)where δL is the perturbation to the wave operator caused by perturbations to the medium and (33)A formal solution to Eq. (32)is obtained in terms of the Green’s function: (34)
4.3. Perturbation to the crosscovariance
With Eq. (34)in hand, the expectation of the perturbation to the crosscovariance is then determined from Eq. (31): (35)where we write the perturbation to the source covariance as a spatially varying change in amplitude, (36)In order to obtain kernels for individual perturbations, we introduce the bilinear operators ℒ_{α} such that for any functions g and h we have (37)These bilinear operators are obtained by integration by parts and are explicitly given by In Eq. (40), denotes the component of the spatial gradient in the k direction. Combining Eqs. (37) and (35), we obtain an explicit linear relationship between the δq_{α} and : (42)Using this expression together with Eqs. (29)and (30), the traveltime sensitivity kernels are (43)for the scatterers δq_{α} ∈ { δc,δρ,δu_{k},δγ } and (44)for the source amplitude δϵ.
5. Convenient source of excitation
5.1. How can we simplify the computation of ?
The kernels from the previous section are well defined once the background medium, attenuation and source functions have been specified. At each ω, the attenuation function γ(r,ω) can be tuned to yield the observed line widths of solar oscillations in the power spectrum. Similarly, at each ω, the source function ϵ(r)Π(ω) can be tuned to give the observed mode amplitudes (Birch et al. 2004; Birch & Gizon 2007). However, the integral relationship between a general ϵ(r) and is far from trivial to compute.
Here we adopt a different strategy. We ask whether a convenient source covariance exists such that the expectation value of the crosscovariance can be written in terms of the Green’s function. As a second step we check whether the resulting power spectrum is solarlike.
It is known in geophysics and acoustics that, under appropriate conditions on the source covariance, the expectation value of the crosscovariance in the frequency domain is related to the imaginary part of the Green’s function (see Snieder et al. 2009, and references therein). If we could write such a simple relationship, our problem would be considerably simpler. For all practical purposes, the Green’s function would be the only truly important quantity.
Let us start by rewriting the equation for the Green’s function (Eq. (16)) for sources at r_{1} and r_{2}, taking the complex conjugate of the first: (45)Multiply the first equation by G(r,r_{2},ω) and the second equation by G^{∗}(r,r_{1},ω), integrate each over r, and then subtract to find (46)The surface term does not vanish unless specific boundary conditions are used, as noted earlier by Snieder (2007). In the solar case (transparent boundary condition), this surface term does not vanish. Reversing the flow u in the above equation, we obtain (47)where the surface integral is the explicit expression for the surface term. By inspection of Eq. (24), we see that the choice of source covariance amplitude (48)implies (49)Thus the expectation value of the crossvariance can be written as a sum of causal and anticausal Green’s functions, when waves are appropriately excited throughout the volume and on the computational boundary S at radius r = R(θ). The volume sources must be proportional to the local attenuation to enforce energy equipartition between the modes (see, e.g., Snieder 2007).
To simplify the problem it would be very nice to assume that equality (49)holds. Does it give a reasonable oscillation power spectrum? The power spectrum can be tuned in frequency space by choosing Π(ω). However, the source of Eq. (48)leaves no freedom for the distribution of power versus radial order at fixed frequency.
A drawback of the choice of source (Eq. (48)) is that the source covariance cannot be updated independently from attenuation and density. Thus, perturbations to the crosscovariance are fully specified by perturbations to c, ρ, γ, and u.
5.2. Sensitivity kernels in terms of G only
As the crosscovariance depends only on the Green’s function via Eq. (49), it is also the case for the kernels. Denoting for simplicity (50)the kernels defined by Eq. (43)can be computed from four Green’s functions (51)For example, the traveltime perturbations induced by a background flow perturbation, (52)can be computed through the kernels (53)
5.3. Special case of no background flow
Following from the previous section and taking the further assumption of no background flow (u = 0), seismic reciprocity implies that G^{†} = G^{∗} and the crosscovariance takes the simple form, (54)This simple form for the crosscovariance simplifies the scattering bilinear operators: Thus, to compute a kernel K_{α} when there is no background flow, all we need is two Green’s function computations, one with a source at r_{1} and the other with a source at r_{2}.
6. Forward solver in the frequency domain
With the theory of the previous sections in hand, we wish to compute the Green’s function at fixed ω, defined by Eq. (16). In order to achieve this we utilize a standard Finite Element Method (FEM), which is further explained in Sect. 6.5.
Ideally we would solve for the Green’s function in a general 3D computational domain. However, as discussed later in Sect. 8.3, 3D computations are very demanding. In the following subsections we consider a reference model that is symmetric about an axis. This setup is appropriate for studying the effects of largescale flows (differential rotation and meridional flow) on solar oscillations. A substantial benefit of this choice is that the solution to the 3D problem can be obtained by decomposition into a set of independent 2D problems, one for each Fourier component in longitude. We refer to problems where the background medium is axisymmetric as 2.5D problems. Below we describe the computational setup and give a weak formulation of the solved equations.
6.1. Montjoie finite element code
We utilize the FEM code Montjoie^{1}, developed at Inria for various wave propagation problems. Montjoie is versatile, welltested, and robust (see Duruflé 2006; Bergot et al. 2010). Montjoie allows us to consider 1.5D, 2.5D, and 3D geometries. The 1.5D problem assumes a spherically symmetric background model and the solution is decomposed into spherical harmonics and 1D finite radial elements. The 2.5D method assumes an axisymmetric background model and the solution is decomposed into azimuthal Fourier modes and 2D finite elements; it is fast and useful for a number of applications. The full 3D method is slow, although it would work for a small fraction of the solar volume. The computational times for the respective geometries are compared in Sect. 8.3.
6.2. Geometrical setup for the 2.5D problem
For an axisymmetric background model, the computational domain Σ is a meridional generating section of the geometry (Fig. 2), i.e. the halfdisk of radius R in the case of the sphere.
Fig. 2 Twodimensional generating section Σ of the volume V, delineated by ∂Σ (thick contour). A point in Σ is specified by coordinates (r,θ), where r is the radius and θ is the colatitude, or, equivalently, by coordinates (ϖ, z), where ϖ = rsinθ is the distance to the axis. 

Open with DEXTER 
The background model is symmetric about the ẑaxis. In 3D space a spatial point is associated with the position vector (59)where belongs to the 2D section Σ and φ is longitude about the axis of symmetry (Fig. 2). The background density, sound speed, and flow are specified at positions in Σ: where is the meridional component of the flow and u_{φ} the rotational velocity component.
6.3. Expansion of solution in longitudinal Fourier modes
As the unknown of the finite element method we choose the function defined by (63)The factor ρc removes the gradients of ρ or c from the weak form of the equation and thus improves convergence (see Sect. 6.5). The equation solved by is (64)In order to use the 2.5D solver, the 3D solution must be expanded as a Fourier series in the longitudinal direction: (65)Inserting the above expansion into Eq. (64), we see that the 3D problem separates into a set of independent 2D problems, one for each m: (66)with (67)and (68)The 2D gradient and divergence operators are defined by (69)On the outer boundary, we apply a Sommerfeldlike boundary condition (70)A boundary condition that would take into account the fast variation of density near the surface (exponential decay) would be preferable, thus removing the need for an extended atmosphere and resultant mesh; this will be studied in a future paper.
6.4. Expansion of Green’s function
For a delta function source at position , (71)and the source coefficients are (72)Thus the Green’s function may be computed using (73)where each solves (74)If is on the zaxis, the only nonzero component is . It is important to note that when there is no flow (u = 0), we have and thus .
6.5. Weak form of the equations
To implement the finite element method, we derive a weak (variational) formulation of the wave equation for each m. Equation (66)is multiplied by a test function Φ and integrated over Σ with integration element ϖdΣ = ϖdϖdz.
Using integration by parts, we obtain (75)where we have used mass conservation and we have assumed that the flow does not cross the computational boundary (). The integral is a line integral over the outer boundary ∂Σ. The boundary condition is used to rewrite the boundary term: (76)To solve this problem by the FEM, the computational domain is subdivided into quadrilateral cells. On each of them, the test function Φ and the solution of the equation are projected on a basis {Φ_{i}} which, in our study, consists of piecewise continuous polynomials. By writing this formulation for all the Φ_{i} basis functions, the problem becomes a linear system that can be solved to obtain . We refer the reader to Zienkiewicz et al. (2005) for an introduction to the method, and to Chabassier & Duruflé (2016) for more specific details.
6.6. Comparison with an exact solution for a piecewise homogeneous layered medium
In order to test the rate of convergence of the finiteelement solver, Chabassier & Duruflé (2016) considered a simple benchmark for which the exact solution is known. In this setup, the background medium consists of a series of concentric spherical shells (with boundaries at radii r_{1} = 0.7R, r_{2} = R, and r_{3} = 2R, where R = 700 Mm) where the sound speed c_{i} and density ρ_{i} are constant within the ith shell. An example of a quadrilateral mesh used for the computation can be seen in Fig. 3. A Sommerfeld boundary condition is applied at the computational boundary r_{3}. The following coefficients were chosen arbitrarily and are not intended to represent the Sun: (77)Here ρ_{i} and c_{i} are in units of ρ_{⊙} and 700 km s^{1} respectively.
Fig. 3 Example of a quadrilateral mesh used to compute the solution for the scattering by spherical layers. 

Open with DEXTER 
Fig. 4 Relative L^{2} error between the FEM solution and the analytical solution for the piecewise homogeneous layered medium (Sect. 6.6 and Fig. 3) for nonaxial incidence of a plane wave at frequency ω/ 2π = 2 mHz. The error is plotted as a function of the maximum mesh size h for different polynomial orders p. The observed rate of convergence is consistent with the theoretical expectation, h^{(p + 1)}. 

Open with DEXTER 
The analytical solution for the scattering of a plane wave traveling from infinity in the + direction can be computed as a linear combination of spherical Hankel functions (Chabassier & Duruflé 2016). A nonaxial incidence has been selected (the wave vector points in the + direction) such that all the modes are excited, not only the mode m = 0. In Fig. 4, the relative L_{2} error between the numerical and the exact solutions is shown. The computation of the analytic solution was performed in multiple precision such that the 16 digits of the reference solution are exact. We obtain optimal convergence in h^{p + 1} where h is the maximum mesh size and p is the order of the polynomial basis (Zienkiewicz et al. 2005). The convergence is exponential in m (spectral accuracy). We note that 75 modes () are sufficient to achieve machine precision accuracy for the above problem.
7. Timedistance helioseismology in a solar model
7.1. Spherically symmetric reference model
We use the sound speed and density profiles from the standard solar model described by ChristensenDalsgaard et al. (1996), which is known as Model S. We interpolate the values of ρ and c on the finiteelement grid using Bsplines.
We implement a Sommerfeldlike radiation boundary condition on the computational boundary. In order to use Eq. (9), which assumes a locally uniform background medium, we match the Model S atmosphere above 500 km with a transition region (Fig. 5). For heights between 500 km and 2100 km, the acoustic cutoff smoothly transitions to zero and the sound speed to a constant value of 6.5 km s^{1}; the derivatives vanish at 2100 km. The density is deduced by integration using as a definition of the acoustic cutoff frequency. Above 2100 km the acoustic cutoff frequency, sound speed, and density are constant. By adding such a layer, waves with frequencies above 5.3 mHz propagate out through this extended region until they are absorbed by the boundary.
To achieve better numerical convergence of the FEM solution, the sound speed and density profiles of Model S are smoothed using Bsplines. Eighthorder Bsplines are defined by fifteen knots selected to capture all the variations of soundspeed and density throughout the Sun, in particular in the nearsurface layers. This ensures that the acoustic cutoff frequency is regular.
Fig. 5 Density (top), sound speed (middle) and acoustic cutoff frequency (bottom) of Model S (blue), the transition region (red), and the constant region (green, required for Sommerfeld BC). Height is computed from the photosphere in Model S. The dashed vertical lines indicate the boundaries between regions. 

Open with DEXTER 
7.2. Meshing the computational domain
Fig. 6 Computational mesh (left) used by the finite element method and the imaginary part of the Green’s function at ω/ 2π = 3 mHz (right). The Dirac source is located at radius 0.8 R_{⊙} along the zaxis. The damping rate was set to γ/ 2π = 30μHz. 

Open with DEXTER 
Fig. 7 Radial wavelength for ω/ 2π = 9 mHz and ℓ = 15 (lines) computed from Eq. (78), with the position of the cell vertices in the radial direction overplotted (crosses). The local wavelength is captured by a tenthorder polynomial within each cell. The wavelength becomes constant in the extended atmosphere described previously. 

Open with DEXTER 
Figure 6 shows the FEM mesh used throughout the remainder of the paper. Meshing the computational domain is performed in two steps. Initially, we mesh the inner part of the domain (r< 0.7 R_{⊙}) with quadrilaterals that are ~60 Mm in size. Above this inner mesh we add concentric mesh layers with a radial thickness equal to the radial wavelength (78)where ℓ = 15 is the minimum angular degree that we want to study. In the above expression we fix the frequency at ω/ 2π = 9 mHz, which is the highest frequency of interest in helioseismology. The number of points in each mesh cell is determined by the order of the polynomials we choose for the finite elements. In practice, we choose polynomials of order ten in the radial direction, corresponding to a spatial sampling of ~λ_{r}/ 10. Figure 7 shows the radial wavelength at ω/ 2π = 9 mHz, ℓ = 15, and selected cell height. In the horizontal direction, subdivisions are performed such that the horizontal length of the cells is not larger than two times their radial height. In this work, we use this mesh for all frequencies below 9 mHz. In future work, constructing meshes that are less refined for lower frequencies could be considered to reduce computational cost.
7.3. Wave attenuation
The full width at half maximum (FWHM) of a peak in the observed power spectrum is proportional to the attenuation of the mode of oscillation. In our framework the FWHM of a single peak is related to the attenuation through γ = FWHM/ 2, where the FWHM is measured in rad s^{1}. Observational studies show that the FWHM of pmode ridges is both dependent upon the harmonic degree ℓ and the frequency (e.g., Korzennik et al. 2004). For this study we restrict ourselves to a frequency dependence approach; this approach is acceptable for a filtered power spectrum that selects a wave packet in a narrow range of phase speeds or a single radial order (ridge filtering) because of the onetoone relationship between frequency and wavenumber. However, when modeling the full power spectrum, this approach can only serve as a reasonable estimate. While a wavenumberdependent damping can be implemented in principle, it is beyond the scope of this paper and we reserve it for a future study.
Fig. 8 Observed full width at half maximum of acoustic modes with radial orders (black dots, Korzennik et al. 2013; Larson & Schou 2015), overplotted with the power law approximation used in the simulations (2γ(ω), orange line). 

Open with DEXTER 
Figure 8 shows the observed values of FWHM reported by Korzennik et al. (2013) (100 <ℓ< 1000) and Larson & Schou (2015) (ℓ< 300) for p modes with radial orders and in the range 1–5.3 mHz. We approximate the attenuation coefficient with a power law in ω of the form, (79)where γ_{0}/ 2π = 4.29μHz, ω_{0}/ 2π = 3 mHz, and β = 5.77. For frequencies above the acoustic cutoff we fix the attenuation to a constant value.
Fig. 9 Snapshots of the inverse temporal Fourier transform of , propagating from a source located on the polar axis at the photosphere, where r belongs to a plane through the source and the center of the Sun. The time t, measured from the source time, is written in the top right of each panel. Each image is multiplied by ρ^{− 1/2}exp(t/ 75min) for display purposes. The values are saturated at a hundredth of the maximum value for the whole time series. The first nine panels have been further saturated (by a factor 2) to improve the visibility of the first arrival wave. 

Open with DEXTER 
7.4. Green’s function
As the behavior of the Green’s function around a Dirac source is difficult to capture accurately, the mesh is refined around the position of the source, which, for numerical stability of our solver, has to coincide with an existing mesh point. This refinement is performed by subdividing the cells neighboring the source point into three quadrilaterals so that the shape of the cells around the source is preserved (see Fig. 6). In practice, we subdivide the cells around the Dirac source five times.
With all the tools in hand, the Green’s function can be computed for any given frequency. Figure 9 shows snapshots in the time domain of the inverse Fourier transform of ImG/ω, for a source located on the zaxis at the solar surface. For this figure we computed 5000 equidistant frequencies (from 0 to 8.33 mHz) in order to cover a time span of about 7 days. In the first six panels the firstarrival wave front is seen propagating away from the source through the core towards the far side of the Sun. At t = 75 min, the secondskip waves become visible. Wave packets with higher skip numbers (which take longer to travel) are also seen at the later times. The firstarrival wave packets reaches the far side of the Sun at around 135 min and is seen at time t = 165 min propagating back towards the source.
7.5. Oscillation power spectrum
In Sect. 5 we demonstrated that under certain assumptions the crosscovariance can be directly related to the Green’s function. But does such a source of excitation produce reasonable oscillation power spectra? This is an important test as the crosscovariance function is directly connected to the oscillation power spectrum (e.g., Sekii & Shibahashi 2003). Here we compute the power spectrum in terms of Im G and compare it with observations, in the case of a spherically symmetric background medium.
We consider the observable measured at radius R and take its spherical harmonics transform: (80)where is the surface element on the unit sphere. The power spectrum is then defined by the expectation value of the squared modulus of the observable, (81)We rewrite in a frame ℛ_{0} with polar axis . In this frame, we denote by the polar angles of . The rotation of Euler angles (α,β,γ) = (π,θ_{0},π−φ_{0}) brings ℛ_{0} to the original frame. Using the rotation formula of spherical harmonics (e.g., Messiah 1960), we have (82)where, for the sake of simplicity, was assumed to depend only on angular distance Θ (horizontal isotropy and u = 0).
Fig. 10 Top: power spectrum computed at the solar surface. The blue crosses indicate the position of the p_{1}–p_{8} modes reported by Korzennik et al. (2013), while the blue dashed line shows the observed fmode ridge which is missing in this simulation owing to the lack of a gravitational term. Bottom: comparisons of the simulated power spectrum (red) and HMI data (black) at ℓ = 500. The small misalignment of the simulated ridges from the observations are due to imperfect modeling of surface layers in Model S (i.e., surface effects, Rosenthal et al. 1999). 

Open with DEXTER 
Using the explicit form of the rotation matrix elements, (83)the expression for the power spectrum simplifies to (84)where the last equality is for the special source without flow. The function P_{ℓ} is the Legendre polynomial of order ℓ.
We have complete freedom in the choice of the frequency dependence of the source power, Π(ω). In the rest of the paper we choose a Lorentzian profile: (85)where ω_{0}/ 2π = 3.3 mHz and Γ/2π = 1.2 mHz. This choice is reasonable for the purposes of this paper.
For a source on the polar axis, only the m = 0 mode needs to be computed. To avoid aliasing, the Green’s function is sampled on a highresolution grid in θ to increase the spatial Nyquist frequency. In these results, we used 20 000 grid points in θ.
Fig. 11 Section of the power spectrum for harmonic degree ℓ = 85 around frequency 3.2 mHz (n = 7–9). Red crosses are mode frequencies of the rotationally split pmodes from the first Montjoie simulation in Table 1. For comparison, the gray scale image shows the GONG observational power spectrum taken from the paper by Hill et al. (1996). Darker shades indicate larger values of the power. We note the side lobes in the observations due to aliasing from observing half of the Sun. The slope (a_{1}) of the frequencies with m is due to the average rotation rate, while the curvature (a_{3}) indicates differential rotation (slower rotation near the poles). 

Open with DEXTER 
Values of the acoefficients for the ℓ = 85 and n = 8 mode near 3.2 mHz computed from Montjoie simulations and ADIPLS eigenvalue calculations, and measured from SDO/HMI observations.
Figure 10 shows the m = 0 power spectrum computed from Eq. (84)with the source located at the photosphere on the zaxis. Here we have computed 8000 frequencies from 0 to 8.3 mHz and harmonic degrees up to 1000. This figure shows a good relationship between the mode frequencies of our simulation and those of MDI/Doppler measured by Korzennik et al. (2013). Unlike the normalmode summation method used in previous work, our power spectrum shows physical ridges above the cutoff frequency. Figure 10 also shows a slice through the power spectrum at ℓ = 500 compared with 72 days of observations from HMI/SDO. We note that at high frequencies the mode frequencies are slightly higher than the observed values. This is due to imperfect modeling of the surface layers in Model S (surface effects), not to numerical issues (the accuracy of the Green’s function is discussed later).
7.6. Frequency splittings due to differential rotation
Having demonstrated the agreement between our simulations and observations in the case of no background flow, we now turn our attention to solar rotation. We wish to check that the differential rotation of the Sun’s convection zone will introduce the correct frequency splittings between the azimuthal modes propagating in the prograde (m> 0) and retrograde (m< 0) directions. We compute the Green’s function for a photospheric source r_{1} located at the equator at longitude 0°, in the presence of a flow . We use a solarlike differential rotation model specified by (86)For comparison with the ℓ = 85 GONG power spectrum near 3.2 mHz as reported by Hill et al. (1996), we compute the Green’s function for frequencies between 2.8 and 3.6 mHz in steps of 0.5 μHz for all azimuthal orders  m  ≤ ℓ. For each value of m and ω, a power spectrum is computed by projecting the crosscovariance onto spherical harmonics as in Eq. (81). The frequencies of the modes with radial orders n = 7–9 were then extracted from each by fitting Lorentzian functions. These mode frequencies are plotted in Fig. 11 over the observational GONG power spectrum from Hill et al. (1996).
In order to quantitatively characterize the frequency splittings due to rotation, we compute the acoefficients as defined by Schou et al. (1994). The mean frequency of the multiplet ℓ = 85 and n = 8 and the first three acoefficients are given in Table 1 in five cases:

Simulation #1: Montjoie simulation using scalar waveEq. (6)with a surface deltafunction source at theequator. The acoefficients are extracted from fits to the modefrequencies measured from the simulated power spectrum.

Simulation #2: Montjoie simulation including the secondorder term u·∇(u·∇ψ) on the lefthand side of Eq. (6). We observe that the a_{2} coefficient (asphericity) vanishes.

Eigenvalue calculation for Eq. (6)with (nonrotating) Model S and a free surface boundary condition at height 0.0007 R_{⊙} above the photosphere, using a modified version of ADIPLS (ChristensenDalsgaard 2008) to compute the eigenfrequencies and the rotational kernels (ChristensenDalsgaard 2003). The modifications to ADIPLS are explained in Appendix A. The odd acoefficients are derived from the firstorder perturbation to the mode frequencies. The even acoefficients are zero to this level of approximation.

Eigenvalue calculation using the standard ADIPLS pulsation code, without neglecting terms in the oscillation equations.

Measurements of acoefficients from 360 days of SDO/HMI observations (Larson & Schou 2015).
As mentioned previously, the mean frequencies of Model S using Montjoie overestimate those of the SDO/HMI observations by ~ 13 μHz. The ADIPLS mean frequencies are also above the SDO/HMI observations by more than 10μHz. This difference comes from imperfect modeling in the nearsurface layers and is often referred to as “the surface effect” (e.g., Rosenthal et al. 1999; Ball et al. 2016). The ~3 μHz frequency difference between the Montjoie and the modified ADIPLS calculations comes from the difference in the atmospheric models.The difference between the modified ADIPLS and the standard ADIPLS frequencies comes from neglecting the buoyancy force in Eq. (2).
Fig. 12 Left: timedistance diagram computed from Eq. (87)at the height of the Dirac source. Right: observational timedistance diagram computed from the Fourier transform of the SOHO/MDI/Doppler mediumdegree power spectrum (Kosovichev et al. 1997). The SOHO/MDI timedistance diagram fades away at large separation distances due to foreshortening. 

Open with DEXTER 
The simulated a_{1} and a_{3} coefficients are of the expected sign and order of magnitude, within a few nHz of each other. The simulated a_{1} coefficients are about 5 nHz smaller than the SDO/HMI observed value even though we did not tune the rotation profile in the simulations. The simulated a_{3} coefficients are also smaller than the observed value, by ~2 nHz.
The value of a_{2} in simulation #1 using Eq. (6)is nonzero, which was (at first) unexpected since our model does not include centrifugal distortion. The nonzero a_{2} is due to the fact that the eigenfunctions are affected by rotation at first order and thus leave a signature in the power spectrum. Adding the term u·∇(u·∇ψ) in simulation #2 restores the eastwest antisymmetry of the advection of the waves by the flow.
Overall, this comparison between simulated and observed mode frequencies is very encouraging.
7.7. Timedistance diagram
For a spherically symmetric solar model, the expectation value of the crosscovariance function in the time domain is (87)where Θ is the angular distance on the surface between the two observation points. The crosscovariance function is also called the timedistance diagram after Duvall et al. (1993). In Fig. 12 we compare the timedistance diagram computed from our power spectrum to an observed timedistance diagram using SOHO/MDI medium degree data (Kosovichev et al. 2000). In order to make this comparison we applied a spatial filter to the simulated power spectrum, F_{ℓ} = [1−tanh^{(}0.03ℓ−3^{)}] /2 for ℓ< 100 and 0 otherwise, to remove highdegree modes.
Comparisons of the two timedistance diagrams is encouraging. However, the amplitude of the backskip ridge at t ~ 250 min is greater in the observations than in the simulations, for which we have no definitive explanation. We think that the most likely explanation is that the damping of the low degree modes is overestimated in our model, resulting in a reduced amplitude of the backskip branch in the timedistance diagram. Further tuning of the power spectrum is required in order to resolve this discrepancy. For further comparison with the observations, Fig. 13 shows time plots at three different travel distances. The widths and relative amplitudes of the first few skips are in general agreement with the observations.
Fig. 13 Temporal crosscovariance function for three angular distances Θ = 30°, 60°, and 90° for the simulations (blue) and the SOHO/MDI Doppler observations (red dashes, Kosovichev et al. 2000). The temporal window functions (w) used in the definitions of travel times are shown in black. 

Open with DEXTER 
Fig. 14 Relative difference between Im G^{(p)}(Θ,ω) and Im G_{ref}(Θ,ω) as defined by Eq. (88)for frequencies ω/ 2π = 3 mHz (left) and ω/ 2π = 7 mHz (right), and different angular distances Θ = 30°, 60° and 90°. Computations were performed without flow (solid lines) and with a meridional flow (dashed lines). 

Open with DEXTER 
8. Validation for helioseismology applications
8.1. Convergence of Green’s function
In order to estimate the accuracy of the forward solver, we first measure the convergence of the Green’s function towards a solution G_{ref} computed for a highly refined mesh (four times more cells) with highorder discretization (order 13). This solution is used as reference since we cannot determine the exact solution to our problem. By choosing a basis of polynomials of order 13 for the finite elements on the refined mesh, the number of degrees of freedom per wavelength is 26 (13 × 2 cells per wavelength), i.e. many more than the 10 points per wavelength used previously.
We compute the Green’s function G^{(p)} by choosing polynomials of order p in the nonrefined mesh containing 10 729 cells. The relative difference between Im G^{(p)} and Im G_{ref}, denoted by ε^{(p)}, is plotted in Fig. 14 as a function of the number of degrees of freedom N^{(p)} = 10729(p + 1)^{2}. Explicitly, (88)In Fig. 14, we plot ε^{(p)} for different angular distances between the source and the receiver Θ (30°, 60°, and 90°) and different frequencies (3 mHz and 7 mHz). We see that ε^{(p)} reaches ~10^{5} for orders of discretization p> 10. We obtain a similar convergence when a meridional flow, as described in Appendix B (with surface velocity U = 20 ms^{1}), is added to the background.
Fig. 15 Travel times δτ^{(p)} as defined in Eq. (89), computed from the difference between crosscovariances C^{(p)}(Θ,ω) computed using polynomials of degrees and a reference crosscovariance C_{ref}(Θ,ω). Computations with and without flow are shown by the dashed and solid lines, respectively. The vertical line labeled “p = 10” indicates the number of degrees of freedom used in all calculations in this paper. 

Open with DEXTER 
Computational times and memory usage for singlefrequency runs of the radial (1.5D), axisymmetric (2.5D), and 3D methods for concentric shells with constant background coefficients and for the solar model.
8.2. Convergence of travel times
Having discussed the accuracy of the Green’s functions at different orders of discretization, we now examine the accuracy of the travel times. Based on Eq. (25), we compute the travel times for waves originating from the pole defined by (89)where C^{(p)} is the crosscovariance computed from the Green’s functions G^{(p)} and C_{ref} is the crosscovariance computed from the G_{ref}, as defined in the previous section. The Green’s functions were computed for a Nyquist frequency of 8.33 mHz with a frequency resolution of 3.3μHz and a constant damping rate of γ/ 2π = 30 μHz. In Fig. 15, we show that our method achieves a traveltime accuracy of 8 ms for the order of discretization p = 10, with or without the presence of a background meridional flow, as in the previous section. This accuracy could be improved by a factor of 10 if we used polynomials of degree 12; however, the CPU time and memory requirements are increased by 100% and 40%.
8.3. Computational times for the 1.5D, 2.5D and 3D problems
Here we compare the computational times of the 2.5D model with two other models (3D and radial 1.5D). Initially, we examine a simple case of constant sound speed and density spherical layers with the Sommerfeld radiation boundary condition described in Sect. 6.6. This is done in order to demonstrate the computational costs in all three model types. Following this, we compute the computational costs for the more demanding solar cases (1.5D and 2.5D) and neglect the full 3D case because of the high cost. However, we note that computational costs must be multiplied by the number of modes and frequencies required to accurately reconstruct the Green’s function.
Table 2 shows CPU time and memory requirements of the three different models for a single frequency. With each additional dimension, the required memory and CPU time become larger, with a dramatic increase in the full 3D case. A comparison of the model requirements shows that even though the axisymmetric method is more demanding than the radial one, the requirements are not unreasonable for most systems.
9. Traveltime sensitivity kernels
Fig. 16 Sensitivity of the mean travel time, [τ(r_{1},r_{2}) + τ(r_{2},r_{1})] /2, to relative soundspeed perturbations δc/c in the interior, with observation point r_{1} located on the polar axis and observation point r_{2} at latitude of 45°. Panel a): cut at r = 0.95 R_{⊙} through the soundspeed kernel. Panel b): slice through a plane containing the two observation points and the center of the Sun. The green line indicates the position of the ray path. Panel c): slice in a plane perpendicular to the ray path connecting the two observation points (green square) at equal distance from the two observation points. In both slices b) and c), the black arc of a circle locates radius r = 0.95 R_{⊙} which corresponds to panel a). 

Open with DEXTER 
In the results discussed thus far, we have focused on obtaining travel times through the direct modeling of waves propagating in an axisymmetric medium. In this section we address the computation and accuracy of the traveltime kernels outlined in Sect. 4 that describe the spatial sensitivity of travel times to local perturbations in the interior.
9.1. Threedimensional kernels
Under the “convenient source” assumption, the computation of 3D kernels requires computing four Green’s functions, as explained in Sect. 6. When the medium does not contain a flow, only two Green’s functions are needed. Even in this case, the computational burden is very demanding for a general 3D background (Sect. 8.3). However, it is feasible to compute 3D kernels when the background model is axially or spherically symmetric, under the 2.5D approach outlined thus far.
In the case of spherical symmetry of the background, the kernels can be built from a single Green’s function where the source is located on the pole at the observation height, after a series of rotations. Only the mode m = 0 needs to be computed in this case. This m = 0 Green’s function is then rotated to the desired source location (r_{1}) and a duplicate is rotated to the receiver location (r_{2}). The construction of the kernels follows Eqs. (55)–(58).
Figure 16 shows slices through a soundspeed kernel with r_{1} at the pole and r_{2} at 45° latitude, computed using the rotation of a m = 0 Green’s function and 800 frequencies between 1.5 mHz and 4.5 mHz. In this kernel we see the traditional bananadoughnut shape reported in geophysics (e.g., Marquering et al. 1999) and helioseismology (Birch et al. 2004). For the choice of observable that we have made, the traveltime sensitivity is very small near the ray path (Fig. 16c). Surrounding the ray path are regions of negative and positive sensitivities corresponding to the consecutive Fresnel zones. The values near the surface are not numerical noise but are due to high spatial frequencies. Figure 16 is a 3D illustration of the soundspeed kernel and also shows the values on the sphere at radius r = 0.95 R_{⊙}. The sensitivity is maximum near the surface around the two observation points. We note that we have plotted the product cK_{c} to better render the deeper layers. The computation of the m = 0 Green’s function took approximately 1 h using 200 cores. The postprocessing consists mostly in computing the rotated Green’s functions and took approximately 3 h.
9.2. Longitudinally averaged kernels
Fig. 17 Left and center panels: kernels ⟨ K_{ur} ⟩ and ⟨ K_{uθ} ⟩ for the r and θ components of the flow. Point r_{1} is at the north pole (photosphere) and point r_{2} is at 45° latitude. The values of the kernels are scaled by the sound speed c and are saturated at 1/400 of the maximum value. The ray path connecting the two points is shown (thick black line) as well as the computational boundary (black half circle). Right panel: comparison between the travel times computed directly from the crosscovariance function (curves, see Eq. (25)) and those computed from the sensitivity kernels (“+” symbols, see Eq. (90)). The blue curve is for the travel times measured from the pole to latitude 45°, the red curve for the travel times measured in the opposite direction. The accuracy of the travel times is approximately 10^{3} s. 

Open with DEXTER 
For an axisymmetric (but not necessary spherically symmetric) background, the Green’s functions must be constructed by summing a sufficient number of m modes. Here, in addition, we only consider axisymmetric perturbations, q_{α} = q_{α}(r,θ), and we wish to determine the 2D spatial sensitivity of the travel times by averaging the kernels over longitude. Useful applications include rotation (u_{φ}) and meridional circulation (u_{θ}). For axisymmetric perturbations, we have (90)where ⟨ K_{α} ⟩ is the longitudinally averaged kernel (91)Let us define the Green’s function and the crosscovariance function in terms of their longitudinal mode components G^{m} and C^{m} as follows: (92)It is important to note the order of the variables in G and in the above notations. Using the formulation of the kernels given by Eq. (43), we obtain (93)Using the explicit expressions for the bilinear operators ℒ_{α} (Eqs. (55)(58)), we can then obtain the longitudinally averaged kernels for all perturbations q_{α}.
As an example, the flow kernels can be written as a sum (94)over azimuthal components: (95)where the operator is either , , or .
In practice, is obtained using generalized seismic reciprocity, , performing a simulation with a source at . Using the convenient source of excitation, the crosscovariance is linked to the Green’s function by Eq. (49). It is possible to obtain a similar relation for the Fourier coefficients: (96)This means that the (3D) kernels can be computed using only the azimuthal modes of the Green’s function obtained from the 2.5D solver.
Several comments can be made:

If one of the observation points is on the rotation axis, then only the m = 0mode of the Green’s function is required. For a measurementbetween two arbitrary points at the surface of the Sun, thecomputation of many modes is required (see next section).

If the background model contains a flow, then the Green’s function for a reversed flow is also required in order to compute the crosscovariance as shown by Eq. (96).

In the case of no background flow, Eq. (66)shows that only the m ≥ 0 solutions need to be computed since G^{− m} = G^{m}.

If the kernels for points (θ_{1},φ_{1}) and (θ_{2},φ_{2}) have been computed, then we obtain for free the kernels for points with different longitudes since the only term in Eq. (95)that depends on the longitudes is the exponential e^{im(φ2−φ1)}.
9.3. Accuracy of traveltime kernels for flows
Fig. 18 Left and center panels: kernels ⟨ K_{ur} ⟩ and ⟨ K_{uθ} ⟩ for the r and θ components of the flow, using all azimuthal components  m  ≤ m_{max} = 35. The separation distance between points r_{1} and r_{2} is 42°; the center point is located at a latitude of 40°. The ray path connecting the two points is shown in black. The right panel shows the convergence of the travel times as a function of m_{max}, for the northsouth (blue) and southnorth (green) travel directions. The red curve shows their difference. 

Open with DEXTER 
Consider a spherically symmetric background reference model with no flow. Figure 17 shows the flow kernels ⟨ K_{ur} ⟩ = ⟨ K_{ur} ⟩ ^{m = 0} and ⟨ K_{uθ} ⟩ = ⟨ K_{uθ} ⟩ ^{m = 0} with a point r_{1} on the pole (at the solar surface) and a point r_{2} located at latitude 45°. The kernel ⟨ K_{uφ} ⟩ ^{m = 0} is zero by construction because ⟨ K_{uφ} ⟩ ^{m} is proportional to m. The 2D kernels for u_{r} and u_{θ} display Fresnel zones as do the 3D kernels, however, the null points along the ray path are absent here because of integration in the longitudinal direction.Additionally, as reported by Birch & Gizon (2007), the kernels exhibit hyperbolalike patterns near point r_{1}. This pattern is due to scattering from distant sources (Gizon & Birch 2002) and is not present in Earth seismology kernels for pointsource earthquakes. We refer the reader to the work by Duvall et al. (2006) for an observational study of 2D horizontal sensitivity kernels. The righthand panel of Fig. 17 compares the travel times computed in two different ways to evaluate the accuracy of the kernels: (1) by multiplying the kernels by a flow model u and integrating (Eq. (90)) and (2) by computing the difference δC = C(u)−C(u = 0) and then measuring the travel time (Eq. (25)). Here we have used the meridional flow model shown in Appendix B with a maximum flow speed of 20 ms^{1} at the surface. For the comparison we have considered three separation distances (30°, 45°, and 60°) and the two directions, δτ(r_{1},r_{2}) and δτ(r_{2},r_{1}). The kernelbased computation of the travel times and the direct computation agree to within 10^{3} s. This is an important result as it demonstrates that our kernels have sufficient accuracy for the interpretation of solar travel times. For meridional circulation measurements, the noise in the travel times is typically in the range 0.1–0.5 s, found by averaging data over four years (Rajaguru & Antia 2015).
For practical applications, the two points should be offaxis since helioseismic observations are limited to a centertolimb distance of ~70°. In such cases the azimuthally averaged kernels require the computation of many components ⟨ K ⟩ ^{m}. Figure 18 shows the components of flow kernel computed from Eq. (95)for all m ≤ 35. In this figure the travel times are measured between points at latitudes 61° and 19°, both along the central meridian. With this separation distance of 42° the ray path reaches a depth of 0.72 R_{⊙}. Like before, the kernels are not symmetric about the center point between the observation points and have features similar those seen in Fig. 17.
To test the convergence of these kernels, we calculate individual kernels including all modes  m  ≤ m_{max} and compute travel times as a function of m_{max} in the presence of the same meridional flow as used previously. The righthand panel of Fig. 18 shows that the travel times converge to an asymptotic value for m_{max}> 25 with an accuracy of ~0.01 s. We note that a larger m_{max} is needed to achieve convergence for shorter separation distances.
9.4. Filtering
Filtering the observations in the ℓω domain is common practice in helioseismology. Several choices of filters have been proposed: filters in ω space, phasespeed filters, and ridge filters (see, e.g., Gizon & Birch 2005). To interpret any particular traveltime measurement, the sensitivity kernels must account for the proper frequency content of the seismic data set by the filtering. This dependency of the kernel on the filter has been studied previously (e.g., Birch et al. 2004; Böning et al. 2016).
Fig. 19 Left: unfiltered ⟨ K_{uθ} ⟩ kernel with r_{1} at the pole and r_{2} at a colatitude of 15.36°. Right: kernel ⟨ K_{uθ} ⟩ for filtered observations, where the Gaussian phasespeed filter F_{ℓ}(ω) is centered at 125.2 km s^{1} with a dispersion of 12.3 km s^{1}, for ℓ up to 1000. The ray path is shown by the black line. 

Open with DEXTER 
Symbolically, the filtered observation, Ψ, is obtained by applying a filtering operator, F, to the original wavefield, ψ: (97)The corresponding kernel is obtained by applying the filtering operator twice to the original kernel: (98)where the filtered Green’s function and crosscovariance are and the W function is computed with the choice C_{ref} =C. Here, F_{1} indicates that the filtering has to be done with respect to the point r_{1}.
The filtered Green’s function may be obtained by filtering the delta source function. To see this, we use generalized seismic reciprocity (101)where (102)is the filtered delta function source. In the above expression F_{ℓ}(ω) can be either a phasespeed filter or a ridge filter. We note that the filtered source function is a function of frequency.
If the background is spherically symmetric with no background flow, it is also possible to compute the nonfiltered Green’s function and to perform the filtering a posteriori. In this case, the Green’s function G(r_{1},r,ω) depends only on the depths r_{1} and r and on the angular distance Θ_{1} between the two points (103)where . The filtered Green’s function takes the form (104)where (105)is the projection of G on the Legendre polynomials. Using seismic reciprocity we obtain G_{ℓ}(r_{1},r,ω) = G_{ℓ}(r,r_{1},ω), and Eq. (104) implies (106)Thus, if the background is spherically symmetric, we also have seismic reciprocity for the filtered Green’s function: the filtering can be seen as a postprocessing operation on r at fixed source position r_{1}. Otherwise, a filtered source should be used, as given by Eq. (102).
As a special case, we computed the azimuthally averaged kernels of Eq. (95)for filtered observations when r_{1} is located at the pole and u = 0. As mentioned previously, only the m = 0 component of the Green’s function is needed. Figure 19 shows the effect of applying a phasespeed filter to the kernel ⟨ K_{uθ} ⟩ with r_{1} at the pole and r_{2} at a colatitude of 15.36°. The phasespeed filter is centered at 125.2 km s^{1} with a width of 12.3 km s^{1}. These values were chosen to be the same as for filter #11 in Duvall & Hanasoge (2013). Once the filter is applied the sensitivity to u_{θ} is predominantly near the base of the ray path. This filtered kernel appears to be similar to the K_{6} case from Böning et al. (2016) indicating general agreement with their work.
10. Conclusion
We have presented a new framework for computational helioseismology by solving the forward problem in the frequency domain. For the sake of simplicity, we considered a simplified scalar acoustic wave equation and assumed spatially uncorrelated sources of excitation distributed through the Sun. Under such conditions, the crosscovariance can be obtained directly from the imaginary part of the frequencydomain Green’s function. This leads to a convenient, flexible, and fast way to compute accurate kernels. The analytical work involved in this framework is less cumbersome than previous work that relies on normalmode expansions of the kernels (e.g., Birch & Gizon 2007; Burston et al. 2015; Böning et al. 2016). The framework can be extended to the vectorial wave equation quite easily using existing Montjoie vectorial setups (e.g. Péron et al. 2016, for Maxwell’s equations). The scalar equation captures the propagation of acoustic waves through a solarlike medium and leads to an oscillation power spectrum that compares well with observations. The present setup will be very useful to test new methods, including instrumental and projection effects, but also to interpret existing traveltime measurements for rotation, meridional circulation, and axisymmetric structures like the average supergranule.
In future work we intend to address the inverse problem. Specifically, we wish to find the parameters δq of the background model (sound speed, density, flows) such that the travel times τ from the model are consistent with the observed travel times τ_{obs}. This is generally done by solving a linear system of the form δτ = K δq + n, where K is a matrix of kernels and n is a vector of traveltime noises, defined through the noise covariance matrix Λ = E [nn^{T}] (see Gizon & Birch 2004; Fournier et al. 2014). This problem can be solved by classical regularization methods such as regularized least square (RLS; e.g., Kosovichev 1996) or optimally localized averaging (OLA, e.g., Haber et al. 2004). Another approach is to solve a nonlinear inverse problem defined in terms of the partial differential equation. Different methods exist such as the Landweber iteration (Hanke et al. 1995) and Newtontype methods such as the iteratively regularized GaussNewton (Bakushinskii 1992) or Newton conjugate gradient methods (Hanke 1997). Most of these methods avoid the explicit computation of sensitivity kernels. All these inverse methods are feasible under the assumption of an axisymmetric background model, thanks to the embarrassingly parallel workload in m and ω. The last iteration of the inversion produces a 3D model of the solar internal properties. Using a full 3D forward solver is currently not practical.
Acknowledgments
L.G. and A.C.B. developed the general theoretical framework. H.B., J.C., and M.D. developed tools to solve the wave equation (Montjoie package). C.S.H., M.L., and D.F. implemented the tools in the solar context, tuned the power spectrum, and computed the kernels. All authors contributed to the writing of the paper. L.G. acknowledges generous financial support from the State of Lower Saxony, Germany, and partial support from the Center for Space Science at the NYU Abu Dhabi Institute, UAE, under grant G1502. The computer infrastructure was provided in part by the German Data Center for SDO funded by the German Aerospace Center (DLR). The research leading to these results has received funding from the EU FP7/20072013 program under grant agreement No. 312844 (SpaceInn). The Global Oscillation Network Group is funded by the National Science Foundation. The HMI project is supported by NASA contract NAS502139. Part of the numerical experiments presented in this paper were carried out using the PLAFRIM experimental testbed, developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI, IMB, Conseil régional d’Aquitaine, Université de Bordeaux, CNRS, and ANR (Programme des investissements d’avenir).
References
 Bakushinskii, A. 1992, Computational Mathematics and Mathematical Physics, 32, 1353 [Google Scholar]
 Ball, W. H., Beeck, B., Cameron, R. H., & Gizon, L. 2016, A&A, 592, A159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birch, A. C., & Gizon, L. 2007, Astron. Nachr., 328, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Bergot, M., Cohen, G., & Duruflé, M. 2010, J. Sci. Comput., 42, 345 [CrossRef] [MathSciNet] [Google Scholar]
 Birch, A. C., Gizon, L., Hindman, B. W., & Haber, D. A. 2007, ApJ, 662, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Birch, A. C., Kosovichev, A. G., & Duvall, Jr., T. L. 2004, ApJ, 608, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Bogdan, T. J., Hindman, B. W., Cally, P. S., & Charbonneau, P. 1996, ApJ, 465, 406 [NASA ADS] [CrossRef] [Google Scholar]
 Böning, V. G. A., Roth, M., Zima, W., Birch, A. C., & Gizon, L. 2016, ApJ, 824, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Bottero, A., Cristini, P., Komatitsch, D., & Asch, M. 2016, J. Acoust. Soc. Am., 140, 3520 [NASA ADS] [CrossRef] [Google Scholar]
 Burston, R., Gizon, L., & Birch, A. C. 2015, Space Sci. Rev., 196, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R., Gizon, L., & Duvall, Jr., T. L. 2008, Sol. Phys., 251, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Chabassier, J., & Duruflé, M. 2016, Research Report 8893, INRIA, https://hal.inria.fr/hal01295077 [Google Scholar]
 ChristensenDalsgaard, J. 2003, Lecture Notes on Stellar Oscillations, 5th edn. (Aarhus University) [Google Scholar]
 ChristensenDalsgaard, J. 2008, Ap&SS, 316, 113 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Dikpati, M., & Choudhuri, A. R. 1995, Sol. Phys., 161, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Duruflé, M. 2006, Doctoral thesis, ENSTA ParisTech, France [Google Scholar]
 Duvall, Jr., T. L., & Hanasoge, S. M. 2013, Sol. Phys., 287, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., Jefferies, S. M., Harvey, J. W., Osaki, Y., & Pomerantz, M. A. 1993, ApJ, 410, 829 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., Birch, A. C., & Gizon, L. 2006, ApJ, 646, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Fournier, D., Gizon, L., Hohage, T., & Birch, A. C. 2014, A&A, 567, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gizon, L., & Birch, A. C. 2002, ApJ, 571, 966 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2005, Liv. Rev. Sol. Phys., 2, 6 [Google Scholar]
 Gizon, L., Birch, A. C., & Spruit, H. C. 2010, ARA&A, 48, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Haber, D., Hindman, B., Toomre, J., & Thompson, M. 2004, Sol. Phys., 220, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., & Duvall, Jr., T. L. 2007, Astron. Nachr., 328, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Birch, A., Gizon, L., & Tromp, J. 2011, ApJ, 738, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Duvall, Jr., T. L., & Sreenivasan, K. R. 2012, Proc. Nat. Acad. Sci., 109, 11928 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Ann. Rev. Fluid Mech., 48, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Hanke, M. 1997, Numerical Functional Analysis and Optimization, 18, 971 [CrossRef] [Google Scholar]
 Hanke, M., Neubauer, A., & Scherzer, O. 1995, Numerische Mathematik, 72, 21 [CrossRef] [MathSciNet] [Google Scholar]
 Hill, F., Stark, P. B., Stebbins, R. T., et al. 1996, Science, 272, 1292 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Jackiewicz, J., Gizon, L., Birch, A. C., & Duvall, Jr., T. L. 2007, ApJ, 671, 1051 [NASA ADS] [CrossRef] [Google Scholar]
 Jackiewicz, J., Gizon, L., & Birch, A. C. 2008, Sol. Phys., 251, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Jouve, L., Brun, A. S., Arlt, R., et al. 2008, A&A, 483, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Korzennik, S. G., RabelloSoares, M. C., & Schou, J. 2004, ApJ, 602, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Korzennik, S. G., RabelloSoares, M. C., Schou, J., & Larson, T. P. 2013, ApJ, 772, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Kosovichev, A. 1996, ApJ, 461, L55 [NASA ADS] [CrossRef] [Google Scholar]
 Kosovichev, A. G., Duvall, Jr., T. L., & Scherrer, P. H. 2000, Sol. Phys., 192, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Kosovichev, A. G., Schou, J., Scherrer, P. H., et al. 1997, Sol. Phys., 170, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, P., Duvall, T. L. Jr., Harvey, J. W., et al. 1990, in Proc. of Conf. on Progress of Seismology of the Sun and Stars, eds. Y. Osaki, & H. Shibahashi (Berlin Springer Verlag), Lect. Notes Phys., 367, 87 [Google Scholar]
 Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Larson, T. P., & Schou, J. 2015, Sol. Phys., 290, 3221 [NASA ADS] [CrossRef] [Google Scholar]
 Liang, Z.C., & Chou, D.Y. 2015, ApJ, 809, 150 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Marquering, H., Dahlen, F. A., & Nolet, G. 1999, Geophys. J. Int., 137, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Messiah, A. 1960, Mécanique quantique, Vol. 2 (Paris: Dunod) [Google Scholar]
 NissenMeyer, T., van Driel, M., Stähler, S. C., et al. 2014, Solid Earth, 5, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Papini, E., Gizon, L., & Birch, A. C. 2014, Sol. Phys., 289, 1919 [NASA ADS] [CrossRef] [Google Scholar]
 Péron, V., Schmidt, K., & Duruflé, M. 2016, SIAM J. Appl. Math., 76, 1031 [CrossRef] [Google Scholar]
 Rajaguru, S. P., & Antia, H. M. 2015, ApJ, 813, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689 [NASA ADS] [Google Scholar]
 Schou, J., ChristensenDalsgaard, J., & Thompson, M. J. 1994, ApJ, 433, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Schunker, H., Cameron, R. H., Gizon, L., & Moradi, H. 2011, Sol. Phys., 271, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Sekii, T., & Shibahashi, H. 2003, in GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. SawayaLacoste, ESA SP, 517, 389 [Google Scholar]
 Snieder, R. 2007, J. Acoust. Soc. Am., 121, 2637 [NASA ADS] [CrossRef] [Google Scholar]
 Snieder, R., & Larose, E. 2013, Annu. Rev. Earth Planet. Sci., 41, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Snieder, R., Miyazawa, M., Slob, E., Vasconcelos, I., & Wapenaar, K. 2009, Surveys in Geophys., 30, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Švanda, M., Gizon, L., Hanasoge, S. M., & Ustyugov, S. D. 2011, A&A, 530, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Ballegooijen, A. A., & Choudhuri, A. R. 1988, ApJ, 333, 965 [NASA ADS] [CrossRef] [Google Scholar]
 van Driel, M., & NissenMeyer, T. 2014, Geophys. J. Int., 199, 880 [NASA ADS] [CrossRef] [Google Scholar]
 van Driel, M., Krischer, L., Stähler, S. C., Hosseini, K., & NissenMeyer, T. 2015, Solid Earth, 6, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, J., Bogart, R. S., Kosovichev, A. G., Duvall, Jr., T. L., & Hartlep, T. 2013, ApJ, 774, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Zienkiewicz, O., Taylor, R., & Zhu, J. 2005, The Finite Element Method Set, sixth edn. (Oxford: ButterworthHeinemann) [Google Scholar]
Appendix A: Solving the eigenvalue problem using ADIPLS
We modified ADIPLS to solve the eigenvalue problem (A.1)where ℋ is given by Eq. (2) and we neglected the gravity terms. Using the same notation as the one used by ChristensenDalsgaard (2008), the corresponding eigenvalue problem for the eigenfrequency and the eigenfunction of a mode takes the form where L^{2} = ℓ(ℓ + 1), and ξ_{r}(r) and ξ_{h}(r) are the radial and horizontal eigenfunctions. The equations can be rewritten as which, in ADIPLS adimensionalized form (see ChristensenDalsgaard 2008), become We modified ADIPLS, by changing the file rhs.n.d.f in such a way that it solves Eqs. (A.6) and (A.7), with a free surface boundary condition.
We calculated the frequency splittings in presence of differential rotation by evaluating the integral (A.8)where the rotation profile Ω(r,θ) is defined by Eq. (86), and K_{nℓm}(r,θ) are the rotational sensitivity kernels (see Chap. 8.4 in ChristensenDalsgaard 2003), which depend on ξ_{r} and ξ_{h}.
Appendix B: Meridional flow model
Various models of meridional flow cells can be found in the literature (e.g., van Ballegooijen & Choudhuri 1988; Dikpati & Choudhuri 1995; Jouve et al. 2008). In these models the flow is expressed through a stream function Ψ as follows: (B.1)This enforces mass conservation, ∇·(ρu) = 0. However, we cannot use directly the various expressions of the meridional flow from previous work, as they were not computed with the density ρ(r) of Model S.
Fig. B.1 u_{r}(r,θ) (top) and u_{θ}(r,θ) (bottom) at θ = 45° of the meridional flow model discussed in Appendix B for a maximum flow velocity U = 20 m s^{1} at the surface. 

Open with DEXTER 
Let us define a stream function Ψ(r,θ) for and , where r_{b} and r_{t} refer to the bottom and the top of the convection zone, respectively. The function Ψ is set to zero outside this region. Let us look for a solution of the form (B.2)where f(r) and g(θ) are dimensionless functions to be determined, U sets the amplitude of the flow, and ρ_{t} = ρ(r_{t}) is the density at the top of the cell. The flow is given by For one cell per hemisphere with r_{b}<r<r_{t}, the functions f and g satisfy the following conditions: (B.5)For the latitudinal dependence, we choose (B.6)We seek the function f in terms of the function (B.7)which controls the radial profile of ρu_{θ}. We have (B.8)The description of the flow cell then relies on the choice of the function h(r). To proceed, we choose the following conditions: (B.9)While the first two conditions are intuitive, the third is arbitrary (other choices would be possible). We choose h as a fourthorder polynomial (B.10)where is implied by h(r_{t}) = 1.
In this paper, we set r_{b} = 0.7 R_{⊙} and r_{t} = R_{⊙}. The depth r_{h} at which the horizontal flow switches sign is such that f(r_{h}) = 0. The value of r_{h} can be obtained by interpolation or by a Newton method. We find r_{h} = 0.859 R_{⊙}. For the amplitude of the flow, we take U = 20 m/s. The radial and colatitudinal flows, given by are plotted as a function of r at colatitude θ = 45° in Fig. B.1.
All Tables
Values of the acoefficients for the ℓ = 85 and n = 8 mode near 3.2 mHz computed from Montjoie simulations and ADIPLS eigenvalue calculations, and measured from SDO/HMI observations.
Computational times and memory usage for singlefrequency runs of the radial (1.5D), axisymmetric (2.5D), and 3D methods for concentric shells with constant background coefficients and for the solar model.
All Figures
Fig. 1 Geometrical setup for a background medium symmetric about an axis ẑ. The thick contour delineates the boundary, S, of the computational domain V. A position vector r in V is specified by sphericalpolar coordinates (r, θ, φ). The seismic observable ψ(r) is measured near the photosphere. 

Open with DEXTER  
In the text 
Fig. 2 Twodimensional generating section Σ of the volume V, delineated by ∂Σ (thick contour). A point in Σ is specified by coordinates (r,θ), where r is the radius and θ is the colatitude, or, equivalently, by coordinates (ϖ, z), where ϖ = rsinθ is the distance to the axis. 

Open with DEXTER  
In the text 
Fig. 3 Example of a quadrilateral mesh used to compute the solution for the scattering by spherical layers. 

Open with DEXTER  
In the text 
Fig. 4 Relative L^{2} error between the FEM solution and the analytical solution for the piecewise homogeneous layered medium (Sect. 6.6 and Fig. 3) for nonaxial incidence of a plane wave at frequency ω/ 2π = 2 mHz. The error is plotted as a function of the maximum mesh size h for different polynomial orders p. The observed rate of convergence is consistent with the theoretical expectation, h^{(p + 1)}. 

Open with DEXTER  
In the text 
Fig. 5 Density (top), sound speed (middle) and acoustic cutoff frequency (bottom) of Model S (blue), the transition region (red), and the constant region (green, required for Sommerfeld BC). Height is computed from the photosphere in Model S. The dashed vertical lines indicate the boundaries between regions. 

Open with DEXTER  
In the text 
Fig. 6 Computational mesh (left) used by the finite element method and the imaginary part of the Green’s function at ω/ 2π = 3 mHz (right). The Dirac source is located at radius 0.8 R_{⊙} along the zaxis. The damping rate was set to γ/ 2π = 30μHz. 

Open with DEXTER  
In the text 
Fig. 7 Radial wavelength for ω/ 2π = 9 mHz and ℓ = 15 (lines) computed from Eq. (78), with the position of the cell vertices in the radial direction overplotted (crosses). The local wavelength is captured by a tenthorder polynomial within each cell. The wavelength becomes constant in the extended atmosphere described previously. 

Open with DEXTER  
In the text 
Fig. 8 Observed full width at half maximum of acoustic modes with radial orders (black dots, Korzennik et al. 2013; Larson & Schou 2015), overplotted with the power law approximation used in the simulations (2γ(ω), orange line). 

Open with DEXTER  
In the text 
Fig. 9 Snapshots of the inverse temporal Fourier transform of , propagating from a source located on the polar axis at the photosphere, where r belongs to a plane through the source and the center of the Sun. The time t, measured from the source time, is written in the top right of each panel. Each image is multiplied by ρ^{− 1/2}exp(t/ 75min) for display purposes. The values are saturated at a hundredth of the maximum value for the whole time series. The first nine panels have been further saturated (by a factor 2) to improve the visibility of the first arrival wave. 

Open with DEXTER  
In the text 
Fig. 10 Top: power spectrum computed at the solar surface. The blue crosses indicate the position of the p_{1}–p_{8} modes reported by Korzennik et al. (2013), while the blue dashed line shows the observed fmode ridge which is missing in this simulation owing to the lack of a gravitational term. Bottom: comparisons of the simulated power spectrum (red) and HMI data (black) at ℓ = 500. The small misalignment of the simulated ridges from the observations are due to imperfect modeling of surface layers in Model S (i.e., surface effects, Rosenthal et al. 1999). 

Open with DEXTER  
In the text 
Fig. 11 Section of the power spectrum for harmonic degree ℓ = 85 around frequency 3.2 mHz (n = 7–9). Red crosses are mode frequencies of the rotationally split pmodes from the first Montjoie simulation in Table 1. For comparison, the gray scale image shows the GONG observational power spectrum taken from the paper by Hill et al. (1996). Darker shades indicate larger values of the power. We note the side lobes in the observations due to aliasing from observing half of the Sun. The slope (a_{1}) of the frequencies with m is due to the average rotation rate, while the curvature (a_{3}) indicates differential rotation (slower rotation near the poles). 

Open with DEXTER  
In the text 
Fig. 12 Left: timedistance diagram computed from Eq. (87)at the height of the Dirac source. Right: observational timedistance diagram computed from the Fourier transform of the SOHO/MDI/Doppler mediumdegree power spectrum (Kosovichev et al. 1997). The SOHO/MDI timedistance diagram fades away at large separation distances due to foreshortening. 

Open with DEXTER  
In the text 
Fig. 13 Temporal crosscovariance function for three angular distances Θ = 30°, 60°, and 90° for the simulations (blue) and the SOHO/MDI Doppler observations (red dashes, Kosovichev et al. 2000). The temporal window functions (w) used in the definitions of travel times are shown in black. 

Open with DEXTER  
In the text 
Fig. 14 Relative difference between Im G^{(p)}(Θ,ω) and Im G_{ref}(Θ,ω) as defined by Eq. (88)for frequencies ω/ 2π = 3 mHz (left) and ω/ 2π = 7 mHz (right), and different angular distances Θ = 30°, 60° and 90°. Computations were performed without flow (solid lines) and with a meridional flow (dashed lines). 

Open with DEXTER  
In the text 
Fig. 15 Travel times δτ^{(p)} as defined in Eq. (89), computed from the difference between crosscovariances C^{(p)}(Θ,ω) computed using polynomials of degrees and a reference crosscovariance C_{ref}(Θ,ω). Computations with and without flow are shown by the dashed and solid lines, respectively. The vertical line labeled “p = 10” indicates the number of degrees of freedom used in all calculations in this paper. 

Open with DEXTER  
In the text 
Fig. 16 Sensitivity of the mean travel time, [τ(r_{1},r_{2}) + τ(r_{2},r_{1})] /2, to relative soundspeed perturbations δc/c in the interior, with observation point r_{1} located on the polar axis and observation point r_{2} at latitude of 45°. Panel a): cut at r = 0.95 R_{⊙} through the soundspeed kernel. Panel b): slice through a plane containing the two observation points and the center of the Sun. The green line indicates the position of the ray path. Panel c): slice in a plane perpendicular to the ray path connecting the two observation points (green square) at equal distance from the two observation points. In both slices b) and c), the black arc of a circle locates radius r = 0.95 R_{⊙} which corresponds to panel a). 

Open with DEXTER  
In the text 
Fig. 17 Left and center panels: kernels ⟨ K_{ur} ⟩ and ⟨ K_{uθ} ⟩ for the r and θ components of the flow. Point r_{1} is at the north pole (photosphere) and point r_{2} is at 45° latitude. The values of the kernels are scaled by the sound speed c and are saturated at 1/400 of the maximum value. The ray path connecting the two points is shown (thick black line) as well as the computational boundary (black half circle). Right panel: comparison between the travel times computed directly from the crosscovariance function (curves, see Eq. (25)) and those computed from the sensitivity kernels (“+” symbols, see Eq. (90)). The blue curve is for the travel times measured from the pole to latitude 45°, the red curve for the travel times measured in the opposite direction. The accuracy of the travel times is approximately 10^{3} s. 

Open with DEXTER  
In the text 
Fig. 18 Left and center panels: kernels ⟨ K_{ur} ⟩ and ⟨ K_{uθ} ⟩ for the r and θ components of the flow, using all azimuthal components  m  ≤ m_{max} = 35. The separation distance between points r_{1} and r_{2} is 42°; the center point is located at a latitude of 40°. The ray path connecting the two points is shown in black. The right panel shows the convergence of the travel times as a function of m_{max}, for the northsouth (blue) and southnorth (green) travel directions. The red curve shows their difference. 

Open with DEXTER  
In the text 
Fig. 19 Left: unfiltered ⟨ K_{uθ} ⟩ kernel with r_{1} at the pole and r_{2} at a colatitude of 15.36°. Right: kernel ⟨ K_{uθ} ⟩ for filtered observations, where the Gaussian phasespeed filter F_{ℓ}(ω) is centered at 125.2 km s^{1} with a dispersion of 12.3 km s^{1}, for ℓ up to 1000. The ray path is shown by the black line. 

Open with DEXTER  
In the text 
Fig. B.1 u_{r}(r,θ) (top) and u_{θ}(r,θ) (bottom) at θ = 45° of the meridional flow model discussed in Appendix B for a maximum flow velocity U = 20 m s^{1} at the surface. 

Open with DEXTER  
In the text 