Influence of metallicity on the near-surface effect on oscillation frequencies

The CoRoT and Kepler missions have provided high-quality measurements of the frequency spectra of solar-like pulsators, enabling us to probe stellar interiors with a very high degree of accuracy by comparing the observed and modeled frequencies. However, the frequencies computed with 1D models suffer from systematic errors related to the poor modeling of the uppermost layers of stars. These biases are what is commonly named the near surface effect. The dominant effect is related to the turbulent pressure that modifies the hydrostatic equilibrium and thus the frequencies. This has already been investigated using grids of 3D RMHD simulations, which also were used to constrain the parameters of the empirical correction models. However, the effect of metallicity has not been considered so far. We study the impact of metallicity on the surface effect across the HR diagram, and provide a method for accounting for it when using the empirical correction models. We computed a grid of patched 1D stellar models with the stellar evolution code CESTAM in which poorly modeled surface layers have been replaced by averaged stratification computed with the 3D RMHD code CO5BOLD. We found that metallicity has a strong impact on the surface effect: keeping T_eff and log g constant, the frequency residuals can vary by up to a factor two. Therefore, the influence of metallicity cannot be neglected. We found that a correct way of accounting for it is to consider the surface Rosseland mean opacity. It allowed us to give a physically-grounded justification as well as a scaling relation for the frequency differences at nu_max as a function of T_eff, log g and kappa. Finally, we provide prescriptions for the fitting parameters of the correction functions. We show that the impact of metallicity through the Rosseland mean opacity must be taken into account when studying and correcting the surface effect.


Introduction
Thespace-bornemissionsCoRoT (Baglin et al.2006;Michel et al. 2008;Auvergne et al. 2009) and Kepler (Borucki et al. 2010) have provided a rich harvest of high-quality seismic data for solar-like pulsators. This has allowed a leap forward in our understanding and modelling of low-mass stars (see the reviews by Chaplin & Miglio 2013;Hekker & Christensen-Dalsgaard 2017). However, for the last three decades (e.g. Dziembowski et al. 1988) it has been known that the comparison between modelled and observed acoustic-mode frequencies suffer from systematic discrepancies. This bias is called the surface effect and has been widely studied in the solar case (Rosenthal et al. 1995;Christensen-Dalsgaard & Thompson 1997;Rosenthal & Christensen-Dalsgaard 1999). They are attributed to our deficient modelling of the uppermost layers of stars with a convective envelope. Indeed, 1D stellar models hardly take into account the complexity of these layers that are subject Present address: Institut d'Astrophysique Spatiale, Université Paris Sud, Orsay, France to highly turbulent flows as well as a complex transition between a convective to a radiative-dominated energy flux transport (e.g. Kupka & Muthsam 2017).
More generally, these frequency residuals prevent a direct comparison between modelled and observed frequencies. Frequency combinations are commonly used to circumvent this problem (e.g. Roxburgh & Vorontsov 2003), but still, an accurate determination of frequencies is highly desirable to take advantage of the full potential of asteroseismology. To reach this goal, a handful of empirical prescriptions with adjustable free parameters have been proposed (Kjeldsen et al. 2008;Ball & Gizon 2014Sonoi et al. 2015) and allow one to apply a posteriori corrections to the modelled frequencies. Such an approach is now widely used (e.g. Silva Aguirre et al. 2017) and has proven to be quite efficient in inferring a stellar model that fits the observed frequencies. However, it suffers from some fundamental drawbacks. The choice of the parameters is not physically motivated. Consequently, there is no guarantee that this optimal model is unique and accurate (i.e. that it properly reproduces the real physical structure of the observed star).
Another complementary approach then consists of investigating the physical nature of the surface effect. This motivated a number of studies to unveil and constrain the physical ingredients responsible for these biases. More precisely, surface effect has been shown to be the result of two distinct effects (e.g. Houdek et al. 2017): structural effects coming mainly from turbulent pressure in the hydrostatic equation which is usually absent in 1D stellar evolution codes, and modal effects gathering modifications of the eigenmodes, mostly due to nonadiabaticity (e.g. Balmforth 1992;Houdek et al. 2017) as well as the perturbation of turbulent pressure induced by the oscillations (Sonoi et al. 2017). Other related processes were also invoked, such as convective backwarming  or magnetic activity (Piau et al. 2014;Magic & Weiss 2016). Nonetheless, as demonstrated by the early work by Rosenthal & Christensen-Dalsgaard (1999) on the Sun using a 3D hydrodynamical simulation, the dominant physical ingredient is thought to be the turbulent pressure that modifies the hydrostatic equilibrium and subsequently introduces an elevation of the star surface. Then, the acoustic cavity is modified and therefore the frequencies are as well.
Based on a grid of 3D numerical simulations, this method was used by Sonoi et al. (2015), Ball et al. (2016),  who investigated the surface effect variations across the Hertzsprung-Russell diagram. These works clearly demonstrated that surface effects sharply depend on effective temperature and surface gravity of star. In addition, Sonoi et al. (2015) presented a way to provide parameters for the empirical surface corrections by fitting them against a physically motivated scaling relation derived by Samadi et al. (2013). However, all these works considered solar metallicity models while the distribution of metallicity for observed solar-like pulsators is quite large (see e.g. Pinsonneault et al. 2014). Our goal is thus to study the influence of metallicity on the surface effects and propose a method to account for it.
The article is organized as follows: in Sect. 2 we describe the method of model patching, which is constructed by replacing the upper layers of a 1D model by horizontally averaged stratification of a 3D model atmosphere, and our set of models. Then in Sect. 3 we show that metallicity has a strong impact on the frequency residuals and therefore its influence cannot be ignored. We also study the variation of the frequency differences with effective temperature, surface gravity and opacity and give a theoretical justification. Finally, in Sect. 4 we provide constraints on the various parameters usually used in the empirical surface effect function across the T eff − log g − log κ space.

