© J. Schou and A. C. Birch 2020

1. Introduction

Soon after frequencies of solar oscillations became available, attempts were made to infer the internal structure of the Sun, typically using inverse methods. It was determined that the frequencies could not be fitted, and it was quickly realized that a major source of the problem was close to the surface (for details see e.g., Christensen-Dalsgaard & Thompson 1997). It was not a surprise that the surface caused problems as the physics in the surface layers is complicated by nonadiabatic effects, high velocities in the granulation, and a rapid variation of the thermodynamic quantities (Brown 1984). Fortunately, as the effect is near the surface, it must have a particular functional form. Given this, a term with a suitably chosen parametrization was added to the expression for the mode frequencies used in the inversions and the parameters were fit as part of the inversion procedure. This term became known as the “surface term”. While this is clearly unsatisfactory from the point of view of a physical understanding, the number of measured frequencies is large enough that the loss of degrees of freedom could be tolerated. As no data were available for other stars, where the loss of degrees of freedom would be more significant, a solution to the problem was not urgent and the approach allowed for significant advances to be made in inferring the deeper solar structure (see, for example Bahcall et al. 2005).

Rosenthal et al. (1999) suggested that part of the problem could be solved by patching averages of 3D hydrodynamic simulations onto 1D evolutionary models. Unfortunately, as they also pointed out, determining how this averaging should be performed is not trivial and the patching approach does not solve the entire problem. In particular, the mode physics is not addressed. A particular problem they point out is that the wave perturbation to the turbulent pressure needs to be modeled. As another example of a dynamical effect, Baldner & Schou (2012) find that the small-scale turbulence causes an apparent net upflow near the solar surface and that this may explain a number of center-to-limb problems observed in helioseismology.

As oscillations started to be observed on other stars, the problem of the surface term became acute. In asteroseismology, few modes are observed and simply fitting a parametrized function results in the loss of a substantial fraction of the available information. Current approaches include scaling the solar surface term and using parametrizations with a small number of terms (e.g., Kjeldsen et al. 2008; Ball & Gizon 2014).

Another approach, employed by Piau et al. (2014), Sonoi et al. (2015, 2017), and Ball et al. (2016), for example, is to use the suggestion of Rosenthal et al. (1999) and to patch averaged 3D simulations onto 1D evolutionary models. Indeed, this has resulted in a significant improvement in the fits with residuals that were reduced from up to 12 μHz below 4000 μHz to about 3 μHz. While this is an important advancement, the residuals are still orders of magnitude larger than the observed errors. Unfortunately, the various dynamical effects are still not properly modeled. Houdek et al. (2017) addresses this problem based on a model of convection. However, while their approach does reduce the residuals down to about 2 μHz, there are still significant residuals, and it is unclear how good the model is and if it is applicable to other stars. Refinements of these approaches, such as including evolutionary effects, were also considered by Jørgensen et al. (2018), Jørgensen & Weiss (2019), and Mosumgaard et al. (2020), among others.

In the method proposed in this paper we take the approach of Rosenthal et al. (1999) one step further by effectively patching the eigenfunctions from a 3D simulation onto those from a 1D model. More precisely, we horizontally average 3D models, as suggested by Houdek et al. (2017) and Ball et al. (2016), for example, and calculate the resonance frequencies as they did, compare those to the frequencies of the modes seen in the simulations, and scale the differences to the solar case, thereby including all of the physical effects that were neglected in the averaging. Doing so allows us to separate the issues pertaining to structural changes and inversions from the mode physics issues; the latter is the subject of interest of this paper.

In Sect. 2 we discuss how to determine the frequency perturbations and apply them to the solar case. In Sect. 3 various results are presented, in Sect. 4 we discuss the results, and in Sect. 5 we conclude.

Apart from the discussion in Sect. 3.2, we do not discuss the behavior of the eigenfunctions near the surface and the physical interpretation of those results. We choose to defer this work to a later paper.

2. Method

The method we propose consists of several steps, which are outlined in the subsections below. In Sect. 2.1 we describe the simulations used, in Sect. 2.2 how the eigenfunctions and frequencies are extracted is explained, in Sect. 2.3 we discuss how to calculate the frequency perturbation relative to the horizontally averaged models using the standard equations, and in Sect. 2.4 how to scale the results to the solar case is explored.

2.1. Models used

The MURaM simulations used here were computed as described in Beeck et al. (2013) and Rempel (2017) for a nonmagnetic Sun. However, they were modified in terms of spatial extent and resolution, as described below.

The simulations were performed in Cartesian boxes with periodic boundary conditions in the horizontal directions. The upper boundary is nonpenetrative, while the flows were allowed to pass through the lower boundary; the entropies of the upflows were selected so as to obtain the desired luminosity. Near the bottom, there is also a diffusive layer. The vertical grid is uniformly spaced. We used a constant gravity of 27 400 cm s⁻¹. No corrections were made for curvature effects. For details on the numerical schemes, the treatment of the radiative transfer, equation of state, etc., we refer to the papers cited above.

