© The Authors 2024

1 Introduction

Spectropolarimetric observations show that the atmospheres of cool main-sequence stars are permeated by magnetic fields covering a wide range of spatio-temporal scales (for reviews, see e.g. Donati & Landstreet 2009; Reiners 2012; Kochukhov et al. 2017 and references therein). The surface of the Sun, the most studied cool main-sequence star, is covered by magnetic fields on spatial scales ranging from the 30 000 km of super-granular cells down to the smallest magnetic features resolvable with current instruments on scales of less than 100 km (see e.g. Solanki et al. 2006; Wedemeyer-Böhm et al. 2009; Fan 2021). The interaction of such magnetic fields with atmospheric plasma flows plays a crucial role in determining the dynamics of particles and heat transfer near the solar surface and the photometric and spectral variability of the Sun. For example, large-scale magnetic fields are responsible for the emergence of dark sunspots, whereas small-scale magnetic flux concentrations appear as bright features. Evidence of magnetic activity leading to the formation of Sun-like spots has also been found for stars other than the Sun (see e.g. Berdyugina 2005; Strassmeier 2009). However, because of the lack of high-resolution observations for cool main-sequence stars other than the Sun, very little is known about magnetic fields on scales much smaller than the radius of these stars. This is a critical issue not only in the context of stellar physics, but also for exoplanet detection. Transit observations, and in particular transmission spectroscopy, the technique used to determine the atmospheric composition of an exoplanet, rely on a precise understanding of the wavelength-dependent brightness fluctuations of the star being occulted. This, in turn, requires a detailed modelling of the magnetic fields and of the turbulent plasma flows present in the proximity of the stellar surface (Rackham et al. 2023).

High-resolution observations of the Sun suggest that small-scale magnetic activity is ubiquitous on the solar surface, and that the magnetic flux in internetwork regions is independent of the solar cycle (Ramelli et al. 2019). A widely accepted explanation for this behaviour is the presence of a small-scale dynamo (SSD) operating in the sub-surface layers of the Sun. In the recent past, radiative magnetohydrodynamic (MHD) numerical simulations have provided many insights into the functioning of the sub-surface solar SSD, showing that a considerable amount of magnetic energy in the solar photosphere could be self-sustained through dynamo action (see e.g. Cattaneo 1999; Vögler & Schüssler 2007; Pietarila Graham et al. 2010; Rempel 2014; Kitiashvili et al. 2015; Thaler & Spruit 2015; Khomenko et al. 2017). On the other hand, numerical studies of small-scale magnetic fields on the surface of cool main-sequence stars other than the Sun mostly focus on how the magnetic field self-organizes and affects the surrounding plasma, neglecting how the magnetic field was initially produced. Seminal work in this respect was done by Beeck and collaborators in 2011, further extended in 2015, who found a substantial variation in the properties of the surface magneto-convection between main-sequence stars of different spectral types as a consequence of including magnetic fields in the simulations (Beeck et al. 2011, 2015). Additional work in this context was done by Steiner and collaborators, who analysed MHD models of spectral types K8V to F5V, showing that the field strength of the strongest flux concentrations was fairly independent of spectral type when measured at the optical depth level τ_R = 1 (Steiner et al. 2014). This work was then extended by Salhab et al. (2018), who found that small-scale magnetism increases the bolometric intensity and flux in all the investigated models, most prominently for the Sun.

The possibility of simulating SSD action on stars other than the Sun has been investigated only very recently. In Bhatia et al. (2022), the authors discuss how small-scale magnetic fields generated by SSD action on the surface of different cool main-sequence stars introduce small but significant changes in thermodynamic stratification due to a reduction in the convective velocities. In a second publication of the series, it is shown how the SSD is affected by the metallicity of the stellar plasma (Witzke et al. 2023). These works provided a first insight into the functioning of SSDs on the surface of stars other than the Sun, but a number of questions related to SSD action in the atmosphere of cool main-sequence stars still remains. This motivates the present work.

In this paper we present simulations and their analysis of the same spectral sequence, that is same effective temperatures and solar metallicity, that was considered in the work of Salhab et al. (2018), which are K8V, K2V, G2V, and F5V, but with a surface gravity more consistent with stars on the main sequence. All simulations were carried out with the CO⁵BOLD code (Freytag et al. 2012), starting from a relaxed non-magnetic state, to which we add a small magnetic field seed, which then is evolved until saturation of the magnetic energy. The aims of our work are (i) to investigate the differences in SSD action among the four models during the kinematic phase; (ii) to study the properties of magnetic fields amplified by SSD action and how they affect the plasma; and (iii) to examine the properties of surface magnetic flux concentrations that originate from the self-organization of magnetic fields produced by SSD action (the magnetic filigree). In Sect. 2, we present the simulation set-up and the parameters used for the computations. Then, in Sect. 3 we provide a general overview of the simulation results. Section 4 discusses the numerical results for the kinematic phase, during which the magnetic field is amplified exponentially in time. The results obtained after saturation of the magnetic energy are discussed in Sect. 5. Finally, in Sect. 6 we report our conclusions. In Appendix A the results of the simulations are discussed in normalized units, whereas the algorithm used to identify magnetic bright features (MBFs) is presented in Appendix B.

2 Methods and stellar models

The simulations discussed in the present paper were carried out with the CO⁵BOLD code using the ‘box-in-a-star’ set-up. This corresponds to simulating a small partial volume near the optical surface of a star, encompassing the layers where energy transport changes from predominantly convective to purely radiative. We considered four different cool stellar models, as discussed in the following.

2.1 The CO⁵BOLD model and implementation

The CO⁵BOLD code is a numerical tool designed to simulate solar and stellar atmospheres, as well as their interiors, by solving the time-dependent ideal MHD equations, in cgs units, $$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho v)=0,$$(1) $$\frac{\partial \rho v}{\partial t}+\nabla \cdot \left[\rho vv+\left(P+\frac{B^2}{8\pi }\right)I-\frac{BB}{4\pi }\right]+\rho \nabla \text{Φ}=0,$$(2) $$\frac{\partial B}{\partial t}+\nabla \times (B\times v)=0,$$(3) $$\frac{\partial \rho e_{\text{tot}}}{\partial t}+\nabla \cdot \left[\left(\rho e_{\text{tot}}+P+\frac{B^2}{8\pi }\right)v-\frac{1}{4\pi }(v\cdot B)B+F_{\text{rad}}\right]=0,$$(4)

in combination with an equation of state (eos), used to relate the mass density ρ and the gas pressure P to the internal energy e_i, and a non-grey radiative transfer equation. Here v is the plasma velocity, B the magnetic field, ρe_tot = ρe_i + ρυ²/2 + B²/(8π) + ρΦ the total energy, Φ the gravitational potential, F_rad the frequency-integrated radiative energy flux vector, I the identity matrix, and a⋅b and ab are the scalar and tensor products of vector a with vector b, respectively. The code can be used to carry out both global (star-in-a-box) and local (box-in-a-star) simulations. In the following, only box-in-a-star simulations are considered.

Equations (1)–(4) are evolved using flux-conserving Riemann-type solvers. While a Harten-Lax-van Leer (HLL) solver (Harten et al. 1983) must be used for magnetic simulations, a Roe solver (Roe 1986) can be used as a less diffusive alternative for the non-magnetic case. In regions of low plasma-β, with β the ratio of thermal to magnetic pressure, the equation of the thermal energy is used instead of the equation of the total energy to avoid negative gas pressures, this at the expense of strict energy conservation. The Alfvén speed, υ_A, can be limited by artificially reducing the strength of the Lorentz force to avoid very small time steps. Moreover, to impose ∇⋅B = 0, the constrained-transport method of Evans & Hawley (1988) is used. Pre-tabulated values as a function of density and internal energy, accounting for the ionisation balance of hydrogen, helium, and a representative metal, are used to solve the eos and to compute the plasma pressure and temperature. To solve the radiative transfer equation, either a short- or a long-characteristic scheme can be used. Opacities are obtained from opacity tables computed using the opacity binning technique, for which spectral lines of similar strength are grouped in bands, and average opacities are computed for each band (Nordlund 1982; Ludwig 1992; Freytag et al. 2012). Additionally, to reduce the computational cost of a simulation, the diffusion approximation can be used to model the radiative transport in the deep, optically thick layers of the numerical domain. The code is written in a modular form using the language Fortran 90 (and some Fortran 77). It is parallelised using a hybrid Message Passing Interface (MPI) and OpenMP parallelization. For more details on CO⁵BOLD, we refer to Freytag et al. (2012).