Model-patching method
In this section we explain the method we used to patch our models and describe our final set of models.

Grid of 3D models
We used a grid of 3D hydrodynamical models from the CIFIST grid of stellar atmosphere including the superadiabatic region to the shallowest layers of the photosphere, computed using the CO 5 BOLD code (see Ludwig et al. 2009;Freytag et al. 2012 for details). The chemical mixture is based on the solar abundances of Grevesse & Sauval (1998) apart from the CNO elements which follow Asplund et al. (2005). We considered 29 models with effective temperature (T eff ) ranging from 4500 K to 6800 K, a surface gravity (log g) ranging from 3.5 to 4.5, and a metallicity [Fe/H] = −1.0, −0.5, +0.0, +0.5. We note that 5000 5500 6000 6500 T  Table 1 summarizes the global parameters of the 3D models. The range of metallicities we considered corresponds to the metallicities of observed solar-like pulsators (Anders et al. 2017;Serenelli et al. 2017). Table 1 exhibits small groups of models (labelled with same first letter) with very similar T eff and log g. Those groups for instance in Fig. 1 at log g = 4.0, help us to investigate the influence of metallicity on the surface effect by keeping other global parameters fixed. However, we pointed out that, whereas within a group the dispersion in log g is rather small (of the order of 0.1%), the dispersion in T eff is much higher (of the order of 1%). Indeed, surface gravity is an input parameter of the hydrodynamical simulations while effective temperature is controlled by adjusting the entropy at the bottom of 3D models. It is therefore difficult to match an accurate effective temperature.

