Does the ν max scaling relation depend on metallicity? Insights from 3D convection simulations

,


INTRODUCTION
Space missions such as CoRoT (Michel et al. 2008), Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2015) have monitored light variations for hundreds of thousands of stars, and revolutionized the study of solar-like oscillations.Solar-like oscillations refer to those detected in Sun-like stars, subgiants, and red giants, which are stochastically excited by near-surface convection.Information about the frequency and amplitude of the oscillations allows us to infer stellar properties (e.g.Silva Aguirre et al. 2017;Yu et al. 2018) and detailed internal structure (e.g.Basu et al. 2009;Bellinger et al. 2019).
Two particularly important asteroseismic observables are the large frequency separation ∆ν and the frequency of maximum mode-power ν max , from which it is possible to obtain estimates of the stellar radii and masses by using the so-called asteroseismic scaling relations that are calibrated against the Sun (Ulrich 1986;Brown et al. 1991;Kjeldsen & Bedding 1995): where R and M stand for the total radius and mass of the star, and the subscript "⊙" represents properties that are associated with the Sun.The large frequency separation ∆ν, defined as the frequency difference between two adjacent radial modes n with the same spherical harmonic (or angular) degree l, is understood to represent the mean density ρ of a star (Ulrich 1986).The frequency of maximum oscillation amplitude ν max has no solid theoretical basis, but is empirically related to the acoustic cut-off frequency ν ac below which oscillations are trapped within the star, which in turn is connected to the stellar surface gravity g and effective temperature T eff in an isothermal atmosphere (Brown et al. 1991).If both ∆ν, ν max and T eff of a star are measured, Eqs. ( 1) and ( 2) can be applied to directly estimate the stellar radius and mass via (Stello et al. 2009;Kallinger et al. 2010): We refer the reader to Hekker (2020) for a thorough review of the asteroseismic scaling relation and its applications.
The asteroseismic scaling relations (1), (2) [equivalently (3)] are widely used for the determination of fundamental stellar parameters due to the simplicity of this method.For example, accurate stellar radii, which can be obtained from the scaling relation, are important for characterizing exoplanet systems (Huber et al. 2013); asteroseismic masses has been applied to quantify the mass loss along the red-giant branch in star clusters (Miglio et al. 2012); stellar ages constrained by asteroseismic radii and masses are crucial to understand the evolution of our Galaxy (Silva Aguirre et al. 2018).
Given their great significance, it is necessary to examine the validity of the scaling relations.By comparing the radii based on the scaling relation (3) against those derived from Gaia parallaxes, Huber et al. (2017) found good agreement (at a 5% level) for stars with radii between 0.8 and 10 R ⊙ (see also Sahlholdt et al. 2018).Gaulme et al. (2016) and Brogaard et al. (2018) tested the asteroseismic mass for three oscillating red giants in eclipsing binaries.It turns out that the scaling relation (3) overestimates masses by about 15% compared to dynamical masses.
Tests on individual scaling relations were also carried out.White et al. (2011), among other works (Guggenberger et al. 2016;Sharma et al. 2016;Rodrigues et al. 2017;Serenelli et al. 2017), calculated ∆ν and ρ from stellar structural models.They highlighted deviations from the ∆ν scaling given in (1) and introduced a correction factor f ∆ν .By applying this f ∆ν -corrected scaling relation, it was found that the radii align more closely with Gaia radii within 2% (Zinn et al. 2019).Moreover, the masses now agree with dynamical masses within observational uncertainties (Brogaard et al. 2018(Brogaard et al. , 2022)).On the other hand, based on the mixing length theory of convection and results from non-adiabatic calculations of mode excitation and damping, Belkacem et al. (2011Belkacem et al. ( , 2013) ) suggested that the ν max scaling relation should depend on Mach number near the stellar surface.Coelho et al. (2015) used model-derived stellar properties to rule out departures from the ν max scaling (2) above 1.5%.In a separate effort, Li et al. (2021) examined the intrinsic scatter of the scaling relations due to hidden variables, which was shown to be nearly negligible.
However, Epstein et al. (2014) studied nine metal-poor halo and thick-disk red giants and found their asteroseismic masses obtained from (3) are systematically higher than expectations by about 15%.The overestimation of stellar masses remains even if the corrected ∆ν scaling is adopted, therefore motivating the examination of the ν max scaling relation in the metal-poor regime.Viani et al. (2017) suggested the ν max scaling relation should depend on mean molecular weight (hence metallicity), based on the relationship between acoustic cut-off frequency and fundamental stellar parameters.Bellinger (2019) and Li et al. (2022) compared stellar properties estimated from the scaling relations with those obtained from individual frequency modeling.They clearly showed that a metallicity term in the scaling relations is needed in order to reach better agreements.Yet, as pointed out by Schonhut-Stasik et al. (2023), the overestimation of masses for metal-poor stars and the trend with metallicity might arise from offsets in the effective temperature scale rather than errors in the scaling relation.
Table 1.Fundamental parameters and basic properties of 3D near-surface convection simulations used in this work.Except for the solar model, model names are constructed from the effective temperature (labeled as "t"), surface gravity ("g"), and metallicity ("p" for positive while "m" for negative metallicity).Since effective temperatures derived from the radiative transfer calculations fluctuate over time, both the mean value and its standard deviation are shown.All models have the same surface gravity, which equals the nominal solar value of Prša et al. (2016).Numerical resolution is not uniform in the vertical direction, therefore both the minimum and the maximum grid spacing are given.Artificial driving simulations are computed with modified bottom boundary conditions.They are particularly designed for the evaluation of damping rates (cf.Sect.4).All artificial driving simulations have identical numerical resolutions in both space and time.In this study, we aim to answer whether the ν max scaling relation depends on metallicity from 3D hydrodynamical simulations.Zhou et al. (2019Zhou et al. ( , 2020) ) demonstrated that 3D simulations of stellar near-surface convection can be used to quantify the excitation, damping, and oscillation amplitude of pressure modes (p-modes) in a self-consistent, parameter-free manner.Here we follow the theoretical formulation and numerical technique developed in previous works to directly estimate ν max (Sect.4) based on six 3D models with almost identical surface temperature and gravity but spanning a wide range of metallicity.The correlation between theoretical ν max and metallicity of the 3D model reflects the role of metallicity in the ν max scaling relation, independent of asteroseismic observations.We also compute the important intermediate component of the ν max scaling relation, the acoustic cut-off frequency, from 3D models and investigate its relationship with metallicity (Sect.3).As a first step, the current work focuses on solar effective temperature and surface gravity.
2. MODELING 2.1.3D stellar atmosphere models 3D model atmospheres used in this study are computed with the Stagger code (Nordlund & Galsgaard 1995;Collet et al. 2018;Stein et al. 2023), a magneto-hydrodynamics (MHD) code that solves the equation of mass, momentum and energy conservation together with the equation of radiative transfer in three spatial dimensions and time.Magnetic fields can be considered by including the induction equation.Partial differential equations are discretized in Cartesian grids, with scalars evaluated at the center of numerical cells while momentum vectors are staggered at cell faces and magnetic fields at edge centers.In this work, we ignore the effect of the magnetic field but consider the radiative transfer processes in detail.The radiative transfer equation is solved at every numerical time step for every mesh point using the long-characteristic Feautrier (1964) method.Spatially, we consider nine directions in the radiative transfer problem, which consists of one vertical plus eight inclined directions representing the combination of two polar angles and four azimuthal angles.The total radiative cooling (heating) rate at a given wavelength is obtained by integrating over the solid angle.The integration over the polar angle is approximated by the Radau quadrature.The numerical technique employed by the Stagger code is explained in detail in Nordlund & Galsgaard (1995) and Collet et al. (2018).
Our simulations are carried out with realistic input physics.We use a customized version of the Mihalas et al. (1988) equation of state (Trampedach et al. 2013) that accounts for all ionization stages of the 17 most abundant elements in the Sun plus the hydrogen molecules.We include a comprehensive collection of continuum absorption and scattering sources as described in Hayek et al. (2010, their Table D.1) andTrampedach et al. (2014, Sect. 3.1).Line opacities are taken from the MARCS opacity sampling data covering wavelength from 91 nm to 20 µm (Gustafsson et al. 2008;Plez 2008).In order to reduce the computational cost, monochromatic opacities at a vast number of wavelengths are grouped into 12 "opacity bins" based on their approximate formation height (reflecting the strength of opacity), and wavelength (Nordlund 1982;Magic et al. 2013;Collet et al. 2018;Zhou et al. 2023).The organization of opacity bins is carefully adjusted such that the predicted radiative cooling rates resemble those from monochromatic calculations.
To study the effect of metallicity on the ν max scaling relation, we consider six models with [Fe/H] ranging from 0.5 to -3, which covers all solar-like oscillating stars discovered to date.We adopt the Asplund et al. (2009) solar metal mixture and apply an abundance enhancement for α-elements of 0.4 dex for metal-poor models with [Fe/H] ≤ −1.In our models, the helium mass fraction Y increases slightly with decreasing [Fe/H], from Y = 0.2422 at [Fe/H] = 0.5 to 0.2526 at [Fe/H] = −3.This trend is at odds with the picture of helium enrichment in the Universe (Jimenez et al. 2003).Future Stagger model atmospheres constructed with updated microphysics need to take into account the consensus that Y increases with metal mass fraction Z.All models are constructed with solar surface gravity and their effective temperatures are close to the solar value.Our 3D model atmospheres cover a small part of the star.Horizontally, the simulation domain has a square shape.The horizontal area is small compared with the surface area of the star but is large enough to enclose at least ten granules at any time of the simulation (Magic et al. 2013).Vertically, the bottom boundary of the simulation is located in the near-surface convection zone, about 2.5-3 Mm below the optical surface.The simulation domain straddles the stellar photosphere and extends up to Rosseland optical depth log τ Ross ∼ −7.Boundaries are periodic in the horizontal directions.The default bottom boundary condition of our simulations is that all inflows (fluid moving toward the surface) have identical, invariant thermal (gas plus radiation) pressure and entropy, whereas thermal pressure of outflows are adjusted to fulfill the equation of motion.For all models constructed with the default boundary conditions, there are 240 grid points in each horizontal direction and 230 points along the vertical.The spacing between two adjacent grid points is constant horizontally, whereas the numerical resolution is not uniform in the vertical direction: the highest vertical resolution is applied near the optical surface in order to sufficiently capture the transition from optically thick to thin regime.Our 3D simulations span a long time sequence for better resolution in the frequency domain.All six models for computing acoustic cutoff frequencies (Sect.3) and mode excitation rates (Sect.4) have the same time duration of 25 hours stellar time with one simulation snapshot saved every 30 seconds, which translates to a frequency resolution of ≈ 11 µHz.Their basic properties are summarized in Table 1.We note that the solar model is nearly identical to that used by Zhou et al. (2019) but extends slightly longer in time, while the five models with non-solar metallicities originate from the Stagger-grid (Magic et al. 2013) with small adjustments in surface gravity to agree with the value used in our solar model.
Computing mode damping rates from near-surface convection simulations for a number of frequencies is challenging.Following Zhou et al. (2019Zhou et al. ( , 2020)), we carried out numerical experiments where we artificially excite radial oscillations by perturbing the bottom boundary at a fixed frequency.Specifically, we add a small sinusoidal term on top of the constant thermal pressure of the incoming flows while keeping the entropy of inflows constant (in first order) at the bottom boundary: where P bot,0 and s bot,0 are constant bottom thermal pressure and entropy per mass used in standard simulations.
The term ϵ is a small, dimensionless number that governs the amplitude of perturbation, ω d is the angular frequency of artificial driving, and t stands for time.The perturbation at the bottom boundary will lead to coherent radial oscillations throughout the simulation domain.Detailed description and validation of this numerical technique has been presented by Zhou et al. (2020).
These artificial driving simulations are specifically designed for evaluating mode damping rates.For the purpose of controlling variables, all artificial driving simulations share the same numerical resolution (the distance between two adjacent grid points) in all dimensions.The solar and metal-poor models cover 6 × 6 Mm 2 area horizontally while the horizontal size of model t58g44p05 is 7 × 7 Mm 2 to ensure that at least ten granules are always included in the simulation domain1 .The metal-rich model therefore has more mesh points in horizontal directions (Table 1).The artificial driving simulations span 50 hours of stellar time.Given the size of the simulation domain, we crudely estimate the linewidths of the fundamental (frequency about 2 mHz) and first overtone simulation mode (frequency about 3 mHz, not to be confused with the driven oscillation) are in the order of 1 − 10 µHz and 10 − 50 µHz, respectively (see also Belkacem et al. 2019).These linewidths correspond to e-folding lifetimes of 8.8-88 hours and 1.76-8.8hours, respectively 2 .Time sequences of the artificial driving simulations thus span at least half of the lifetime of the fundamental simulation mode and several lifetimes for the first overtone mode.

