Issue 
A&A
Volume 638, June 2020



Article Number  A51  
Number of page(s)  9  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/201936530  
Published online  10 June 2020 
Estimating the nonstructural component of the helioseismic surface term using hydrodynamic simulations
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: schou@mps.mpg.de
Received:
20
August
2019
Accepted:
11
April
2020
Context. As the amount of asteroseismic data available continues to grow, the inability to accurately model observed oscillation frequencies is becoming a critical problem for interpreting these frequencies. A major component of this problem is the modeling of the nearsurface layers.
Aims. Our aim is to develop a method to estimate the effect of the nearsurface layers on oscillation frequencies.
Methods. In the proposed method we numerically estimate eigenfunctions in 3D hydrodynamic simulations. We match those to the eigenfunctions calculated from the classic equations applied to the horizontal averages of the structure variables. We use this procedure to calculate the frequency perturbation resulting from the dynamical part of the interaction of the oscillations with nearsurface convection. As the last step we scale the numbers to the Sun. To provide a qualitative test of our method we performed a series of simulations, calculated the perturbations using our procedure, and compared them to previously reported residuals relative to solar models.
Results. We find that we can largely reproduce the observed frequency residuals without resorting to poorly justified theoretical models. We find that, while the calculations of Houdek et al. (2017, MNRAS, 464, L124) produce similar frequency perturbations, the densitypressure phase differences computed here do not match those of that work.
Key words: Sun: helioseismology / asteroseismology / stars: oscillations / convection / waves
© J. Schou and A. C. Birch 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
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., ChristensenDalsgaard & 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 smallscale turbulence causes an apparent net upflow near the solar surface and that this may explain a number of centertolimb 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^{−1}. 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 signaltonoise 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.
Simulations used.
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 Γ_{1} 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^{−7} g cm^{−3}. 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}, ρ, Γ_{1}, 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 horizontallyaveraged quantities are averaged in time to obtain the mean thermodynamic quantities. We denote the resulting quantities with an overbar, for example is the horizontally and timeaveraged gas pressure.
The residuals with respect to the time averages are due to both the radial acoustic oscillations and the timedependent 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 (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 k_{h} ≠ 0. 

Open with DEXTER 
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
where 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 and f_{n} are complex. The singular value decomposition (SVD) decomposes a signal into a sum of such terms:
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(ω)^{2}, using a maximum likelihood fit (Anderson et al. 1990). Given the mode frequency, we also calculate the displacement eigenfunction .
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 horizontallyaveraged 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. 

Open with DEXTER 
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 nearsurface convection for an extended MURaM model. In the second step, which is described in the next subsection, we use a modemass 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 Γ_{1}. We then extrapolate Γ_{1} 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
again to roughly match the behavior of the deepest simulation. Finally we set . For convenience, the averaged quantities are indicated by , 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 Γ_{1} should be used in the oscillation equations near the surface. Rosenthal et al. (1999) considered two models. The gasgamma model (GGM) assumes that the turbulent and gas pressures behave similarly, that is that the effective Γ_{1} is
where is the horizontally averaged Γ_{1}. This is the approach taken by Ball et al. (2016).
The reducedgamma model (RGM) assumes that
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 Γ_{1}. 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:
and
where
is the inverse of the pressure scale height, and
is the square of the buoyancy frequency, r_{c} is the local radius of curvature,
is the squared sound speed, and we use 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^{2} and c^{2} from the mean pressure, Γ_{1}, and density, rather than using their averages from the simulations. The quantity 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 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. 

Open with DEXTER 
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:
where the (realvalued) 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 modemass scaling to correct for the truncation. To see the intuitive motivation for our approach, consider the eigenvalue problem for radial modes Fξ = ω^{2}ξ where F is the wave equation operator and ω is the normalmode frequency. The firstorder perturbation theory relationship between a perturbation to the waveequation operator δF and the resulting perturbation to the square of the mode frequency δω^{2} is
where
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^{2} with in the integral in the numerator of Eq. (12). For the extended model, we similarly get
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:
where we have assumed the approximation δω^{2} ≈ 2ωδω, which is valid for δω ≪ ω. The mode masses M_{extend} and M_{patch} have different units. This is because one is a modemass for a planeparallel 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 (ChristensenDalsgaard 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 modemass 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 k_{h} start to appear within the fitting interval, as is discussed in the main text. 

Open with DEXTER 
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 (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}), Γ_{1} 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 pmode 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 dashdotted 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). 