Computation of patched models
For each 3D model, both a patched (hereafter PM) and an unpatched model (hereafter UPM) have been constructed. A patched model is a model computed using a 1D stellar evolutionary code in which we replaced the surface layers with the stratification obtained by horizontally averaging a 3D model computed with a R-MHD code. The fully 1D model is called an unpatched model. The construction of PM and UPM has been widely described in Trampedach (1997), Samadi et al. (2007Samadi et al. ( , 2008, Sonoi et al. (2015), Jørgensen et al. (2017). The 1D counterparts of 3D hydrodynamical models have been obtained using the 1D stellar evolutionary code CESTAM (Morel 1997;Marques et al. 2013) by tuning the age (or the central temperature for advanced stages), the total stellar mass M, and the mixing length parameter α MLT in order to match the effective temperature, the surface gravity and the temperature at the bottom of the 3D model, located just below the superadiabatic region. We chose to remove the first four bottom layers and the last top layer of the 3D hydrodynamical model to be sure to remove any numerically induced errors and that the patching point is deeply inside the Model , respectively. T b is the mean temperature at the bottom of the 3D model and ν max is the frequency with the largest amplitude in the oscillation power spectrum (ν max = 3050.0(M/M )(R 2 /R 2 )(5777/T eff ) 1/2 , Kjeldsen & Bedding 1995) and M is the mass of the PM which differs by a fraction 10 −7 from the one of UPM. The initial helium and metal abundances are close to the ones at the surface. We recall that the metal abundance is different from the iron abundance [Fe/H] imposed in our models. The evolutionary stages PMS, MS, and SG stand for pre-main-sequence, main-sequence and sub-giant.
adiabatic region, which has been shown to be a condition for a obtaining reliable PM (Jørgensen et al. 2017). The 1D models use the equation of state, and opacities given by OPAL2005 (Rogers & Nayfonov 2002;Iglesias & Rogers 1996) and implement standard mixing-length theory (Böhm-Vitense 1958) with no overshoot. We ignore diffusion processes, rotation and turbulent pressure. The atmosphere is computed using the Eddington approximation. The helium abundance in 1D models is set to the one used in 3D models.
Finally, we note that for some 3D models one can find a degenerate solution for the corresponding 1D model: we could patch either a PMS or a sub-giant model. We opted for evolved models since they corresponds to stars in which solar-like oscillations are observed so far. However, when the evolved models are too old (older than the age of the Universe) we kept the PMS model, except if lying on the Hayashi track. Table 1 also summarizes the stellar parameters of both UPM and PM together with relative radius differences R PM /R UPM − 1. Our set of models covers a wide portion of the Hertzsprung-Russel Diagram for intermediate mass stars. We note that our patched models with metallicity [Fe/H] = +0.5 only have log g = 4.0. Indeed, 3D models from the CIFIST grid with [Fe/H] = +0.5 were only available for log g ≥ 4.0. In addition, 3D models with log g 4.5 are located below the main sequence diagonal, and therefore it is impossible to find a 1D model matching their characteristics (with the physical ingredients we used). Thus, a large portion of our initial [Fe/H] = +0.5 3D hydrodynamical models were not suitable for our purposes.

Computation of oscillation frequencies
In this work, we consider only structural effects and an adiabatic treatment of the oscillations. The frequencies are computed using the ADIPLS code (Christensen-Dalsgaard 2011) for both UPM and PM by assuming the gas Γ 1 approximation, which assumes that the relative Lagrangian perturbations of gas pressure and turbulent pressure are equal (Rosenthal & Christensen-Dalsgaard 1999;Sonoi et al. 2017). Besides this distinction in the treatment of Γ 1 entering the calculation of the model frequencies, we emphasize that the frequency differences studied in this work are only emerging from structural effects. Therefore, it should be emphasized that the frequency differences studied in this work concern only purely structural effects. We have checked that we recovered the previous results of Sonoi et al. (2015) for the solar metallicity. For the sake of simplicity, we mainly focussed on the surface effect affecting radial modes: non-radial modes exhibit a mixed behaviour that would make our analysis more complex (however, see Sect. 4.2.2 for a discussion).

Influence of metallicity
Until now, surface effects have always been studied assuming a solar metallicity. Corrections depend only on T eff and log g such as the power law proposed by Kjeldsen et al. (2008), cubic and combined inverse-cubic laws (Ball & Gizon 2014), or a modified Lorentzian (Sonoi et al. 2015). This section is intended to motivate the investigation of the dependence of the surface effect on metallicity.

Qualitative influence of metallicity on frequency differences
We begin this section by quickly describing the effects of a change of metallicity on the frequency residuals. Frequency differences are induced by the surface elevation between PM and UPM due to turbulent pressure, which extends the size of the resonant cavities and therefore decreases the mode frequencies for PM, leading to negative frequency differences δν = ν PM − ν UPM . Up to now, only the influence of effective temperature and surface gravity on surface effects have been studied. However, the abundance of heavy elements has a strong impact on opacity and hence on the convective flux imposed by a change in the radiative flux. In turn, a change in the convective flux leads to a change of convective velocity and therefore a change of turbulent pressure and finally it changes the location of the surface. We mention here that metallicity also has an effect on gas pressure, through the mean molecular weight µ, which varies in the opposite direction of the turbulent pressure and therefore counteracts its effect. Finally, while mechanisms by which a change of metallicity can act on the surface effect are known, those mechanisms are too intricate to identify the resulting effect on the variations of surface term without a deeper analysis as will be demonstrated in the following (see Sect. 4). Figure 2 shows the (purely structural) frequency differences for three groups of models that have approximately the same effective temperature and surface gravity. The discrepancies in ν n between two models appear at relatively low frequencies and generally increase towards high frequencies. As for finding a general trend of the evolution of the surface effect against the metallicity, it seems from Fig. 2 no such trend exists: in the top panel, frequency differences, at ν max for instance, slightly decrease from [Fe/H] = −0.5 to 0.0 and then are much higher for the [Fe/H] = +0.5 model. In the middle panel, the frequency residual at ν max significantly increases from [Fe/H] = −0.5 to 0.0. Finally, in the bottom panel, very little variations at ν max can be noticed from one composition to an other. However, the variation of the frequency differences seems to follow closely the variations of the elevation of the stellar surface between UPM and PM: The

Effect of the elevation on the frequency differences
To gain some insight into the influence of metallicity on surface effect, we tried to scale the normalized frequency differences at ν max for our set of models. This is a necessary step to allow an estimate of the surface effect correction parameters (see Sect. 4). Thus, let us start with the perturbative approach as adopted by (Christensen-Dalsgaard & Thompson 1997; see also Goldreich et al. 1991;Balmforth et al. 1996). The authors show that the frequency difference can be well approximated by where c is the adiabatic sound speed, the variable v is defined by v = Γ 1 /c,K n c 2 ,v andK n v,c 2 are the kernels that can be determined from eigenfunctions, δ m c 2 and δ m v are the Lagrangian differences of c 2 and v, respectively, at fixed mass. Rosenthal & Christensen-Dalsgaard (1999) further approximated the frequency differences for radial modes, based on the expression ofK n v,c 2 and using a first-order asymptotic expansion for the eigen-function, by where ∆ν is the asymptotic large frequency separation, ∆r is the previously defined elevation, and c ph the photospheric sound speed (see Appendix A for a demonstration of this relation).
This relation has been previously tested by Sonoi et al. (2015) at solar metallicity using surface effect derived from a grid of 3D numerical simulations. It turns out that Eq. (4) reproduces the overall scale of the surface effect (such as in Fig. 3 were the surface effect is considered at ν max ) for a set of models. It is thus necessary to determine whether this relation holds for models with a non-solar metallicity. To this end, we have compared frequency residuals at ν = ν max given by Eq. (4) as shown in Fig. 3 (top panel). There is still a good agreement between the frequency differences and the approximated expression given by Eq. (4). Moreover, it appears that the frequency differences are dominated by the surface elevation ∆r. To understand the link to metallicity, it is thus necessary to go a step further and to investigate the relation between surface elevation and metallicity.

Scaling law for the frequency differences
In this section, we aim to determine a relation between frequency differences at ν max and global parameters of the models. First, as shown in the previous section, there is no clear trend between the surface term and metallicity. Indeed, at constant metallicity and considering our rather large range of effective temperatures and surface gravities, the dominant opacity mechanisms are not the same from a model to an other for instance, the opacity at the surface is dominated by the negative hydrogen ions for T eff 5000 K. Therefore, the relation between δν/ν and Z is non-trivial. To overcome this problem, we directly consider the Rosseland mean opacity at the photosphere instead of the metallicity as a global parameter in addition to the effective temperature and to the surface gravity (in the following, the photosphere is defined as the radius at which T = T eff ).
Let us begin by considering the elevation in Eq.
(1) which must be expressed as a function of these global parameters. Using the hydrostatic equilibrium equation, it reads where H PM p and H UPM p are the pressure scale heights at the photosphere associated with the patched and unpatched models, p tot is the total pressure such as p tot = p turb + p g with p turb and p g the turbulent and gas pressure, respectively. Further assuming that Finally, since the pressure scale-height scales as T eff /g, the elevation scales as ∆r ∝ (T eff p turb )/(gp g ).
To go further, we need to find an expression for p turb /p g . Near the photosphere, the turbulent pressure can be written as where v conv is the vertical component of the convective velocity. We now need an expression for this velocity and for the density. Assuming a standard Eddington grey atmosphere, the optical depth is approximated by τ = H p ρκ, and in the Eddington approximation, we have τ = 2 /3 at the bottom of the photosphere. Then, and accordingly: As for finding an expression for v conv , we note that F tot = F rad +F conv , with F rad and F conv the radiative and convective component of the total energy flux respectively. The convective flux is proportional to the kinetic energy flux (as shown for instance within the MLT framework). Then, The ratio F rad /F conv is assumed to remain nearly constant from one model to an other. Therefore, v conv finally reads, Inserting the expressions of Eqs. (8) and (10) into Eq. (7) leads to: where κ = 0.415 cm 2 g −1 . From the perfect gas law for, p g ∝ ρT eff and using Eq. (8), we can rewrite ∆r as Replacing ∆r into Eq. (4) one finally obtains the following estimate: This expression provides us with a simple relation between the frequency differences and the global parameters. The dependence on the metallicity is embedded into the Rosseland mean opacity. We note that it is possible to go further and to explicitly introduce the metallicity. For instance, in the vicinity of the solar effective temperature and gravity, the opacity is dominated by the H − so that κ ∝ ρ − 1 /2 T 9 eff Z. However, given the wide range of effective temperatures and surface gravities of our grid of models, it is more relevant to keep the Rosseland mean opacity at T = T eff (surface opacity) as a global parameter. Indeed, the Rosseland mean opacity is a quantity available in any 1D stellar evolutionary code.
Then, using Eq. (13) as a guideline, we performed a fit where the powers of the temperature (p), gravity (q), and opacity (s) have been adjusted at ν = ν max for each model. Figure 3, bottom panel, displays the result. This figure shows a very good agreement between exponents derived in Eq. (13) and the one actually obtained using our simulations. Consequently this scaling can be used to provide a physically-grounded values for the parameters of the empirical correction function of the surface effect. Finally, we note that using the opacity instead of the metallicity allows us to take a detailed mixture into account.
In addition to our crude approximations, a possible source of discrepancies between values predicted by Eq. (13) and the one calculated can be that we did not fix the helium abundance from one model to the other when varying the metallicity. The changing helium abundances have an impact both on the evolution of the model and on its opacity at the surface. However, the helium abundances [He/H] range between −5.8 × 10 −3 and +1.2 × 10 −2 and should be a negligible source of uncertainty. A final source of error comes from the method we used to average the 3D stratifications. Indeed, since the Rosseland opacity is involved Eq. (13), it would be more precise to patch the models using a stratification averaged against the Rosseland optical depth instead of the actual geometrically averaged stratification, but this is beyond the scope of this paper and will be investigated in a forthcoming work.

Surface-effect corrections
A handful of empirical functions have been suggested to perform a posteriori corrections on the modelled frequencies. After having given a theoretical background that explains variations of δν/ν, we considered the most commonly used correction models to study the evolution of the related free parameters as a function of effective temperature, surface gravity, and surface opacity. This is intended to provide constraints on those parameters and thus to provide physically-grounded values for use on seismic observations.

Empirical functions for correcting modelled frequencies
4.1.1. Kjeldsen et al. (2008) power law Kjeldsen et al. (2008) proposed a power law which was found to match the frequency differences obtained between the observed and modelled solar frequencies: where a and b are the parameter to be adjusted. They found a = −4.73 and b = 4.9 for their model of the Sun by matching a subset of nine radial modes centred on ν max . Kjeldsen et al. (2008) provided a method to correct the frequency for a star similar to the Sun without having to calibrate b. Let us assume we want to model a star with near solar global parameters and we want to constrain our model using the individual frequencies. The radial mode frequencies spectrum of our best model which include a surface term are denoted ν i,best and the frequencies of solar radial modes for the same order are denoted ν i,ref . Then, Kjeldsen et al. (2008) proposed that the frequencies can be linked, to a good approximation, by ν i,best rν i,ref , using the proportionality factor r between mean densities of both models:ρ best = r 2ρ ref .
Using this relation and the large separations of both models, they provided a way to obtain a and b. Further assuming b constant (the value of which depends of the physical ingredients used in the model), they derived a value for a for a set of theoretical models close to the Sun.
This power law has been widely used since and many authors (e.g. Metcalfe et al. 2009;Bedding et al. 2010) have used a constant value for b (not necessarily 4.90 though) derived from solar frequency measurements. Keeping b constant is often necessary in the case for which observations do not provide enough constraints to adjust it. However, using the solar value leads to a bad correction if the modelled star is too different from the Sun (e.g. Kallinger et al. 2010). Furthermore, b depends on the input physics. Otherwise, b can be considered as a variable parameter in the modelling and therefore significantly improve the correction. Different models of the star HD 52265 have been compared by  using various input physics and found approximatively the same predicted age models when either frequency ratios (Roxburgh & Vorontsov 2003) or individual corrected frequencies were used as constraints. The age dispersion was slightly higher with models constrained by individual corrected frequencies (∼±9.5%) and using uncorrected individual frequencies lead to ages 40% larger .
In the following, we have studied two versions of this parametric function. The first, adjusted on the whole radial mode frequency spectrum for frequency less than the acoustic cut-off frequency, will be referred to as K08. The second, adjusted on a reduced frequency interval 0 < ν/ν max < 1.05 is refered to as K08r (see Fig. 6 and Appendixes B and C).

Ball & Gizon (2014) cubic and combined inverse-cubic laws
Ball & Gizon (2014) suggested a new function to correct frequency differences. It is partially based on the early work by Gough (1990; for the cubic part). They accounted for two leading effects introducing systematic errors in the theoretical computation of the frequency spectrum: the modification of the sound speed caused by a magnetic field concentrated into a filament by convective motions, causing a frequency shift scaling as ν 3 /E (Libbrecht & Woodard 1990), E being the normalized mode inertia; and the modification of the pressure scale height caused by a poor description of convection, inducing a frequency shift scaling as ν −1 /E. This correction funtional has the advantages of being independent of a solar calibration and including a dependence on the normalized mode inertia which allows us to correct non-radial modes, without the need of re-scaling their frequency differences. Because of this, they suggested a cubic correction taking only into account the dominant effect and a combined inverse-cubic correction including the perturbation. The cubic correction (in the following BG1) is defined by and the combined inverse-cubic correction (in the following BG2) is where E is the normalized mode mass: where R, M, and ρ are respectively the photospheric radius, mass and density of the star, and ξ r and ξ h are the radial and the horizontal component of the displacement of an eigenmode of degree . a 3,BG1 , a −1,BG2 , and a 3,BG2 are the parameters to be adjusted. They used the acoustic cut-off frequency ν c instead of ν max in order to normalize their fitting parameters: it only results in a modification of a −1 and a 3 and does not change the law itself.

Sonoi et al. (2015) modified Lorentzian
The final function to be introduced was a modified Lorentzian (Sonoi et al. 2015) that was found to better correct the surface effect derived from the 3D simulations at high frequency. It reads where α and β parameters are to be determined. When ν PM /ν max 1 we get back to Kjeldsen et al. (2008) law. When ν PM = ν max , δν/ν max = α/2. Therefore, a and α are directly linked to δν/ν given by Eq. (13), which gives physical justification for its variations. In the following, we will refer to this correction law as S15. A recent comparison of the above correction laws has been performed by Ball & Gizon (2017) on six sub-and red giants from the Kepler Input Catalog. We note that since these are evolved stars, they display mixed-modes which have their frequency residuals off the general trend of radial p-modes frequency residuals. This should have consequences on the quality of the correction. They computed stellar models matching their six stars constrained by the effective temperature, the metallicity and the individual frequencies. They tested five correction relations: BG1, BG2, S15, and K08 with power b = 5.0 calibrated with a solar model computed with their own input physics and K08 with power b left free. Logarithms of K08 (first row), S15 (second row) and BG2 (third row; log |a 3,BG1 | behaves almost identically to log |a 3,BG2 |) fitting coefficients as a function of log z 1 (Eq. 13): log f = slope × log z 1 + offset. Purple dots (resp. yellow, green and red) corresponds to [Fe/H] = −1.0 (resp. −0.5, 0.0 and +0.5). Ball & Gizon (2017) found no correction to be clearly superior than the others for all stars. However BG2 and then BG1 performed slightly better than the others, followed by the free power law, S15 and finally K08. S15 was shown to poorly correct high frequencies and the K08 with b = 5.0 correction gave worse results than no correction for their of their stars. We present very similar conclusions in the following.

.1. Prescriptions for radial modes
In order to fit the parameters of the correction functions and to determine the fitting parameters, we used a least square minimization algorithm implementing a Levenberg-Marquardt A107, page 8 of 14 method which minimizes the squared deviation defined as where i corresponds to the eigenmode index, N to the total number of radial modes and δν the correction computed from the considered correction relation. Tables B.1 and C.1 summarize the coefficients and their squared deviations from our computations. We also define √ D/N as the root-mean frequency differences after correction.
Coefficients a and b involved in Eq. (14) and a 3,BG2 from Eq. (16) are presented in T eff − log g plane in Fig. 4. Furthermore, all coefficients are represented as a function of log z 1 in Fig. 5. a, α and a 3,BG2 (and therefore a 3,BG1 ) show similar trends. Indeed, they are related to the amplitude of the surface effect at ν max and Eq. (13) allows one to understand their variations. However, this theoretical justification for the variations of a, α, a 3,BG1 and a 3,BG2 does not provide a way to favour one correction law over an other. Coefficient a −1,BG2 exhibits the same behaviour as above. However, we cannot offer the same explanation for its trend because the inverse term in BG2 is a second correction to the cubic term and is not related to the amplitude of the surface effect at ν max .
The trends followed by b and β are related to the slope of the frequency differences. As shown in Sonoi et al. (2015) and in Fig. 4, the coefficients b (whatever the metallicity) increase significantly towards cooler stars, which again contradicts the assumption of a constant b. Regarding the relevance of giving a prescription for log b and log β thanks to the linear relationship with log z 1 , we can see in Fig. 5 that log b and log β are affected by a high dispersion compared to the grey line. This could mean either that we omitted a physical dependency in Eq. (13) that only affects the agreement with log b and log β, or a prescription based on same other physical basis should be investigated. Table 2 shows the prescriptions for the variations in the T eff − log g − κ space of all coefficients c 0 studied in this article in the form log c 0 = c 1 log ∆ν/∆ν + c 2 log T eff /T eff + c 3 log g/g + c 4 log κ/κ + c 5 .
We note that the opacity has a strong impact on each of the coefficients and must, therefore be taken into account when correcting the surface effect.
The top panel of Fig. 6, top panel also shows the value of the root-mean frequency differences after correction for each model and each correction law. From this, we see that BG1 is the worst performer followed by K08. Those laws provide a correction that leaves frequency residuals comprised between 1 and 10 µHz which are still higher than the frequency resolution provided by CoRoT and Kepler. The better performance of K08 over BG1 can be explained by the fact that K08 have two degrees of freedom whereas BG1 has only one. For radial modes, the inclusion of the normalized mode mass E n in BG1 does not compensate the loss of a degree of freedom.
Then, the remaining laws K08r, S15 and BG2 provide correction almost as good as the resolution of CoRoT and Kepler. K08r and BG2 are slightly better than S15, yet K08r is applied only on the frequency range 0 < ν/ν max < 1.05.

Mixed-modes case
We also performed the same test as in Ball & Gizon (2017) on evolved models that present mixed-modes in their frequency spectrum. In Sect. 4.1.2 we see that, thanks to the dependence in the normalized mode inertia, BG1 and BG2 can be applied to non-radial modes without any change of the law. However, in order to be able to compare all empirical corrections on non-radial modes, one had to rescale the frequency differences on which K08, K08r and S15 by mean of the inertia ratio Q n for a mode of frequency ν n defined as the ratio of the inertia of this mode by the inertia of a radial mode interpolated at the frequency ν n : Q n = E n /E n0 (ν n ) (e.g. Rosenthal & Christensen-Dalsgaard 1999). Furthermore, we added one last empirical relation by modifying the expression given for S15 in Eq. (18) similarly to BG1 and BG2 in which we replaced α by α/E where E is defined in Eq. (17) (the new function is denoted S15E). This allows S15E to be applied directly on non-radial modes frequency differences. The empirical relations K08, K08r and S15 were then adjusted on Q n δν n , with 0 ≤ ≤ 2 and S15E, BG1 and BG2 were adjusted directly on δν n , with 0 ≤ ≤ 2. For 9 of the 16 evolved models considered, the least-square algorithm converge to a solution of K08, K08r or S15 very remote from the general trend of frequency differences, whereas for the second group of relations (S15E, BG1 and BG2), the residual root-mean frequency differences after correction are greatly improved to a value between 0.1 and 1 µHz.
As a third test, we performed the same fits excluding the quadrupolar modes (i.e. we fit modes with 0 ≤ ≤ 1). This time, corrections laws accuracies are similar to the one presented in Sect. 4.2.1. As for the newly introduced S15E, it performs slightly worse than S15 but still better than K08. This third test suggests that the failure of K08, K08r and S15 in fitting Q n δν n , with 0 ≤ ≤ 2 is due to quadripolar modes. There are two reasons for it.
First, the p and g cavity are less coupled for = 2 than for = 1 mixed-modes which induces more important changes on the behaviour of a mode when the surface layers are changed between UPM and PM. Indeed, modifying the surface layers changes the frequency of pure p modes that couple with different g modes for PM and UPM (Ball & Gizon 2017). As a consequence, when computing Q n,2 δν n,2 , it so happens that we deal with mixed-modes from PM and UPM that have different properties. Second, due to the presence of mixed-modes, Q n,2 sometimes becomes higher than ten, while it is normally of the order of unity (see example of Cm10 in Fig. 7). It over-scales the corresponding quadrupolar mixed-modes and gives much weight UC K08 K08r S15 BG1 BG2 Fig. 6. Root-mean frequency differences ν max √ D/N of radial modes after correction for each empirical law and for each model. It shall be noticed that the deviation for K08r is computed only on the range 0 < ν/ν max < 1.05. Black dots corresponds to uncorrected frequencies (UC), that is root mean frequency differences between PM and UPM frequencies.
to those modes, which in turns has strong impact on the quality of the fit. The peaks in the value of Q n,2 arise for mixed-modes having most of their amplitude in the g mode cavity, contrary to dipolar mixed-modes which have their Q n,1 staying close to unity. On the other hand, fitting directly δν n with S15E, BG1 and BG2 does not amplify the frequency differences affecting mixed-modes, providing a much better correction.
For all these reasons, we recommend the use of BG2 or S15E when correcting sets of radial and non-radial modes and we recommend either S15 or BG2 when correcting only radial modes (K08r can be used for frequencies ν max ). The remaining advantages of S15 over BG2 is that the coefficients α and β are better described as a function of Eq. (13) than coefficients of BG2 which present a bigger dispersion on Fig. 5. This being said, many of the correction laws considered in this paper gives a rootmean frequency difference of the order of 0.1 µHz (at least for few models), similar to the frequency uncertainties of CoRoT and Kepler. Furthermore, the search for the best a posteriori correction law should not set aside the need of a theoretical understanding of the surface effect.

Conclusion
We have computed a grid of 29 couples of one dimensional models using the method of patched models consisting in replacing poorly modelled surface layers of a 1D model by the stratification, averaged over the geometrical depth and time, computed from 3D hydrodynamical models. The grid includes models with effective temperature ranging from T eff = 5000 K to 6800 K, surface gravity ranging from log g = 3.5 to 4.5 and iron abundance ranging from [Fe/H] = −1.0 to +0.5.
Our aim was to estimate and understand the impact of varying metallicities on the surface effect. Our main result is that, in the considered range of metallicities (i.e. [Fe/H] = −1.0 to +0.5) the amplitude of the surface effect computed at ν max , and for models with same effective temperature and same surface gravity, can be up to a factor of three between the model with the Each colour corresponds to a degree . Yellow (resp. green) dots breaking from the general trend correspond to dipolar (resp. quadrupolar) mixed modes.
lowest amplitude and the model with the highest one. However, it appears that studying the amplitude as a function of the metallicity does not lead to a clear trend, whereas the Rosseland mean opacity κ turned out to be the adapted quantity for understanding the variation of the surface effect. Based on relatively simple physical arguments, consolidated using the grid of 3D models, we found a scaling relation between the amplitude of the surface effect and the global parameters T eff , log g and the opacity κ computed at the photosphere.
We also tested the accuracy of existing surface effect empirical corrections of radial modes frequency differences on each model of our grid in order to obtain a prescription for the A107, page 10 of 14 coefficients. Then, we tested those laws on radial and non-radial modes for evolved models exhibiting mixed-modes, in order to test how the empirical corrections perform when mixed-modes are involved. Overall, the combined correction law proposed by Ball & Gizon (2014) is found to give the best performer, closely followed by the law proposed by Sonoi et al. (2015). These two laws leave frequency differences that are less than 1 µHz on average, even reaching 0.1 µHz for the coolest stars of our set of model, which is of the order the frequency resolution provided by CoRoT and Kepler. We note that, on a low frequency range (0 < ν/ν max < 1.05), the Kjeldsen et al. (2008) power law (calibrated on this reduced range) gives equivalent results. Then the Kjeldsen et al. (2008) power law calibrated on the whole range of frequency and the purely cubic correction proposed by Ball & Gizon (2014) are the worst performer with remaining mean frequency differences of the order of few µHz. When applying those corrections on frequency spectra including mixed-modes, only the empirical corrections BG1 and BG2 proposed by Ball & Gizon (2014) and the modified S15E where we added a factor of 1/E improve the mean frequency dispersion. Only S15E and BG2 leave a satisfying root-mean dispersion of the order of the CoRoT and Kepler frequency resolution.
Therefore, we derived prescriptions for the fitting parameters of those radial modes correction empirical models as functions of log ∆ν, log T eff , log g and log κ which are quantities easily computed by 1D stellar evolution model. The next step will be to test our prescriptions against observed frequency spectrum in order to determine their degree of accuracy. We will focus on this in a furture work.
Finally, we only considered in this article the issue of structural effects. However, other effects such as non-adiabaticity effects may also play a non-negligible role in the propagation of acoustic waves in the surface layers. This will be studied in a forthcoming paper. correction In Table C.1 we gather the values of the root-mean frequency differences ν max √ D/N of radial modes after correction for each empirical law shown in Fig. 6.