Stellar structural models
Models of stellar structure are necessary for the calculation of mode eigenfunctions and mode masses, which are important components for evaluating excitation and damping rates (Sect.4).In the solar case, we use the standard solar model of Christensen-Dalsgaard et al. (1996, model S).At other metallicities, stellar interior models are computed with the Garching Stellar Evolution Code (GARSTEC; Weiss & Schlattl 2008).However, the interior model is only constrained by atmosphere parameters that define 3D models, namely T eff , log g, and [Fe/H].Therefore, different combinations of input parameters such as mass and mixing length parameter α MLT might give rise to multiple models that fit the three constraints well.In order to break this degeneracy, we construct stellar structural models based on a novel method developed by Jørgensen et al. (2018) and Mosumgaard et al. (2020), which couples the (interpolated) mean 3D stratification from the Stagger-grid with the 1D interior model at every step of stellar evolution.In practice, we replace the near-surface regime of the 1D interior model with the (interpolated) mean 3D model and set the outer boundary condition of the stellar evolution calculation from the mean 3D model in every timestep of stellar evolution (cf.Jørgensen et al. 2018 andMosumgaard et al. 2020 for a detailed description of the methodology).As the outer boundary of the interior model is located in the convection zone where the temperature gradient is close to adiabatic, the evolution calculation is not sensitive to α MLT (Mosumgaard et al. 2020 Sect. 3.2).To this end, we fix α MLT to the solar-calibrated value and only adjust the stellar mass to obtain the best-fitting model at non-solar metallicities.We also note that stellar structural models generated with the 3D-1D coupling method are consistent with 3D model atmospheres used in this study in terms of mean structures near the stellar surface.

