# Computational helioseismology in the frequency domain: acoustic waves in axisymmetric solar models with flows

Laurent Gizon, Hélène Barucq, Marc Duruflé, Chris Hanson, Michael

Leguèbe, Aaron Birch, Juliette Chabassier, Damien Fournier, Thorsten

Hohage, Emanuele Papini

► To cite this version:

Laurent Gizon, Hélène Barucq, Marc Duruflé, Chris Hanson, Michael Leguèbe, et al. Computational helioseismology in the frequency domain: acoustic waves in axisymmetric solar models with flows. Astronomy and Astrophysics - A&A, 2017, 600, pp.A35. 10.1051/0004-6361/201629470. hal-01403332

HAL Id: hal-01403332

https://hal.science/hal-01403332

Submitted on 2 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Computational helioseismology in the frequency domain: acoustic waves in axisymmetric solar models with flows

Laurent Gizon<fnref n="1" />,<fnref n="2" />,*, Hélène Barucq<fnref n="3" />,**, Marc Duruflé<fnref n="4" />, Chris S. Hanson<fnref n="1" />, Michael Leguèbe<fnref n="1" />, Aaron C. Birch<fnref n="1" />, Juliette Chabassier<fnref n="3" />, Damien Fournier<fnref n="5" />, Thorsten Hohage<fnref n="5" />, Emanuele Papini<fnref n="1" />

Received 15 June 2016; accepted XXX

## ABSTRACT

**Context.** Local helioseismology has so far relied on semi-analytical 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 time-distance helioseismology in the frequency domain. The fundamental quantity to be computed is the cross-covariance of the seismic wavefield.

**Methods.** We choose sources of wave excitation that enable us to relate the cross-covariance 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 high-order finite-element 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 travel-time 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 travel-time measurements.

**Conclusions.** The method presented here is computationally efficient and will be used to interpret travel-time measurements in order to infer, e.g., the large-scale meridional flow in the solar convection zone. It allows the implementation of (full-waveform) iterative inversions, whereby the axisymmetric background model is updated at each iteration.

**Key words.** helioseismology – solar physics – numerical methods

## 1. Introduction

Time-distance 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 two-point cross-covariance function of the random solar oscillations. The cross-covariance function tells us about the travel time of wave packets between any two locations, in either direction. Flows break the time-symmetry of the cross-covariance 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 travel-time 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 cross-covariance function caused by weak heterogeneities. This framework has proven to be useful for local helioseismology (Birch et al. 2004; Birch & Gizon 2007; Birch et al. 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 PDE-constrained 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

<footnote n="1">¹ Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany</footnote>

<footnote n="2">² Institut für Astrophysik, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany</footnote>

<footnote n="3">³ Magique-3D, Inria Bordeaux Sud-Ouest, Université de Pau et des Pays de l'Adour, 64013 Pau, France</footnote>

<footnote n="4">⁴ Magique-3D, Inria Bordeaux Sud-Ouest, Université de Bordeaux, 33400 Talence, France</footnote>

<footnote n="5">⁵ Institut für Numerische und Angewandte Mathematik, Georg-August-Universität Göttingen, Lotzestraße 16, 37083 Göttingen, Germany</footnote>

<footnote>* Corresponding author: gizon@mps.mpg.de</footnote>

<footnote>** Corresponding author: helene.barucq@inria.fr</footnote>
framework of Gizon & Birch (2002). Of particular relevance to the present work is the proposal by Nissen-Meyer 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 & Nissen-Meyer 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 cross-covariance 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:

C1. Oscillation power spectrum. The purpose of helioseismology is to infer a model for the solar interior (density, sound speed, flows, etc., as functions of position) that is consistent with the seismic data. In global mode helioseismology, the emphasis is on the comparison between the model and observed mode frequencies. In contrast, local helioseismology requires models for the mode amplitudes and line profiles in addition to the mode frequencies. Thus a solar model must also describe wave excitation and damping. Most reference solar models are 1D and their quality can be assessed by comparing the model and observed oscillation power spectra. In linear inversions for local helioseismology, the reference oscillation power spectrum has to be very solar-like.

C2. Flexibility. The semi-analytical 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 cross-covariance 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 point-spread function) and geometrical effects (e.g., foreshortening, line-of-sight 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.

C3. Representation of the Green's function. Often Green's functions are computed using a normal-mode 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).

C4. Time-domain 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).

C5. Computational challenge and inverse problem. Modern helioseismic observations consist of large data sets (long time series of 16 megapixel images). Cross-correlations 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:

1. Rewriting the perturbation to the cross-covariance function in terms of the Green's function, $G$, and the expectation value of the cross-covariance function, $\bar{C}$, in the reference model. In this formalism $G$ and $\bar{C}$ are the fundamental quantities that enter any forward calculation in local helioseismology. The problem then becomes deterministic. Additionally, many systematic effects can be accounted for by treating them as numerical operations on $G$ and $\bar{C}$.

2. Computing $G$ and $\bar{C}$ in the frequency-domain. 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 multi-core environment. Furthermore, the computation can easily be restricted to the frequency range of interest, at the necessary frequency resolution.

3. Adopting a flexible geometrical setup using a finite-element 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.

4. 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.

5. Demonstrating that there is a choice of source covariance such that the expectation value of the cross-covariance 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, time-distance diagrams, and travel-time sensitivity kernels. An application including a large-scale 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 $\xi(r, t)$ of small amplitude waves obeys an equation of the form

$$ \left( \frac{\partial}{\partial t} + \gamma \otimes +u \cdot \nabla \right)^2 \xi + \mathcal{H}[\xi] = f, \qquad (1) $$
where $\mathcal{H}$ is a spatial operator and $\mathbf{u}(\mathbf{r})$ is a background flow. Wave attenuation is accounted for by temporal convolution (⊗) with the function $\gamma(\mathbf{r}, t)$. The function $f(\mathbf{r}, t)$ represents forcing by turbulent convection and is thus a realization drawn from a random process. The spatial differential operator $\mathcal{H}$ is Hermitian for appropriate choices of the boundary condition and inner product (see Lynden-Bell & Ostriker 1967). An approximation of $\mathcal{H}$ that captures the physics of acoustic oscillations is:

$$ \mathcal{H}[\xi] = -\frac{1}{\rho}\nabla(\rho c^2 \nabla \cdot \xi) + \text{terms involving gravity}, \quad (2) $$

where $\rho$ 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

$$ -(\omega + i\hat{\gamma} + i\mathbf{u} \cdot \nabla)^2 \hat{\xi} + \mathcal{H}[\hat{\xi}] = \hat{f}, \quad (3) $$

using the Fourier convention

$$ \hat{\xi}(\mathbf{r}, \omega) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \xi(\mathbf{r}, t) e^{i\omega t} dt. \quad (4) $$

If, in addition, the source of excitation is statistically stationary, then the frequency components $\hat{f}(\mathbf{r}, \omega)$ are independent random variables and the frequency components of the wave displacement, $\hat{\xi}(\mathbf{r}, \omega)$, are also statistically independent of each other. Thus, in frequency space, the problem separates into a set of independent boundary-value 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 $\psi = c\nabla \cdot \xi$. Neglecting the gravity terms, the second-order terms in $\gamma$ and $\mathbf{u}$, and the cross term in $\gamma\mathbf{u}$ in Eq. (3), and then taking the divergence, we obtain

$$ \begin{aligned} & - (\omega^2 + 2i\omega\gamma)\psi - 2i\omega c \nabla \cdot (\mathbf{u} \cdot \nabla \xi) - 2i\omega c \xi \cdot \nabla \gamma \\ & - c\nabla \cdot \left(\frac{1}{\rho}\nabla(\rho c\psi)\right) = c\nabla \cdot \mathbf{f}. \end{aligned} \quad (5) $$

For slow variations of $\mathbf{u}$, $c$ and $\gamma$ compared to the wavelength, the wave equation for $\psi$ simplifies to

$$ L[\psi] := -\sigma^2\psi - 2i\omega\mathbf{u} \cdot \nabla\psi + H[\psi] = s \quad (6) $$

with

$$ \begin{aligned} \sigma^2 &:= \omega^2 + 2i\omega\gamma \\ H[\psi] &:= -c\nabla \cdot \left(\frac{1}{\rho}\nabla(\rho c\psi)\right). \end{aligned} \quad (7) $$

The random source of excitation, $s(\mathbf{r}, \omega) = c(\mathbf{r})\nabla \cdot f(\mathbf{r}, \omega)$, and the attenuation, $\gamma(\mathbf{r}, \omega) > 0$, depend on frequency.

Another possible choice in order to obtain a scalar equation is to use the ansatz $\xi = \rho^{-1}\nabla(\rho c\psi)$. Under the assumption that the forcing term $\rho f$ is curl free, this leads to the same scalar equation as above.

When $\rho$ and $c$ are solar-like, the scalar equation (6) captures most of the interesting physics of $\rho$ modes. For example, a solar-like density gives the correct acoustic cutoff frequency of

<figure><img src="image_2.png" /><figcaption>Fig. 1: Geometrical setup for a background medium symmetric about an axis $\hat{z}$. The thick contour delineates the boundary, $S$, of the computational domain $V$. A position vector $\mathbf{r}$ in $V$ is specified by spherical-polar coordinates $(r, \theta, \phi)$. The seismic observable $\psi(\mathbf{r})$ is measured near the photosphere.</figcaption></figure>

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