Ideally, the simulations would be deep, as this would allow for a large number of radial orders. Unfortunately, while it is not important for the present use to have an accurate stratification deep in the box, it is nonetheless best to avoid constraining the convection by having a box that is too narrow; the boxes should be at least twice, or preferably three times, as wide as their extent below the surface. Using deep boxes also implies that the maximum sound speed is larger, requiring smaller time steps, and this also increases the thermal relaxation time. Given these problems and the need for very long simulation times in order to resolve the modes and improve the signal-to-noise ratio (S/N), we thus had to compromise and instead chose to use several simulations with the same width, but different depths, in order to obtain a reasonable frequency coverage. Table 1 summarizes the sizes, resolutions, time span, and other properties of the simulations, as well as the parameters used in the analysis. The main simulations used here are cases 10, 12, and 13. These are all very close in temperature and only slightly warmer than the Sun. Cases 11 and 15 are higher and lower resolution versions of case 10, intended to show that an adequate resolution was used. In order to save space, the snapshots were generally deleted and only averaged quantities were kept.

Throughout this paper, we use the symbols p_gas for the gas pressure, p_turb for the turbulent pressure, p_tot = p_gas + p_turb for the total pressure, ρ for the density, v_z for the vertical velocity, and Γ₁ for the first adiabatic exponent. To define a height coordinate that is consistent with the 1D models, where τ is not given, the height z = 0 is given as the location where horizontally averaged density is 2.3 × 10⁻⁷ g cm⁻³. This was chosen so as to make the average height of the τ_Rosseland surface coincide with z = 0 for the regular resolution cases, and to be consistent with Beeck et al. (2013). The distance from the center of the Sun is given by r = z + R_⊙, where R_⊙ is the solar radius.

2.2. Extraction of eigenfunctions

One way to characterize the acoustic oscillations in the simulations is to measure their eigenfunctions and eigenfrequencies. In order to extract the eigenfunctions from a simulation, we compute the unweighted horizontal averages of p_gas, ρ, Γ₁, and v_z at each saved time step and height in the simulation domain. At each time and height, we subtract the corresponding horizontal average to obtain a set of residuals from which the horizontally averaged p_turb is computed. From that, we compute the horizontally averaged total pressure p_tot = p_gas + p_turb.

These horizontally-averaged quantities are averaged in time to obtain the mean thermodynamic quantities. We denote the resulting quantities with an overbar, for example ${\overline{p}}_{\mathrm{gas}}$ is the horizontally- and time-averaged gas pressure.

The residuals with respect to the time averages are due to both the radial acoustic oscillations and the time-dependent convection. As described later in this section, we use frequency filtering to mostly separate these two contributions.

The wave component for each quantity is then linearly interpolated at each height to a uniform grid in time with twice the average resolution of the simulation grid and Fourier transformed. Figure 1 shows an example of a resulting power spectrum. The power spectrum shows several resonances; these are the radial oscillations. The power spectrum also shows a number of weak peaks in between the main modes. An examination of the spectra of the modes with a horizontal wavenumber of k_h ≠ 0 shows that these peaks are due to leaks from modes with k_h = 1 and $k_{\mathrm{h}}=\sqrt{2}$ (with k_h in units of wavelengths across the box width), but they appear wider in the k_h = 0 spectra. A likely explanation for their appearance is that the convection has large components at those horizontal wavenumbers, which interact with the corresponding modes, leading to power with k_h = 0. As the convection is time variable, it also leads to widening of the peaks.

Fig. 1. Power spectrum of the pressure perturbation for case 10, averaged over the depth range used for the SVD. The power was smoothed to reduce the effect of the stochastic excitation. The large peaks at 1856 μHz, 2812 μHz, and 3858 μHz represent the modes with n = 0, 1, and 2, respectively, where n is the number of nodes in the displacement eigenfunction. The n = 3 peak around 5000 μHz is not analyzed further, as it is very close to the acoustic cutoff frequency. The small peaks around 2100 μHz and 2300 μHz are leaks from modes with kh ≠ 0.

Based on the stochastic nature of the modes, we expect each physical quantity to be separable in time (or equivalently frequency) and height near a resonance, with the former showing a Lorentzian power profile in frequency. For example, in the neighborhood of the nth resonance, the Fourier transform of a physical variable q has the form

$$\begin{array}{c}\hfill q(z,\omega )\approx q_n^{\prime }(z)f_n(\omega ),\end{array}$$(1)

where $q_n^{\prime }$ is the eigenfunction (denoted here by a superscript ′) in height, f_n is the Fourier transform of the time dependence, and ω is the frequency. In general both $q_n^{\prime }$ and f_n are complex. The singular value decomposition (SVD) decomposes a signal into a sum of such terms:

$$\begin{array}{c}\hfill q(z,\omega )={\displaystyle \sum _i{w}_i{q}_i^{\prime }(z)f_i(\omega ),}\end{array}$$(2)

where the sum is taken over the singular vectors q′ and f, which are generally complex, and the singular values w, which are real and sorted in decreasing value. In particular, the first term represents the best fit of Eq. (1). Performing the SVD on the entire dataset turns out not to be desirable for a number of reasons. The convective signal is very large near the photosphere and would introduce noise in the frequency dependence. There is significant convective power at low frequencies, which is not interesting for the purposes of this work. Also, as the eigenfunctions may be similar over a limited fitting range, they are not cleanly separated. A joint SVD using all physical quantities would be another option, but this approach is not effective as the S/Ns for different physical quantities can be very different.

