Estimating the nonstructural component of the helioseismic surface term using hydrodynamic simulations

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 near-surface layers. Our aim is to develop a method to estimate the effect of the near-surface layers on oscillation frequencies. 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 near-surface 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. 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 density-pressure phase differences computed here do not match those of that work.


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), Ball et al. (2016), and Sonoi et al. (2017), 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 Section 2 we discuss how to determine the frequency perturbations and apply them to the solar case. In Section 3 various results are presented, in Section 4 we discuss the results, and in Section 5 we conclude.
Apart from the discussion in Section 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.

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

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 entropy 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 27400 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 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 Table 1. List of simulations used. The quantity ∆t save indicates the average spacing between the saved (analyzed) time steps. As simulations that are started from a horizontally uniform state take time to develop convection and to reach a statistically steady state, a section at the beginning of each simulation was discarded based on a visual inspection. These sections are not included in the length shown in the table. We note that T eff was calculated from the average of the energy flux over the last 1000 points of the simulation.  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 Γ 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.

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 horizontally-averaged quantities are averaged in time to obtain the mean thermodynamic quantities. We denote the resulting quantities with an overbar, for example p 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 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.
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 h = √ 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.
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 n th resonance, the Fourier transform of a physical variable q has the form where q ′ n 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 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 ξ = −v ′ z /(iω 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 horizontallyaveraged vertical flow near the surface.  and  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.

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 Sec. 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 M 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 p gas = p tot − p turb . 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 gas-gamma model (GGM) assumes that the turbulent and gas pressures behave similarly, that is that the effective Γ 1 is where Γ 1 is the horizontally averaged Γ 1 . This is the approach taken by Ball et al. (2016).
The reduced-gamma 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 Section 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 M 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 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 2 and c 2 from the mean pressure, Γ 1 , and density, rather than using their averages from the simulations. The quantity ω 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πGρ for spherical models. The symbol G denotes the gravitational constant. In equation (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 M extend , we solve equations (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.
We compute the mode frequencies for the M extend stratification, neglecting the dynamical effects of the convection, by solving the eigenvalue problem of equations (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 M 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 M extend are: 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.

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ξ = ω 2 ξ where F is the wave equation operator and ω is the normal-mode frequency. The firstorder perturbation theory relationship between a perturbation to the wave-equation operator δF and the resulting perturbation to the square of the mode frequency δω 2 is where M patch = ξ 2 ρr 2 dr (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 2 with R 2 ⊙ 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 mode-mass for a planeparallel model (M 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 equation (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 M patch . Based on this, we then solve the oscillation equations again (eq. 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 the model M patch . We then calculate the scaled frequency perturbation using equation (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. Figure 4 shows the frequency perturbations computed using equation (15) for the various simulations. These frequency perturbations are due to all of the physics not captured in the patched model M patch . As cases 10, 12, and 13 have the same resolution, the results have been connected for clarity. The error bars were estimated from equation (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.

Results
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 h = √ 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 ), Γ 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.  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).

Comparison with previously published residuals
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 im-  Ball et al. (2016). Regarding the latter, we note that we were unable to reproduce the frequencies from the model (obtained from Ball, private communications) 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

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  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).
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, private communications). 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.

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

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.

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