For the present paper we used the HLL solver in combination with a second-order reconstruction scheme (of FRweno type; see Freytag 2013). The gravitational field, g = −∇Φ, was assumed vertical and constant. The thermal energy equation was activated for β < 0.1 and υ_A was limited to a maximum of 50 km s⁻¹. The equations were then evolved in time with a Hancock predictor step. Moreover, we solved the radiative transfer equation with the long-characteristic radiation transport module, using Rosseland mean opacities (grey opacity) and employing the diffusion approximation for depths with min(τ_R) > 100, where the minimum was evaluated over horizontal planes at constant geometrical height. The Rosseland mean opacities entering the opacity tables were adapted from the MARCS stellar atmosphere package (Gustafsson et al. 2008).

2.2 Stellar parameters

The main governing parameters in these simulations are the effective temperature, T_eff; the gravitational acceleration, g; and the chemical composition. Solar metallicity is used for all models. For the effective temperatures, we chose to approximately replicate the nominal temperatures of the four models in Salhab et al. (2018), that is T_eff = 4000 K (K8V), 5000 K (K2V), 5770 K (G2V), and 6500 K (F5V). Given these parameters, we used the PARSEC code (Nguyen et al. 2022) to generate two different isochrones, considering 2.0 Gy and 4.6 Gy old stars. The first isochrone was then used to estimate the surface gravity of the F5V model (log 𝑔 = 4.24), whereas the second isochrone was used for the K models (log 𝑔 = 4.66 and log 𝑔 = 4.59 for K8V and K2V models, respectively). Finally, for the G2V model, we considered the solar surface value, log 𝑔 = 4.44. These parameters are summarized in Table 1, where values of four Gaia benchmark stars are also reported to give an idea of where the present models lie on the Hertzsprung-Russell diagram (first two and last three rows, respectively).

2.3 Numerical set-up

The numerical set-up of the simulations in this paper is similar to that used in Riva & Steiner (2022) for the case dh8. More precisely, we consider a Cartesian computational domain with horizontal extension such that the box widths, L_x = L_y, are approximately six mean granular scales. The vertical extension, L_z, is chosen such that the bottom boundary is at a depth where the entropy of the downflowing plasma approaches the essentially constant entropy flow of the upflowing plasma (cf. Fig. 11 of Salhab et al. 2018), whereas the top boundary is located well above the entropy minimum in the photosphere. The exact physical scales and sizes of the four models are reported in Table 1, which also displays the height extent of the boxes in terms of pressure scale heights, with N^below = ln(P^bottom/P⁰) and N^above = ln(P⁰/P^top) the number of pressure scale heights below and above the optical surface, and P⁰, P^bottom, and P^top the horizontal and time-averaged gas pressures at the geometrical height z = 0, at the bottom, and at the top of the domain, respectively. Here, z = 0 denotes the space and time-averaged optical depth surface τ_r = 1.

The granular size is defined as L_gran ≡ 〈L_int(z = 0, t)〉_t, where 〈−〉, denotes time averaging and with $$L_{\text{int }}(z,t)=2\pi {{{\displaystyle \int }}^{\text{}}}_0^{\infty }{k}_{\text{h}}^{-1}{E}_{vz}\left(k_{\text{h}},z,t\right)\text{d}{k}_{\text{h}}/{{{\displaystyle \int }}^{\text{}}}_0^{\infty }{E}_{vz}\left(k_{\text{h}},z,t\right)\text{d}{k}_{\text{h}}$$(5)

the integral length scale for the turbulent motions, k_h the horizontal wave number, and $E_{vz}\left(k_{\text{h}},z,t\right)={\widehat{v}}_z\left(k_{\text{h}},z,t\right){\widehat{v}}_z^{\ast }\left(k_{\text{h}},z,t\right)/2$ the vertical velocity spectrum. This definition is motivated by the strong correlation between υ_z(τ_R = 1) and the emerging vertical bolometric intensity. The pressure scale height is computed as H_P = −〈1/∂_z ln(〈P〉_z)〉_t, with 〈−〉_z denoting horizontal averaging computed at constant geometrical height. For all models we consider a uniform grid spacing, ∆x = ∆y = ∆z, and n_x × n_y × n_z = 750 × 750 × n_z grid cells, where n_x, n_y, and n_z are the number of grid cells in the horizontal directions x and y and in the vertical direction z, respectively, such that we have approximately the same number of grid cells per granule and per pressure scale height, as detailed in Table 1. This is to guarantee similar Reynolds and magnetic Reynolds numbers, Re and Re_m, in the four simulations and, consequently, to disentangle the impact of the numerical set-up from the impact of the stellar parameters when drawing conclusions in the following.

To verify that the flow parameters are similar in the four models, we computed the effective viscosity and magnetic diffusivity, v and η, following the methodology detailed in Riva & Steiner (2022) and we used Re = 〈υ_rms〉_t〈L_int〉_t/〈v〉_t and Re_m = 〈υ_rms〉_t〈L_int〉_t/〈η〉_t, vertically averaged between ln(〈P〉_z,t/P⁰) = 3 and ln(〈P〉_z,t/P⁰) = 1, which is the region where the drive amplitude of the dynamo is largest, as discussed in Sect. 4. Here υ_rms is the root mean squared (rms) velocity evaluated at constant geometrical height. The magnetic Prandtl number is defined as Pr_m = Re_m/Re. The results are displayed in Table 2. Although Re and Re_m show a mild but nearly monotonic increase from model K8V to F5V, we see that the change in Pr_m is rather modest, smaller than 25% for the four models.

From the values reported in Table 1, we see that both the granular size and the pressure scale height increase with effective temperature. However, when comparing L_gran and H_P(z = 0) to the results reported in Table 1 of Salhab et al. (2018), we see that g has a significant impact on stratification properties and on the size of the granules. Larger gravitational accelerations lead to smaller granules and shallower atmospheres, and vice versa for smaller g. This is an expression of the fact that H_P ∝ T_eff/g.

We use periodic boundaries in both horizontal directions. The top boundary is open for fluid flows and outward radiation, with the density decreasing exponentially in the boundary (ghost) cells outside the domain, whereas magnetic fields are forced to be vertical. Moreover, we impose ∂_tv_z = −c_s∂_zv_z and ∂_zv_h = 0, with c_s being the local sound speed and v_h the horizontal velocity. At the bottom we set an open boundary condition for fluid variables, enforcing ∂_zv = 0, vanishing horizontally averaged vertical mass flux, and a control on the specific entropy of the inflowing material to fine tune the resulting effective temperature. The parameters C_sChange = 0.1, C_pChange = 0.3, and C_lin = C_sqrt = 0.003 were used to reduce deviations of the entropy, the pressure, and the vertical velocity from the horizontal means, as detailed in Freytag et al. (2012) and Riva & Steiner (2022). For the magnetic field at the bottom boundary, we set ∂_zB = 0 in locations of outflowing plasma, whereas we imposed B_x = C_B, B_y = 0, and ∇ ⋅ B = 0 for locations of inflowing plasma, with C_B a constant input value discussed in Sect. 3.1. More details on the boundary conditions used in CO⁵BOLD are given in Freytag (2017).