Open with DEXTER 
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 solarlike oscillations in Sunlike 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 waveconvection 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 centertolimb 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 pmode 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
 Anderson, E. R., Duvall, T. L., Jr., & Jefferies, S. M. 1990, ApJ, 364, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., Serenelli, A. M., & Basu, S. 2005, ApJ, 621, L85 [NASA ADS] [CrossRef] [Google Scholar]
 Baldner, C. S., & Schou, J. 2012, ApJ, 760, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Ball, W. H., & Gizon, L. 2014, A&A, 568, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ball, W. H., Beeck, B., Cameron, R. H., & Gizon, L. 2016, A&A, 592, A159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Balmforth, N. J. 1992, MNRAS, 255, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013, A&A, 558, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, T. M. 1984, Science, 226, 687 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 ChristensenDalsgaard, J., & Thompson, M. J. 1997, MNRAS, 284, 527 [NASA ADS] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hanasoge, S. M., Gizon, L., & Bal, G. 2013, ApJ, 773, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Houdek, G., Trampedach, R., Aarslev, M. J., & ChristensenDalsgaard, J. 2017, MNRAS, 464, L124 [NASA ADS] [CrossRef] [Google Scholar]
 Jørgensen, A. C. S., & Weiss, A. 2019, MNRAS, 488, 3463 [CrossRef] [Google Scholar]
 Jørgensen, A. C. S., Mosumgaard, J. R., Weiss, A., Silva Aguirre, V., & ChristensenDalsgaard, J. 2018, MNRAS, 481, L35 [CrossRef] [Google Scholar]
 Kjeldsen, H., Bedding, T. R., & ChristensenDalsgaard, J. 2008, ApJ, 683, L175 [NASA ADS] [CrossRef] [Google Scholar]
 Libbrecht, K. G. 1992, ApJ, 387, 712 [NASA ADS] [CrossRef] [Google Scholar]
 Mosumgaard, J. R., Jørgensen, A. C. S., Weiss, A., Silva Aguirre, V., & ChristensenDalsgaard, J. 2020, MNRAS, 491, 1160 [CrossRef] [Google Scholar]
 Nordlund, Å., & Stein, R. F. 2001, ApJ, 546, 576 [NASA ADS] [CrossRef] [Google Scholar]
 Piau, L., Collet, R., Stein, R. F., et al. 2014, MNRAS, 437, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2017, ApJ, 834, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689 [NASA ADS] [Google Scholar]
 Schou, J. 2018, A&A, 617, A111 [CrossRef] [EDP Sciences] [Google Scholar]
 Sonoi, T., Samadi, R., Belkacem, K., et al. 2015, A&A, 583, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sonoi, T., Belkacem, K., Dupret, M. A., et al. 2017, A&A, 600, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2001, ApJ, 546, 585 [NASA ADS] [CrossRef] [Google Scholar]
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.
Unscaled and scaled frequency changes resulting from different choices for which zero crossing, z_{zero}, was used to define the lower boundary condition.
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 Γ_{1} 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 Γ_{1} 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.
Unscaled and scaled frequency changes for the standard Γ_{1} extrapolation (top half) and constant Γ_{1} (bottom half).
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.
Unscaled and scaled frequency changes for different choices of where to match the simulation and analytical eigenfunctions. In all cases z_{fit, max} is 1 Mm above z_{fit, min}.
All Tables
Unscaled and scaled frequency changes resulting from different choices for which zero crossing, z_{zero}, was used to define the lower boundary condition.
Unscaled and scaled frequency changes for the standard Γ_{1} extrapolation (top half) and constant Γ_{1} (bottom half).
Unscaled and scaled frequency changes for different choices of where to match the simulation and analytical eigenfunctions. In all cases z_{fit, max} is 1 Mm above z_{fit, min}.
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 k_{h} ≠ 0. 

Open with DEXTER  
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. 

Open with DEXTER  
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. 

Open with DEXTER  
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 k_{h} start to appear within the fitting interval, as is discussed in the main text. 

Open with DEXTER  
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 dashdotted 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). 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.