To overcome these issues, we use the pressure perturbation over the interior part of the simulations (excluding the lower 10% and from 0.5 Mm below the surface and above) for the SVD, as this height range has a good S/N. Also, we only use a small range in frequency around each observed peak, which is ± four half width at half maximum (HWHM), but in no case is it smaller than ±20 μHz or larger than ±200 μHz. In all cases, it was found that a single component, by far, dominates the SVD, capturing above 99% of the variance for modes below 4500 μHz and in some cases over 99.999%. We thus only study the dominant component and suppress the subscript i in the following. In other words, we simply use the SVD as a way to fit the data to the model given by Eq. (1). As this procedure uses almost all of the power in a peak, it gives substantially less noisy results than using the spatial dependency at the peak frequency, for instance. The complex conjugate of the dominant frequency dependence vector (f in Eq. (2)) is then multiplied on all variables at all heights to obtain the vertical eigenfunctions for these variables.

Given the stochastic nature of the excitation of the modes, the mode frequency ω_n was then computed by fitting a Lorentzian to |f(ω)|², using a maximum likelihood fit (Anderson et al. 1990). Given the mode frequency, we also calculate the displacement eigenfunction $\xi =-v_z^{\prime }/(i{\omega }_n)$.

The SVD singular vectors are degenerate up to a phase factor. For convenience, we shift the phases at each radial order to make the displacement eigenfunction real 50 gridpoints above the bottom.

Figure 2 shows examples of eigenfunctions. Beyond the classic behavior, which is defined by Eqs. (6) through (10) below, the most notable feature is probably the large imaginary component near the surface, which was previously ascribed by Balmforth (1992) to nonadiabatic effects and nonlocal convection effects and later by Baldner & Schou (2012) to the large horizontally-averaged vertical flow near the surface. Nordlund & Stein (2001) and Stein & Nordlund (2001) also extracted eigenfunctions from simulations, but unfortunately they did not show the imaginary component. The former paper also discusses some of the formalism needed to analyze the waves. While this is clearly interesting in terms of understanding the relevant physics, we defer the study of the eigenfunctions to a future paper, except for a brief discussion that appears in Sect. 3.2.

Fig. 2. Displacement eigenfunctions corresponding to the three largest peaks in Fig. 1. The eigenfunctions were multiplied by the square root of the background density to improve the visibility and as the square of this quantity is proportional to the energy density. Colors indicate the different modes. Solid lines indicate the real part, and dashed lines illustrate the imaginary part. For ease of display, the scaled eigenfunctions were normalized to unity at the bottom boundary.

2.3. Determination of the frequency perturbation

Our goal is to compute the frequency perturbation caused by the physics not captured by the simple patched models in a solar model that extends from the center of the star to the surface. We tackle this calculation in two steps. In the first step, we compute the change in the resonance frequencies caused by near-surface convection for an extended MURaM model. In the second step, which is described in the next subsection, we use a mode-mass correction to compute the frequency changes that would be expected for a solar model that extends down to the center of the star. An advantage of this approach is that it avoids most issues resulting from different physics in the 3D box model and the 1D evolutionary models. In particular, issues relating to the curvature and variable gravity are avoided to lowest order, as are issues related to large scale structural problems. Alternatives to this approach are discussed in Sect. 4 and tests of various assumptions appear in the appendix.

The background model, as well as the eigenfunctions computed in the previous section, have artifacts near the bottom boundary. These artifacts are due to the effect of the boundary on the stratification and to the diffusive layer near the bottom. If a layer at the bottom is simply removed, the boundary condition is ill defined. To mitigate these problems, and as the n = 0 mode in any case does not have a zero crossing of the displacement eigenfunction within the simulation domain, we extend the MURaM model downward to −50 Mm to capture at least one zero crossing of the displacement eigenfunction for each mode. This extended model, which is denoted with ℳ_extend, is constructed as follows: First the part of the model below z_extend is removed. In all cases, z_extend is chosen to be 0.5 to 0.6 Mm above the bottom of the simulation, which ensures that the diffusive layer is not included and that the stratification at z_extend is very close to adiabatic at the averaged Γ₁. We then extrapolate Γ₁ linearly downward to a value of 1.6 at −9 Mm and keep it at that value below −9 Mm to correspond roughly to the results of deeper simulations. For the extension, we use a grid with the same vertical spacing as in the simulation. We note that the details of the extension have a negligible effect on the final result, as is demonstrated in Appendix A. Assuming hydrostatic equilibrium, constant gravity, and perfect adiabatic stratification, we can extend the model downward to obtain p_tot and ρ. From that, we calculate

$$\begin{array}{c}\hfill \frac{{\overline{p}}_{\mathrm{turb}}(z)}{{\overline{p}}_{\mathrm{tot}}(z)}=exp\left[\frac{z-z_{\mathrm{extend}}}{3.03\mathrm{Mm}}\right]\frac{{\overline{p}}_{\mathrm{turb}}(z_{\mathrm{extend}})}{{\overline{p}}_{\mathrm{tot}}(z_{\mathrm{extend}})},\end{array}$$(3)

again to roughly match the behavior of the deepest simulation. Finally we set ${\overline{p}}_{\mathrm{gas}}={\overline{p}}_{\mathrm{tot}}-{\overline{p}}_{\mathrm{turb}}$. For convenience, the averaged quantities are indicated by $\stackrel{̅}{}$, even over the extended region, to make it clear as to which quantities are averaged and which are derived.