All simulations were started from the non-magnetic runs of Salhab et al. (2018). The log g values and grids were adapted to match the parameters displayed in Table 1 and the simulations were evolved for several thousand seconds to damp out fluctuations introduced by the interpolations. A seed magnetic field B_z = 1 mG was then added to these relaxed non-magnetic states. From the resulting initial states, we ran four simulations, one for each model, until saturation of the magnetic energy (first saturation), imposing C_B = 0. However, to investigate the impact of the bottom boundary condition on the results, after the first saturation stage we continued the simulations with C_B ≠ 0, until reaching another saturated state (second saturation). We consider C_B ≠ 0 a less stringent condition that matches reality more closely as it allows for some replenishment of the magnetic energy from outside the computational domain. We note that using C_B ≠ 0 implies that the magnetic energy is not produced by SSD action alone, but also by amplification of the magnetic field advected into the domain across the bottom boundary. However, although magnetic energy enters the domain in upflow regions, the magnetic energy leaving the domain in downflow regions overcompensates this contribution, thus ensuring a net Poynting flux leaving the domain at the bottom boundary. A more detailed discussion of the impact of these boundary conditions on the results is given in Sect. 5.

Fig. 1 Temporal profiles of magnetic to kinetic energy ratios averaged over the convection zone of the four models K8V, K2V, G2V, and F5V. The round, square, and diamond markers denote phase 1, phase 2, and phase 3, respectively.

3 Synopsis of the numerical runs

3.1 Evolution

The time-traces of the magnetic to kinetic energy ratio, E_m/E_k, are displayed in Fig. 1 for the four models, where 〈−〉_c denotes averaging over the convection zone (z < 0). We can identify three main phases. After a few minutes of quick magnetic amplification during which the magnetic flux is expelled from the granules and is transported to the intergranular space, phase 1 starts, in which the magnetic energy is amplified exponentially by the SSD. We consider the end of phase 1 to be the time when 〈E_m/E_k〉_c = 0.2%, denoted by the round markers in Fig. 1. This low value of 〈E_m/E_k〉_c corresponds to a still negligible Lorentz-force feedback on the fluid motion. After phase 1, a balance is reached between the amplification of the magnetic field and losses of magnetic energy due to magnetic energy dissipation and to pumping out of the magnetic field across the bottom boundary. Therefore, the simulations reach a statistical steady-state in phase 2 (the first saturated state), whose initial and final points are denoted by the square markers in Fig. 1. Finally, setting C_B = 110, 110, 155, and 140 G for models K8V, K2V, G2V, and F5V, respectively, a second saturation state is reached, phase 3, whose initial points are denoted by the diamond markers in Fig. 1. This corresponds to setting $C_B=0.08B_{\text{eq\hspace{0.17em}dyn}}^{\text{bottom }}$, with $$B_{\text{eq\hspace{0.17em}dyn }}={(4\pi \rho )}^{1/2}{\upsilon }_{\text{rms }}$$(6)

the dynamic equipartition field strength and ρ and v_rms evaluated at the bottom of the computational domain. The condition $C_B=0.08B_{\text{eq\hspace{0.17em}dyn}}^{\text{bottom }}$ is motivated by the fact that max(B/B_{eq dyn}) ≃ 8% in the convection zone during phase 2; the maximum value is attained in the deep convection zone, before decreasing to smaller values when approaching the bottom boundary. Therefore, this boundary condition mimics the presence of a SSD operating below the bottom boundary, which generates a magnetic field of strength B ≃ 0.08B_{eq dyn}. We note that the initial points of phases 2 and 3 are determined by eye inspection of the curves displayed in Fig. 1.

Table 3 gives more details on simulation times and output frequencies for each model, with N_output denoting the number of outputs for the full simulation box (the mass density, internal energy, and velocity and magnetic fields for every grid cell). The emerging vertical bolometric intensity (I_bol), some horizontally averaged quantities, as well as mass density, internal energy, and velocity and magnetic fields at z = 0 were saved every 10 s. Here phase 1 starts at t = 0 and not after the initial expulsion of the magnetic flux from the granules. This is justified by the fact that phase 1 is used in the following as a non-magnetic proxy, and we do not expect the magnetic field to have any noticeable impact on the plasma flow during these first few minutes of the simulations.

3.2 Morphology

Figure 2 provides an overview of the basic properties of the simulation results, from the coolest model in the top (K8V) to the warmest in the bottom (F5V). The first three columns display representative time instances of the vertically directed bolometric intensity, normalized to its mean, during phases 1, 2, and 3, respectively; the last column shows the vertical field strength at z = 0 corresponding to the time instances of the third column. The standard picture of bright granules, where the plasma raises to the surface, and dark intergranular lanes, where the plasma flows back into the convection zone, is clearly visible for all models. Moreover, consistent with findings in Beeck et al. (2011), Freytag et al. (2012), and Salhab et al. (2018), among others, the visible continuum intensity contrast of the granular pattern increases from the coolest to warmer models. This increasing intensity contrast with increasing T_eff is related to the increase in the relative temperature fluctuations with increasing T_eff, as discussed in Beeck et al. (2013) and Salhab et al. (2018), among others. During phase 3 (third column), bright features appear in locations of strong magnetic field. We discuss these features further and the impact that magnetic fields have on plasma flows in Sect. 5.

3.3 Notation

In the following, all the numerical results are evaluated during one of the three main phases discussed above; the initial and final times are reported in Table 3. Uncertainties on the results are computed as the 1σ deviation in temporal fluctuations of the quantities. Horizontal averages are computed either at constant geometrical height (i.e. constant z, denoted in the following as 〈⋅〉_z) or over a surface of constant Rosseland opacity (denoted in the following as 〈⋅〉_τ). The former is adopted for displaying quantities as a function of height and for analysis in the convection zone, where stratification is nearly adiabatic and largely independent of the radiation field, whereas the latter is favoured for analysis near the optical surface, consistently with the considerations in the Appendix of Beeck et al. (2013). We note that the time t in the index of the angle brackets has the meaning of averaging in time, for example 〈⋅〉_zt.

To have a meaningful comparison between models of different spectral types, we used ln(〈P〉_z,t/P⁰) as vertical coordinate instead of z, with the time average of the gas pressure evaluated during the kinematic phase. Quantities on surfaces of constant optical depth were obtained with linear interpolations in ln(τ_R). To obtain quantities on constant τ_R we used the full simulation boxes. On the other hand, quantities averaged at constant geometrical height, as well as quantities related to the bolometric intensity, were computed from the data saved each 10 s.

4 Kinematic phase

From Fig. 1, it is clear that the SSD growth rate is different in the four models. By fitting 〈E_m/E_k〉_c during the kinematic phase with an exponential function, we obtain the energy e-folding times τ_E ≃ 2700, 2200, 2600, and 3500 s for models K8V, K2V, G2V, and F5V, respectively. These correspond to the magnetic energy growth rates, γ_E = 1/τ_E, displayed in the fifth column of Table 2.

To better understand these results, we use the ideal induction equation, Eq. (3), to express the time evolution of the magnetic energy density as $${\partial }_t{E}_{\text{m}}=\frac{1}{4\pi }\left[-\upsilon \cdot \nabla \frac{B^2}{2}+B\cdot (B\cdot \nabla v)-B^2\nabla \cdot \upsilon \right].$$(7)

Figure 3 displays 〈∂_tE_m/E_m〉_t for the four models (solid lines, left axis). This provides an indication of the drive amplitude of the dynamo and where it is most efficient in amplifying the magnetic energy. It results that ∂_tE_m/E_m peaks just beneath the optical surface, at approximately ln(〈P〉_z,t/P⁰) = 2 for all models, with values that are about a factor two higher for the K2V model than for the F5V model. The values for the K8V and G2V models lie between those of the other two models, in agreement with the magnetic energy e-folding times reported above.

The large differences among the four models of the SSD growth rates (with the K2V and F5V models having the lowest and highest γ_E, respectively) nearly vanish when normalizing the timescales in the system to the turnover time, τ, reported in Table 1. This can be seen from the last column of Table 2, where γ_Eτ is almost independent of the spectral type of the model. More details on the temporal evolution of the models in normalized units are given in Appendix A, where the growth of magnetic to kinetic energy for each model is shown with time units normalized by the corresponding turnover times. Here, we define the turnover time as τ = 1/ max_z(〈v_z,rms〉_t/〈L_int〉_t), where max_z denotes taking the maximum value over the vertical coordinate.