THE EFFECT OF METALLICITY ON THE ACOUSTIC CUT-OFF FREQUENCY
The acoustic cut-off frequency, which can be understood as a barrier to acoustic waves, is an important quantity that motivates the link of ν max to g/ √ T eff .Assuming the stellar atmosphere is isothermal, the kinetic energy of pmodes with frequency lower than ν ac,1 will decay exponentially in the stellar atmosphere, hence trapped in the stellar interior.Higher frequency modes with ν > ν ac,1 , also called pseudo-modes (Roxburgh & Vorontsov 1995;Karoff 2007), are able to propagate through the atmosphere.Viani et al. (2017) proposed that under the assumption of isothermal atmosphere and ideal gas equation of state, where µ represents the mean molecular weight, which originates from the ideal gas equation of state and was ignored in the derivation of the canonical ν max scaling relation.Eq. ( 5) suggests that at fixed effective temperature and surface gravity, the acoustic cut-off frequency increases with metallicity (see Fig. 13 of Viani et al. 2017).A more general expression of the acoustic cut-off frequency was derived in Deubner & Gough (1984).For adiabatic oscillations whose wavelength is much less than the stellar radius, the acoustic cut-off frequency reads Figure 1.Left panel: First-order acoustic cut-off frequency computed via Eq. 5 as a function of spatially and temporally averaged Rosseland optical depth for all models included in this work.The gray shaded region outlines standard deviations of the temporal fluctuations of νac,1 for the solar model.We only show the distribution of acoustic cut-off frequencies up to log τRoss = −4, as our models are likely less reliable in regions higher up in the atmosphere due to the neglect of magnetic fields.
Right panel: The distribution of acoustic cut-off frequencies computed based on the general definition (Eq.( 6)).  which reduces to ν ac,1 in the isothermal scenario.
We compute acoustic cut-off frequencies in the isothermal limit (Eq.( 5)) as well as from the general definition (Eq.( 6)) for all 3D models.All components in Eqs. ( 5) and ( 6) are mean quantities averaged over constant geometric depth (horizontal average) and the entire time series.The distribution of acoustic cut-off frequencies near the optical surface and in the lower stellar atmosphere is shown in Fig. 1.The violent oscillation of ν ac just below the optical surface (right panel of Fig. 1) is associated with the large super-adiabatic temperature gradient at the top of the hydrogen ionization zone (Stothers & Chin 1995).This gives rise to a "density inflection region" (cf.Fig. 1 of Mosumgaard et al. 2020) and a rapid change in density scale height.Because acoustic cut-off frequencies change with optical depth in the stellar atmosphere, it is difficult to determine an exact ν ac value from our models.Nevertheless, we can see that the maximum of ν ac,1 and ν ac within −4 ≤ log τ Ross ≤ −1 anti-correlate with metallicity when −1 ≤ [Fe/H] ≤ 0.5.This Excitation rate [erg/s] t58g44p05 solar t58g44m05 t58g44m10 t58g44m20 t58g44m30 Figure 3. Excitation rates of radial oscillations, Pexc, as a function of cyclic frequency computed via Eq.( 7) for six models at solar T eff and log g.All curves shown here are smoothed from the original results with Gaussian kernels with a full width at half maximum (FWHM) of 0.47 mHz.
trend is in tension with the conclusion of Viani et al. (2017), who suggested at given T eff and log g, acoustic cut-off frequency increases with metallicity because of the dependence on the mean molecular weight (see their Fig.13).We remark that, unlike the present work, Viani et al. (2017) assumed the helium enrichment law of ∆Y /∆Z = 1 or 2 in their calculations of µ, whereas helium anti-correlates with metal mass fraction in our 3D models (Sect.2).With this in mind, we computed mean molecular weights for all 3D models and found that µ in stellar atmospheres decreases with metallicity in the range −1 ≤ [Fe/H] ≤ 0.5.Therefore, different helium mass fractions used in Viani et al. (2017) and this work is not the main cause of the tension in the acoustic cut-off frequency.The disagreement is associated with the assumption of isothermal atmosphere when deriving Eq. ( 5).Although isothermal is a good approximation for gray atmospheres, temperature profiles predicted by 3D models are far from isothermal and show substantial departure from the solar relation between T and T eff (Fig. 2).To this end, the full expression of the acoustic cut-off frequency is more appropriate in our context because it does not rely on the assumption of an isothermal atmosphere.
On the other hand, the physical meaning of acoustic cut-off frequency becomes quite ambiguous in the case of very metal-poor models t58g44m20 and t58g44m30, where their ν ac fluctuate considerably with optical depth.This is to a large extent due to their steep temperature gradient in the atmosphere (Fig. 2).We note that temperature stratification in the optically thin layers is largely dictated by the relative strength of adiabatic cooling in the upflows and radiative heating due to the absorption of photons (Asplund et al. 1999;Collet et al. 2007).For the [Fe/H]= −3 simulation, with the weakest line opacity, the balance between radiative heating and adiabatic cooling is dominated by the latter, [Fe/H]= −2 is a transition case, and [Fe/H]≥ −1 is dominated by radiative heating.In short, 3D model atmospheres suggest that in the range of −1 ≤ [Fe/H] ≤ 0.5, the acoustic cut-off frequency increases with decreasing metallicity.