$$ M(\mathbf{r}, \mathbf{r}', \omega) = \mathbb{E}[s^*(\mathbf{r}, \omega)s(\mathbf{r}', \omega)] = \epsilon(\mathbf{r})\Pi(\omega)\delta(\mathbf{r} - \mathbf{r}'), \quad (8) $$

where $*$ takes the complex conjugate and the functions $\epsilon(\mathbf{r})$ and $\Pi(\omega)$ 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 spherical-polar coordinate system $(r, \theta, \phi)$. The solar (or stellar) model is cylindrically symmetric around the $\hat{z}$-axis. The computational boundary $S$ encloses the solar photosphere. The wave equation is supplemented by a boundary condition at $S$ on $\psi$. We denote by $\hat{n}$ the unit normal vector to $S$ and by $\partial_n\psi$ the normal derivative of $\psi$.

A free-surface boundary condition ($\psi \propto \nabla \cdot \xi = 0$, Dirichlet) is often used in helioseismology (Bogdan et al. 1996; Christensen-Dalsgaard 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 high-frequency 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,

$$ \partial_n \psi = ik_n \psi \quad \text{on } S, \quad (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 equation (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 non-attenuating waves in the Sun is the Hermiticity of the wave operator, which implies that the eigenvalues are real. Lynden-Bell & Ostriker (1967) proved Hermiticity of the wave operator for a self-gravitating fluid bounded by a free surface.

Here we ask if this condition is satisfied for the time-harmonic scalar wave equation defined by Eq. (6). The scalar wave field belongs to a complex Hilbert space. For any two functions $\psi$ and $\varphi$ in this space we define the inner product

$$ \langle \psi, \varphi \rangle = \int_V \psi^* \varphi \rho \, dr, \qquad (10) $$

where $V$ is the computational domain. Performing integration by parts twice implies

$$ \langle \psi, H[\varphi] \rangle = \langle H[\psi], \varphi \rangle - \int_S [\psi^* \partial_n(\rho c \varphi) - \varphi \partial_n(\rho c \psi^*)] c \, dS. \quad (11) $$

We see that the spatial operator $H$ is Hermitian if $\psi$ and $\varphi$ 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 equation (6) are the damping and advection terms. Using integration by parts once, we ascertain that the damping operator is anti-Hermitian and the advection operator is Hermitian:

$$ \langle \psi, 2i\omega\gamma\varphi \rangle = -\langle 2i\omega\gamma\psi, \varphi \rangle, \quad (12) $$

$$ \langle \psi, 2i\omega u \cdot \nabla \varphi \rangle = \langle 2i\omega u \cdot \nabla \psi, \varphi \rangle. \quad (13) $$

The last result is a consequence of mass conservation, $\nabla \cdot (\rho u) = 0$, for a steady flow that does not cross the computational boundary ($u \cdot \hat{n} = 0$ on $S$).

Combining Eqs. (6), (11), (12), and (13) yields

$$ \langle \psi, L[\varphi] \rangle = \langle L^{\dagger}[\psi], \varphi \rangle - \int_S [\psi^* \partial_n(\rho c \varphi) - \varphi \partial_n(\rho c \psi^*)] c \, dS, \quad (14) $$

where

$$ L^{\dagger} = L_{(-\gamma)} = L_{(-u)}^{*} \quad (15) $$

is the operator obtained by switching the sign of $\gamma$ in $L$ or, equivalently, by switching the sign of $u$ and taking the complex conjugate. We note that for a radiation boundary condition, $L^\dagger$ is not the adjoint of $L$ because the surface integral remains.

Article number, page 4 of 23page.23

## 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

$$ L[G(\mathbf{r}, \mathbf{r}', \omega)]=\frac{1}{\rho(\mathbf{r})} \delta(\mathbf{r}-\mathbf{r}'), \quad (16) $$

where $\delta$ is the Dirac delta function and $G$ satisfies the same boundary condition as the wave field $\psi$. Consider a source at $\mathbf{r}_1$ in the physical domain and a source at $\mathbf{r}_2$ in the domain in which the flow has opposite sign:

$$ \begin{cases} \rho L[G(\mathbf{r}, \mathbf{r}_1, \omega)] = \delta(\mathbf{r} - \mathbf{r}_1) \\ \rho L_{(-\omega)}[G(\mathbf{r}, \mathbf{r}_2, \omega; -\mathbf{u})] = \delta(\mathbf{r} - \mathbf{r}_2), \end{cases} \quad (17) $$

where $G(\cdot; -\mathbf{u})$ is the Green's function associated with the operator $L_{(-\omega)}$, with the same boundary condition as $G$. Multiplying the first equation by $G(\mathbf{r}, \mathbf{r}_2, \omega; -\mathbf{u})$ and the second by $G(\mathbf{r}, \mathbf{r}_1, \omega)$, integrating each over position $\mathbf{r}$, and subtracting the two equations, we obtain

$$ \begin{aligned} G(\mathbf{r}_1, \mathbf{r}_2, \omega; -\mathbf{u}) - G(\mathbf{r}_2, \mathbf{r}_1, \omega) = & \langle G^*(\mathbf{r}, \mathbf{r}_2, \omega; -\mathbf{u}), L[G(\mathbf{r}, \mathbf{r}_1, \omega)] \rangle \\ & - \langle L^{\dagger}[G^*(\mathbf{r}, \mathbf{r}_2, \omega; -\mathbf{u})], G(\mathbf{r}, \mathbf{r}_1, \omega) \rangle. \end{aligned} \quad (18) $$

Following the same logic as earlier (integration by parts),

$$ \langle \psi^*, L[\varphi] \rangle = \langle L^{\dagger}[\psi^*], \varphi \rangle - \int_S [\psi \partial_n(\rho c \varphi) - \varphi \partial_n(\rho c \psi)] c dS \quad (19) $$

for any two functions $\psi$ and $\varphi$ in the space of solutions. This time the surface integral vanishes for either type of boundary condition (Dirichlet or radiation). In particular,

$$ \langle L^{\dagger}[G^*], \varphi \rangle = \langle G^*, L[\varphi] \rangle. \quad (20) $$

This relationship shows that the right-hand side of Eq. (18) vanishes and so we find

$$ G(\mathbf{r}_2, \mathbf{r}_1, \omega) = G(\mathbf{r}_1, \mathbf{r}_2, \omega; -\mathbf{u}). \quad (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. Cross-covariance function and travel times

For the sake of simplicity, let us suppose that $\psi$ can be directly measured on the solar surface. This choice is not unreasonable as $\psi$ 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 $\psi$, for example a proxy for the line-of-sight velocity.

At frequency $\omega$, consider the cross-covariance in Fourier space as the product of the wave field at two locations of measurement,

$$ C(\mathbf{r}_1, \mathbf{r}_2, \omega) = \psi^*(\mathbf{r}_1, \omega) \psi(\mathbf{r}_2, \omega). \quad (22) $$

In terms of the Green's function, we have

$$ \psi(\mathbf{r}_j, \omega) = \int_V G(\mathbf{r}_j, \mathbf{r}, \omega) s(\mathbf{r}, \omega) \rho \, d\mathbf{r}, \quad (23) $$
where the source $s(\mathbf{r}, \omega)$ is a realization of a random process. Under the assumption of spatially uncorrelated sources (see Eq. (8)), the cross-covariance can be written as one volume integral:

$$ 
\begin{aligned}
\bar{C}(\mathbf{r}_1, \mathbf{r}_2, \omega) &= \mathbb{E}[C(\mathbf{r}_1, \mathbf{r}_2, \omega)] \\
&= \int_V \rho \mathrm{d}\mathbf{r} \int_V \rho' \mathrm{d}\mathbf{r}' G^*(\mathbf{r}_1, \mathbf{r}, \omega) G(\mathbf{r}_2, \mathbf{r}', \omega) M(\mathbf{r}, \mathbf{r}', \omega) \\
&= \Pi(\omega) \int_V G^*(\mathbf{r}_1, \mathbf{r}, \omega) G(\mathbf{r}_2, \mathbf{r}, \omega) \epsilon(\mathbf{r}) \rho^2 \mathrm{d}\mathbf{r}.
\end{aligned} 
\quad (24) $$

Following Gizon & Birch (2002), we define the perturbation to the travel time $\delta\tau$ between points $\mathbf{r}_1$ and $\mathbf{r}_2$ as

$$ \delta\tau(\mathbf{r}_1, \mathbf{r}_2) = \int_{-\infty}^{\infty} W^*(\omega) [C(\mathbf{r}_1, \mathbf{r}_2, \omega) - C_{\text{ref}}(\mathbf{r}_1, \mathbf{r}_2, \omega)] \mathrm{d}\omega, \quad (25) $$

where

$$ W(\omega) = - \frac{\int_{-\infty}^{\infty} w(t) [\partial_t C_{\text{ref}}(\mathbf{r}_1, \mathbf{r}_2, t)] e^{i\omega t} \mathrm{d}t}{\int_{-\infty}^{\infty} w(t') [\partial_t' C_{\text{ref}}(\mathbf{r}_1, \mathbf{r}_2, t')]^2 \mathrm{d}t'}, \quad (26) $$

$C_{\text{ref}}$ is a reference cross-covariance 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\sigma^2]$ where $t_g > 0$ is a group time and $\sigma$ 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 travel-time sensitivity kernels in terms of the Green's function and the cross-covariance 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(\mathbf{r})$, $\rho(\mathbf{r})$, and $u(\mathbf{r})$, as well as spatial changes to the attenuation $\gamma(\mathbf{r})$ and the source amplitude $\epsilon(\mathbf{r})$. The vector flow

$$ \mathbf{u}(\mathbf{r}) = \sum_k u_k(\mathbf{r}) \hat{\mathbf{e}}_k(\mathbf{r}) = u_r(\mathbf{r}) \hat{\mathbf{r}} + u_\theta(\mathbf{r}) \hat{\mathbf{\theta}} + u_\phi(\mathbf{r}) \hat{\mathbf{\phi}}, \quad (27) $$

is specified by the three components $u_k$ on the basis $\{\hat{\mathbf{e}}_k\}$ of unit vectors $\hat{\mathbf{r}}$, $\hat{\mathbf{\theta}}$ and $\hat{\mathbf{\phi}}$ for the spherical-polar coordinate system. In short, the physical variables of interest can be combined into a set

$$ \{q_{\alpha}\} = \{c, \rho, u_r, u_\theta, u_\phi, \gamma, \epsilon\}. \quad (28) $$

We are looking for travel-time sensitivity kernels $K_\alpha$ such that the travel-time perturbation $\delta\tau$ can be written as

$$ \mathbb{E}[\delta\tau] = \sum_{\alpha} \int_V \delta q_{\alpha}(\mathbf{r}) K_{\alpha}(\mathbf{r}) \mathrm{d}\mathbf{r} + \text{surface terms}, \quad (29) $$

for infinitesimal perturbations $\delta q_\alpha$. For some variables there is a surface integral on the computational boundary, because integration by parts is needed to write travel-time 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 travel-time perturbation

$$ \delta\tau(\mathbf{r}_1, \mathbf{r}_2) = \int_{-\infty}^{\infty} W^*(\omega) \delta C(\mathbf{r}_1, \mathbf{r}_2, \omega) \mathrm{d}\omega \quad (30) $$

is written in terms of the first-order perturbation to the cross-covariance:

$$ \delta C(\mathbf{r}_1, \mathbf{r}_2, \omega) = \psi^*(\mathbf{r}_1, \omega) \delta\psi(\mathbf{r}_2, \omega) + \delta\psi^*(\mathbf{r}_1, \omega) \psi(\mathbf{r}_2, \omega). \quad (31) $$

The next step is then to express $\delta\psi$ in terms of the $\delta q_\alpha$.

## 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,

$$ 
\begin{aligned}
L[\delta\psi] &= -\delta L[\psi] + \delta s \\
&= 2i\omega \delta\gamma\psi + 2i\omega \delta\mathbf{u} \cdot \nabla\psi - \delta H[\psi] + \delta s,
\end{aligned} 
\quad (32) $$

where $\delta L$ is the perturbation to the wave operator caused by perturbations to the medium and

$$ \delta H[\psi] = \frac{\delta c}{c} H[\psi] + H\left[\frac{\delta c}{c}\psi\right] - c\nabla \cdot \left(\frac{\delta \rho}{\rho^2} \nabla(\rho c \psi)\right) + H\left[\frac{\delta \rho}{\rho}\psi\right]. \quad (33) $$

A formal solution to Eq. (32) is obtained in terms of the Green's function:

$$ \delta\psi(\mathbf{r}_j, \omega) = - \int_V G(\mathbf{r}_j, \mathbf{r}, \omega) \delta L[\psi(\mathbf{r}, \omega)] \rho \mathrm{d}\mathbf{r} \\ + \int_V G(\mathbf{r}_j, \mathbf{r}, \omega) \delta s(\mathbf{r}, \omega) \rho \mathrm{d}\mathbf{r}. \quad (34) $$

## 4.3. Perturbation to the cross-covariance

With Eq. (34) in hand, the expectation of the perturbation to the cross-covariance is then determined from Eq. (31):

$$ 
\begin{aligned}
\delta \bar{C}(\mathbf{r}_1, \mathbf{r}_2, \omega) = & - \int_V G(\mathbf{r}_2, \mathbf{r}, \omega) \delta L[\bar{C}(\mathbf{r}_1, \mathbf{r}, \omega)] \rho \mathrm{d}\mathbf{r} \\
& - \int_V G^*(\mathbf{r}_1, \mathbf{r}, \omega) \delta L^*[\bar{C}^*(\mathbf{r}_2, \mathbf{r}, \omega)] \rho \mathrm{d}\mathbf{r} \\
& + \Pi(\omega) \int_V G^*(\mathbf{r}_1, \mathbf{r}, \omega) G(\mathbf{r}_2, \mathbf{r}, \omega) \delta\epsilon(\mathbf{r}) \rho^2 \mathrm{d}\mathbf{r},
\end{aligned} 
\quad (35) $$

where we write the perturbation to the source covariance as a spatially varying change in amplitude,

$$ \mathbb{E}[s^*(\mathbf{r}, \omega)\delta s(\mathbf{r}', \omega) + \delta s^*(\mathbf{r}, \omega)s(\mathbf{r}', \omega)] = \delta\epsilon(\mathbf{r})\Pi(\omega)\delta(\mathbf{r}-\mathbf{r}'). \quad (36) $$

In order to obtain kernels for individual perturbations, we introduce the bilinear operators $\mathcal{L}_\alpha$ such that for any functions $g$ and $h$ we have

$$ \int_V g \, \delta L[h] \rho \mathrm{d}\mathbf{r} = \sum_{\alpha} \int_V \delta q_{\alpha}(\mathbf{r}) \mathcal{L}_{\alpha}[g, h] \mathrm{d}\mathbf{r} + \text{surface terms.} \quad (37) $$
These bilinear operators are obtained by integration by parts and are explicitly given by

$$ \mathcal{L}_c[g, h] = (g H[h] + h H[g]) \rho/c, \quad (38) $$

$$ \mathcal{L}_{\rho}[g, h] = -\nabla(\rho c g) \cdot \nabla(\rho c h) / \rho^2 + h H[g], \quad (39) $$

$$ \mathcal{L}_{u_k}[g, h] = -2i\rho\omega g (\partial_k h), \quad (40) $$

$$ \mathcal{L}_{\gamma}[g, h] = -2i\rho\omega g h. \quad (41) $$

In Eq. (40), $\partial_k = \hat{\epsilon}_k(\mathbf{r}) \cdot \nabla_{\mathbf{r}}$ denotes the component of the spatial gradient in the $k$ direction. Combining Eqs. (37) and (35), we obtain an explicit linear relationship between the $\delta q_{\alpha}$ and $\delta \bar{C}$:

$$ \begin{aligned} \delta \bar{C}(\mathbf{r}_1, \mathbf{r}_2, \omega) = & -\sum_{\alpha} \int_V \delta q_{\alpha}(\mathbf{r}) \mathcal{L}_{\alpha} [G(\mathbf{r}_2, \mathbf{r}, \omega), \bar{C}(\mathbf{r}_1, \mathbf{r}, \omega)] d\mathbf{r} \\ & -\sum_{\alpha} \int_V \delta q_{\alpha}(\mathbf{r}) \mathcal{L}_{\alpha}^{*} [G^{*}(\mathbf{r}_1, \mathbf{r}, \omega), \bar{C}^{*}(\mathbf{r}_2, \mathbf{r}, \omega)] d\mathbf{r} \\ & + \Pi(\omega) \int_V \delta\epsilon(\mathbf{r}) G^{*}(\mathbf{r}_1, \mathbf{r}, \omega) G(\mathbf{r}_2, \mathbf{r}, \omega) \rho^2 d\mathbf{r} \\ & + \text{surface terms.} \end{aligned} \quad (42) $$

Using this expression together with Eqs. (29) and (30), the travel-time sensitivity kernels are

$$ \begin{aligned} K_{\alpha}(\mathbf{r}; \mathbf{r}_1, \mathbf{r}_2) = & - \int_{-\infty}^{\infty} W^{*}(\omega) \mathcal{L}_{\alpha} [G(\mathbf{r}_2, \mathbf{r}, \omega), \bar{C}(\mathbf{r}_1, \mathbf{r}, \omega)] d\omega \\ & - \int_{-\infty}^{\infty} W^{*}(\omega) \mathcal{L}_{\alpha}^{*} [G^{*}(\mathbf{r}_1, \mathbf{r}, \omega), \bar{C}^{*}(\mathbf{r}_2, \mathbf{r}, \omega)] d\omega \end{aligned} \quad (43) $$

for the scatterers $\delta q_{\alpha} \in \{\delta c, \delta \rho, \delta u_k, \delta \gamma\}$ and

$$ K_{\epsilon}(\mathbf{r}; \mathbf{r}_1, \mathbf{r}_2) = \rho^2(\mathbf{r}) \int_{-\infty}^{\infty} \Pi(\omega) W^*(\omega) G^*(\mathbf{r}_1, \mathbf{r}, \omega) G(\mathbf{r}_2, \mathbf{r}, \omega) d\omega \quad (44) $$

for the source amplitude $\delta \epsilon$.

# 5. Convenient source of excitation

## 5.1. How can we simplify the computation of $\bar{C}$?

The kernels from the previous section are well defined once the background medium, attenuation and source functions have been specified. At each $\omega$, the attenuation function $\gamma(\mathbf{r}, \omega)$ can be tuned to yield the observed line widths of solar oscillations in the power spectrum. Similarly, at each $\omega$, the source function $\epsilon(\mathbf{r})\Pi(\omega)$ can be tuned to give the observed mode amplitudes (Birch et al. 2004; Birch & Gizon 2007). However, the integral relationship between a general $\epsilon(\mathbf{r})$ and $\bar{C}$ 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 cross-covariance can be written in terms of the Green's function. As a second step we check whether the resulting power spectrum is solar-like.

It is known in geophysics and acoustics that, under appropriate conditions on the source covariance, the expectation value of the cross-covariance 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 $\mathbf{r}_1$ and $\mathbf{r}_2$, taking the complex conjugate of the first:

$$ \left\{ \begin{aligned} \delta(\mathbf{r} - \mathbf{r}_1) &= \rho(\mathbf{r}) L^*[G^*(\mathbf{r}, \mathbf{r}_1, \omega)], \\ \delta(\mathbf{r} - \mathbf{r}_2) &= \rho(\mathbf{r}) L[G(\mathbf{r}, \mathbf{r}_2, \omega)]. \end{aligned} \right. \quad (45) $$

Multiply the first equation by $G(\mathbf{r}, \mathbf{r}_2, \omega)$ and the second equation by $G^*(\mathbf{r}, \mathbf{r}_1, \omega)$, integrate each over $\mathbf{r}$, and then subtract to find

$$ \begin{aligned} G(\mathbf{r}_1, \mathbf{r}_2, \omega) - G^*(\mathbf{r}_2, \mathbf{r}_1, \omega) = & \langle L[G(\mathbf{r}, \mathbf{r}_1, \omega)], G(\mathbf{r}, \mathbf{r}_2, \omega) \rangle - \langle G(\mathbf{r}, \mathbf{r}_1, \omega), L[G(\mathbf{r}, \mathbf{r}_2, \omega)] \rangle \\ = & \langle L[G(\mathbf{r}, \mathbf{r}_1, \omega)], G(\mathbf{r}, \mathbf{r}_2, \omega) \rangle - \langle L^\dagger[G(\mathbf{r}, \mathbf{r}_1, \omega)], G(\mathbf{r}, \mathbf{r}_2, \omega) \rangle \\ & + \text{surface term} \\ = & 4i\omega \langle \gamma(\mathbf{r}, \omega) G(\mathbf{r}, \mathbf{r}_1, \omega), G(\mathbf{r}, \mathbf{r}_2, \omega) \rangle + \text{surface term} \\ = & 4i\omega \langle \gamma(\mathbf{r}, \omega) G(\mathbf{r}_1, \mathbf{r}, \omega; -\mathbf{u}), G(\mathbf{r}_2, \mathbf{r}, \omega; -\mathbf{u}) \rangle + \text{surface term}. \end{aligned} \quad (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 $\mathbf{u}$ in the above equation, we obtain

$$ \begin{aligned} G(\mathbf{r}_2, \mathbf{r}_1, \omega) - G^*(\mathbf{r}_2, \mathbf{r}_1, \omega; -\mathbf{u}) = & \\ = & 4i\omega \int_V \gamma(\mathbf{r}, \omega) G^*(\mathbf{r}_1, \mathbf{r}, \omega) G(\mathbf{r}_2, \mathbf{r}, \omega) \rho dV \\ & + 2i\omega \int_S c(\mathbf{r}) G^*(\mathbf{r}_1, \mathbf{r}, \omega) G(\mathbf{r}_2, \mathbf{r}, \omega) \rho dS, \end{aligned} \quad (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

$$ \epsilon := \frac{\gamma(\mathbf{r}, \omega)}{\rho(\mathbf{r})} + \frac{c}{2\rho} \delta(r - R(\theta)) \quad (48) $$

implies

$$ \mathbb{E}[C(\mathbf{r}_1, \mathbf{r}_2, \omega)] = \frac{\Pi(\omega)}{4i\omega} [G(\mathbf{r}_2, \mathbf{r}_1, \omega) - G^*(\mathbf{r}_2, \mathbf{r}_1, \omega; -\mathbf{u})]. \quad (49) $$

Thus the expectation value of the cross-correlation can be written as a sum of causal and anti-causal Green's functions, when waves are appropriately excited throughout the volume and on the computational boundary $S$ at radius $r = R(\theta)$. 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 $\Pi(\omega)$. 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 cross-covariance are fully specified by perturbations to $c$, $\rho$, $\gamma$, and $\mathbf{u}$.
## 5.2. Sensitivity kernels in terms of G only

As the cross-covariance depends only on the Green's function via Eq. (49), it is also the case for the kernels. Denoting for simplicity

$$ 
\begin{aligned}
G_i(\mathbf{r}) &:= G(\mathbf{r}, \mathbf{r}_i, \omega), \\
G_i^\dagger(\mathbf{r}) &:= G^*(\mathbf{r}, \mathbf{r}_i, \omega; -\mathbf{u}),
\end{aligned} 
\quad (50) $$

the kernels defined by Eq. (43) can be computed from four Green's functions

$$ G_1(\mathbf{r}), G_2(\mathbf{r}) \text{ and } G_1^\dagger(\mathbf{r}), G_2^\dagger(\mathbf{r}). \quad (51) $$

For example, the travel-time perturbations induced by a background flow perturbation,

$$ \delta\tau(\mathbf{r}_1, \mathbf{r}_2) = \sum_k \int_V \delta u_k(\mathbf{r}) K_{u_k}(\mathbf{r}; \mathbf{r}_1, \mathbf{r}_2) \, d\mathbf{r}, \quad (52) $$

can be computed through the kernels

$$ K_{u_k} = \frac{\rho}{2} \int_{-\infty}^{\infty} \Pi W^* \left[ G_2^{\dagger *} \partial_k (G_1 - G_1^\dagger) + G_1^\dagger \partial_k (G_2^* - G_2^{\dagger *}) \right] d\omega. \quad (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 = 0$), seismic reciprocity implies that $G^\dagger = G^*$ and the cross-covariance takes the simple form,

$$ \bar{C}_i(\mathbf{r}) := \bar{C}(\mathbf{r}_i, \mathbf{r}, \omega) = \frac{\Pi(\omega)}{2\omega} \operatorname{Im} G_i(\mathbf{r}). \quad (54) $$

This simple form for the cross-covariance simplifies the scattering bilinear operators:

$$ \mathcal{L}_c[G_i, \bar{C}_j] = \frac{\rho \Pi}{2 \omega c} [G_i \operatorname{Im}(H[G_j]) + H[G_i] \operatorname{Im}G_j], \quad (55) $$

$$ \mathcal{L}_{\rho}[G_i, \bar{C}_j] = -\frac{\Pi}{2\omega} \left( \frac{1}{\rho^2} \nabla(\rho c G_i) \cdot \nabla(\rho c \operatorname{Im}G_j) - H[G_i] \operatorname{Im}G_j \right), \quad (56) $$

$$ \mathcal{L}_{u_k}[G_i, \bar{C}_j] = -i\rho\Pi G_i \partial_k(\operatorname{Im}G_j), \quad (57) $$

$$ \mathcal{L}_{\gamma}[G_i, \bar{C}_j] = -i\rho\Pi G_i \operatorname{Im}G_j. \quad (58) $$

Thus, to compute a kernel $K_\alpha$ when there is no background flow, all we need is two Green's function computations, one with a source at $\mathbf{r}_1$ and the other with a source at $\mathbf{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 $\omega$, 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 large-scale 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<fnref n="1" />, developed at Inria for various wave propagation problems. Montjoie is versatile, well-tested, 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 section 8.3.

## 6.2. Geometrical setup for the 2.5D problem

For an axisymmetric background model, the computational domain $\Sigma$ is a meridional generating section of the geometry (Fig. 2), i.e. the half-disk of radius $R$ in the case of the sphere.

The background model is symmetric about the $\hat{z}$-axis. In 3D space a spatial point is associated with the position vector

$$ \mathbf{r} = (r, \theta, \phi) := (\tilde{r}, \phi), \quad (59) $$

where $\tilde{r}$ belongs to the 2D section $\Sigma$ and $\phi$ is longitude about the axis of symmetry (Fig. 2). The background density, sound speed, and flow are specified at positions $\tilde{r}$ in $\Sigma$:

$$ \rho = \rho(\tilde{r}), \quad (60) $$

$$ c = c(\tilde{r}), \quad (61) $$

$$ \mathbf{u} = \tilde{u}(\tilde{r}) + u_{\phi}(\tilde{r})\hat{\phi}, \quad (62) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2: Two-dimensional generating section $\Sigma$ of the volume $V$, delineated by $\partial\Sigma$ (thick contour). A point $\tilde{r}$ in $\Sigma$ is specified by coordinates $(r, \theta)$, where $r$ is the radius and $\theta$ is the co-latitude, or, equivalently, by coordinates $(\varpi, z)$, where $\varpi = r \sin\theta$ is the distance to the axis.</figcaption></figure>

<footnote n="1">¹ http://montjoie.gforge.inria.fr/</footnote>
where $\tilde{u}$ is the meridional component of the flow and $u_\phi$ the ro-
tational velocity component.

6.3. Expansion of solution in longitudinal Fourier modes

As the unknown of the finite element method we choose the
function $\tilde{\psi}$ defined by

$$
\tilde{\psi} = \rho c \psi. \tag{63}
$$

The factor ρc removes the gradients of ρ or c from the weak form
of the equation and thus improves convergence (see Section 6.5).
The equation solved by $\tilde{\psi}$ is

$$
- \frac{\sigma^2}{\rho c^2} \tilde{\psi} - \frac{2i\omega}{c} \mathbf{u} \cdot \nabla \left( \frac{\tilde{\psi}}{\rho c} \right) - \nabla \cdot \left( \frac{1}{\rho} \nabla \tilde{\psi} \right) = \frac{s}{c}. \quad (64)
$$

In order to use the 2.5D solver, the 3D solution $\tilde{\psi}$ must be ex-
panded as a Fourier series in the longitudinal direction:

$$
\tilde{\psi}(\mathbf{r}, \omega) = \sum_{m=-\infty}^{\infty} \tilde{\psi}^m(\tilde{\mathbf{r}}, \omega) e^{im\phi}. \quad (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*:

$$
-\sigma_m^2 \frac{\tilde{\psi}^m}{\rho c^2} - \frac{2i\omega}{c} \tilde{u} \cdot \tilde{\nabla} \left( \frac{\tilde{\psi}^m}{\rho c} \right) - \tilde{\nabla} \cdot \left( \frac{1}{\rho} \tilde{\nabla} \tilde{\psi}^m \right) = \frac{s^m}{c}, \quad (66)
$$

with

$$
\sigma_m^2 = \omega^2 + 2i\omega\gamma - \frac{2m\omega u_\phi}{\omega} - \frac{m^2 c^2}{\omega^2} \quad (67)
$$

and

$$
s^m(\tilde{\mathbf{r}}, \omega) = \frac{1}{2\pi} \int_0^{2\pi} s(\mathbf{r}, \omega) e^{-im\phi} d\phi. \quad (68)
$$

The 2D gradient and divergence operators are defined by

$$
\begin{align}
\tilde{\nabla} F &= \hat{\omega} \partial_{\omega} F + \hat{z} \partial_z F, \\
\tilde{\nabla} \cdot \mathbf{F} &= \frac{1}{\omega} \partial_{\omega} (\omega F \omega) + \partial_z F_z. \tag{69}
\end{align}
$$

On the outer boundary, we apply a Sommerfeld-like boundary
condition

$$
\partial_n \tilde{\psi}^m = i \frac{\sigma_m}{c} \tilde{\psi}^m. \quad (70)
$$

A boundary condition that would take into account the fast vari-
ation 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 $\mathbf{r}_s = (\tilde{\mathbf{r}}_s, \phi_s)$,

$$
s(\mathbf{r}; \mathbf{r}_s) = \delta(\phi - \phi_s) \delta(\tilde{\mathbf{r}} - \tilde{\mathbf{r}}_s), \quad (71)
$$

and the source coefficients are

$$
s^m(\tilde{\mathbf{r}}; \mathbf{r}_s) = \begin{cases} \frac{1}{2\pi} e^{-im\phi_s} \delta(\tilde{\mathbf{r}} - \tilde{\mathbf{r}}_s) & \text{if } \mathbf{r}_s \text{ is not on the } \hat{z} \text{ axis,} \\ \frac{1}{2\pi} \delta_{m,0} \delta(\tilde{\mathbf{r}} - \tilde{\mathbf{r}}_s) & \text{otherwise.} \end{cases} \quad (72)
$$

Article number, page 8 of 23page.23

Thus the Green’s function may be computed using

$$
\tilde{G}(\mathbf{r}, \mathbf{r}_s, \omega) = \sum_{m=-\infty}^{\infty} \tilde{G}^m(\tilde{\mathbf{r}}, \tilde{\mathbf{r}}_s, \omega) e^{im(\phi - \phi_s)} = \rho(\tilde{\mathbf{r}}) c(\tilde{\mathbf{r}}) G(\mathbf{r}, \mathbf{r}_s, \omega), \quad (73)
$$

where each $\tilde{G}^m$ solves

$$
-\sigma_m^2 \frac{\tilde{G}^m}{\rho c^2} - \frac{2i\omega}{c} \tilde{u} \cdot \tilde{\nabla} \left( \frac{\tilde{G}^m}{\rho c} \right) - \tilde{\nabla} \cdot \left( \frac{1}{\rho} \tilde{\nabla} \tilde{G}^m \right) = \frac{1}{2\pi c} \delta(\tilde{\mathbf{r}} - \tilde{\mathbf{r}}_s). \quad (74)
$$

If $\tilde{\mathbf{r}}_s$ is on the z-axis, the only non-zero component is $\tilde{G}^{m=0}$. It is
important to note that when there is no flow ($\mathbf{u} = \mathbf{0}$), we have
$\sigma_{-m}^2 = \sigma_m^2$ and thus $\tilde{G}^{-m} = \tilde{G}^m$.

6.5. Weak form of the equations

To implement the finite element method, we derive a weak (vari-
tational) formulation of the wave equation for each *m*. Equa-
tion (66) is multiplied by a test function $\Phi$ and integrated over $\Sigma$
with integration element $\omega d\Sigma = \omega d\omega dz$.

Using integration by parts, we obtain

$$
\begin{align}
& - \int_{\Sigma} \Phi \frac{\sigma_m^2}{\rho c^2} \tilde{\psi}^m \, \omega \, d\Sigma \nonumber \\
& - i\omega \int_{\Sigma} \frac{1}{\rho c^2} \tilde{u} \cdot [\tilde{\nabla} \tilde{\psi}^m) \Phi - (\tilde{\nabla} \Phi) \tilde{\psi}^m] \, \omega \, d\Sigma \nonumber \\
& + \int_{\Sigma} \frac{1}{\rho} \tilde{\nabla} \Phi \cdot \tilde{\nabla} \tilde{\psi}^m \, \omega \, d\Sigma - \oint_{\partial \Sigma} \frac{1}{\rho} \Phi (\partial_n \tilde{\psi}^m) \, \omega \, dl \tag{75} \\
& = \int_{\Sigma} \Phi \frac{s^m}{c} \, \omega \, d\Sigma. \nonumber
\end{align}
$$

where we have used mass conservation $\tilde{\nabla} \cdot (\rho \tilde{u}) = 0$ and we have assumed that the flow does not cross the computational boundary $(\tilde{u} \cdot \hat{n} = 0)$ on $\partial \Sigma$. The integral $\oint dl$ is a line integral over the outer boundary $\partial \Sigma$. The boundary condition is used to rewrite the boundary term:

$$
\oint_{\partial \Sigma} \frac{1}{\rho} \Phi (\partial_n \tilde{\psi}^m) \omega \, dl = i \oint_{\partial \Sigma} \Phi \frac{\sigma_m}{\rho c} \tilde{\psi}^m \omega \, dl. \quad (76)
$$

To solve this problem by the FEM, the computational do-
main is subdivided into quadrilateral cells. On each of them, the
test function Φ and the solution of the equation ψ̃m are projected
on a basis {Φi} which, in our study, consists of piecewise con-
tinuous polynomials. By writing this formulation for all the Φi
basis functions, the problem becomes a linear system that can
be solved to obtain ψ̃m. 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 finite-element
solver, Chabassier & Duruflé (2016) considered a simple bench-
mark 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₁ = 0.7R, r₂ = R, and r₃ = 2R,
where R = 700 Mm) where the sound speed cᵢ and density ρᵢ are
constant within the *i*-th 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₃.
<figure><img src="image_4.png" /><figcaption>Fig. 3: Example of a quadrilateral mesh used to compute the solution for the scattering by spherical layers.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Relative L² error between the FEM solution and the analytical solution for the piecewise homogeneous layered medium (Sect. 6.6 and Fig. 3) for non-axial 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).</figcaption></figure>

The following coefficients were chosen arbitrarily and are not intended to represent the Sun:

$$ \begin{cases} \omega/2\pi = 2 \text{ mHz}, \\ 1/\rho_1 = 1.5, & 1/\rho_2 = 2, \quad 1/\rho_3 = 1, \\ 1/(\rho_2 c_1^2) = 0.8, & 1/(\rho_2 c_2^2) = 0.2, \quad 1/(\rho_3 c_3^2) = 1. \end{cases} \quad (77) $$

Here ρᵢ and cᵢ are in units of ρ☉ and 700 km/s respectively.

The analytical solution for the scattering of a plane wave traveling from infinity in the +x̂ direction can be computed as a linear combination of spherical Hankel functions (Chabassier & Duruflé 2016). A non-axial incidence has been selected (the wave vector points in the +x̂ direction) such that all the modes are excited, not only the mode m = 0. In Fig. 4, the relative L₂ 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 (|m| ≤ 37) are sufficient to achieve machine precision accuracy for the above problem.

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

# 7. Time-distance 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 Christensen-Dalsgaard et al. (1996), which is known as Model S. We interpolate the values of ρ and c on the finite-element grid using B-splines.

We implement a Sommerfeld-like 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; the derivatives vanish at 2100 km. The density is deduced by integration using ω²c = ρ¹/²c²∂r[r²∂r(ρ⁻¹/²)]/r² 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 attenuated by the boundary.

To achieve better numerical convergence of the FEM solution, the sound speed and density profiles of Model S are smoothed using B-splines. Eighth-order B-splines are defined by fifteen knots selected to capture all the variations of sound-speed and density throughout the Sun, in particular in the near-surface layers. This ensures that the acoustic cutoff frequency is regular.

## 7.2. Meshing the computational domain

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.7R☉) 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

$$ \lambda_r = 2\pi \left( \frac{\omega^2}{c^2} - \frac{\ell(\ell+1)}{R_{\odot}^2} \right)^{-1/2}, \quad (78) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6: Computational mesh (left) used by the finite element method and the imaginary part of the Green's function at $\omega/2\pi = 3$ mHz (right). The Dirac source is located at radius $0.8R_\odot$ along the z-axis. The damping rate was set to $\gamma/2\pi = 30$ $\mu$Hz.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Radial wavelength for $\omega/2\pi = 9$ mHz and $\ell = 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 tenth-order polynomial within each cell. The wavelength becomes constant in the extended atmosphere described previously.</figcaption></figure>

where $\ell = 15$ is the minimum angular degree that we want to study. In the above expression we fix the frequency at $\omega/2\pi = 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 $\sim \lambda_r/10$. Figure 7 shows the radial wavelength at $\omega/2\pi = 9$ mHz, $\ell = 15$,

<figure><img src="image_9.png" /><figcaption>Fig. 8: Observed full width at half maximum of acoustic modes with radial orders $1 \le n \le 22$ (black dots, Korzennik et al. 2013; Larson & Schou 2015), overplotted with the power law approximation used in the simulations ($2\gamma(\omega)$, orange line).</figcaption></figure>

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 $\gamma = \text{FWHM}/2$, where the FWHM is measured in rad s$^{-1}$. Observational studies show that the FWHM of p-mode ridges is both dependent upon the harmonic degree $\ell$ 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 one-to-one relationship between frequency and wavenumber. However, when modeling the full power spectrum, this approach can only serve as a reasonable estimate. While a wavenumber-dependent damping can be implemented in principle, it is beyond the scope of this paper and we reserve it for a future study.

Figure 8 shows the observed values of FWHM reported by Korzennik et al. (2013) ($100 < \ell < 1000$) and Larson & Schou (2015) ($\ell < 300$) for p modes with radial orders $1 \le n \le 22$ and in the range 1-5.3 mHz. We approximate the attenuation coefficient with a power law in $\omega$ of the form,

$$ \gamma(\omega) = \gamma_0 \left| \frac{\omega}{\omega_0} \right|^{\beta}, \qquad (79) $$

where $\gamma_0/2\pi = 4.29$ $\mu$Hz, $\omega_0/2\pi = 3$ mHz, and $\beta = 5.77$. For frequencies above the acoustic cutoff we fix the attenuation to a constant value of 250 $\mu$Hz.
<figure><img src="image_10.png" /><figcaption>Fig. 9: Snapshots of the inverse temporal Fourier transform of ImG̃(r, ω)/ω, 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 ρ⁻¹/² exp(t/75 min) 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.</figcaption></figure>

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 po-
sition 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 com-
puted 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 z-axis 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 pan-
els the first-arrival wave front is seen propagating away from the source through the core towards the far side of the Sun. At $t = 75$ min, the second-skip waves become visible. Wave packets with higher skip numbers (which take longer to travel) are also seen at the later times. The first-arrival 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 cross-covariance 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 cross-covariance 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:

$$ \psi_{\ell}^{m}(\omega) = \oint d\hat{\mathbf{r}}_{0} Y_{\ell}^{m*}(\hat{\mathbf{r}}_{0}) \int_{V} G(R\hat{\mathbf{r}}_{0}, \mathbf{r}, \omega) s(\mathbf{r}, \omega) \rho d\mathbf{r}, \quad (80) $$

where $d\hat{\mathbf{r}}_{0} = \sin \theta_{0} d\theta_{0} d\phi_{0}$ 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,

$$ \begin{aligned} \mathcal{P}_{\ell}^{m}(\omega) &= \mathbb{E}[|\psi_{\ell}^{m}(\omega)|^2] \\ &= \oint d\hat{\mathbf{r}}_{0} Y_{\ell}^{m}(\hat{\mathbf{r}}_{0}) \oint d\hat{\mathbf{r}}_{0}' Y_{\ell}^{m*}(\hat{\mathbf{r}}_{0}') \bar{C}(R\hat{\mathbf{r}}_{0}, R\hat{\mathbf{r}}_{0}', \omega). \end{aligned} \quad (81) $$

We rewrite $Y_{\ell}^{m}(\hat{\mathbf{r}}_{0}')$ in a frame $\mathcal{R}_{0}$ with polar axis $\hat{\mathbf{r}}_{0} = (\theta_{0}, \phi_{0})$. In this frame, we denote by $\hat{\mathbf{r}} = (\Theta, \Phi)$ the polar angles of $\hat{\mathbf{r}}_{0}$. The rotation of Euler angles $(\alpha, \beta, \gamma) = (\pi, \theta_{0}, \pi - \phi_{0})$ brings $\mathcal{R}_{0}$ to the original frame. Using the rotation formula of spherical harmonics (e.g., Messiah 1960), we have

$$ \begin{aligned} \mathcal{P}_{\ell}^{m}(\omega) &= \oint d\hat{\mathbf{r}}_{0} Y_{\ell}^{m}(\hat{\mathbf{r}}_{0}) \oint d\hat{\mathbf{r}} \left( \sum_{m'=-\ell}^{\ell} Y_{\ell}^{m'}(\hat{\mathbf{r}}) \mathcal{D}_{m'm}^{(\ell)}(\alpha, \beta, \gamma) \right)^{*} \bar{C}(R\hat{\mathbf{r}}, R\hat{\mathbf{r}}, \omega) \\ &= \left( \oint d\hat{\mathbf{r}}_{0} Y_{\ell}^{m}(\hat{\mathbf{r}}_{0}) \mathcal{D}_{0m}^{(\ell)*}(\alpha, \beta, \gamma) \right) \oint d\hat{\mathbf{r}} Y_{\ell}^{0}(\hat{\mathbf{r}}) \bar{C}(R\hat{\mathbf{r}}, R\hat{\mathbf{r}}, \omega), \end{aligned} \quad (82) $$

where, for the sake of simplicity, $\bar{C}$ was assumed to depend only on angular distance $\Theta$ (horizontal isotropy and $\mathbf{u} = 0$).

Using the explicit form of the rotation matrix elements,

$$ \mathcal{D}_{0m}^{(\ell)}(\alpha, \beta, \gamma) = (-)^m \sqrt{\frac{4\pi}{2\ell+1}} Y_{\ell}^{m*}(\beta, \gamma) = \sqrt{\frac{4\pi}{2\ell+1}} Y_{\ell}^{m}(\hat{\mathbf{r}}_{0}), \quad (83) $$

the expression for the power spectrum simplifies to

$$ \begin{aligned} \mathcal{P}_{\ell}^{m}(\omega) &= \sqrt{\frac{4\pi}{2\ell+1}} \oint d\hat{\mathbf{r}} Y_{\ell}^{0}(\hat{\mathbf{r}}) \bar{C}(R\hat{\mathbf{r}}, R\hat{\mathbf{r}}, \omega) \\ &= 2\pi \int_{0}^{\pi} P_{\ell}(\cos \Theta) \bar{C}(\Theta, \omega) \sin \Theta d\Theta \\ &= \frac{\Pi(\omega)}{\omega} \int_{0}^{\pi} \mathrm{Im}G(\Theta, \omega) P_{\ell}(\cos \Theta) \sin \Theta d\Theta, \end{aligned} \quad (84) $$

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>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 f-mode 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 $\ell = 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).</figcaption></figure>

where the last equality is for the special source without flow. The function $P_{\ell}$ is the Legendre polynomial of order $\ell$.

We have complete freedom in the choice of the frequency dependence of the source power, $\Pi(\omega)$. In the rest of the paper we choose a Lorentzian profile:

$$ \Pi(\omega) = \left[ 1 + \left( \frac{|\omega| - \omega_0}{\Gamma/2} \right)^2 \right]^{-1}, \quad (85) $$

where $\omega_0/2\pi = 3.3$ mHz and $\Gamma/2\pi = 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 high-resolution grid in $\Theta$ to increase the spatial Nyquist frequency. In these results, we used 20000 grid points in $\Theta$.
<table><caption>Table 1: Values of the *a*-coefficients for the *l* = 85 and *n* = 8 mode near 3.2 mHz computed from Montjoie simulations and ADIPLS eigenvalue calculations, and measured from SDO/HMI observations. The frequency resolution of the Montjoie simulations is Δν = 0.5 μHz. Simulation #1 was performed using Eq. (6) and simulation #2 includes the second-order term **u** · ∇(**u** · ∇**ψ**). For the ADIPLS calculations, the odd *a*-coefficients are obtained from the rotational sensitivity kernels. The modified ADIPLS calculation ignores gravity terms. See main text.</caption><thead><tr><th>a-coefficients<br>(ℓ, n) = (85, 8)</th><th>Simulation #1<br>(nHz)</th><th>Simulation #2<br>(nHz)</th><th>Modified ADIPLS<br>(nHz)</th><th>ADIPLS<br>(nHz)</th><th>SDO/HMI<br>(nHz)</th></tr></thead><tbody><tr><td>ν&#773;<sub>ln</sub></td><td>3 218 437.9</td><td>3 218 150.0</td><td>3 221 934.5</td><td>3 215 796.6</td><td>3 205 271.5 ± 7.2</td></tr><tr><td>a<sub>1</sub></td><td>437.8</td><td>437.0</td><td>436.8</td><td>435.8</td><td>442.9 ± 0.1</td></tr><tr><td>a<sub>2</sub></td><td>1.8</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.6 ± 0.2</td></tr><tr><td>a<sub>3</sub></td><td>17.9</td><td>19.0</td><td>20.7</td><td>20.7</td><td>22.1 ± 0.2</td></tr></tbody></table>

Figure 10 shows the *m* = 0 power spectrum $\mathcal{P}_{\ell}^{0}(\omega)$ computed from Eq. (84) with the source located at the photosphere on the z-axis. 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 normal-mode 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 $\ell$ = 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 $\mathbf{u} = \Omega(r, \theta)\omega\hat{\phi}$. We use a solar-like differential rotation model specified by

$$ \frac{\Omega(r, \theta)}{2\pi} = \begin{cases} 454 - 55 \cos^2 \theta - 76 \cos^4 \theta \text{ nHz} & r > 0.7 R_{\odot}, \\ 435 \text{ nHz} & r < 0.7 R_{\odot}. \end{cases} \qquad (86) $$

For comparison with the $\ell$ = 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| \le \ell$. For each value of *m* and *ω*, a power spectrum is computed by projecting the cross-covariance $\bar{C}(r, \omega)$ onto spherical harmonics as in Eq. (81). The frequencies of the modes with radial orders *n* = 7, 8, 9 were then extracted from each $\mathcal{P}_{\ell}^{m}(\omega)$ 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 *a*-coefficients as defined by Schou et al. (1994). The mean frequency of the multiplet $\ell$ = 85 and *n* = 8 and the first three *a*-coefficients are given in Table 1 in five cases:

<figure><img src="image_13.png" /><figcaption>Fig. 11: Section of the power spectrum $\mathcal{P}_{\ell}^{m}(\omega)$ for harmonic degree $\ell$ = 85 around frequency 3.2 mHz ($n$ = 7, 8, 9). Red crosses are mode frequencies of the rotationally split p-modes 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).</figcaption></figure>

1. Montjoie simulation using scalar wave equation (6) with a surface delta-function source at the equator. The *a*-coefficients are extracted from fits to the mode frequencies measured from the simulated power spectrum.

2. Montjoie simulation including the second-order term $\mathbf{u} \cdot \nabla(\mathbf{u} \cdot \nabla\mathbf{u})$ on the left-hand side of Eq. (6). We observe that the $a_2$ coefficient (asphericity) vanishes.

3. Eigenvalue calculation for Eq. (6) with (non-rotating) Model S and a free surface boundary condition at height 0.0007$R_{\odot}$ above the photosphere, using a modified version of ADIPLS (Christensen-Dalsgaard 2008) to compute the eigenfrequencies $\bar{\nu}_{\ell n}$ and the rotational kernels (Christensen-Dalsgaard 2003). The modifications to ADIPLS are explained in Appendix A. The odd *a*-coefficients are derived from the first-order perturbation to the mode frequencies. The even *a*-coefficients are zero to this level of approximation.

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

5. Measurements of *a*-coefficients 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
<figure><img src="image_14.png" /><figcaption>Fig. 12: Left: Time-distance diagram computed from Eq. (87) at the height of the Dirac source. The right panel is the observational time-distance diagram computed from the Fourier transform of the SOHO/MDI/Doppler medium-degree power spectrum (Kosovichev et al. 1997). The SOHO/MDI time-distance diagram fades away at large separation distances due to foreshortening.</figcaption></figure>

difference comes from imperfect modeling in the near-surface 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).

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 non-zero, which was (at first) unexpected since our model does not include centrifugal distortion. The non-zero $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 \cdot \nabla(u \cdot \nabla\psi)$ in simulation #2 restores the east-west antisymmetry of the advection of the waves by the flow.

Overall, this comparison between simulated and observed mode frequencies is very encouraging.

## 7.7. Time-distance diagram

For a spherically symmetric solar model, the expectation value of the cross-covariance function in the time domain is

$$ \begin{aligned} \bar{C}(\Theta, t) &= \int_{-\infty}^{\infty} \bar{C}(\Theta, \omega) e^{-i\omega t} d\omega \\ &= \int_{-\infty}^{\infty} \frac{\Pi(\omega)}{2\omega} \mathrm{Im}G(\Theta, \omega) e^{-i\omega t} d\omega, \end{aligned} \qquad (87) $$

Article number, page 14 of 23page.23

where $\Theta$ is the angular distance on the surface between the two observation points. The cross-covariance function is also called the time-distance diagram after Duvall et al. (1993). In Fig. 12 we compare the time-distance diagram computed from our power spectrum to an observed time-distance 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_\ell = [1 - \tanh(0.03\ell - 3)]/2$ for $\ell < 100$ and 0 otherwise, to remove high-degree modes.

Comparisons of the two time-distance diagrams is encouraging. However, the amplitude of the back-skip ridge at $t \sim 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 back-skip branch in the time-distance 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.

# 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_{\text{ref}}$ computed for a highly refined mesh (four times more cells) with high-order 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
<figure><img src="image_15.png" /><figcaption>Fig. 13: Temporal cross-covariance function for three angular distances $\Theta = 30^{\circ}$, $60^{\circ}$, and $90^{\circ}$ 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.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14: Relative difference between $\mathrm{Im} G^{(p)}(\Theta, \omega)$ and $\mathrm{Im} G_{\mathrm{ref}}(\Theta, \omega)$ as defined by Eq. (88) for frequencies $\omega/2\pi = 3$ mHz (left) and $\omega/2\pi = 7$ mHz (right), and different angular distances $\Theta = 30^{\circ}$, $60^{\circ}$ and $90^{\circ}$. Computations were performed without flow (solid lines) and with a meridional flow (dashed lines).</figcaption></figure>

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 non-refined mesh containing 10729 cells. The relative difference between $\mathrm{Im} G^{(p)}$ and $\mathrm{Im} G_{\mathrm{ref}}$, denoted by $\varepsilon^{(p)}$, is plotted in Fig. 14 as a function of the number of degrees of freedom $N^{(p)} = 10729(p+1)^2$. Explicitly,

In Fig. 14, we plot $\varepsilon^{(p)}$ for different angular distances between the source and the receiver $\Theta$ ($30^{\circ}$, $60^{\circ}$, and $90^{\circ}$) and different frequencies ($3$ mHz and $7$ mHz). We see that $\varepsilon^{(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.

$$ \varepsilon^{(p)}(\Theta, \omega) = \frac{|\mathrm{Im} G^{(p)}(\Theta, \omega) - \mathrm{Im} G_{\mathrm{ref}}(\Theta, \omega)|}{\|\mathrm{Im} G_{\mathrm{ref}}(\cdot, \omega)\|_{L^2}}. \qquad (88) $$
<table><caption>Table 2: Computational times and memory usage for single-frequency runs of the radial (1.5D), axisymmetric (2.5D), and 3D methods for concentric shells with constant background coefficients. The number of modes (number of *l* times number of *m* for the radial case and number of *m* for the axisymmetric case) is chosen such that the error is smaller than $10^{-4}$. Computational costs are also given for a single-mode run of the solar case in 1.5D and 2.5D. All computational times are given for single-core computations.</caption><thead><tr><th>Simulation</th><th>Radial</th><th>Axisymmetric</th><th>3D</th></tr></thead><tbody><tr><td>Concentric shells (Sect. 6.6)</td><td></td><td></td><td></td></tr><tr><td>Number of degrees of freedom</td><td>41</td><td>4101</td><td>453001</td></tr><tr><td>Number of modes</td><td>1369</td><td>75</td><td>—</td></tr><tr><td>CPU time</td><td>6.5 s</td><td>42.4 s</td><td>1753.2 s</td></tr><tr><td>Memory usage</td><td>128 MB</td><td>173 MB</td><td>26.3 GB</td></tr><tr><td>Solar model (Sect. 7)</td><td></td><td></td><td></td></tr><tr><td>Number of degrees of freedom</td><td>381</td><td>1072900</td><td></td></tr><tr><td>Number of modes</td><td>1</td><td>1</td><td></td></tr><tr><td>CPU time</td><td>0.48 s</td><td>141 s</td><td></td></tr><tr><td>Memory usage</td><td>744 kB</td><td>6.6 GB</td><td></td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Fig. 15: Travel times $\delta\tau^{(p)}$ as defined in Eq. (89), computed from the difference between cross-covariances $C^{(p)}(\Theta, \omega)$ computed using polynomials of degrees $2 \le p \le 12$ and a reference cross-covariance $C_{\text{ref}}(\Theta, \omega)$. 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.</figcaption></figure>

## 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

$$ \delta\tau^{(p)} = \int_{-\infty}^{\infty} W^*(C^{(p)} - C_{\text{ref}}) d\omega, \qquad (89) $$

where $C^{(p)}$ is the cross-covariance computed from the Green’s functions $G^{(p)}$ and $C_{\text{ref}}$ is the cross-covariance computed from the $G_{\text{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 $\gamma/2\pi = 30$ µHz. In Fig. 15, we show that our method achieves

a travel-time 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 section 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 2D) 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. Travel-time sensitivity kernels

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 travel-time kernels outlined in section 4 that describe the spatial sensitivity of travel times to local perturbations in the interior.

## 9.1. Three-dimensional 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 back-
<figure><img src="image_18.png" /><figcaption>Fig. 16: Sensitivity of the mean travel time, $[\tau(\mathbf{r}_1, \mathbf{r}_2) + \tau(\mathbf{r}_2, \mathbf{r}_1)]/2$, to relative sound-speed perturbations $\delta c/c$ in the interior, with observation point $\mathbf{r}_1$ located on the polar axis and observation point $\mathbf{r}_2$ at latitude of 45°. Panel (a): Cut at $r = 0.95R_\odot$ through the sound-speed 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.95R_\odot$ which corresponds to panel (a).</figcaption></figure>

ground (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 ($\mathbf{r}_1$) and a duplicate is rotated to the receiver location ($\mathbf{r}_2$). The construction of the kernels follows Eqs. (55)-(58).

Figure 16 shows slices through a sound-speed kernel with $\mathbf{r}_1$ at the pole and $\mathbf{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 banana-doughnut 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 travel-time 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. Fig. 16 is a 3D illustration of the sound-speed kernel and also shows the values on the sphere at radius $r = 0.95R_\odot$. 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 hr using 200 cores. The post-processing consists mostly in computing the rotated Green's functions and took approximately 3 hr.

## 9.2. Longitudinally averaged kernels

For an axisymmetric (but not necessary spherically symmetric) background, the Green's functions must be constructed by sum-

ming a sufficient number of $m$ modes. Here, in addition, we only consider axisymmetric perturbations, $q_\alpha = q_\alpha(r, \theta)$, and we wish to determine the 2D spatial sensitivity of the travel times by averaging the kernels over longitude. Useful applications include rotation ($u_\phi$) and meridional circulation ($u_\theta$). For axisymmetric perturbations, we have

$$ \delta\tau(\mathbf{r}_1, \mathbf{r}_2) = 2\pi \sum_{\alpha} \int_{0}^{\pi} d\theta \int_{0}^{R} dr \delta q_{\alpha}(r, \theta) \langle K_{\alpha} \rangle(r, \theta; \mathbf{r}_1, \mathbf{r}_2) \quad (90) $$

where $\langle K_\alpha \rangle$ is the longitudinally averaged kernel

$$ \langle K_{\alpha} \rangle(r, \theta; \mathbf{r}_1, \mathbf{r}_2) = \frac{r^2 \sin \theta}{2\pi} \int_{0}^{2\pi} K_{\alpha}(r, \theta, \phi; \mathbf{r}_1, \mathbf{r}_2) d\phi. \quad (91) $$

Let us define the Green's function and the cross-covariance function in terms of their longitudinal mode components $G^m$ and $C^m$ as follows:

$$ G(\mathbf{r}, \mathbf{r}_1, \omega) = \sum_{m=-\infty}^{\infty} G^m(\tilde{\mathbf{r}}, \tilde{\mathbf{r}}_1, \omega) e^{im(\phi_1 - \phi)}, \\ \bar{C}(\mathbf{r}_1, \mathbf{r}, \omega) = \sum_{m=-\infty}^{\infty} C^m(\tilde{\mathbf{r}}_1, \tilde{\mathbf{r}}, \omega) e^{im(\phi_1 - \phi)}. \quad (92) $$

It is important to note the order of the variables in $G$ and $\bar{C}$ in the above notations. Using the formulation of the kernels given by Eq. (43), we obtain

$$ K_{\alpha}(r, \theta, \phi; \mathbf{r}_1, \mathbf{r}_2) = - \int_{-\infty}^{\infty} d\omega W^* \sum_{m,m'} \mathcal{L}_{\alpha} [G^m(\tilde{\mathbf{r}}_2, \tilde{\mathbf{r}})e^{im(\phi_2 - \phi)}, C^{m'}(\tilde{\mathbf{r}}_1, \tilde{\mathbf{r}})e^{-im'(\phi_1 - \phi)}] \\ - \int_{-\infty}^{\infty} d\omega W^* \sum_{m,m'} \mathcal{L}_{\alpha} [G^{m*}(\tilde{\mathbf{r}}_1, \tilde{\mathbf{r}})e^{-im(\phi_1 - \phi)}, C^{m'*}(\tilde{\mathbf{r}}_2, \tilde{\mathbf{r}})e^{im'(\phi_2 - \phi)}]. $$
<figure><img src="image_19.png" /><figcaption>Fig. 17: Left and center panels: Kernels $\langle K_{u_r} \rangle$ and $\langle K_{u_\theta} \rangle$ for the $r$ and $\theta$ 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 cross-covariance 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.</figcaption></figure>

(93) 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.

Using the explicit expressions for the bilinear operators $L_\alpha$ (Eqs. (55)-(58)), we can then obtain the longitudinally averaged kernels for all perturbations $q_\alpha$.

As an example, the flow kernels can be written as a sum

$$ \langle K_{u_k} \rangle(r, \theta; \mathbf{r}_1, \mathbf{r}_2) = \sum_{m=-\infty}^{\infty} \langle K_{u_k} \rangle^m(r, \theta; \mathbf{r}_1, \mathbf{r}_2) \quad (94) $$

over azimuthal components:

$$ \langle K_{u_k} \rangle^m = 2i\rho r^2 \sin\theta \int_{-\infty}^{\infty} d\omega \, \omega W^* e^{im(\phi_2 - \phi_1)} \\ \times [G^m(\tilde{\mathbf{r}}_2, \tilde{\mathbf{r}}) \tilde{\partial}_k C^{-m}(\tilde{\mathbf{r}}_1, \tilde{\mathbf{r}}) - G^{m*}(\tilde{\mathbf{r}}_1, \tilde{\mathbf{r}}) \tilde{\partial}_k C^{-m*}(\tilde{\mathbf{r}}_2, \tilde{\mathbf{r}})], \quad (95) $$

where the operator $\tilde{\partial}_k$ is either $\tilde{\partial}_r = \partial_r$, $\tilde{\partial}_\theta = \partial_\theta/r$, or $\tilde{\partial}_\phi = im/(r \sin\theta)$.

In practice, $G^m(\tilde{\mathbf{r}}_j, \tilde{\mathbf{r}}, \omega)$ is obtained using generalized seismic reciprocity, $G^m(\tilde{\mathbf{r}}_j, \tilde{\mathbf{r}}, \omega) = G^{-m}(\tilde{\mathbf{r}}, \tilde{\mathbf{r}}_j, \omega; -\mathbf{u})$, performing a simulation with a source at $\tilde{\mathbf{r}}_j$. Using the convenient source of excitation, the cross-covariance is linked to the Green's function by Eq. (49). It is possible to obtain a similar relation for the Fourier coefficients:

$$ C^m(\tilde{\mathbf{r}}_1, \tilde{\mathbf{r}}, \omega) = \frac{\Pi(\omega)}{4i\omega} [G^m(\tilde{\mathbf{r}}, \tilde{\mathbf{r}}_1, \omega) - G^{-m*}(\tilde{\mathbf{r}}, \tilde{\mathbf{r}}_1, \omega; -\mathbf{u})]. \quad (96) $$

Several comments can be made:

- If one of the observation points is on the rotation axis, then only the $m=0$ mode of the Green's function is required. For a measurement between two arbitrary points at the surface of the Sun, the computation 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 cross-covariance as shown by Eq. (96).

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

- If the kernels for points $(\theta_1, \phi_1)$ and $(\theta_2, \phi_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(\phi_2 - \phi_1)}$.

### 9.3. Accuracy of travel-time kernels for flows

Consider a spherically symmetric background reference model with no flow. Figure 17 shows the flow kernels $\langle K_{u_r} \rangle = \langle K_{u_\theta} \rangle^{m=0}$ and $\langle K_{u_\theta} \rangle = \langle K_{u_\phi} \rangle^{m=0}$ with a point $\mathbf{r}_1$ on the pole (at the solar surface) and a point $\mathbf{r}_2$ located at latitude 45°. The kernel $\langle K_{u_\phi} \rangle^{m=0}$ is zero by construction because $\langle K_{u_\phi} \rangle^m$ is proportional to $m$. The 2D kernels for $u_r$ and $u_\theta$ display Fresnel zones as
<figure><img src="image_20.png" /><figcaption>Fig. 18: Left and center panels: Kernels $\langle K_{u_r} \rangle$ and $\langle K_{u_\theta} \rangle$ for the $r$ and $\theta$ components of the flow, using all azimuthal components $|m| \le 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 north-south (blue) and south-north (green) travel directions. The red curve shows their difference.</figcaption></figure>

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 hyperbola-like 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 point-source earthquakes. We refer the reader to the work by Duvall et al. (2006) for an observational study of 2D horizontal sensitivity kernels. The right-hand 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 $\delta 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\text{ms}^{-1}$ at the surface. For the comparison we have considered three separation distances ($30^\circ, 45^\circ$, and $60^\circ$) and the two directions, $\delta\tau(r_1, r_2)$ and $\delta\tau(r_2, r_1)$. The kernel-based 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 off-axis since helioseismic observations are limited to a center-to-limb distance of ~70°. In such cases the azimuthally averaged kernels require the computation of many components $\langle K \rangle^m$. Figure 18 shows the components of flow kernel computed from Eq. (95) for all $m \le 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_\odot$. 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| \le m_{max}$ and compute travel times as a function of $m_{max}$ in the presence of the same meridional flow as used previously. The right-hand 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 $\ell-\omega$ domain is common practice in helioseismology. Several choices of filters have been proposed: filters in $\omega$ space, phase-speed filters, and ridge filters (see, e.g., Gizon & Birch 2005). To interpret any particular travel-time 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).

Symbolically, the filtered observation, $\Psi$, is obtained by applying a filtering operator, $\mathcal{F}$, to the original wavefield, $\psi$:

$$ \Psi(\mathbf{r}, \omega) = \mathcal{F}[\psi(\mathbf{r}, \omega)]. \qquad (97) $$
<figure><img src="image_21.png" /><figcaption>Fig. 19: Left: Unfiltered $\langle K_{u\theta} \rangle$ kernel with $\mathbf{r}_1$ at the pole and $\mathbf{r}_2$ at a co-latitude of 15.36°. Right: Kernel $\langle \mathcal{H}_{u\theta} \rangle$ for filtered observations, where the Gaussian phase-speed filter $F_\ell(\omega)$ is centered at 125.2 km/s with a dispersion of 12.3 km/s, for $\ell$ up to 1000. The ray path is shown by the black line.</figcaption></figure>

The corresponding kernel is obtained by applying the filtering operator twice to the original kernel:

$$ \mathcal{H}_{\alpha}(\mathbf{r}, \mathbf{r}_1, \mathbf{r}_2) = - \int_{-\infty}^{\infty} W^* \mathcal{L}_{\alpha}[\mathcal{G}(\mathbf{r}_2, \mathbf{r}, \omega), \mathcal{C}(\mathbf{r}_1, \mathbf{r}, \omega)] d\omega \\ - \int_{-\infty}^{\infty} W^* \mathcal{L}_{\alpha}^*[\mathcal{G}^*(\mathbf{r}_1, \mathbf{r}, \omega), \mathcal{C}^*(\mathbf{r}_2, \mathbf{r}, \omega)] d\omega, \quad (98) $$

where the filtered Green's function and cross-covariance are

$$ \mathcal{G}(\mathbf{r}_1, \mathbf{r}, \omega) = \mathcal{F}_1[G(\mathbf{r}_1, \mathbf{r}, \omega)], \quad (99) $$

$$ \begin{aligned} \mathcal{C}(\mathbf{r}_1, \mathbf{r}, \omega) &= \mathcal{F}_1[\overline{\mathcal{C}(\mathbf{r}_1, \mathbf{r}, \omega)}] \\ &= \frac{\Pi(\omega)}{4i\omega} [\mathcal{G}(\mathbf{r}_1, \mathbf{r}, \omega; -\mathbf{u}) - \mathcal{G}^*(\mathbf{r}_1, \mathbf{r}, \omega)], \quad (100) \end{aligned} $$

and the W function is computed with the choice $C_{\text{ref}} = \mathcal{C}$. Here, $\mathcal{F}_1$ indicates that the filtering has to be done with respect to the point $\mathbf{r}_1$.

The filtered Green's function may be obtained by filtering the delta source function. To see this, we use generalized seismic reciprocity

$$ \begin{aligned} \mathcal{G}(\mathbf{r}_1, \mathbf{r}, \omega) &= \mathcal{F}_1[G(\mathbf{r}, \mathbf{r}_1, \omega; -\mathbf{u})] \\ &= \int_V G(\mathbf{r}, \mathbf{r}', \omega; -\mathbf{u}) \mathcal{F}_1[\delta(\mathbf{r}_1 - \mathbf{r}')] d\mathbf{r}', \quad (101) \end{aligned} $$

where

$$ \begin{aligned} \mathcal{F}_1[\delta(\mathbf{r}_1 - \mathbf{r}')] &= \delta(\mathbf{r}_1 - \mathbf{r}') \mathcal{F}_1 \left[ \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} Y_{\ell}^m(\hat{\mathbf{r}}_1) Y_{\ell}^{m*}(\hat{\mathbf{r}}') \right] \\ &= \delta(\mathbf{r}_1 - \mathbf{r}') \sum_{\ell=0}^{\infty} F_{\ell}(\omega) \frac{2\ell + 1}{4\pi} P_{\ell}(\hat{\mathbf{r}}_1 \cdot \hat{\mathbf{r}}') \quad (102) \end{aligned} $$

is the filtered delta function source. In the above expression $F_\ell(\omega)$ can be either a phase-speed 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 non-filtered Green's function and to perform the filtering a posteriori. In this case, the Green's function $G(\mathbf{r}_1, \mathbf{r}, \omega)$ depends only on the depths $\mathbf{r}_1$ and $\mathbf{r}$ and on the angular distance $\Theta_1$ between the two points

$$ G(\mathbf{r}_1, \mathbf{r}, \omega) = G(\mathbf{r}_1, \mathbf{r}, \Theta_1, \omega), \quad (103) $$

where $\cos \Theta_1 = \hat{\mathbf{r}}_1 \cdot \hat{\mathbf{r}}$. The filtered Green's function takes the form

$$ \mathcal{G}(\mathbf{r}_1, \mathbf{r}, \omega) = \sum_{\ell=-\infty}^{\infty} F_{\ell}(\omega) G_{\ell}(\mathbf{r}_1, \mathbf{r}, \omega) P_{\ell}(\cos \Theta_1), \quad (104) $$

where

$$ G_{\ell}(\mathbf{r}_1, \mathbf{r}, \omega) = \frac{2\ell + 1}{2} \int_{0}^{\pi} G(\mathbf{r}_1, \mathbf{r}, \Theta, \omega) P_{\ell}(\cos \Theta) \sin \Theta \, d\Theta \quad (105) $$

is the projection of $G$ on the Legendre polynomials. Using seismic reciprocity we obtain $G_\ell(\mathbf{r}_1, \mathbf{r}, \omega) = G_\ell(\mathbf{r}, \mathbf{r}_1, \omega)$, and Eq. (104) implies

$$ \mathcal{G}(\mathbf{r}_1, \mathbf{r}, \omega) = \mathcal{G}(\mathbf{r}, \mathbf{r}_1, \omega). \quad (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 post-processing operation on $\mathbf{r}$ at fixed source position $\mathbf{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 $\mathbf{r}_1$ is located at the pole and $\mathbf{u} = \mathbf{0}$. As mentioned previously, only the $m = 0$
component of the Green's function is needed. Figure 19 shows
the effect of applying a phase-speed filter to the kernel $\langle K_{\theta\theta} \rangle$ with
$r_1$ at the pole and $r_2$ at a co-latitude of 15.36°. The phase-speed
filter is centered at 125.2 km/s with a width of 12.3 km/s. 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_\theta$ 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 helio-
seismology by solving the forward problem in the frequency
domain. For the sake of simplicity, we considered a simpli-
fied scalar acoustic wave equation and assumed spatially uncor-
related sources of excitation distributed through the Sun. Un-
der such conditions, the cross-covariance can be obtained di-
rectly from the imaginary part of the frequency-domain 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 normal-mode 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 us-
ing existing Montjoie vectorial setups (e.g. Péron et al. 2016, for
Maxwell's equations). The scalar equation captures the propaga-
tion of acoustic waves through a solar-like medium and leads to
an oscillation power spectrum that compares well with observa-
tions. The present setup will be very useful to test new methods,
including instrumental and projection effects, but also to inter-
pret existing travel-time measurements for rotation, meridional
circulation, and axisymmetric structures like the average super-
granule.

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 travel-time noises, defined through the noise covariance ma-
trix Λ = E[n n^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 non-linear inverse prob-
lem defined in terms of the partial differential equation. Differ-
ent methods exist such as the Landweber iteration (Hanke et al.
1995) and Newton-type methods such as the iteratively regular-
ized Gauss-Newton (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 work-
load 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.

Acknowledgements. LG and ACB developed the general theoretical framework. HB, JC, and MD developed tools to solve the wave equation (Montjoie package). CSH, ML, and DF implemented the tools in the solar context, tuned the power spectrum, and computed the kernels. All authors contributed to the writing of the paper. LG 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 Global Oscillation Network Group
is funded by the National Science Foundation. The HMI project is supported
by NASA contract NAS5-02139. 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

Ball, W. H., Beeck, B., Cameron, R. H., & Gizon, L. 2016, A&A, 592, A159

Bergot, M., Cohen, G., & Duruflé, M. 2010, Journal of Scientific Computing, 42,
345

Birch, A. C. & Gizon, L. 2007, Astron. Nachr., 328, 228

Birch, A. C., Gizon, L., Hindman, B. W., & Haber, D. A. 2007, ApJ, 662, 730

Birch, A. C., Kosovichev, A. G., & Duvall, Jr., T. L. 2004, ApJ, 608, 580

Bogdan, T. J., Hindman, B. W., Cally, P. S., & Charbonneau, P. 1996, ApJ, 465,
406

Böning, V. G. A., Roth, M., Zima, W., Birch, A. C., & Gizon, L. 2016, ApJ, 824,
49

Bottero, A., Cristini, P., Komatitsch, D., & Asch, M. 2016, ArXiv:1604.04237

Burston, R., Gizon, L., & Birch, A. C. 2015, Space Sci. Rev., 196, 201

Cameron, R., Gizon, L., & Duvall, Jr., T. L. 2008, Sol. Phys., 251, 291

Chabassier, J. & Duruflé, M. 2016, High-order finite element method for solving convected Helmholtz equation in radial and axisymmetric domains: application to helioseismology, Research Report 8893, INRIA, https://hal.inria.fr/hal-01295077

Christensen-Dalsgaard, J. 2003, Lecture Notes on Stellar Oscillations, 5th edn.
(Aarhus University)

Christensen-Dalsgaard, J. 2008, Ap&SS, 316, 113

Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272,
1286

Dikpati, M. & Choudhuri, A. R. 1995, Sol. Phys., 161, 9

Duruflé, M. 2006, Doctoral thesis, Numerical integration and high-order finite element methods applied to time-harmonic Maxwell equations, ENSTA Paris-Tech

Duvall, Jr., T. L., Birch, A. C., & Gizon, L. 2006, ApJ, 646, 553

Duvall, Jr., T. L. & Hanasoge, S. M. 2013, Sol. Phys., 287, 71

Duvall, Jr., T. L., Jefferies, S. M., Harvey, J. W., Osaki, Y., & Pomerantz, M. A. 1993, ApJ, 410, 829

Fournier, D., Gizon, L., Hohage, T., & Birch, A. C. 2014, A&A, 567, A137

Gizon, L. & Birch, A. C. 2002, ApJ, 571, 966

Gizon, L. & Birch, A. C. 2004, ApJ, 614, 472

Gizon, L. & Birch, A. C. 2005, Living Reviews in Solar Physics, 2, 6

Gizon, L., Birch, A. C., & Spruit, H. C. 2010, Annual Review of Astronomy and Astrophysics, 48, 289

Haber, D., Hindman, B., Toomre, J., & Thompson, M. 2004, Solar Physics, 220,
371

Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annual Review of Fluid
Mechanics, 48, 191

Hanasoge, S. M., Birch, A., Gizon, L., & Tromp, J. 2011, ApJ, 738, 100

Hanasoge, S. M. & Duvall, Jr., T. L. 2007, Astron. Nachr., 328, 319

Hanasoge, S. M., Duvall, Jr., T. L., & Sreenivasan, K. R. 2012, Proceedings of the National Academy of Sciences, 109, 11928

Hanke, M. 1997, Numerical Functional Analysis and Optimization, 18, 971

Hanke, M., Neubauer, A., & Scherzer, O. 1995, Numerische Mathematik, 72, 21

Hill, F., Stark, P. B., Stebbins, R. T., et al. 1996, Science, 272, 1292

Jackiewicz, J., Gizon, L., & Birch, A. C. 2008, Sol. Phys., 251, 381

Jackiewicz, J., Gizon, L., Birch, A. C., & Duvall, Jr., T. L. 2007, ApJ, 671, 1051

Jouve, L., Brun, A. S., Arlt, R., et al. 2008, A&A, 483, 949

Korzennik, S. G., Rabello-Soares, M. C., & Schou, J. 2004, ApJ, 602, 481

Korzennik, S. G., Rabello-Soares, M. C., Schou, J., & Larson, T. P. 2013, ApJ,
772, 87

Kosovichev, A. 1996, ApJ, 461, L55

Kosovichev, A. G., Duvall, Jr., T. L., & Scherrer, P. H. 2000, Sol. Phys., 192, 159

Kosovichev, A. G., Schou, J., Scherrer, P. H., et al. 1997, Sol. Phys., 170, 43

Kumar, P., Duvall, Jr., T. L., Harvey, J. W., et al. 1990, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 367, Progress of Seismology of the Sun and Stars, ed. Y. Osaki & H. Shibahashi, 87

Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67

Larson, T. P. & Schou, J. 2015, Sol. Phys., 290, 3221

Liang, Z.-C. & Chou, D.-Y. 2015, ApJ, 809, 150

Lynden-Bell, D. & Ostriker, J. P. 1967, MNRAS, 136, 293

Marquering, H., Dahlen, F. A., & Nolet, G. 1999, Geophysical Journal International, 137, 805

Messiah, A. 1960, Mécanique quantique, Vol. 2 (Paris: Dunod)
Nissen-Meyer, T., van Driel, M., Stähler, S. C., et al. 2014, Solid Earth, 5, 425

Papini, E., Gizon, L., & Birch, A. C. 2014, Sol. Phys., 289, 1919

Péron, V., Schmidt, K., & Duruflé, M. 2016, SIAM J. Appl. Math., 76, 1031

Rajaguru, S. P. & Antia, H. M. 2015, ApJ, 813, 114

Rosenthal, C. S., Christensen-Dalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689

Schou, J., Christensen-Dalsgaard, J., & Thompson, M. J. 1994, ApJ, 433, 389

Schunker, H., Cameron, R. H., Gizon, L., & Moradi, H. 2011, Sol. Phys., 271, 1

Sekii, T. & Shibahashi, H. 2003, in ESA Special Publication, Vol. 517, GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. Sawaya-Lacoste, 389–392

Snieder, R. 2007, J. Acoust. Soc. Am., 121, 2637

Snieder, R. & Larose, E. 2013, Annual Review of Earth and Planetary Sciences, 41, 183

Snieder, R., Miyazawa, M., Slob, E., Vasconcelos, I., & Wapenaar, K. 2009, Surveys in Geophysics, 30, 503

Švanda, M., Gizon, L., Hanasoge, S. M., & Ustyugov, S. D. 2011, A&A, 530, A148

van Ballegooijen, A. A. & Choudhuri, A. R. 1988, ApJ, 333, 965

van Driel, M., Krischer, L., Stähler, S. C., Hosseini, K., & Nissen-Meyer, T. 2015, Solid Earth, 6, 701

van Driel, M. & Nissen-Meyer, T. 2014, Geophysical Journal International, 199, 880

Zhao, J., Bogart, R. S., Kosovichev, A. G., Duvall, Jr., T. L., & Hartlep, T. 2013, ApJ, 774, L29

Zienkiewicz, O., Taylor, R., & Zhu, J. 2005, The Finite Element Method Set, sixth edn. (Oxford: Butterworth-Heinemann)

## Appendix A: Solving the eigenvalue problem using ADIPLS

We modified ADIPLS to solve the eigenvalue problem

$$ \mathcal{H}[\xi] = \omega^2 \xi \qquad (A.1) $$

where $\mathcal{H}$ is given by Eq. (2) and we neglected the gravity terms. Using the same notation as the one used by Christensen-Dalsgaard (2008), the corresponding eigenvalue problem for the eigenfrequency and the eigenfunction of a mode takes the form

$$ \omega^2 \xi_r = -\frac{1}{\rho} \frac{d}{dr} \left[ \rho c^2 \left( \frac{d\xi_r}{dr} + 2 \frac{\xi_r}{r} - L^2 \frac{\xi_h}{r} \right) \right], \qquad (A.2) $$

$$ \omega^2 \xi_h = -c^2 \left( \frac{1}{r} \frac{d\xi_r}{dr} + 2 \frac{\xi_r}{r^2} - L^2 \frac{\xi_h}{r^2} \right), \qquad (A.3) $$

where $L^2 = \ell(\ell + 1)$, and $\xi_r(r)$ and $\xi_h(r)$ are the radial and hori- zontal eigenfunctions. The equations can be rewritten as

$$ r \frac{d\xi_r}{dr} = -2\xi_r + \left( L^2 - \frac{\omega^2 r^2}{c^2} \right) \xi_h, \qquad (A.4) $$

$$ r \frac{d\xi_h}{dr} = \xi_r - \left( 1 + \frac{d \ln \rho}{d \ln r} \right) \xi_h, \qquad (A.5) $$

which, in ADIPLS adimensionalized form (see Christensen- Dalsgaard 2008), become

$$ x \frac{dy_1}{dx} = -2y_1 + \left(1 - \frac{V_g}{\eta}\right) y_2, \qquad (A.6) $$

$$ x \frac{dy_2}{dx} = L^2 y_1 + (A + V_g - 1) y_2. \qquad (A.7) $$

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 differ- ential rotation by evaluating the integral

$$ \delta\omega_{nlm} = m \int_{0}^{R} dr \int_{0}^{\pi} K_{nlm}(r, \theta) \Omega(r, \theta) r d\theta \qquad (A.8) $$

where the rotation profile $\Omega(r, \theta)$ is defined by Eq. (86), and $K_{nlm}(r, \theta)$ are the rotational sensitivity kernels (see Chapter 8.4 in Christensen-Dalsgaard 2003), which depend on $\xi_r$ and $\xi_h$.

Article number, page 22 of 23page.23

<figure><img src="image_22.png" /><figcaption>Fig. B.1: Plot of $u_r(r, \theta)$ (left) and $u_\theta(r, \theta)$ (right) at $\theta = 45^\circ$ of the meridional flow model discussed in the Appendix for a maximum flow velocity $U = 20$ m/s at the surface.</figcaption></figure>

## Appendix B: Meridional flow model

Various models of meridional flow cells can be found in the lit- erature (e.g., van Ballegooijen & Choudhuri 1988; Dikpati & Choudhuri 1995; Jouve et al. 2008). In these models the flow is expressed through a stream function $\Psi$ as follows:

$$ u(r, \theta) = \frac{1}{\rho} \nabla \times (\Psi \hat{\phi}). \qquad (B.1) $$

This enforces mass conservation, $\nabla \cdot (\rho u) = 0$. However, we can- not use directly the various expressions of the meridional flow from previous work, as they were not computed with the density $\rho(r)$ of model S.

Let us define a stream function $\Psi(r, \theta)$ for $r_b \le r \le r_t$ and $0 \le \theta \le \pi$, where $r_b$ and $r_t$ refer to the bottom and the top of the convection zone, respectively. The function $\Psi$ is set to zero outside this region. Let us look for a solution of the form

$$ \Psi(r, \theta) = \rho_t r_t U f(r) g(\theta), \qquad (B.2) $$

where $f(r)$ and $g(\theta)$ are dimensionless functions to be deter- mined, $U$ sets the amplitude of the flow, and $\rho_t = \rho(r_t)$ is the density at the top of the cell. The flow is given by

$$ u_r(r, \theta) = U \frac{\rho_t r_t}{\rho r} \frac{f(r)}{\sin \theta} \frac{d}{d\theta} [g(\theta) \sin \theta], \qquad (B.3) $$

$$ u_\theta(r, \theta) = -U \frac{\rho_t r_t}{\rho r} g(\theta) \frac{d}{dr} [r f(r)]. \qquad (B.4) $$
For one cell per hemisphere with $r_b < r < r_t$, the functions $f$ and $g$ satisfy the following conditions:

$$f(r_b) = f(r_t) = 0 \text{ and } g(0) = g(\pi) = 0. \qquad (B.5)$$

For the latitudinal dependence, we choose

$$g(\theta) = \sin(2\theta). \qquad (B.6)$$

We seek the function $f$ in terms of the function

$$h(r) = -\frac{\rho_t r_t}{\rho r} \frac{d}{dr}(rf), \qquad (B.7)$$

which controls the radial profile of $\rho u_\theta$. We have

$$f(r) = -\frac{1}{\rho_t r_t} \int_{r_b}^{r} h(r') \rho(r') r' dr'. \qquad (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:

$$h(r_t) = 1, h(r_b) = 0, \frac{dh}{dr}(r_b) = 0 \qquad (B.9)$$

While the first two conditions are intuitive, the third is arbitrary (other choices would be possible). We choose $h$ as a fourth-order polynomial

$$h(r) = (r - r_b)^2 (r - r_h)(r - a)/r_t^4, \qquad (B.10)$$

where $a = r_t - r_t^4(r_t - r_b)^{-2}(r_t - r_h)^{-1}$ is implied by $h(r_t) = 1$. In this paper, we set $r_b = 0.7R_\odot$ and $r_t = R_\odot$. 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.859R_\odot$. For the amplitude of the flow, we take $U = 20$ m/s. The radial and co-latitudinal flows, given by

$$u_r(r, \theta) = U \frac{\rho_t r_t}{\rho r} f(r) 2[\cos^2 \theta + \cos(2\theta)], \qquad (B.11)$$

$$u_\theta(r, \theta) = U h(r) \sin(2\theta), \qquad (B.12)$$

are plotted as a function of $r$ at co-latitude $\theta = 45^\circ$ in Fig. B.1.