The small variation in γ_Eτ can be understood as follows. It is well known that magnetic energy amplification by SSD is mostly due to the second term on the right-hand side of Eq. (7), B ⋅ (B ⋅ ∇v), called the stretching term, whereas SSD action is suppressed by compression and magnetic diffusivity (see e.g. Batchelor & Taylor 1950). The stretching term is expected to be proportional to the velocity shear term, |∇v|, which is also displayed in Fig. 3 (dashed lines, right axis). Since |∇v| ≈ 1/τ, we have a nearly constant SSD driving term in normalized units.

The impact of compression, which corresponds to the last term of Eq. (7), gets larger when the Mach number, M, of the simulation increases (see e.g. Martins Afonso et al. 2019; Seta & Federrath 2021). At ln(〈P〉_z,t/P⁰) = 2, where the efficiency of the SSD is greatest, we have M = 〈v_rms/c_s,rms〉_t = 0.18, 0.27, 0.34, and 0.48 for models K8V, K2V, G2V, and F5V, respectively, with c_s,rms the rms value of the sound speed. Hence, the effect of compression is smallest for K8V and largest for F5V. However, the impact of compression on the SSD operation is marginal for the present models, as discussed in the following, and therefore it is not expected to noticeably affect the value of γ_Eτ. Finally, concerning the effect of magnetic diffusivity, in normalized units one would expect higher growth rates for larger Re_m and larger Pr_m, and vice versa for smaller Re_m and smaller Pr_m (Brandenburg & Subramanian 2005; Käpylä et al. 2018). However, the changes in Re_m and in Pr_m are rather mild (see the third and fourth columns of Table 2). Consequently, γ_Eτ is almost constant in the present models.

Energy transfer analysis. To gain a deeper insight into the mechanism driving SSD action, we performed an energy transfer analysis similar to Pietarila Graham et al. (2010), Rempel (2014), and Riva & Steiner (2022). In particular, we write the spectral magnetic energy as $$E_{\text{m}}\left(k_{\text{h}}\right)=\widehat{B}\left(k_{\text{h}}\right)\cdot {\widehat{B}}^{\ast }\left(k_{\text{h}}\right)/(8\pi ),$$(8)

whose time evolution can be decomposed into the three components $${\partial }_t{E}_{\text{m}}\left(k_{\text{h}}\right)=T_{\text{ma}}\left(k_{\text{h}}\right)+T_{\text{ms}}\left(k_{\text{h}}\right)+T_{\text{mc}}\left(k_{\text{h}}\right),$$(9)

where $$T_{\text{ma}}\left(k_{\text{h}}\right)=-{\widehat{B}}^{\ast }\left(k_{\text{h}}\right)\cdot \widehat{\upsilon \cdot \nabla B}\left(k_{\text{h}}\right)/(8\pi )+\text{ c.c. },$$(10) $$T_{\text{ms}}\left(k_{\text{h}}\right)={\widehat{B}}^{\ast }\left(k_{\text{h}}\right)\cdot \widehat{B\cdot \nabla \upsilon }\left(k_{\text{h}}\right)/(8\pi )+\text{ c.c. },$$(11) $$T_{\text{mc}}\left(k_{\text{h}}\right)=-{\widehat{B}}^{\ast }\left(k_{\text{h}}\right)\cdot \widehat{B\nabla \cdot \upsilon }\left(k_{\text{h}}\right)/(8\pi )+\text{ c.c}.$$(12)

are the energy transfers to and from the magnetic energy due to advection, stretching, and compression of the magnetic field, respectively. Here $\widehat{\text{-}}$ denotes the Fourier transform, * the complex conjugate, and c.c. the respective complex conjugate expression. A detailed analysis of these terms close to ln(〈P〉_z,t/P⁰) = 2 (not shown here for conciseness) reveals that the term T_ms is responsible for amplifying magnetic energy at all scales, with a peak at λ ≃ 0.1L_gran − 0.2L_gran. At the same time, magnetic field advection (T_ma) removes magnetic energy from large scales, breaking it down into smaller scale magnetic fields. In particular, T_ma + T_ms almost exactly cancels out at large scales, thus resulting in magnetic field amplification at small scales only. The contribution of T_mc is almost negligible, being significantly smaller than the other two terms at all scales. This behaviour is observed for all four of the models, and is consistent with magnetic energy being amplified by SSD action, as discussed in Riva & Steiner (2022).

From our analysis, we conclude that the stretching of magnetic fields by velocity shear is the main drive of SSD action, in agreement with the literature. Moreover, the behaviour observed for the growth time of the magnetic energy, τ_E, can to a great degree be explained by the amplitude of the velocity shear, |∇v|, which is inversely proportional to the turnover time given in Table 1. Thus, in the convection zone, |∇v| is highest for the K2V model and lowest for the F5V model, as can also be seen in Fig. 3. This behaviour is not due to a different operation of the SSD, but simply to the different turnover times in SI units for the different models. The energy transfer analysis proves that the magnetic field amplification occurs on small scales for all models, as expected for a small-scale dynamo. For completeness, the spatial spectrum of the kinetic and magnetic energy of all models are shown in Appendix A with the wave number normalized to the integral length scale of each individual model and temporally averaged over the kinematic phase.

Fig. 2 Maps of bolometric intensity and vertical magnetic field representative of the three simulation phases. The time instants of the bolometric intensity during phases 1 (first column), 2 (second column), and 3 (third column) are normalized to their respective mean intensities. The fourth column displays the vertical magnetic field at z = 0 during phase 3 at the same time instant as for the third column. The rows (from top to bottom) correspond to the spectral types K8V, K2V, G2V, and F5V. The red arrows in the K2V model indicate non-magnetic bright points, whereas the red contours in the third column denote the MBFs detected by the algorithm (see Sect. 5.3 and Appendix B).

Fig. 3 Time-averaged profiles of 〈∂tEm〉z/〈Em〉z (solid lines, left axis, smoothed for clarity with a moving average of width 10∆z) and computed during the kinematic phase, and of 〈|∇v|〉z,t (dashed lines, right axis) for the four models K8V, K2V, G2V, and F5V. The black dotted line denotes the geometrical height z = 0. To the left of this line is the convection zone; the photosphere is to the right.

5 Saturation phases

As discussed in Sect. 2.3, the bottom boundary condition C_B = 0 prevents replenishment from outside the domain of the magnetic field pumped out of the domain across its bottom boundary. Consequently, the saturation level reached during phase 2 can be considered to be a lower limit for the magnetic energy reached in our models. As is discussed in the following, the amount of magnetic energy generated with this bottom boundary condition does not have a noticeable impact on the plasma flows and on the radiative field, consistently with the solar results in Rempel (2014). For this reason, we only briefly discuss the results of phase 2, while putting more emphasis on the results of phase 3, which displays MBFs similar to those observed on the Sun (see Sect. 5.3).

5.1 Characteristics of the magnetic field

From Fig. 1 it is clear that the simulations attain a magnetic energy that is significantly lower in phase 2 than in phase 3 for all spectral types. In addition, during both saturation phases, the magnetic-to-kinetic energy ratio, averaged over the convection zone, is independent of the stellar model, with 〈〈E_m〉_c/〈E_k〉_c〉_t = 1.3–2.5% during phase 2 and 〈〈E_m〉_c/〈E_k〉_c〉_t = 17–21% during phase 3.