ESTIMATING MODE EXCITATION, DAMPING RATES, AND ν max
In solar-type stars, turbulent convection near the stellar surface stochastically excites numerous p-modes at different frequencies, appearing as peaks in the observed power spectrum.Observationally, oscillation amplitudes are measured by first smoothing the raw power spectrum and then subtracting the background noise (e.g.Kjeldsen et al. 2005;Michel et al. 2009).The frequency corresponding to the peak of the amplitude spectrum is ν max .However, because the extent of our simulation domain is very small compared to the whole star, only 3-4 radial modes are spontaneously excited in the simulation box3 , and their amplitudes are not comparable to the measured amplitudes of stellar p-modes.In order to obtain realistic estimations of mode amplitudes and ν max , instead of analyzing the simulation modes, we quantify the energy injection (excitation) and energy dissipation (damping) rates for radial oscillations at different frequencies and derive mode amplitudes according to the energy balance between excitation and damping.Figure 4. Left panel: Linear damping rates computed from artificial driving solar simulations at different cyclic frequencies (blue triangles) are divided by π to compare with observed line widths from BiSON l = 0 data (black dots).Below 4 mHz, the observational data, with uncertainties, is taken from Chaplin et al. (2005).BiSON line widths above 4 mHz are derived from Table 1 of Chaplin et al. (1998), therefore uncertainties are not accessible.The cyan line is obtained by smoothing the raw simulation data with a Gaussian kernel with an FWHM equal to 0.2 mHz.Right panel: Distribution of smoothed damping rates for six models with different metallicities.An identical smoothing kernel is used for all models.The black dashed line corresponds to the cyan line in the left panel.Figure 5. Predicted photospheric velocity amplitudes for models with different metallicities.The range of theoretical νmax, 2.88 ≤ νmax ≤ 3.01 mHz, estimated for the solar and all sub-solar metallicity models is shaded in gray.The vertical dotted gray line marks the measured νmax for the Sun (Kjeldsen et al. 2008).
Excitation rates of radial oscillations are computed as where F represents the Fourier transform from time to frequency domain, Re and Im represent real and imaginary part, respectively.The terms ω = 2πν, A box , and T tot are the angular frequency, horizontal area, and total time coverage of the simulation, respectively.Numerically, the integration is carried out from the radius corresponding to the bottom boundary of the simulation domain, r 3D bot , to the uppermost point of the coupled stellar structure model introduced in Sect.2.2, r top .δ Pnad is the horizontally averaged non-adiabatic pressure fluctuation caused by convection.It is obtained by subtracting the adiabatic part from the total pressure fluctuation in 3D simulations, therefore including contributions from entropy fluctuation and fluctuations in turbulent pressure (Reynold stress).How pressure fluctuations are evaluated is detailed in Zhou et al. (2019, their Sect. 3.2).The adiabatic eigenfunction of radial modes (also called the amplitude function) ξ r and mode kinetic energy per unit surface area are calculated using the Aarhus adiabatic oscillation package (ADIPLS, Christensen-Dalsgaard 2008) based on the coupled stellar structure model.R phot is the photospheric radius.We note that Eq. ( 7) is only valid for adiabatic eigenfunctions (the "quasi-adiabatic" approximation).The derivation of Eq. ( 7) and the physical mechanism of p-mode excitation are explained in detail in Nordlund & Stein (2001), Stein & Nordlund (2001) and Zhou et al. (2019).
Excitation rates for all models included in this study are depicted in Fig. 3. Following previous works, we smooth P exc calculated from simulation data in order to eliminate noise associated with stochastic convection.Below 2.5 mHz, it is evident that mode excitation positively correlates with metallicity.However, between 2.5 and 4 mHz where P exc exhibits a plateau feature, there is no clear trend between P exc and [Fe/H].
Damping rates of radial oscillations, which are related to the observed mode line widths Γ through η = πΓ, are computed in frequency space according to Here δ ρ * is the complex conjugate Fourier component of the horizontally averaged density fluctuation caused by oscillations and ρ0 is the equilibrium density obtained by averaging ρ over all simulation snapshots.δ ρ * /ρ 0 is related to the mode eigenfunction via the perturbed fluid continuity equation (cf.Aerts et al. 2010 Chapter 3).In principle, the work integral in Eq. ( 9) should be computed for the entire propagation cavity of the oscillation mode.Practically, the integration range is limited by the vertical coverage of the simulation, with the lower and upper bound y bot and y top being the geometric depth at the bottom and top simulation boundary, respectively.Nevertheless, for n ≳ 15 solar radial modes, damping processes take place around the photosphere and the super-adiabatic layers just below the optical surface according to theoretical predictions (Balmforth 1992 Sect. 7.2 andHoudek 1996 Sect. 3.4).Therefore, truncating the work integral will not affect the evaluation of damping rates for medium-and high-frequency p-modes.For all six models, our numerical results also show that contribution from layers near the bottom boundary to the work integral is negligible for oscillations with ν ≳ 2.6 mHz, indicating that the vertical extent of our simulations is sufficient for quantifying damping rates for these higher frequency oscillations.The denominator of Eq. ( 9) is proportional to the mode kinetic energy, where m mode is the mode mass per unit surface area and V is the vertical velocity amplitude.We note that m mode is computed based on the coupled stellar structured model, whereas V is extracted from the 3D simulation.
Quantifying damping rates from Eq. ( 9) is challenging because it requires the knowledge of density and pressure fluctuation as well as their phase difference.Therefore, density fluctuations associated with mode eigenfunctions need to be extracted from simulations.For standard surface convection simulations, this can be done only at the frequency of the simulation mode (cf.Belkacem et al. 2019), that is, p-mode oscillations naturally excited in the simulation box.Given the extent of our simulation domain, there are only 3-4 radial simulation modes from which we can compute damping rates.This difficulty was overcome by carrying out multiple artificially driven simulations (see Sect. 2.1) at different driving frequencies, which gives η as a function of frequency4 (Fig. 4).In the left panel, predictions from solar simulations are compared with measured radial mode line widths from BiSON (Chaplin et al. 1998(Chaplin et al. , 2005)).In the solar case, the agreement between modeling and observation is encouraging above ≈ 2.6 mHz.The location of the "dip" (i.e.local minimum) in damping rates is also well-reproduced.Our simulations systematically underestimate damping rates below 2.6 mHz.The problem at low frequencies is a known deficiency of our approach, which can be partly attributed to the limited vertical extent of the simulation domain compared to the whole star.This cuts off significant excitation and damping, interior of the simulations, at low frequency.Also, the numerically predicted damping rates demonstrate stronger fluctuations across frequencies than the measured values due to the relatively short time coverage of our simulations compared to BiSON observations, which are in the order of decades (Chaplin et al. 2005).Nevertheless, by extending the time span of artificial driving simulations from 10 hours in Zhou et al. (2019) to 50 hours in the present investigation, fluctuations in η are slightly reduced, especially in the "plateau" region (see Fig. 4 and Zhou et al. 2019 Fig. 7).As further increasing the simulation time sequence is difficult because of practical issues such as data storage, we smooth the damping rates quantified from the simulations to make them more comparable with observations.The smoothed damping rates for all models investigated in this work are shown in the right panel of Fig. 4.Although the magnitude of η is apparently affected by how smoothing is performed, the location of the damping rate dip, which is closely related to ν max , is not sensitive to the smoothing method adopted.For all models with solar and sub-solar metallicities, the depression of their damping rates is located at approximately the same frequency.For model t58g44p05, the local minimum of η occurs at about 2.7 mHz, being 0.2 mHz lower than other models.We emphasize that the result of the metal-rich model is likely less reliable, as the depression in η resides at the frequency edge where our numerical results could be trusted.
Once both excitation and damping rates are available, the mean kinetic velocity amplitude of the mode at the stellar surface can be quantified from the energy balance between driving and damping: where M mode is the total mode mass (Aerts et al. 2010 Eq. 3.140, not to be confused with the mode mass per unit surface area m mode ) calculated from the coupled stellar structure model.Smoothed excitation and damping rates are used to compute the velocity amplitude, which is shown in Fig. 5. From definition, the theoretical ν max is the frequency corresponding to the peak of the velocity amplitude.The thus estimated ν max is about 2.7 mHz for model t58g44p05.However, due to larger uncertainties in η at lower frequencies, the ν max estimation for this model is less reliable therefore we will exclude it in subsequent discussions.The frequency of maximum power predicted from the solar simulation is in reasonable agreement with the measured value (ν max,⊙ = 3.1 mHz, Kjeldsen et al. 2008).The remaining discrepancy is likely due to intrinsic errors in our determination of excitation and damping rates, such as the limitation of box-in-a-star simulations or the artificial driving method.Owing to the same reason and also being aware that the smoothing kernel applied to P exc and η affects the exact value of theoretical ν max , we elect to provide a range rather than exact values for ν max as the latter is ambiguous within our approach.The solar and metal-poor models all have well-defined velocity peaks located at nearly identical frequencies.Their ν max are estimated to be in the range of [2.88, 3.01] mHz, coinciding with the dip in damping rates.