Given that the turbulent pressure is important for the structure, it is important to ask which Γ₁ should be used in the oscillation equations near the surface. Rosenthal et al. (1999) considered two models. The gas-gamma model (GGM) assumes that the turbulent and gas pressures behave similarly, that is that the effective Γ₁ is

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1^{\mathrm{GGM}}={\overline{\mathrm{\Gamma }}}_1,\end{array}$$(4)

where ${\overline{\mathrm{\Gamma }}}_1$ is the horizontally averaged Γ₁. This is the approach taken by Ball et al. (2016).

The reduced-gamma model (RGM) assumes that

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1^{\mathrm{RGM}}=\frac{{\overline{p}}_{\mathrm{gas}}}{{\overline{p}}_{\mathrm{tot}}}{\overline{\mathrm{\Gamma }}}_1,\end{array}$$(5)

in other words the turbulent pressure is not perturbed by the oscillations. This approach is one of the cases investigated by Houdek et al. (2017), who also consider more elaborate models based on theories of convection. As we have the eigenfunctions for all of the variables, we are able to address the accuracy of these approximations, as we briefly discuss in Sect. 3.2. In the equations below, we refer to a generic Γ₁. In the subsequent calculations, we consider both choices and compare each to the respective results shown in Ball et al. (2016) and Houdek et al. (2017).

To calculate the computed eigenfunctions in the extended MURaM model ℳ_extend, we solve the standard (“classical”) oscillation equations for radial (k_h = 0) oscillations:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\xi }{\mathrm{d}z}=-(\frac{2}{r_{\mathrm{c}}}-\frac{1}{{\mathrm{\Gamma }}_1{H}_{\mathrm{p}}})\xi -\frac{1}{\overline{\rho }{c}^2}{p}^{\prime }\end{array}$$(6)

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}{p}^{\prime }}{\mathrm{d}z}=\overline{\rho }({\omega }^2-N^2+{\omega }_{\mathrm{\Phi }}^2)\xi -\frac{1}{{\mathrm{\Gamma }}_1{H}_{\mathrm{p}}}{p}^{\prime },\end{array}$$(7)

where

$$\begin{array}{c}\hfill H_{\mathrm{p}}^{-1}=-\frac{\mathrm{d}ln\overline{p}}{\mathrm{d}z}\end{array}$$(8)

is the inverse of the pressure scale height, and

$$\begin{array}{c}\hfill N^2=g(\frac{1}{{\mathrm{\Gamma }}_1\overline{p}}\frac{\mathrm{d}\overline{p}}{\mathrm{d}z}-\frac{1}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}z})\end{array}$$(9)

is the square of the buoyancy frequency, r_c is the local radius of curvature,

$$\begin{array}{c}\hfill c^2={\mathrm{\Gamma }}_1\frac{\overline{p}}{\overline{\rho }}\end{array}$$(10)

is the squared sound speed, and we use $\overline{p}$ for the total pressure. As the equations can be written in different forms by substituting variables, it is important to use a minimal number of variables, such that the results are insensitive to the rewrites. We therefore compute N² and c² from the mean pressure, Γ₁, and density, rather than using their averages from the simulations. The quantity ${\omega }_{\mathrm{\Phi }}^2$ is due to the perturbation to the gravitational potential; it is zero for the Cartesian box case (the acceleration due to gravity is constant in this case) and is given by $4\pi G\overline{\rho }$ for spherical models. The symbol G denotes the gravitational constant. In Eq. (6), the radius of curvature r_c is infinity for the Cartesian box case and r_c = r for the spherical case.

To extend the eigenfunctions of the MURaM simulation down through the extended model ℳ_extend, we solve Eqs. (6) and (7) for the frequencies of each of the modes observed in the simulation. As no boundary conditions are imposed, this yields two solutions. We then compute the linear combination of these two solutions that best fits, in the least squares sense, the real part of the observed displacement eigenfunction over the fitting range given by z_fit, min and z_fit, max in Table 1. In other words, we effectively use the observed eigenfunction as an upper boundary condition.

Figure 3 shows examples of the fitted eigenfunctions computed using the above procedure. As can be seen, the fits are very good over the fitting range.

Fig. 3. Examples of the real parts of the numerical eigenfunctions (solid) and fitted eigenfunctions (dashed, barely visible) for case 10 using the GGM approximation. Dotted curves show 100 times the difference. Vertical dotted lines show the fitting range.

We compute the mode frequencies for the ℳ_extend stratification, neglecting the dynamical effects of the convection, by solving the eigenvalue problem of Eqs. (6) and (7). As the bottom boundary condition, we use a zero in the displacement eigenfunction at the location of a, deeply located, zero crossing in the fitted eigenfunction. At the top, we match the eigenfunction to the solution for an isothermal atmosphere. The matching is carried out at 0.5 Mm above the photosphere to roughly match the boundary conditions used by the results we compare to. For the sake of comparison, we chose this boundary condition rather than a boundary condition that is higher in the atmosphere.

As the mean stratification is the same in the two cases, the difference between a mode frequency in the MURaM simulation and in ℳ_extend is the frequency perturbation caused by the impact of the effects that are not taken into account when averaging the background state and patching it onto a 1D model. We denote this frequency shift as δω_extend.