A deeper insight into the properties of the magnetic fields generated by SSD action provides the first four rows of Table 4, which display the rms of B_z and B (rows 1 and 2, respectively) and the horizontal averages of |B_z| and B (rows 3 and 4, respectively) at τ_R = 1 and time averaged over phase 3, with B = |B|. Additionally, the left panel of Fig. 4 shows 〈B〉_z,t as a function of ln(〈P〉_z,t/P⁰) (solid lines). Interestingly, the magnetic field strengths are independent of the spectral type. These amplitudes are similar to the solar value found in Rempel (2014), which for model O16bM is 〈B〉_τ=1,t = 160 G. Moreover, we obtain magnetic field strengths similar to those found in Bhatia et al. (2022), which are in the range 〈B〉_τ=1,t = 103–188 G. The main difference with the results in Bhatia et al. (2022) are the magnetic field strengths for stars of spectral type F, which in our case are slightly lower than for the G2V model and not about 40% higher as in Bhatia et al. (2022). The stronger magnetic field for the F3V model with respect to the cooler models is linked to a larger fraction of kinetic energy converted into magnetic energy, which we do not observe for our F5V model (as discussed later in this section). This could be related to the different log ɡ and T_eff values used for the F5V and F3V models. Moreover, comparing 〈B_rms(τ_R = 1)〉_t to the values reported in Beeck et al. (2015, their Table 2), we see that our results are similar to the results of the magnetic simulations that started with a uniform vertical seed field of B_z = 20 G.

To investigate the reason for the magnetic field strength and for the magnetic pressure, P_m = B²/(8π) (row 8 of Table 4), being approximately independent of the spectral type, the left panel of Fig. 4 also displays the amplitudes of the dynamic equipartition field, Eq. (6), as a function of height for the four models during the kinematic phase (dashed lines). We see that the equipartition fields are of similar strengths for all models at all heights. This means that higher densities, and hence higher gas pressures in cooler models are compensated for by higher velocities in hotter models (see rows 7 and 11 of Table 4), with the two effects approximately balancing each other.

In addition, at the optical surface, B_{eq dyn} decreases by approximately 6% from phase 1 to phase 3 for all models (see row 5 of Table 4), which corresponds to a decrease in the kinetic energy density of approximately 11–12% (see rows 9 and 10 of Table 4, with P_dyn = ρv²/2). This is mainly related to a decrease in the rms velocity from phase 1 to phase 3 (row 11 of Table 4), and it means that the saturation mechanism of the SSD is similar for all spectral types and that it is independent of the growth rates discussed in Sect. 4. The magnetic to dynamic equipartition field ratios show very similar amplitudes for all models, with B/B_{eq dyn} ≃ 0.2 near the surface and increasing up to 0.3 near the bottom boundary.

The right panel of Fig. 4 displays the ratio of the horizontal, B_h, to the absolute vertical magnetic field strength for the four models, time averaged over phase 3. For an isotropically distributed magnetic field, one would expect this ratio to be π/2. For the present simulations, the magnetic fields are slightly more vertical than horizontal in the convection zone, in agreement with the findings reported in Bhatia et al. (2022). On the other hand, magnetic fields are close to isotropy at the optical surface, and become more horizontal in the photosphere, in particular for the coolest models. A peak in the solar photosphere of the horizontal-to-vertical magnetic field ratio is consistent with the findings in Steiner et al. (2008), where it is shown that the predominance of horizontal fields in the photosphere can be considered a consequence of the flux expulsion process by overshooting convection. This is also consistent with the results in Sect. 4.6 of Rempel (2014). We note that the strong decrease in the B_h-to-|B_z| ratios near the top boundary is due to the boundary conditions, which impose a vertical magnetic field at the top of the domain.

Finally, Fig. 5 displays the flux-based probability density function (PDF), PDF_Φ(|B_z|) = |B_z|PDF(|B_z|)/Φ_tot (solid lines, left axis), and the corresponding empirical distribution function, ${\text{EDF}}_{\text{Φ}}\left(\left|B_z\right|\right)={{{\displaystyle \int }}^{\text{}}}_0^{\left|B_z\right|}{\text{PDF}}_{\text{Φ}}\left(\left|B_z^{\prime }\right|\right)\text{d}{B}_z^{\prime }$ (dashed lines, right axis), for the four models at z = 0 (left panel) and τ_R = 1 (right panel) during phase 3, where Φ_tot is the absolute magnetic flux trough the area L_xL_y and PDF(|B_z|) the probability density function of |B_z|. PDF_Φ(|B_z|)dB_z represents the probability that a given flux element, dΦ, has a flux density in the interval [|B_z|, |B_z + dB_z|], whereas EDF_Φ(|B_z|) gives the contribution to the total absolute magnetic flux from magnetic fields weaker than |B_z| (Steiner 2003). The peaks of PDF_Φ(|B_z|) are located at approximately the same flux density, about |B_z| = 15–25 G, for all four models, and both at z = 0 and τ_R = 1. At z = 0, the four PDF_Φ(|B_z|) display similar behaviour up to |B_z| = 300 G, while they extend to stronger magnetic fields for cooler models than for hotter models. On the other hand, at τ_R = 1, the profiles of PDF_Φ(|B_z|) show similar behaviour also at high flux densities, although the peaks of the two coolest models at |B_z| ≃ 20 G are higher than for models G2V and F5V. Overall, these results are similar to the findings in Beeck et al. (2015) for the simulations that started with B₀ = 20 G. We note that our simulations do not display a second peak of PDF_Φ(|B_z|) at |B_z| ≃ 1000–2000 G, as observed for simulations that started with an initial magnetic field higher than or equal to 50 G (see e.g. Beeck et al. 2015; Salhab et al. 2018). Likewise, the EDF_Φ(|B_z|) also shows similar profiles for all models, except for large |B_z| at z = 0. The EDF_Φ(|B_z|) in Fig. 5 show that about half of the absolute magnetic flux has a field strength |B_z| < 200 G for all models and that kG fields contribute at z = 0 to Φ_tot by 9%, 6%, 7%, and 1% for models K8V, K2V, G2V, and F5V, respectively. Thus, more than 90% of the magnetic flux through the stellar surface has a field strength of less than 1 kG. More about strong magnetic flux concentrations can be found in Sect. 5.3.

At the origin of the transfer of energy from the kinetic to the magnetic energy reservoir is the work done by the Lorentz force on the plasma. Considering the second saturation phase and integrated over the entire convection zone of the model, this work is actually substantial compared to the radiative flux from the stellar surface. Indeed, ${{{\displaystyle \int }}^{\text{}}}_{z_{\text{bottom }}}^{z=0}v\cdot (j\times B)\text{d}z/F_{\text{rad}}=$ −0.57, −0.36, −0.35, and −0.34 for model K8V, K2V, G2V, and F5V, respectively, where j = ∇ × B/(4π) and F_rad is the radiative flux from the stellar surface. Thus, the energy flux transferred by the Lorentz force of the SSD is of the order of the stellar luminosity, although the exact number strongly depends of the numerical set-up.

Fig. 4 Magnetic field structure for the four models. Left: temporal and horizontal averages of the magnetic (solid lines) and dynamic equipartition (dashed lines) fields. Right: horizontal to vertical magnetic field ratio. Temporal averages of the dynamic equipartition and of the magnetic field strengths are evaluated during phases 1 and 3, respectively. The black dotted lines denote the geometrical height z = 0; to the left of the line is the convection zone, and to the right the photosphere.

Fig. 5 Magnetic flux distributions. Flux-based probability density functions (solid lines, left axis) and corresponding empirical distribution functions (dashed lines, right axis) as a function of |Bz| at z = 0 (left panel) and τR = 1 (right panel) for the four models during phase 3. The insets in the top left corners of the larger panels display zoomed-in images of PDFΦ at high magnetic field strengths.

5.2 Impact of SSD action

Row 12 of Table 4 displays the visible continuum rms intensity contrast of the granular pattern, defined as $$c_{\text{rms}}=\sqrt{{\langle {\langle {\left(\frac{I_{\text{bol}}-{\langle I_{\text{bol}}\rangle }_{\text{h}}}{{\langle I_{\text{bol}}\rangle }_{\text{h}}}\right)}^2\rangle }_{\text{h}}}_{}}{\rangle }_t,$$(13)