DISCUSSION
Results from Sect. 4 lead to the conclusion that for [Fe/H] ≤ 0, the theoretical ν max does not demonstrate any correlation with metallicity.In order words, the ν max scaling relation does not depend on metallicity at solar T eff and log g.It is worth noting that slight differences in T eff among our 3D models (Table 1) have a negligible impact on our conclusion.The mean effective temperature ranges from 5768 K for model t58g44m10 to 5794 K for t58g44m05.Such differences will result in 0.23% difference in ν max according to the scaling relation, translating to ≈ 7 µHz absolute difference in frequency.Considering the intrinsic systematics of our approach, the uncertainty associated with T eff is minor.On the other hand, while the predicted velocity amplitude spectra match observations in the bulk part (Zhou et al. 2019 Fig. 8), our numerical results are not sufficiently accurate to draw any conclusions about the relationship between the oscillation amplitude and [Fe/H]5 .5.1.Why is the ν max scaling relation not sensitive to metallicity?When [Fe/H] ≤ 0, the effect of metallicity on ν max is less than 4.5%.The upper limit arises from uncertainties in our numerical method rather than suggesting any correlation between ν max and [Fe/H].Given the fact that stars at the same location in the HR diagram but differing in chemical composition have different stratification (e.g., T − τ relation and pressure scale height, see Fig. 2) near the photosphere, this conclusion is not trivial.
In Sect.4, we found that the maximum velocity amplitude is linked to the depression in damping rates, in accordance with the conclusions of Balmforth (1992) and Belkacem et al. (2011).As suggested by Balmforth (1992), the dip in η is caused by the negative peak (destabilizing, or exciting) in the thermal contribution to the total damping, which counteracts the convective turbulence that stabilizes (dampens) solar-like oscillations (cf.Balmforth 1992 Fig. 13).In 3D simulations, however, the destabilizing effect from thermal pressure fluctuations is strongest at 3.5-4 mHz (Zhou The turnover frequency is computed from the time-averaged rms of vertical velocity fluctuation and the horizontal-and timeaveraged pressure scale height.The relative difference in νturn between the metal-poor and solar model is shown in the lower panel.Shaded areas in gray, light green and light yellow represent standard deviations of the time variation of νturn (and their relative differences) for the solar model, t58g44m10 and t58g44m30, respectively.et al. 2019, Fig. 7).The location of the dip around 2.95 mHz results from subtle cancelations between the positive turbulence and the negative thermal contributions to mode damping.Instead of explaining the cause of the dampingrate depression with convection theory (see Belkacem et al. 2011 for efforts in this direction), we opt to examine if the characteristic timescale (frequency) near the stellar surface depends on [Fe/H].A frequently used timescale in the convective region is the convective turnover time, defined as the time spent for materials to travel one pressure scale height in the vertical direction (Freytag et al. 2012), H P /v y,rms , where v y,rms is the root-mean-square (rms) of vertical velocity fluctuation.Here we compare the convective turnover frequency from different models.Although our approach cannot provide a solid theoretical basis for the ν max scaling relation since the detailed relationship between convective turnover frequency and ν max is unknown, it offers useful insights into the problem in a straightforward manner.
The distribution of convective turnover frequencies, ν turn = v y,rms /H P , within −2 ≤ log τ Ross ≤ 5 is depicted in Fig. 6 for all 3D models studied except the metal-rich one, together with the relative difference in ν turn with respect to the solar model.We note that layers with lower or higher optical depths are irrelevant in this context because for radial oscillations near 3 mHz, contributions to the work integral (Eq.( 9)) mostly originate from the near-surface region with −2 ≲ log τ Ross ≲ 5 (approximately from 0.2 Mm above the optical surface to 1 Mm below it in the solar case) as revealed by our artificial driving simulations.Fig. 6 shows that models with different metallicity have similar turnover frequencies near the stellar surface, with fractional differences less than 5% when τ Ross > 0.1.As indicated in Magic et al. (2013), metal-poor stars have smaller pressure scale height and granule size.Their convective velocity is also smaller than that predicted by the solar-metallicity model in the super-adiabatic region.Here we demonstrate that these trends with [Fe/H] nearly cancel out when considering frequencies or timescales.The lack of metallicity dependence for ν max is thus reinforced by the insensitivity of convective turnover frequency to [Fe/H].Furthermore, it is worth mentioning that Rodríguez Díaz et al. ( 2022) quantified the characteristic granulation timescale for four metallicities at solar T eff and log g and found less than 4% relative difference across −1 ≤ [Fe/H] ≤ 0.5 (cf. Rodríguez Díaz et al. 2022 Table 1).To summarize, frequencies that characterize the near-surface region, including ν max , are insensitive to metallicity in general, at least at solar surface temperature and gravity.In the original derivation of the ν max scaling relation (Brown et al. 1991), ν max was related to fundamental stellar parameters via the acoustic cut-off frequency.Results from 3D simulations indicate that ν max is insensitive to metallicity, whereas a clear anti-correlation between ν ac and [Fe/H] is seen within −1 ≤ [Fe/H] ≤ 0.5 (Fig. 1).Putting these in the context of the ν max − ν ac − g/ √ T eff relation, it is then implied that the ν max − ν ac relation depends on [Fe/H] in an opposite way as the ν ac − g/ √ T eff one such that effects of metallicity cancels out for the scaling relation overall.To assess whether this implication is reasonable, we build our analysis on the detailed physical explanation of the ν max − ν ac relation proposed by Belkacem et al. (2011Belkacem et al. ( , 2013)), which is derived from the mixing length theory of convection, meanwhile drawing on results from non-adiabatic calculations of stellar oscillations.Belkacem et al. (2011) suggested that ν max is inversely proportional to the thermal timescale in the super-adiabatic region, which in turn relates to ν ac through thermodynamic derivatives, the Mach number, and the mixing length parameter (Belkacem et al. 2011, Eq. (15)).Moreover, most of the dispersion in the ν max − ν ac relation is due to the Mach number.With inputs from 3D simulations, Belkacem et al. (2013) proposed that The distribution of Mach number, M = v y,rms /c s , in the near-surface layers predicted from our simulations is shown in Fig. 7.It is evident that below the optical surface, the Mach number decreases with metallicity.Relation (11) therefore suggests a positive correlation between ν max and [Fe/H] at fixed ν ac , qualitatively in line with implications from our numerical results.Another prominent feature of M is that it varies significantly in the region of interest.In combination with the fact that ν ac is not a constant in stellar atmosphere, these factors make quantitative interpretation of the ν max − ν ac scaling (11) difficult.Here we crudely estimate the effect of metallicity on the ν max − ν ac relation by arbitrarily focusing on a single depth τ Ross = 10.At this location, the Mach number of the solar model is 1.08 times larger than model t58g44m10.Inserting this ratio into relation (11) gives a factor of 1.24 deviation of the ν max − ν ac scaling from [Fe/H] = 0 to −1.Considering the difference in ν ac between the solar model and t58g44m10 is roughly 3%, the relationship (11) appears to be far too sensitive to metallicity from the perspective of our 3D simulations.
In short, the ν max − ν ac relation proposed by Belkacem et al. (2011Belkacem et al. ( , 2013) ) suggests a positive correlation between ν max and [Fe/H] at given ν ac , which qualitatively agrees with what implied from our numerical results.However, it is challenging to examine the relationship (11) quantitatively because both the Mach number and acoustic cut-off frequency vary within the near-surface region.This illustrates the complexity of the ν max − ν ac − g/ √ T eff relation from another angle, that ν max is a constant, whereas the excitation and damping processes that determine its value primarily occur in the near-surface, super-adiabatic region where physical quantities vary considerably in the radial direction.