The mode masses in ℳ_extend are:

$$\begin{array}{c}\hfill M_{\mathrm{extend}}={\displaystyle \int \overline{\rho }(z){\xi }^2(z)\mathrm{d}z,}\end{array}$$(11)

where the (real-valued) eigenfunction ξ is normalized to be one at the surface and the integral runs from the chosen zero crossing to the top of the box. In the next section, we show that these mode masses are key to scaling the frequency shifts δω_extend to the solar case.

2.4. Scaling to the Sun

The frequency perturbation δω_extend is, however, not the one expected for the Sun as the model is truncated in depth. We use a simple mode-mass scaling to correct for the truncation. To see the intuitive motivation for our approach, consider the eigenvalue problem for radial modes Fξ = ω²ξ where F is the wave equation operator and ω is the normal-mode frequency. The first-order perturbation theory relationship between a perturbation to the wave-equation operator δF and the resulting perturbation to the square of the mode frequency δω² is

$$\begin{array}{c}\hfill \delta {\omega }^2=\frac{\int \xi \delta F(\xi )\rho r^2\mathrm{d}r}{M_{\mathrm{patch}}}\approx R_{\odot }^2\frac{\int \xi \delta F(\xi )\rho \mathrm{d}r}{M_{\mathrm{patch}}},\end{array}$$(12)

where

$$\begin{array}{c}\hfill M_{\mathrm{patch}}={\displaystyle \int {\xi }^2\rho r^2\mathrm{d}r}\end{array}$$(13)

is the mode mass of the patched (solar) model, which is defined later in this subsection, and where we have assumed that the perturbation is confined close to the surface, allowing us to replace r² with $R_{\odot }^2$ in the integral in the numerator of Eq. (12). For the extended model, we similarly get

$$\begin{array}{c}\hfill \delta {\omega }_{\mathrm{extend}}^2=\frac{\int \xi \delta F(\xi )\rho \mathrm{d}z}{M_{\mathrm{extend}}}·\end{array}$$(14)

As the perturbations we consider (analogous to δF) are confined very close to the surface and well within the extended model, the numerators in Eqs. (12) and (14) should be identical for the Sun and for the MURaM case. From this, it follows that the frequency perturbation in the solar case is given by:

$$\begin{array}{c}\hfill \delta \omega =\frac{R_{\odot }^2{M}_{\mathrm{extend}}}{M_{\mathrm{patch}}}\delta {\omega }_{\mathrm{extend}},\end{array}$$(15)

where we have assumed the approximation δω² ≈ 2ωδω, which is valid for δω ≪ ω. The mode masses M_extend and M_patch have different units. This is because one is a mode-mass for a plane-parallel model (ℳ_extend) and the other is for a spherical model.

It is essential that the eigenfunctions are identically normalized in the two cases in order to use the scaling of Eq. (15). To achieve this, we first patch the outermost part of the averaged MURaM stratification onto Model S (Christensen-Dalsgaard et al. 1996). This is done by shifting the averaged simulation in height, so as to, on average, match the pressure and density at a depth of 1.5 Mm in Model S. Model S and the averaged simulation are then averaged between depths of 1 Mm and 2 Mm by a weight that varies linearly from 100% Model S at 2 Mm to 100% simulation average at 1 Mm, where MURaM and Model S match well. Using a gradual transition, as opposed to a step change, smooths out any remaining discontinuities. We denote this model ℳ_patch. Based on this, we then solve the oscillation equations again (Eqs. (6) and (7)), together with the boundary conditions of regularity at the center of the star and the same top boundary condition that is used for the Cartesian box in order to obtain consistent mode masses. For each resonant mode, we use the eigenfunction that results from this calculation to compute the mode mass M_patch for the model ℳ_patch. We then calculate the scaled frequency perturbation using Eq. (15), having interpolated the mode masses M_patch to the frequency grid of the δω_extend.

An important effect of the mode-mass scaling of the simulation results is that it eliminates the sensitivity linked to how the model is extended below the simulation domain and that the results are independent of which zero crossing is selected, as demonstrated numerically in Appendix A. We also emphasize that the only use of the patched model is to calculate the mode masses, which are used in the scaling. As such, small errors are of little consequence.

3. Results

Figure 4 shows the frequency perturbations computed using Eq. (15) for the various simulations. These frequency perturbations are due to all of the physics not captured in the patched model ℳ_patch. As cases 10, 12, and 13 have the same resolution, the results have been connected for clarity. The error bars were estimated from Eq. (2) of Libbrecht (1992) using the estimated linewidths, zero background noise, and then scaled based on the mode mass as was done for the frequency perturbations. They do not include the contribution from the noise in the measured eigenfunctions. However, the smooth variation of the frequency shift between simulations with different depths indicate that other effects do not contribute significant noise relative to the observed effect. More importantly, they do not include any systematic contributions, such as those from the imperfect match between the observed and the analytical eigenfunctions.