with 〈I_bol〉_h denoting the average of the bolometric intensity over the horizontal plane. Consistently with the discussion of Fig. 2 in Sect. 3.2, c_rms increases from the coolest to warmer models. Row 12 of Table 4 also shows that c_rms decreases from phase 1 to phase 3 for all models, and that this decrease is greater with increasing T_eff. To investigate this behaviour, we display in Fig. 6 the PDF of the relative intensity fluctuations, $\delta {\tilde{I}}_{\text{bol }}=\left(I_{\text{bol}}-{\langle I_{\text{bol}}\rangle }_{\text{h},t}\right)/{\langle I_{\text{bol}}\rangle }_{\text{h},t}$, for the four models during phases 1 (solid lines) and 3 (dashed lines). The appearance of bright features in the intergranular lanes in phase 3 (see red contours in the third column of Fig. 2) manifests in Fig. 6 as a decrease in the PDFs at the low end, except for the K8V model, and the emergence of an extra tail at the high end of the intensity distributions (solid vs. dashed lines). However, the decrease in c_rms from phases 1 to 3 is not related to changes in the tails of the PDFs, which have much lower realization rates than the core of the distribution, but rather to changes near the centre of the PDFs (see inset in the top right corner of Fig. 6).

The lack of MBFs in Fig. 2 during phase 2 (cf. second and third columns) is due to the low saturation level reached by the magnetic energy with C_B = 0. Indeed, 〈B〉_τ=1,t ≃ 29–53 G during phase 2, which is about a factor of three smaller than when using C_B ≠ 0 (cf. row 4 of Table 4). As also shown by Rempel (2014) for the solar case, this is not enough magnetic flux to form kG flux concentrations. We note that tiny bright dots within the intergranular lanes occur in all spectral types during phases 1 and 2, and are particularly conspicuous in the K2V model. Two of them are indicated by red arrows in Fig. 2. These dots are not of magnetic origin, but they form from swirling downdrafts due to the low pressure at the centre of the vortices, giving rise to non-magnetic bright points (Vögler et al. 2005; Calvo et al. 2016).

The impact of the magnetic fields generated by SSD action on the granulation is demonstrated in Fig. 7, which displays the integral length scale of vertical velocity fluctuations, Eq. (5), averaged over phases 1 and 3 (solid and dashed lines, respectively). The convection zone (left of the black dotted line) shows a 5%–10% reduction in L_int from phase 1 to 3 for all models. As explained in Bhatia et al. (2022), a reduction in kinetic energy because of SSD action, and a corresponding reduction in con-vective velocities (rows 9–11 of Table 4), leads to a reduction in the pressure scale height. This, in turn, explains the decrease in L_int from phase 1 to 3. On the other hand, differences at the optical surface are much smaller. Above the optical surface (to the right of the black dotted line), models K8V, K2V, and G2V also display a reduction in L_int, whereas the hottest model, F5V, shows an opposite trend.

Fig. 6 Probability density functions of the relative intensity fluctuations $\delta {\tilde{I}}_{\text{bol }}=\left(I_{\text{bol}}-{\langle I_{\text{bol}}\rangle }_{\text{h},t}\right)/{\langle I_{\text{bol}}\rangle }_{\text{h},t}$. The solid and dashed lines denote phases 1 and 3, respectively. The inset in the top right corner shows a zoomed-in image near the maximum of the PDFs.

5.3 Characteristics of magnetic bright features

5.3.1 Properties of magnetic flux concentrations

During phase 3, all the stellar models show bright features located in the intergranualar space, where strong magnetic fields form magnetic flux concentrations (see the red contours in the third column of Fig. 2). These correspond to the MBFs, also known as magnetic filigree, magnetic bright points, or G-band bright points, observed with high-resolution instruments on the solar disk (Dunn & Zirker 1973; Mehltretter 1974; Berger & Title 1996). We now characterize the MBFs observed in our models during phase 3.

As the name suggests, MBFs harbour strong magnetic fields and are brighter than their surrounding space. Therefore, to identify MBFs in the simulations and perform a statistical analysis, we developed an algorithm based on the following conditions: pixels corresponding to an individual MBF have I_bol and |B_z|(z = 0) larger than some given thresholds and they are located within a closed contour of sufficiently large area. For the intensity threshold we choose I_bol > 〈I_bol〉_h,t; for the absolute vertical magnetic field we use |B_z|(z = 0) > 875.7, 704.2, 700.3, and 584.2 G for model K8V, K2V, G2V, and F5V, respectively, and we discard contours enclosing fewer than 20 pixels, with one pixel corresponding to an area of ∆x∆y. A discussion of these criteria and more details on the algorithm are given in Appendix B.

The first row of Table 5 gives the number of MBFs per snapshot found with the algorithm. We see that the K2V and G2V models display the highest number of MBFs, while the K8V model counts less than half their MBFs. Additionally, Fig. 8 displays the PDFs of the area of MBFs. The distributions are very similar for the four models and approximately follow an exponential function, with more MBFs of small areas and fewer MBFs of large areas. Since the four models have similar MBF area PDFs, but different numbers of MBFs per snapshot, the fraction of total area covered by MBFs changes among the four models, with the K2V and G2V models having the largest fraction (see second row of Table 5). This is also in qualitative agreement with the third column of Fig. 2, which displays more pixels covered by MBFs for the K2V and G2V models than for the other two models.

Rows 3–10 of Table 5 give the magnetic properties of the MBFs, where 〈q〉_MBF denotes first computing the quantity q = B_z,rms, max(|B_z|), and B_h|B_Z| over the area of each MBF individually, and then averaging over the sample of all the identified MBFs. We note that averaging over each MBF individually weights the smaller MBFs more than the larger ones (see Fig. 8). Nevertheless, this does not affect the trends discussed in the following¹. Rows 3 and 5 of Table 5 show that the rms and the maximum of the absolute vertical magnetic field decreases as T_eff increases at z = 0. Here, max(|B_z|) is evaluated as the average of |B_z| over the nine pixels closest to the maximum value of |B_z| of each MBF. At τ_R = 1 the differences in the rms and maximum values of |B_z| for the four models are smaller than at z = 0, but the decreasing trend with increasing T_eff is still significant (rows 4 and 6 of Table 5). This is different from what is found for magnetic simulations that started from a homogeneous vertical field of 50 G in Steiner et al. (2014) and Salhab et al. (2018), where the field strength at τ_R = 1 is independent of the spectral type.

Rows 7 and 8 of Table 5 show the fraction of vertical magnetic flux and magnetic energy contained by the MBFs with respect to the full horizontal plane at z = 0. Even if their area is less than 1% of the horizontal domain, MBFs contain a significant fraction of the magnetic flux and magnetic energy. Moreover, as is the case for the number density of MBFs, their contribution to the magnetic flux and magnetic energy is also larger for the K2V and G2V models than for K8V and F5V.

For solar simulations carried out with the MURaM code under similar conditions, it is found that kG flux concentrations contribute about 9–16% to the total magnetic energy, depending on the resolution (Rempel 2014). This is consistent with the present results.

Rows 9 and 10 of Table 5 display the ratio of B_h to |B_z| in the MBFs at z = 0 and τ_R = 1, respectively. At z = 0, |B_z| is much larger than B_h, indicating that the magnetic field is preferentially vertical in the MBFs. This is not true at τ_R = 1, in particular for the MBFs of the K2V and G2V models, for which the magnetic field is instead isotropically distributed.

Fig. 7 Integral length scale of vertical velocity fluctuations for the four models. The solid and dashed coloured lines denote phases 1 and 3, respectively. The black dotted line denotes the geometrical height z = 0; to the left of the line is the convection zone, and to the right the photosphere.

Fig. 8 Probability density functions of MBF areas for the four models. The black dotted line denotes the 20-pixel threshold, below which MBFs are discarded. The dashed lines are the extension of the PDF as if small MBFs were retained.

5.3.2 Discussion

We now investigate the reason for the decreasing field strength of MBFs with increasing effective temperature, showing that it is the result of two competing effects: the increasing near-surface opacities and increasing evacuation efficiency with increasing T_eff. As discussed, for example, in Salhab et al. (2018), in the present range of temperatures and gas pressures an increase in effective temperature leads to an increase in opacity, which in turn gives higher gas pressures at the optical surface for cooler in comparison to hotter models (row 7 of Table 4). Consequently, the thermal equipartition field strength, $B_{\text{eq,th }}=\sqrt{8\pi P_{\text{gas}}}$, decreases with increasing T_eff (row 6 of Table 4), thus leading to flux concentrations of lower magnetic field strength at higher effective temperatures. However, as the decreasing trend of B_eq,th with increasing T_eff is much more pronounced than for the magnetic field strength (cf. row 6 of Table 4 with rows 3 and 5 of Table 5), the opacity effect alone is not able to explain the results in Table 5.