Influence of magnetic fields on 3D surface convection simulations
In this study, we ignore the effect of the magnetic field.Magnetic fields can be included in the simulation by imposing a fixed field strength and orientation at the boundary of the simulation domain (e.g.vertical magnetic fields as adopted in Vögler 2005 andLudwig et al. 2023) or setting up an initial seed field with zero net flux and a small field strength that allow it to grow and converge to a natural state of the simulation (small-scale dynamo [SSD], e.g.Thaler &Spruit 2015 andBhatia et al. 2022).In MHD simulations of stellar surface convection, magnetic fields tend to concentrate at intergranular lanes.The average size of granules is smaller in simulations with a strong field strength (Vögler 2005;Jacoutot et al. 2008).For the solar case, MHD models predict flatter temperature gradients near the photosphere compared with hydrodynamical models, hence have an impact on the modeled center-to-limb variations (Pereira et al. 2013;Ludwig et al. 2023).Also, including magnetic fields will slightly reduce the pressure scale height at the optical surface.Dynamically, magnetic fields inhibit surface convection, resulting in smaller convective velocities and velocity fluctuations for solar-type stars (Bhatia et al. 2022).The amplitude of simulation modes decreases with increasing field strength, indicating magnetic fields suppress p-mode oscillations in 3D simulations (Jacoutot et al. 2008;Kitiashvili et al. 2011).This trend is in line with the effect of solar and stellar activity on oscillation amplitude as observed by BiSON and Kepler (Chaplin et al. 2000(Chaplin et al. , 2011)).Furthermore, helioseismic observations uncovered a positive correlation between solar activity and measured ν max , which shifts approximately 25 µHz between low and high activity (Howe et al. 2020).Since the effect of magnetic fields on ν max is much smaller than the intrinsic systematics of ν max estimated from numerical simulations (Sect.4), our neglect of magnet fields will unlikely affect the result of this work.