Fig. 4. Frequency perturbations scaled to the Sun, plotted as a function of the simulation frequencies. Squares are for cases 10, 12, and 13, diamonds are for case 11, and crosses are for case 15. The solid lines connect the result of cases 10, 12, and 13. The top curve shows the RGM cases; the bottom curve shows the GGM results. The results for the first zero crossing below 10 Mm were used (as opposed to a shallower crossing) in order to ensure that the physics is well represented by the classic equations. The effects of using different crossings are discussed in Appendix A. For the sake of comparison, the upper dashed black line shows the residuals (Sun minus model) from the RGM patched model of Houdek et al. (2017, model “B”) and the lower dashed black line shows the residuals from the GGM patched model of Ball et al. (2016). Regarding the latter, we note that we were unable to reproduce the frequencies from the model (obtained from Ball, priv. comm.) exactly, due to what appears to be numerical problems in the model. However, the differences were small enough so as to not affect the overall conclusion. The red dashed lines show the results from Fig. 5 of Jørgensen et al. (2018), while the blue lines show the RGM numbers shown in blue in Fig. 6 of Jørgensen & Weiss (2019). The corresponding GGM numbers are courtesy of Jørgensen (priv. comm.). The numbers from Sonoi et al. (2017) are not shown as they have deviations at low frequencies, which are indicative of an overall structural error. The vertical dotted line indicates roughly where leaks from other kh start to appear within the fitting interval, as is discussed in the main text.

As noted earlier, we see that modes with other k_h leak weakly into k_h = 0. To ensure that these do not bias the results, we inspected the spectra of k_h = 1 and $k_{\mathrm{h}}=\sqrt{2}$ (with k_h in units of wavelengths across the box width) to determine if any peaks fall within the fitting interval. As it turns out, none are close to the fitting window below 3400 μHz. Above this frequency, some of the leaks are within fitting window. This is largely due to the rapidly increasing linewidth. We therefore do not trust the results above 3400 μHz, as indicated in Fig. 4.

As part of the procedure, a number of choices were made, including which zero crossing to use as the lower boundary condition, the depth from which to extend the model analytically (given by z_extend), Γ₁ in the extension, and the fitting interval (given by z_fit, min and z_fit, max). While these would not be expected to have a large effect if chosen reasonably, we nonetheless study the effect of changing various parameters in Appendix A. As expected, most parameters have a negligible effect, except for the fitting range where, not surprisingly, small changes are seen. However, these are small relative to the remaining discrepancies when not close to the photosphere, so they are of little importance for the present study.

The results from case 11 agree well with those of case 10, while those of simulation 15, which has poorer resolution, show significant differences. In other words, it appears that the resolution of case 10 (and thus 12 and 13) is adequate for the purpose of this study.

3.1. Comparison with previously published residuals

Figure 4 also shows the frequency residuals (Sun minus model) for the patched model with turbulent pressure and the RGM approximation from Houdek et al. (2017, model “B”, dashed line in their Fig. 1). In the middle of the p-mode band (around 3000 μHz), we explain almost all of their residuals with our RGM calculation. At lower frequencies, we clearly overestimate the residuals, and we underestimate them at high frequencies. Given the results in the next subsection, we did not compare our model with the more elaborate model used by Houdek et al. (2017) directly. For a comparison, the RGM numbers of Jørgensen et al. (2018) are somewhat smaller, while the numbers from Jørgensen & Weiss (2019) are a bit larger than those of Houdek et al. (2017).

As another test of the model presented here, Fig. 4 shows the residuals from the GGM patched model from Ball et al. (2016, Fig. 1, open circles). Around 3000 μHz, we explain almost all of these residuals. At low frequencies, where the residuals are small and the errors are relatively large, there is no clear improvement. At frequencies above about 3400 μHz, our results provide a poor match; however, as noted earlier, this may be due to leaks from other modes. For a comparison, the GGM numbers of Jørgensen et al. (2018) are somewhat smaller, while the GGM numbers of Jørgensen & Weiss (2019) are larger than the numbers from Ball et al. (2016).

3.2. Comparison with Houdek et al. (2017) eigenfunctions

As discussed earlier, we can study the approximations suggested by Houdek et al. (2017) since we have access to the eigenfunctions of the relevant quantities. Figure 5 shows the quantities displayed in Fig. 2 of Houdek et al. (2017) for a mode with a frequency that is close to the mode used for that figure; the frequency is quite close and the quantities do not strongly depend on frequency. The correspondence is rather poor. The phases of the pressure perturbation eigenfunctions agree quite well over the log pressure range between 5.5 and 5.9, but they are very different closer to the surface. Also, the norm of the pressure perturbation eigenfunction has a substantially different behavior from what Houdek et al. (2017) predict. We note that the curves agree better if they are shifted horizontally (i.e., shifted in depth). This may be due to the fact that the background state used is somewhat inaccurate (G. Houdek, priv. comm.).

Fig. 5. Quantities from Fig. 2 of Houdek et al. (2017) for the mode at 2815 μHz from case 11. The dashed line shows the phase difference between the Lagrangian density perturbation and the Lagrangian turbulent pressure perturbation. The dash-dotted line shows the phase difference relative to the Lagrangian gas pressure perturbation. The solid line shows the norm of the ratio of the Lagrangian turbulent pressure perturbation to the total pressure, which was multiplied by an arbitrary factor to match the range of the other quantities. The dotted lines show the same three quantities from Houdek et al. (2017).