To further investigate the trend displayed by the magnetic field strength of MBFs with increasing T_eff, we compute the mass-density, temperature, and pressure contrasts of the MBFs at z = 0 and τ_R = 1 as c_q(z = 0) = (〈q〉_MBF − 〈q〉_z=0)/〈q〉_z=0 and c_q(τ_R = 1) = (〈q〉_MBF − 〈q〉_τ=1)/〈q〉_τ=1, respectively, with q = ρ, T, or P_gas. These are displayed in rows 11–16 of Table 5. The mass density contrast at z = 0 is significantly lower than zero for all models, and it strongly decreases with increasing T_eff. This indicates that the degree of evacuation of MBFs strongly increases with increasing T_eff. As discussed in Salhab et al. (2018), in the convective collapse framework (Parker 1978; Webb & Roberts 1978; Spruit 1979), a more efficient evacuation is associated with a more efficient MBF formation process, hence leading to a higher ratio of the magnetic to thermal equipartition field strengths (cf. rows 5 and 6 of Table 5 with row 6 of Table 4). The higher efficiency of MBF formation at high T_eff is due to the larger superadiabaticity of hot models with respect to the cooler ones (see e.g. Salhab et al. 2018; Beeck et al. 2013; Magic et al. 2013). At τ_R = 1 the MBF evacuation is much lower than at z = 0 (cf. rows 11 and 12 of Table 5). This is due to the fact that the strong magnetic field in MBFs, in combination with horizontal pressure balance, leads to a downward shift of the optical surface with respect to the surrounding atmosphere. Therefore, the τ_R = 1 surface is highly corrugated and it is deeply depressed within MBFs, thus reaching higher density values. The compensating effect of the evacuation process explains the milder dependence of the magnetic field strength of MBFs on T_eff, in contrast to the strong dependence of B_eq,th on T_eff. The two competing effects are probably also at the origin of the increasing trend for A_MBF/A_tot, Φ_MBF/Φ_tot, and E_m,MBF/E_m,tot for models K8V, K2V, and G2V, followed by a decrease in these three quantities when further increasing T_eff from 5766 K to 6506 K (see rows 2, 7, and 8 of Table 5).

The strong magnetic field in MBFs mitigates the horizontal heat exchange responsible for thermal equilibrium, resulting in colder temperatures at 𝓏 = 0 in MBFs with respect to the surrounding plasma at the same geometrical height (see row 13 of Table 5). On the other hand, as the visible surface within MBFs is located below the 𝓏 = 0 level, the MBF temperature at τ_R = 1 is higher than the normal surface temperature (see row 14 of Table 5), although it is lower than the averaged temperature at the geometrical height level of the depressed surface. As temperature deviations are much smaller than mass density deviations in MBFs (cf. rows 13–14 with rows 11–12 of Table 5), c_P ≈ c_ρ both at 𝓏 = 0 and τ_R = 1 (cf. rows 15–16 with rows 11–12 of Table 5).

Finally, we report the average shift of the τ_R = 1 surface in MBFs with respect to the 〈τ_R〉𝓏 = 1 surface in row 17 of Table 5, denoted WD in reminiscence of the Wilson depression of sunspots. It results that WD strongly increases with increasing T_eff. Although the increase in the pressure scale height with increasing T_eff (see Table 1) could explain this trend, WD increases more rapidly than H_P, consistently with results in Steiner et al. (2014) and Salhab et al. (2018). This, again, is due to the increased evacuation efficiency of hotter spectral models with respect to the cooler ones.

6 Conclusions

In the present paper we numerically investigated the action of a small-scale dynamo operating in the near-surface layers of four cool main-sequence stars of spectral types K8V, K2V, G2V, and F5V. The numerical domains were set up to ensure similar Reynolds and magnetic Reynolds numbers for all models. Starting from a seed field of 1 mG, all simulations display the operation of a SSD close beneath the stellar surface, which transfers energy from the kinetic to the magnetic reservoir, mostly through stretching of magnetic fields by velocity shear. As the four models have different turnover times, the SSD growth rates are different for the four stellar models; the highest and lowest growth rates are for the K2V and F5V models, respectively. When normalizing the time in the system by the turnover time, the growth rates are very similar, in agreement with the fact that the magnetic Reynolds and Prandtl numbers are similar for the four models.

When a balance between amplification of the magnetic field and losses of magnetic energy is attained, the simulations reach a statistical steady-state. The SSD field strength reached in the saturation phase is similar in all models, irrespective of the kinematic SSD growth rate. However, the bottom boundary condition strongly affects the amount of magnetic energy that can be sustained during the saturation phase. In particular, the magnetic field amplified by the SSD has a noticeable impact on plasma flows only if a partial replenishment of the magnetic energy that is pumped out across the bottom boundary is allowed to enter from outside the computational domain. In this case, 1721% of the energy is transferred from the kinetic energy to the magnetic energy reservoir. This leads to a reduction in the convective velocities that is enough to reduce the convective horizontal length scale in the convection zone by 5%–10%, although not at the visible surface. In this situation, the SSD magnetic fields display properties very similar to the simulations that started with a uniform initial field of B_𝓏 = 20 G. The rms of the absolute magnetic field strength at the visible surface is approximately 200 G, independent of the spectral type. Moreover, while the magnetic field is close to isotropic near the τ_R = 1 level, it is predominantly horizontal further up in the photosphere for all models.

In the saturation stage including replenishment of pumped out magnetic fields, the SSD action also results in the formation of strong kG magnetic flux concentrations at the surface, which appear as MBFs on the intensity map of all stellar models. Strong magnetic flux concentrations are twice as numerous in the K2V and G2V models as in the K8V and F5V models. It results in an absolute vertical magnetic flux located in the MBFs, normalized to the total absolute vertical magnetic flux, that is about twice in the K2V and G2V models with respect to the K8V and F5V models. However, more than 90% of the magnetic flux through any of these stellar surfaces has a field strength less than 1 kG. The vertical component of the field strength of the MBFs at optical depth unity increases from approximately 1.0 to 1.6 kG with decreasing temperature from F5V to K8V due to the competing effects of opacity and formation efficiency.

Overall, these results give a deeper insight into the operation of the SSD in the near-surface layers of cool main-sequence stars and the self-organization of magnetic fields into kG flux concentrations. In addition, the present simulations provide the basis for investigating the impact of the SSD self-generated MBFs on interferometric and planet transit observations, thus extending the work in Chiavassa et al. (2010) and Chiavassa et al. (2017) to magnetic phenomena. In the future, it will be important to broaden the parameter space of the numerical models by including more stars, where the SSD magnetic fields could have stronger effects, such as stars hotter than those in the F5V model (see e.g. Bhatia et al. 2022). An ideal case could be Atype stars with shallow convection zones. This would mitigate the impact of the bottom boundary condition on the simulation results as the convection zones can be fully contained in the numerical domain.

Acknowledgements

This work was supported by the Swiss National Science Foundation under grant ID 200020_182094. Part of the numerical simulations were carried out on Piz Daint at CSCS under project IDs s1059, s1172, and u14, and the rest was carried out on the HPC ICS cluster at USI. We sincerely thank the anonymous referee for the careful reading of the manuscript and his/her comments, which helped us to improve the quality of the paper.

Appendix A Results in normalized units

Fig. A.1 Temporal profiles of magnetic-to-kinetic energy ratios averaged over the convection zone, with the time coordinate normalized to the turnover times of the four models.

Fig. A.2 Kinetic and magnetic energy spectra (solid and dashed lines, respectively) time averaged over the kinematic phase and normalized to their surface area for the four models K8V, K2V, G2V, and F5V. The x-axis is the wave number normalized to the integral length scale. It is also visible that the magnetic energy resides mostly at small scales, as expected for a SSD.