CONCLUSIONS
In this paper, we investigated whether the frequency of maximum power depends on metallicity at solar effective temperature and surface gravity by quantifying velocity amplitudes for radial modes based on ab initio simulations of near-surface convection covering a wide metallicity range.Velocity amplitudes were determined from the balance between energy supply and dissipation, both of which were calculated from first principles following methods described in Zhou et al. (2019Zhou et al. ( , 2020)).Our formulation and numerical approach enable purely theoretical predictions of ν max independent from observational data.
The main conclusion of this work is that when [Fe/H] ≤ 0, ν max does not demonstrate any trend with [Fe/H].Furthermore, the convective turnover frequencies in the super-adiabatic layers just below the optical surface are nearly invariant across [Fe/H] as both the length scale and convective velocity are smaller in metal-poor stars, implying that for solar-type stars, frequencies characterizing the near-surface region are not sensitive to [Fe/H] in general.We also calculated the acoustic cut-off frequency for all 3D models for better insights into the ν max − ν ac − g/ √ T eff scaling.In the metallicity range −1 ≤ [Fe/H] ≤ 0.5, an anti-correlation between ν ac and [Fe/H] is unambiguously seen, which is in tension with the conclusion of Viani et al. (2017).The disagreement is attributed to the assumption of isothermal atmosphere in the derivation of the ν ac − g/ √ T eff scaling.On the other hand, ν max and ν ac predicted from 3D simulations suggest a positive correlation between ν max and [Fe/H] at fixed ν ac .This is qualitatively in line with the ν max − ν ac scaling (11) proposed by Belkacem et al. (2013).Nevertheless, relation ( 11) is likely to be over-sensitive to metallicity owing to its strong nonlinear dependence on Mach number.
For solar-type stars, the effect of metallicity on the ν max scaling relation is estimated to be less than 4.5%.The upper limit here represents the intrinsic uncertainties of our numerical approach rather than indicating any possible correlation between ν max and [Fe/H].It is worth noting that the exact value of observed ν max is affected by the instrument, the data analysis process, and by the stellar activity cycle to a lesser extent (Howe et al. 2020).For example, the measured solar ν max ranges from 3060 (measured by the VIRGO instrument aboard the Solar and Heliospheric Observatory, Hekker et al. 2013) to 3229 µHz (measured by GOLF, Breton et al. 2020).Unless homogeneous analyses are carried out for both the Sun and other stars, the systematics in ν max associated with observations is comparable to the upper limit of the metallicity dependence suggested by this work.To this end, we argue that the original ν max scaling relation (2) can be applied to metal-poor solar-type stars with confidence.
The present analysis is confined to the parameter space of the Sun.In the next step, it will be meaningful to investigate the effect of metallicity on the ν max scaling relation for solar-like oscillators across the HR diagram, especially for red-giant stars.Thousands of oscillating giants from solar metallicity down to [Fe/H] ≃ −2.5 were detected by Kepler and TESS (Yu et al. 2018;Schonhut-Stasik et al. 2023).Their masses and radii inferred from the asteroseismic scaling relation have important applications in Galactic archaeology and the study of globular clusters (e.g.Casagrande et al. 2016;Miglio et al. 2016;Tailo et al. 2022;Howell et al. 2023).
Meanwhile, we are cautious about the limitations of the adopted numerical approach.The artificial driving simulations failed to predict reasonable damping rates at low frequencies because of their small spatial and temporal coverage.The same shortcoming gives rise to fluctuations of η across frequencies, making the final velocity amplitude affected by how smoothing is performed.These deficiencies prevent us from investigating the relationship between oscillation amplitude and metallicity, where a positive correlation was concluded from asteroseismic observations (Yu et al. 2018), or studying how the width (FWHM) of the envelope of oscillation power excess changes with [Fe/H].Resolving these problems from a theoretical angle would require accurate oscillation amplitudes for a broad range of frequencies predicted from more advanced simulations that cover the whole stellar convective envelope (e.g., Hotta et al. 2019;Popovas et al. 2022) and possibly includes magnetic fields (Käpylä 2022), a more realistic, analytical theory of the interaction between convective turbulence and oscillations (Philidet et al. 2022), or the combination of both.

Figure 2 .
Figure 2. The relationship between horizontal-and time-averaged temperature and Rosseland optical depth (T − τ relation) predicted by 3D model atmospheres at different metallicities.The Eddington T − τ relation, which is an approximation to the gray atmosphere, is shown in solid black line.

Figure 6 .
Figure6.Distribution of the convective turnover frequency in the near-surface region for the solar and metal-poor models.The turnover frequency is computed from the time-averaged rms of vertical velocity fluctuation and the horizontal-and timeaveraged pressure scale height.The relative difference in νturn between the metal-poor and solar model is shown in the lower panel.Shaded areas in gray, light green and light yellow represent standard deviations of the time variation of νturn (and their relative differences) for the solar model, t58g44m10 and t58g44m30, respectively.

Figure 7 .
Figure 7. Mach number, computed as the ratio between the rms vertical velocity fluctuation and the spatial and temporal mean of sound speed, as a function of Rosseland optical depth for five models with different metallicities.Standard deviations of the time fluctuation of Mach number are shaded for the solar model, model t58g44m10 and t58g44m30.5.2.The ν max − ν ac − g/ √ T eff relation from the perspective of surface convection simulations