Given this situation, it is perhaps somewhat surprising that the frequency perturbations calculated by Houdek et al. (2017) are quite good. The discrepancy suggests that the results of applying the theory to stars of a different spectral type should be used with caution.

4. Discussion

As we have demonstrated, our approach allowed us to determine the surface term quite well. In particular, we were able to estimate the surface term without relying on the accuracy of the GGM or RGM approximations, both of which lack a good physical justification.

Having said that, the fits are not perfect, and this is particularly the case at high frequencies. There are several potential reasons for this. While the errors on the observed solar frequencies are negligible relative to the residuals over most of the frequency range, the effects of line asymmetry may provide some contribution. The residuals shown in Fig. 4 are based on the differences to a solar model and may thus contain model errors. While they are difficult to assess, these errors are also unlikely to explain the discrepancies.

While we implemented GGM and RGM corrections, as used by Ball et al. (2016), Houdek et al. (2017), and others, and while Ball et al. (2016) essentially used the same MURaM models as those employed for the purposes of this study, we did not use the exact same models and we did not implement the exact same boundary conditions. While these differences are expected to be small, they are both located very close to the surface and the resulting frequency changes will thus likely appear as a surface term. Given this they may explain all or much of the discrepancies, especially at high frequencies.

Finally, it is possible that either our models or those we compare to have inaccuracies in the physics. Possible problems include incorrect opacities, equation of state, or approximations in the radiative transfer.

While our results are quite promising and show that the suggested approach to modeling the surface term works, there are many possibilities for improvements, ranging from technical ones to a better understanding of the underlying physics. On the technical side, the lack of a dense grid of frequencies is clearly a problem. Running simulations with different depth ranges, as done here, is one way to mitigate this problem. Unfortunately, this is very computationally expensive. More fundamentally, it involves extending the simulations to depths where the physics is simple and well described by the classic equations and is thus wasteful.

An alternative approach would be to change the bottom boundary condition or to put a partly reflecting layer inside the simulation. Both of these options would shift the resonance frequencies of the computational domain. These solutions would, of course, need to be carried out with extreme care in order to not affect the background state by applying a change that only affects the dynamics at high frequencies, for example. A related technical problem is the long simulation times needed to reduce the effect of the stochastic excitation and rather large damping on our ability to measure the frequencies.

An intriguing way around these problems would be to add a source of excitation of the acoustic waves well below the part of the simulation with nonstandard physics. In this way, one or more well defined frequencies could be excited, possibly at a large amplitude and covering the desired range of frequencies. This approach would eliminate the problem of estimating the frequencies. Unfortunately, simply driving at a frequency that is far from a resonance is, at best, problematic; the presence of nodes in the eigenfunctions of the perturbed quantity requires careful placement of the excitation and the simulations still have naturally excited modes. A way to circumvent these problems would be to add a damping layer near the bottom of the simulations. Even if imperfect, this would allow for driving at any desired frequency. Of course, the concerns regarding a possible modification of the background need to be kept in mind. Also, it is important to realize that the stochastic nature of the convection still affects the interaction of the modes with the convection.

Another thing to study would be nonradial (k_h ≠ 0) modes. However, for the global modes used for inversions of solar-like oscillations in Sun-like stars, the results are not expected to be significantly different.

While we have chosen to divide the procedure into two parts, one to estimate the frequency perturbation in a shallow box and one to scale it to the Sun, it may be possible to combine these by actually patching the eigenfunctions or by deriving an effective boundary condition at some suitable point, for instance. While this may seem attractive and it avoids the mixing of the patching approach and the perturbation approach, there are significant complications. If the matching point is chosen close to the surface, the classic equations do not work well there. If chosen deeper down, curvature and variable gravity effects have to be considered. The approach we have taken carefully avoids these issues, at the price of, perhaps, being conceptually more complicated.

Toward the more theoretical side, studying the numerical eigenfunctions may allow us to obtain significant insights, as already shown in Sect. 3.2. A physical understanding of the wave-convection interaction, perhaps along the lines of Houdek et al. (2017), would clearly be highly desirable. However, this would also probably be difficult. The approach of Hanasoge et al. (2013) is another promising possibility.

A better understanding of the interactions is also important for other problems in helio- and asteroseismology, such as the center-to-limb effects present in helioseismology and the visibilities of oscillation modes in other stars (Schou 2018). Indeed the former was the original motivation for Baldner & Schou (2012). Having said that, the most important application is likely to be in asteroseismology, where the need to fit for the surface term or make questionable assumptions is currently impeding progress.

5. Conclusion

We have shown that it is possible to estimate the surface term to within about 1 μHz in the middle of the p-mode band using hydrodynamic simulations.

However, it is also clear that further work is needed before the results can be applied in helio- and asteroseismic inversions. In particular, the remaining discrepancies need to be understood, a way needs to be found to obtain a denser frequency grid, and computations have to be made for a range of stars.

Acknowledgments

The authors would like to thank Robert Cameron for assistance with running the MURaM simulations and for general discussions, Warrick Ball for help with understanding details of Ball et al. (2016) and for providing values from plots, Günter Houdek for explaining details of Houdek et al. (2017) and various physics issues, and Regner Trampedach for many useful discussions. The authors also thank Andreas Christ Sølvsten Jørgensen and Jakob Mosumgaard for providing us with the data shown in Fig. 4 and for many useful discussions. We acknowledge partial support from the European Research Council Synergy Grant WHOLE SUN #810218.