When comparing the physical properties of the four simulations, it can be useful to normalize the time and length scales in the models to the characteristic scales in the system, which are the turnover time and the integral length scale, respectively. Since the efficiency of the SSD is greatest close to ln(〈P〉_𝓏,t/P⁰) = 2 (see Sect. 4), the integral length scales discussed in the following of this Appendix are computed at this position (see also Fig. 7), whereas the turnover times τ are given in Table 1.

In Fig. A.l we display the temporal profiles of the magnetic-to-kinetic energy ratio averaged over the convection zone, with the time of each model normalized to the corresponding turnover time. During the kinematic phase, the four curves are very close to each other. This is an expression of the small deviations of the normalized SSD growth times among the four models (see the last column of Table 2). It is also in agreement with the mild variations in the flow parameters of the four models (see second, third, and fourth columns of Table 2).

Figure A.2 displays the kinetic (solid lines) and magnetic (dashed lines) energy spectra time averaged over the kinematic phase. At each time instant, each spectrum is normalized to its integral area before performing the time average. The length scales are normalized to the integral length scale of each model. It results that, in normalized units, the spectra of the four models are very close to each other. This is again an expression of the small variations in the flow parameters among the four models (see second, third, and fourth columns of Table 2).

Appendix B Magnetic bright feature recognition

Fig. B.1 Bivariate distribution, in logarithmic scale, of Ibol and |B𝓏 | during phase 3. The rows (from top to bottom) correspond to the spectral types K8V, K2V, G2V, and F5 V, respectively. The black vertical and red horizontal lines denote 〈Ibol〉h,t and B𝓏,tbreshoid, respectively.

Magnetic flux concentrations that form during phase 3 lead to the emergence of bright features in the intergranular space. Thus, it is natural to set a condition on |B_𝓏| and on I_bol in order to identify magnetic bright features (MBFs). More precisely, combining the two-dimensional maps of |B_𝓏|(𝓏 = 0) with those of I_bol, we define an individual MBF as the area with |B_𝓏| and I_bol above a given threshold and within a closed contour.

The |B_𝓏|(z = 0) and I_bol thresholds are defined as follows. Figure B.1 displays the logarithm of the bivariate distribution of I_bol and |B_𝓏|(𝓏 = 0) during phase 3 for the four models. An extended tail of bright magnetic pixels can be clearly seen in the top right part of each panel. These are the MBFs. Thus, for the I_bol threshold, we simply use the averaged bolometric intensity, 〈I_bol〉_h,t, indicated by the black vertical lines in Fig. B.1. Furthermore, we set the |B_𝓏| threshold, B_{𝓏,threshold}, as the local minimum of the distribution of |B_𝓏|(I_bol > 〈I_bol〉_h,t), indicated by the red horizontal lines in Fig. B.1. This gives B_{𝓏,threshold} = 875.7, 704.2, 700.3, and 584.2 G for models K8V, K2V, G2V, and F5V, respectively.

This approach for the MBF identification is different from that used by Salhab et al. (2018), where one threshold on the magnetic field alone is used, and with the same B_{𝓏,threshold} for all models. This for two reasons. First, the PDF_Φ(|B_𝓏|) values in Fig. 5 do not display a secondary peak at high magnetic flux densities as in Beeck et al. (2015) or Salhab et al. (2018). On the other hand, when considering |B_𝓏|(I_bol > 〈I_bol〉_h,t) alone, one can clearly discern a strong field component in Fig. B.1. This justifies the use of a threshold on I_bol for identifying an MBF. Second, the strong field component at τ_R = 1 of the present SSD simulations shows a stronger dependence on T_eff than for simulations that started with an initial magnetic field higher than or equal to 50 G (see e.g. Beeck et al. 2015; Salhab et al. 2018). This justifies the use of a B_{𝓏,threshold} that depends on T_eff.

As an additional constraint, we also require MBFs to occupy an area equal to or larger than 20 pixels, with the area of a pixel being ∆x∆y. We do this to discard small transient features. Tests with a reasonable range of thresholds (between 0 and 30 pixels) proved that the precise choice of the threshold does not impact the trends discussed in Sect. 5.3.

References

All Tables

All Figures

	Fig. 1 Temporal profiles of magnetic to kinetic energy ratios averaged over the convection zone of the four models K8V, K2V, G2V, and F5V. The round, square, and diamond markers denote phase 1, phase 2, and phase 3, respectively.
In the text

Fig. 2 Maps of bolometric intensity and vertical magnetic field representative of the three simulation phases. The time instants of the bolometric intensity during phases 1 (first column), 2 (second column), and 3 (third column) are normalized to their respective mean intensities. The fourth column displays the vertical magnetic field at z = 0 during phase 3 at the same time instant as for the third column. The rows (from top to bottom) correspond to the spectral types K8V, K2V, G2V, and F5V. The red arrows in the K2V model indicate non-magnetic bright points, whereas the red contours in the third column denote the MBFs detected by the algorithm (see Sect. 5.3 and Appendix B).

In the text

Fig. 3 Time-averaged profiles of 〈∂tEm〉z/〈Em〉z (solid lines, left axis, smoothed for clarity with a moving average of width 10∆z) and computed during the kinematic phase, and of 〈|∇v|〉z,t (dashed lines, right axis) for the four models K8V, K2V, G2V, and F5V. The black dotted line denotes the geometrical height z = 0. To the left of this line is the convection zone; the photosphere is to the right.

In the text

Fig. 4 Magnetic field structure for the four models. Left: temporal and horizontal averages of the magnetic (solid lines) and dynamic equipartition (dashed lines) fields. Right: horizontal to vertical magnetic field ratio. Temporal averages of the dynamic equipartition and of the magnetic field strengths are evaluated during phases 1 and 3, respectively. The black dotted lines denote the geometrical height z = 0; to the left of the line is the convection zone, and to the right the photosphere.

In the text

Fig. 5 Magnetic flux distributions. Flux-based probability density functions (solid lines, left axis) and corresponding empirical distribution functions (dashed lines, right axis) as a function of |Bz| at z = 0 (left panel) and τR = 1 (right panel) for the four models during phase 3. The insets in the top left corners of the larger panels display zoomed-in images of PDFΦ at high magnetic field strengths.

In the text

	Fig. 6 Probability density functions of the relative intensity fluctuations $\delta {\tilde{I}}_{\text{bol }}=\left(I_{\text{bol}}-{\langle I_{\text{bol}}\rangle }_{\text{h},t}\right)/{\langle I_{\text{bol}}\rangle }_{\text{h},t}$. The solid and dashed lines denote phases 1 and 3, respectively. The inset in the top right corner shows a zoomed-in image near the maximum of the PDFs.
In the text

	Fig. 7 Integral length scale of vertical velocity fluctuations for the four models. The solid and dashed coloured lines denote phases 1 and 3, respectively. The black dotted line denotes the geometrical height z = 0; to the left of the line is the convection zone, and to the right the photosphere.
In the text

	Fig. 8 Probability density functions of MBF areas for the four models. The black dotted line denotes the 20-pixel threshold, below which MBFs are discarded. The dashed lines are the extension of the PDF as if small MBFs were retained.
In the text

	Fig. A.1 Temporal profiles of magnetic-to-kinetic energy ratios averaged over the convection zone, with the time coordinate normalized to the turnover times of the four models.
In the text

	Fig. A.2 Kinetic and magnetic energy spectra (solid and dashed lines, respectively) time averaged over the kinematic phase and normalized to their surface area for the four models K8V, K2V, G2V, and F5V. The x-axis is the wave number normalized to the integral length scale. It is also visible that the magnetic energy resides mostly at small scales, as expected for a SSD.
In the text

	Fig. B.1 Bivariate distribution, in logarithmic scale, of Ibol and |B𝓏 | during phase 3. The rows (from top to bottom) correspond to the spectral types K8V, K2V, G2V, and F5 V, respectively. The black vertical and red horizontal lines denote 〈Ibol〉h,t and B𝓏,tbreshoid, respectively.
In the text