References

Appendix A: Tests of various approximations

In the main part of this paper, we have made a number of approximations and somewhat arbitrary choices (such as the matching depth). While the approximations would be expected to be reasonable and the overall approach should be insensitive to the choices, it is nonetheless instructive to examine some of them in detail, and we do so in the following subsections. We note that the smooth variation of the scaled perturbations with the simulation depth already validates some of the choices, and we have already discussed the effect of the spatial resolution. For the tests, we used the GGM version of case 12 as this has more frequencies and can be truncated to correspond to cases 10 and 13.

A.1. Zero crossing choice

Table A.1 shows the results of using different zero crossings. As expected, the unscaled frequency changes vary substantially, but this is compensated by the change in the mode mass, such that the resulting differences are negligible, even when using crossings that are very close to the surface.

A.2. Model extrapolation

Another somewhat arbitrary choice is how to extrapolate the models downward. Based on the discussion in Sect. 2.4, the deeper stratification should not affect the final results as the resulting changes in the resonant frequencies would be compensated for by the corresponding changes in mode masses. The most obvious choice we made was how to extrapolate downward. In the main part of the paper, we chose to linearly increase Γ₁ with depth to roughly match deeper simulations. We then assumed adiabatic stratification, constant gravity, and hydrostatic equilibrium to obtain the structure. As an alternative, we kept Γ₁ constant, resulting in a large change in the sound speed. The results are given in Table A.2. While the positions of the zero crossings change significantly, the resulting frequency changes are negligible, giving us further confidence in the insensitivity to the deeper structure.

A.3. Fitting interval

As shown in Fig. 3, the eigenfunctions from the simulation and the analytical ones do not match perfectly. Given this, some effect of the fitting range may be expected; therefore, we explore the effect in Table A.3. As expected, there is some variation, especially when fitting close to the surface, but over the range used for the different cases, the changes are small relative to the remaining discrepancies.

All Tables

All Figures

Fig. 1. Power spectrum of the pressure perturbation for case 10, averaged over the depth range used for the SVD. The power was smoothed to reduce the effect of the stochastic excitation. The large peaks at 1856 μHz, 2812 μHz, and 3858 μHz represent the modes with n = 0, 1, and 2, respectively, where n is the number of nodes in the displacement eigenfunction. The n = 3 peak around 5000 μHz is not analyzed further, as it is very close to the acoustic cutoff frequency. The small peaks around 2100 μHz and 2300 μHz are leaks from modes with kh ≠ 0.

In the text

Fig. 2. Displacement eigenfunctions corresponding to the three largest peaks in Fig. 1. The eigenfunctions were multiplied by the square root of the background density to improve the visibility and as the square of this quantity is proportional to the energy density. Colors indicate the different modes. Solid lines indicate the real part, and dashed lines illustrate the imaginary part. For ease of display, the scaled eigenfunctions were normalized to unity at the bottom boundary.

In the text

	Fig. 3. Examples of the real parts of the numerical eigenfunctions (solid) and fitted eigenfunctions (dashed, barely visible) for case 10 using the GGM approximation. Dotted curves show 100 times the difference. Vertical dotted lines show the fitting range.
In the text

Fig. 4. Frequency perturbations scaled to the Sun, plotted as a function of the simulation frequencies. Squares are for cases 10, 12, and 13, diamonds are for case 11, and crosses are for case 15. The solid lines connect the result of cases 10, 12, and 13. The top curve shows the RGM cases; the bottom curve shows the GGM results. The results for the first zero crossing below 10 Mm were used (as opposed to a shallower crossing) in order to ensure that the physics is well represented by the classic equations. The effects of using different crossings are discussed in Appendix A. For the sake of comparison, the upper dashed black line shows the residuals (Sun minus model) from the RGM patched model of Houdek et al. (2017, model “B”) and the lower dashed black line shows the residuals from the GGM patched model of Ball et al. (2016). Regarding the latter, we note that we were unable to reproduce the frequencies from the model (obtained from Ball, priv. comm.) exactly, due to what appears to be numerical problems in the model. However, the differences were small enough so as to not affect the overall conclusion. The red dashed lines show the results from Fig. 5 of Jørgensen et al. (2018), while the blue lines show the RGM numbers shown in blue in Fig. 6 of Jørgensen & Weiss (2019). The corresponding GGM numbers are courtesy of Jørgensen (priv. comm.). The numbers from Sonoi et al. (2017) are not shown as they have deviations at low frequencies, which are indicative of an overall structural error. The vertical dotted line indicates roughly where leaks from other kh start to appear within the fitting interval, as is discussed in the main text.

In the text

Fig. 5. Quantities from Fig. 2 of Houdek et al. (2017) for the mode at 2815 μHz from case 11. The dashed line shows the phase difference between the Lagrangian density perturbation and the Lagrangian turbulent pressure perturbation. The dash-dotted line shows the phase difference relative to the Lagrangian gas pressure perturbation. The solid line shows the norm of the ratio of the Lagrangian turbulent pressure perturbation to the total pressure, which was multiplied by an arbitrary factor to match the range of the other quantities. The dotted lines show the same three quantities from Houdek et al. (2017).

In the text