© The Authors 2026

1 Introduction

Magnetism is ubiquitous in outer convection zones of cool stars. Evidence of magnetism appears not just in indicators of chromospheric activity (Baliunas & Vaughan 1985), stellar flares (Vasilyev et al. 2024) and photometry (Strassmeier 2009), but also in polarization (Kochukhov 2016). From an observational point of view, magnetic fields can have a wavelength-dependent signature in the emergent stellar spectrum and limb darkening (Ludwig et al. 2023; Kostogryz et al. 2024), which affects the interpretation of stellar light curves. This is especially important when it comes to characterizing exoplanets and their atmospheres around magnetically active stars (de Wit & Seager 2013). Accounting for stellar magnetism is of prime importance in the field of exoplanet detection and characterization. In the post-Kepler era, the detections of exoplanets from photometry has shot up to thousands of candidates. These candidates are then targeted by ground-based radial velocity (RV) campaigns to accurately characterize the properties of the exoplanet candidate. The biggest obstacle to precise and accurate characterization is stellar variability (Crass et al. 2021). In addition, transmission spectroscopy, which is used to characterize exoplanetary atmospheres and the wavelength-dependent radius of exoplanets, is also affected by stellar activity (Rackham et al. 2023). Hence, there is a requirement for models that capture the radiative MHD properties of faculae and spots as accurately as possible.

The most widely used models for interpreting stellar spectra are simplified, hydrodynamic, 1D models that either do not account for magnetism at all (Kurucz 1970; Husser et al. 2013; Gustafsson et al. 2008) or implement only the surface radiative effects of starspots (i.e., flux blocking) in a parameterized manner (Somers et al. 2020). A more consistent approach to modeling magnetic fields in a mixing length framework was carried out by Feiden & Chaboyer (2013).

Early studies of non-radiative compressible 3D MHD convection have focused on the types of convective patterns that form over a range of imposed vertical field strengths and the transitions between these patterns with varying field strength. Chapter 8 of Weiss & Proctor (2014) gives an overview of these studies in a dimensionless setup. To summarize, it was found that the convective pattern depends strongly on the strength of the imposed vertical magnetic field. The magnetic structures ranged from point- and ribbon-like concentrations at weaker field strengths to large-scale flux-separated weak-field convecting regions and strong-field regions with suppressed convection (Weiss et al. 2002) at stronger field strengths. In the presence of the strongest fields, regular convection gets replaced by over-stable small-scale oscillations.

In the last four decades or so, the sophistication of numerical models has increased in tandem with available computing resources, allowing realistic 3D radiative MHD simulations of near-surface convection with magnetic fields. Recently, simulations of faculae-like magnetic fields (Beeck et al. 2015a; Salhab et al. 2018) as well as spots (Panja et al. 2020; Bhatia et al. 2025) on cool main-sequence stars have been conducted with a box-in-a-star approach to study the effect of magnetism on structure as well as the wavelength dependence of intensity (Norris et al. 2023; Smitha et al. 2025; Witzke et al. 2022).

The present work is an update as well as an extension to the facular studies of Beeck et al. (2015a). We focus more on the effects of field on convection and structure in this paper, in addition to the effects on stellar photospheric velocities, intensity, and field distribution as in Beeck et al. (2015a). To our knowledge, no other study has probed this aspect for realistic stellar surface simulations, except for the analysis in Bhatia et al. (2024), which was restricted to SSD simulations. In addition, our reference cubes are saturated SSD atmospheres already containing a basal level of magnetism, instead of being purely hydrodynamic, and the simulations are run with a new equation of state (Irwin 2012). We note that both the equation of state as well as the opacities for radiative transport use abundances from Asplund et al. (2009), which is not the case for Beeck et al. (2015a).

The paper is organized as follows: in Sect. 2, we describe the setup as well as the data analysis choices, in Sect. 3, we describe the changes in thermodynamic and convective structure, and in Sect. 4, we provide explanations for these changes and their implications for our understanding of magnetized near-surface stellar convection. Finally, in Sect. 5 we summarize the main results.

2 Methods

In this paper, we considered the stellar types F3V, G2V, K4V, and M0V and conducted near-surface simulations using a box-in-a-star approach. We used the MURaM radiative MHD code (Vögler et al. 2005; Rempel 2014) to conduct all the simulations. We employed the small-scale dynamo (SSD) models described in Bhatia et al. (2022, 2024) as initial conditions. Briefly, the simulation boxes are periodic in the horizontal directions and have an open (symmetric) bottom boundary for flows and magnetic fields. Two parameters determine the stellar type: the effective temperature, T_eff, and the surface gravity, g_surf. The latter was set as a constant in these simulations, whereas the former was determined indirectly by specifying the entropy of inflows at the bottom boundary. In addition, specifying the pressure at the bottom boundary sets the height of the photosphere in the simulation domain. For more details of the models and setup, we refer to the publications cited above.

We then imposed uniform vertical fields of 100 G, 200 G, 300 G, and 500 G on the existing SSD simulations. After inserting the fields, the simulations were run for a few hours up to the point where variations in root-mean-square velocity throughout the box and the effective temperature had stabilized. This allowed the field to organize itself naturally.

Up to this point, the simulations were run with gray opacity (no wavelength dependence). The next step was to introduce a more realistic treatment of radiative transfer. The basic idea is to compute the average opacity for specific “bins” in wavelength and a reference optical depth (τ_ref), instead of solving the radiative transport equation for a huge number of wavelength points. For the present set of simulations, 12 opacity bins were used, with thresholds in wavelength and τ₅₀₀ nm formation height similar to those used in Table 2 of Beeck et al. (2012) for the STAGGER simulations. The bins themselves were averaged from an opacity distribution function (Kurucz 1979; Witzke et al. 2021) following the approach of Vögler et al. (2004). The choice of 12 bins was made based on the excellent match obtained between solar surface simulations using MURaM and solar observations of the irradiance spectrum, photospheric line strengths, shapes (including bisectors), as well as limb darkening (Witzke et al. 2024). A description of the binning procedure is provided in Appendix B. For the imposed field simulations, we used the final SSD opacity bins. At this stage, the F-star simulations were doubled in resolution in the vertical direction to better resolve the relatively sharp gradients near the optical surface, since we plan to use these cubes for future spectral synthesis.

All simulations were run for 60 minutes after converging to a statistical steady state. Only the output produced during this extra runtime was analyzed. The analysis was conducted over snapshots with a cadence of 90 seconds. We used the same notation as the one described in Bhatia et al. (2024) for temporal and horizontal means: overline $\overline{q}$ denotes temporal averages and angular brackets ⟨q⟩ denote spatial averages of the analyzed quantities. The transparent error bars, unless otherwise specified, are 1σ standard deviation in averages over time. The horizontal averages were computed on geometric slices (unless otherwise stated) and are plotted against the average number of pressure scale heights computed for the SSD case (log(p/p₀), where p₀ is the mean pressure at z_⟨τ⟩=1). Unless otherwise specified, τ refers to the Rosseland mean optical depth. We provide a color palette in Table 1 to show the list of colors used for plots showing multiple simulations.

3 Results

The effect of magnetic fields on the near-surface structure, convection, and bolometric intensity is varied and depends on stellar type as well as field strength. Stronger fields (500 G simulations) result in the formation of field concentrations that appear dark in intensity (micro-pores). This is due to the stabilizing effects of strong vertical magnetic field on convection (Gough & Tayler 1966, Eq. (1.1)). Weaker fields (300 G and less) result in the formation of bright points and fibrils for the F, G, and K stars. Magnetic fields evacuate intergranular lanes, reducing the density and, accordingly, the opacity, and allowing radiation to escape from surrounding hotter upflows in deeper layers (also termed the “hot-wall” effect) (Spruit 1976).

Cooler stars (M0V) have higher surface densities and, for the same field strength, are more likely to form darker concentrations. Cooler stars have a lower contrast between upflows and downflows in density and temperature (Beeck et al. 2013, Fig. 12). Since the corresponding optical depression in the τ = 1 layer in downflows is also reduced, the hot-wall effect is not so effective. These results are quite consistent with previous studies of faculae-like simulations of near-surface convection in cool dwarfs (Beeck et al. 2015a; Salhab et al. 2018). Fig. 1 shows a snapshot of intensity, field, and velocity for the 300 G cases for all investigated spectral types, illustrating the aforementioned points. We also show similar snapshots for the 100 G case and the 500 G case in Figs. C.2 and C.4, respectively, for comparison. Plots for the other field strength cases are available in Appendix C. We note that despite imposing a moderate to strong vertical magnetic field in the various simulations, there are indications of weak, opposite polarity fields, especially at granular edges near the strong field concentrations (slight red areas in middle columns of Fig. 1 and subsequent snapshots). This is consistent with observations of solar plage regions (Buehler et al. 2015) and is a result of near-surface convection-driven recirculation of magnetic field, where the strong shear between up- and downflows churns and drags up field to the surface. On the other hand, there are also regions of mixed polarity further away from the strong flux concentrations, which could either be from recycled fields or from a SSD.

Despite the aforementioned effects on bolometric intensity in the presence of magnetic fields, the variation in the effective temperature calculated from the angle-averaged bolometric radiative flux is within 1% of the reference SSD simulations (see Fig. 2). The trends in the changes of effective temperature with field strength are consistent with previous studies; for example, Beeck et al. (2015a, Table 1) and Norris et al. (2023, Table 2). We note that the effective temperature calculated here assumes that the whole stellar surface has the same field strength. It is not meant to represent the actual effective temperature of a facular region on a star’s surface, which would have a center-to-limb variation.

In addition to the previously established results in Beeck et al. (2015a), we report here the effect of faculae-like fields on thermodynamic structure and convection, and offer probable explanations for the changes observed compared to the reference SSD simulations. We take up the 300 G case simulations as a typical strong facular region for reference for these comparisons. Other field strengths show similar behavior. The magnitude of changes compared to SSD simulations depends on the field strength.

Fig. 1 Snapshot of the bolometric intensity I (in 1010 ergs per centimeter squared per second) for the 300 G case (left column) and the corresponding vertical magnetic field, Bz (in kilogauss) (middle column), and the vertical velocity, vz (in kilometers per second) (right column), for spectral types (from top to bottom) F, G, K, and M, respectively. The field and velocity plots correspond to the τ = 1 optical surface.

Fig. 2 Teff as a function of field strength, normalized by the SSD Teff. The shaded regions represent changes of ±1%. The vertical bars at each point are 1σ average variation in Teff. The legend contains the average SSD Teff in Kelvin for each stellar type.

3.1 Changes in the thermodynamic structure

The horizontally averaged plots of density, temperature, and pressure show changes relative to the SSD case, with the magnitude of change depending on the field strength (Fig. 3). To examine the changes in detail, in Fig. 4, we plot the relative difference between the 300 G cases and their corresponding SSD runs for a given horizontally averaged quantity, q, as (q − q_SSD)/q_SSD against a pressure scale axis, log(p/p₀), as described in Sect. 2. We computed the deviations on a geometric height axis and used the SSD log(p/p₀) axis to plot them. We see that both the density, ρ, and gas pressure, p_gas, are reduced around the stellar surface in the facular simulation for all the stars, with the magnitude and depth of change somewhat proportional to T_eff. Near the bottom of the box (left side of the plots), the changes are negligible. The temperature profile shows a similar dip below the surface (log(p/p₀) ~ 0.5–1), but at the surface (dotted vertical lines), the change becomes almost negligible. Possible explanations for this trend are discussed in Sect. 4.

We note that the bottom panel of Fig. 4 seems to indicate that the average temperature in faculae is always lower than the corresponding SSD temperature, even above the optical surface (log(p/p₀) < 0), which seems to contradict solar observations (Solanki 1986; Solanki & Brigljevic 1992; Buehler et al. 2019, Table 2). This is a consequence of horizontal averaging over geometric slices. The observed temperatures are averages over τ surfaces. In the top panel of Fig. 5, we take the full 3D temperature cubes from single snapshots for each star and field case, and plot the average of surfaces of equal Rosseland optical depth, τ. The change in the temperature profiles (Fig. 5, bottom panel) shows that the F, G, and K stars have a higher average temperature above log₁₀ τ ~ −1, with the magnitude of change somewhat proportional to field strength and T_eff. We also note that the enhancement in temperature relative to the SSD case seems to begin deeper in the atmosphere for cooler stars. The M-star simulations do not show any such trend and the difference in the τ stratification of temperature between the different simulations is rather small. For the hotter stars, the trend is consistent with solar observations.

Fig. 3 Density, ρ (top), gas pressure, pgas (middle), and temperature, T (bottom), for all magnetic field strengths and stellar types. The horizontal axis is the number of pressure scale heights below the surface computed for the SSD case for each star (the left side of the plot is toward the bottom boundary and the right side is toward the top). The vertical dotted line marks the height at which ⟨τ⟩ = 1 in each case.

Fig. 4 Changes in density, ρ (top), gas pressure, pgas (middle), and temperature, T (bottom), for the 300 G case relative to the SSD case for all four stellar types in percent (%). The error bands represent the 1σ variation in time for the horizontal means. All subsequent error bands have the same meaning.

Fig. 5 Temperature, T, averaged over equal Rosseland tau τ surfaces instead of geometric horizontal slices for each stellar and magnetic field case. Absolute profiles (top) and difference in the imposed field temperature profiles relative to the SSD cases (bottom).

Fig. 6 Changes in square of vertical velocity, $v_z^2$ (top), and in the ratio of horizontal to vertical velocities, vh,rms/vz,rms (bottom), for the 300 G case relative to the SSD case for all four stellar types.

3.2 Changes in the velocities

The presence of a magnetic field is known to suppress convective velocities in general (Weiss et al. 2002). This is usually attributed to the Lorentz force acting on plasma motions through the (v × B) × B term in the MHD momentum equation and quenching the overturning motion. In Fig. 6, we see the expected decrease in convective velocities (considered here in terms of $v_z^2$) for all the stars, with the decrease roughly around 20%, all the way down to the bottom of the box. Near the surface, the trends are reversed. Here, since the mode of energy transport changes from convective to radiative, the interpretation of the increase in $v_z^2$ above the surface is not so straightforward. The ratio of horizontal to vertical velocities shows a small decrease in the deeper convection, indicating a corresponding decrease in horizontal velocities as well. This is well correlated with a decrease in the horizontal extent of convection in the presence of magnetic fields (Bhatia et al. 2024, Sect. 4.1).

Fig. 7 PDF of logarithmic bolometric intensity for all cases (top) and the difference in PDFs between the imposed field cases and the SSD case (bottom).

3.3 Bolometric intensity

The PDF of disk-center bolometric intensity in Fig. 7 shows trends consistent with previous simulations of near-surface stellar convection with faculae-like fields (Salhab et al. 2018; Beeck et al. 2015b). Briefly, the distributions can be interpreted in terms of two peaks, corresponding to granules and intergranular lanes. The presence of magnetic fields adds to the dark (left) flank through formation of dark pore-like structures. This effect increases with increasing field strength and is strongest for the M star. For the K-star (and, to a lesser degree, the other stars), the bright (right) flank is also affected by the formation of bright points and bright filigree in the intergranular lanes. More generally, the shape of the two peaks also seems to be slightly affected, with the average granular brightness slightly decreasing (peak shifting to the left, visible in the residuals plot below as an increase in density on the left side of the peak granule intensity and a decrease on the right side). We note that these results hold only for the bolometric intensity at the disk center and the effect of magnetic fields on intensity distribution at different viewing angles as well as in different spectral ranges is expected to be different (Norris et al. 2017, 2023). We plan to investigate this more generally in a subsequent paper.

3.4 Morphology of magnetic field concentrations and associated flows at and below the stellar surface

The simulations exhibit a diversity of flux tube and flux sheet sizes as well as magnetic field strengths. To examine these further, we plot the horizontal and vertical cuts of B_z and v_z in Fig. C.1 for all the simulations with 300 G imposed magnetic field. We note that, for the F star, most of the magnetic field seems to concentrate into a wide and vertically contiguous intergranular lane going all the way to the bottom boundary (for the z = 8.53 Mm slice, at x = 10 Mm and y = 13.5 Mm, roughly), whereas for the cooler stars, the lanes become narrower and more spread out. For all the stars, most of the magnetic field tends to collect in the intergranular lanes that are contiguous almost all the way to the bottom of the simulation domain. The red contours in the field maps in Fig. C.1 enclose regions where B_z > B_eqp, with $B_{\text{eqp}}=\sqrt{8\pi ⟨p_{\text{gas}}⟩}$ (with averaging done over the horizontal simulation domain). The blue contours enclose regions where the magnetic energy in B_z is stronger than the average kinetic energy at a given height. We note that the magnetic field for the F star stays relatively strong compared to the gas pressure for the upper part of the box, whereas for the M star this is true only near the surface. For all stars, however, the magnetic energy remains stronger than the kinetic energy in magnetic concentrations throughout the simulated convection zone. We also note that the magnetic field concentrations become quite fragmented in the lower half of the simulation boxes for all stars.

In the velocity cuts, downflows are associated with the edges of field concentrations. We note that, even in the large field concentrations, the v_z flow maps show a network of small-scale motions that seem convective in nature. This is especially apparent for the F star. The corresponding bolometric intensity map for the F star (Fig. 1, top left panel) also shows intensity variations somewhat correlated with these velocity structures. These flows are a version of flows associated with umbral dots and are, in general, a signature of convection in strong vertical fields (Narayan & Scharmer 2010). A zoom in of the F-star pore-like feature is presented in Fig. 8, showing multiple clear cusplike structures in upflows embedded in a network of relatively weak downflows. The magnetic field structure shows a weakening of field in the cusp, compared to the surroundings. These features are correlated with upward bumps in the optical surface and bright features in the bolometric intensity, as was seen in previous simulations (Schüssler & Vögler 2006), as well as observations (Riethmüller et al. 2008) of umbral dots.

Fig. 8 Zoom-in of the F-star pore-like structure in Fig. 1 showing umbral-dot-like features. Bolometric intensity of the zoomed-in region (left) with vertical cuts in convective velocity, vz (middle), and the vertical magnetic field, Bz (right), along the y and x axes. The dashed colored lines represent the locations of the cuts and the solid green line represents the τ = 1 surface. The bottom of the simulation box corresponds to z = 0 for reference.

3.5 PDFs of surface magnetic fields

The PDFs of the magnitude of magnetic field at the optical (τ = 1) surface (Fig. 9, top panel) show a rather similar distribution for the facular simulations in the hecto-gauss (hG) range, with a slight trend with T_eff (hotter stars have slightly higher fraction of stronger hG fields). For comparison, the PDF of SSD simulations are plotted as solid lines. They show the same trend until roughly 500 G. For higher values they diverge from the facular fields. In the kilogauss (kG) range, the distribution of the 300 G simulations is similar till roughly 2.5 kG. The strongest fields seem to follow a trend (hotter stars have higher fraction of stronger fields) with stellar type, except for the M star. All stellar types with an imposed field seem to show field strengths in excess of the corresponding gas plus hydrodynamic turbulent pressure equipartition value (Bhatia et al. 2024). We note that inferring anything about the relation between field strength and pressure on a τ surface is difficult since the surface itself is corrugated, with stronger fields lying geometrically deeper in intergranular lanes due to evacuation of plasma and associated opacity changes.

To compare against gas pressure, we must look at the PDF of fields at the geometric surface corresponding to ⟨τ⟩ = 1 (bottom panel), the idea being that at any given height, the simulation box is expected to be in horizontal pressure balance on average. Here we see that the kG fields show a distinct peak for the imposed field cases, with the peak at the pressure equipartition field strength for the F and G star. For the K and M star, the peaks are less pronounced and lie at decidedly sub-equipartition values. The peak field strength is related to the efficiency of the convective collapse mechanism (Spruit & Zweibel 1979). Cooler dwarfs with higher surface gravity tend to be more stable against convective collapse and are, accordingly, less likely to develop field concentrations with equipartition field strength (Rajaguru et al. 2002). Coming back to the top plot, we note that stronger fields are observed in this plot because they are observed at a deeper geometrical depth, where p_gas is higher than at z_⟨τ⟩=1.

Fig. 9 Top: PDFs of the strength of magnetic field B for the 300 G (dashed) and the SSD (solid) cases, calculated for the τ = 1 surface. Bottom: same as above, but calculated for the geometric surface z⟨τ⟩=1 corresponding to the height at which ⟨τ⟩ = 1. The vertical lines correspond to the pressure equipartition field, Beqp, for both field cases for all stars, calculated as described in Bhatia et al. (2024, Fig. 2 and Sect. 4.1).

Fig. 10 Plasma β as a function of log(p/p0) for the 300 G (dashed) and SSD (solid) cases for all stars.

4 Discussion

4.1 Factors affecting the temperature stratification

The presence of magnetic fields primarily affects the thermodynamics by introducing an additional component to the horizontal pressure balance from magnetic pressure p_mag = B²/8π. The main consequence of the first effect is a reduction in gas pressure at the surface. In the regions where β = p_gas/p_mag ~ 1–10, this effect becomes significant. In Fig. 10, we see that for the 300 G case for the F star, this is the case near the optical surface. The decrease in p_gas is simply balanced by extra p_mag to maintain horizontal pressure balance. For the cooler stars, since β is significantly higher at the same depths, this effect is quite weak. The changes in density largely follow the changes in pressure.

In addition to the pressure changes, magnetic field also restricts convective velocities, as is noted in Sect. 3.2. Consequently, this is expected to reduce the convective flux reaching the surface. The picture is not so simple, however. The total flux has four components: convective (F_conv), kinetic (F_kin), Poynting (F_poy), and radiative (F_rad) flux (see Appendix A for details). Below the optical surface, F_conv is always positive, F_kin and F_poy are negative on average, and F_rad is negligibly small (Beeck et al. 2013). Near the surface, the mode of energy transport changes from convective to radiative, and the sum of first three terms goes to zero. As the magnetic field evacuates plasma, the drop in density also causes a drop in opacity. This allows radiation to escape from deeper layers. In Fig. 11, it can be seen that this effect is the strongest for the hottest star in the 300 G case, with a noticeable decrease in the total non-radiative flux (relative to the SSD case) already visible up to 2 pressure scale heights below the surface. As the deeper regions lose energy, they become cooler compared to the SSD case (by up to 15% in the F-300 G case). At the surface, however, this difference becomes small again.

The fact that the mode of energy transport starts transitioning from convective to radiative at a larger depth for the 300 G case implies that the plasma in strong field regions becomes stably stratified, i.e., sub-adiabatic, with a positive entropy gradient. In Fig. 12, we plot the area fraction of regions where magnetic pressure is stronger than the local gas pressure and the corresponding strong-field average entropy in these regions. We contrast these profiles to the horizontally averaged entropy profiles for the entire time series in Fig. 13. It is quite apparent that the regions with strong fields occupy a significant area fraction near the surface for the F-300 G case and the corresponding entropy profile has the entropy minimum located deeper in the convection zone compared to the average profile. For the F-SSD case, the area fraction is quite small so the contribution from these regions is negligible. Other stars show a similar trend, but to a lesser degree. This effect is the same as that seen in average entropy profiles of the umbral region in starspot simulations (Bhatia et al. 2025, Fig. 5).

The complete picture then is as follows: the entropy minimum of the strong-field regions in the facular simulation is located geometrically deeper compared to the weak-field regions. This results in a decrease in the horizontally averaged temperature below the surface. Above the optical surface of these strong-field regions, the sub-adiabatic (stable) stratification causes the temperature profile to be flatter than expected when compared to the surrounding weak-field convecting medium. In addition, there is a probable contribution from the hot walls of the surrounding convecting plasma. When compared to the respective SSD simulation, these effects together result in a temperature profile for the facular simulation that is cooler below the optical surface (by up to 15% for the F-300 G simulation) but rather similar above the surface (a change of <5%).

Fig. 11 Fractional non-radiative flux, F; that is, sum of the convective Fconv, kinetic Fkin, and Poynting Fpoy flux normalized by the angleaveraged outgoing bolometric flux for the 300 G (dashed) and SSD (solid) cases.

Fig. 12 Area fraction of regions with $B>\sqrt{8\pi p_{\text{gas}}}(top)$ and the corresponding entropy normalized by the upflow entropy at the bottom boundary, s/⟨sbot⟩up, in the strong-field regions calculated from single cubes for the 300 G (dashed) and SSD (solid) cases.

Fig. 13 Horizontally averaged entropy normalized by the upflow entropy at the bottom boundary, s/⟨sbot⟩up, for the 300 G (dashed) and SSD (solid) cases.

Fig. 14 Kinetic energy equipartition field strength as a function of pressure scale height for the 300 G (dashed) and SSD (solid) cases.

4.2 Factors affecting the surface magnetic field distribution

The magnetic field is, in general, concentrated in downflows, as flux expulsion from convective flows (Parker 1963; Weiss 1966) drives the magnetic field into the intergranular lanes. At the fixed geometrical depth corresponding to the average optical surface, the characteristic field strength of these concentrations is lower for hotter stars (peaks in Fig. 9, bottom panel). This is a consequence of horizontal balance against ambient gas pressure, which decreases with T_eff. Since the total magnetic flux must remain conserved, this is consistent with strong fields at the same depth occupying a larger area fraction for hotter stars (see Fig. C.1). Deeper in the simulation boxes, the field strength falls significantly below the pressure equipartition value and the spatial distribution no longer follows the trend with effective temperature. The blue outlines delineating contours of fields comparable to kinetic equipartition look rather similar between the different stars. This is due to the fact that the mean kinetic energy equipartition value of magnetic fields is rather similar with depth, independent of stellar type or imposed field (Fig. 14).

We note that the field strengths achievable are limited by equipartition with the plasma pressure, which is lower for hotter stars. For the F and (to some extent) the G stars, pressure equipartition values are actually achieved in the strongest flux concentrations (red contours). In the cooler K and M stars, however, the field at the surface is sub-equipartition. This is consistent with the result of Rajaguru et al. (2002), who find that the minimum critical plasma beta (β = β_c) for which convective collapse can be triggered is higher for cooler stars. The typical maximum β can be estimated from the constraint that a flux tube will only form if the initial field is weak compared to plasma flows. Following the same kinetic energy equipartition arguments as in Sect. 4.2 of Spruit (1979), the maximum β in a flux tube must be less than 2/(γM²), where γ is the adiabatic index and M is the Mach number. For a rough estimate, we used the values from Fig. 6 in Beeck et al. (2013) for M (~0.22 for the M star and ~1 for the F star) and took γ as 1.3 for the F star and 1.66 for the M star (Rajaguru et al. 2002, see Table 1). With these values, maximum β varies from ~1.5 for the F star to ~25 for the M star. The critical β_c values were estimated from a bi-variate spline fit of Table 2 of Rajaguru et al. (2002), giving β_c ~ 0.026 for the F star and β_c ~ 8.756 for the M star.

To estimate the saturated state of a flux tube undergoing convective collapse on different stars, we assumed that the relation between the initial and the final β for such a flux tube is relatively constant and has the same slope as in Fig. 2 of Spruit (1979); that is, ${\beta }_{\text{sat.}}\sim {\beta }_{\text{initial}}^a$, with a ~ −1/3. Since β_c is different for different stars, the relation can be written as $$\log {\beta }_{\text{sat.}}\approx (1-a)\log {\beta }_c+a\log {\beta }_{\text{initial}}.$$(1)

In the extreme case in which the flux tube is completely evacuated (β_sat. → 0), the maximum field strength is roughly the gas pressure equipartition value. Hence, the observed field strength can be estimated as $B_{\text{eqp}}/\sqrt{1+{\beta }_{\text{sat}}}$. The minimum achieved value of β_sat. determines the observed field strength. Combining the estimated β_sat. values computed using the above expression and the B_eqp values from Fig. 9 bottom panel, we obtain the maximum field strength of ~1 kG for the F star and ~1.3 kG for the M star. Despite the surface gas pressure of the M star being almost an order of magnitude higher than that of the F star, the field strengths achieved are quite similar and in line with Fig. 9. This rough analysis highlights the importance of considering the efficiency of the convective collapse mechanism for estimating small-scale magnetic field strength on stellar surfaces. We note that we have neglected the effect of the depth of the Wilson depression on the observed field strength, since it is a nontrivial function of opacity, density, and temperature stratification in the flux tube.

The trends in disk-center bolometric intensity and surface field strength distribution (see Sect. 3.3) are largely similar to those obtained in previous such studies, such as those of Beeck et al. (2015a) and Salhab et al. (2018). We note the remarkable similarity of field strength distribution at the optical surface for all the stellar types considered. The hG field strength distribution is expected to depend on the kinetic energy distribution because of a rough balance between kinetic and magnetic energy in the plasma flow. Since the surface kinetic energy is roughly similar at all scales for the stellar types considered (Bhatia et al. 2024, see Fig. 7), it is to be expected that the magnetic field distribution is also rather similar in this range.

The observed kG field values at the surface also depend on their Wilson depression. The depth of the Wilson depression depends on the opacity and density of the plasma in the intergranular lanes (where the opacity itself depends largely on temperature). If we assume that a magnetic flux tube is in temperature equilibrium with its surroundings, we can take β to be constant. In this case, a deeper Wilson depression corresponds to a higher observed field strength. For an F star, since the saturated field strength is basically at pressure equipartition, the degree of evacuation is quite high, and the observed field is in equipartition with gas pressure at a larger depth. This is not the case for the M star, where the degree of evacuation is sub-equipartition, even for the strongest possible field from the discussion in the previous paragraph. The field strengths actually observed end up being quite similar between the stellar cases considered in this paper. We note that this is a qualitative explanation. A possible analysis in the framework of thin flux tube and/or sheet approximation was done previously by Yelles Chaouche et al. (2009).

5 Conclusion

In this study, we investigate the effects of facular magnetic fields on thermodynamic structure, the emergent bolometric intensity, and the distribution of field strengths at the optical surface of cool main-sequence stars ranging from F to M spectral type. We find that, as was expected, strong magnetic fields reduce the density and pressure near the surface due to plasma evacuation in the presence of magnetic fields (Stix 1989, see Chapter 8). We also see a change in the entropy (temperature) structure, indicating a transition to radiative energy transport in strong-field regions deeper in the convection zone as compared to surrounding weak-field regions. While there is a decrease in average temperature just below the surface, the temperature profile above the surface is rather similar to the corresponding SSD simulation.

The morphology of magnetic features throughout the simulated convection zone show a larger and more contiguous area fraction of strong fields for the hotter stars for a variety of imposed average field strengths (Figs. 1, C.2, C.4). This is a consequence of lower ambient pressure on hotter stars. The distributions of field strengths show peaks for strong fields at different values for different stellar types on a constant geometric slice corresponding to the optical surface. When the corrugation of the optical surface is accounted for, the distributions appear rather similar for strong field (>1 kG) values, with the F star showing a slightly higher fraction of strong fields. The largest concentrations on the F and G stars show similarities to solar pores. The convective velocity and the magnetic field structure in the pore-like regions in the F star shows similarities to the convective velocities and the magnetic field structure of solar umbra dot simulations.

The distributions of intensities show clear extended bright and dark flanks in the I_bol PDF for almost all simulations with magnetic fields, compared to the reference SSD case (Fig. 7). The bright flank remains rather similar for different field strengths for a given stellar type, whereas the dark flank becomes stronger with increasing field strength. The simulations exhibit a preference for darker features with decreasing T_eff.

Magnetic fields affect not just the emergent intensity but also the structure and convective properties of plasma in the convection zone. The rapidly expanding field of exoplanet characterization requires precise and accurate modeling of the stellar signal, which itself depends on detailed and accurate simulations of near-surface stellar magnetism. This study aims to fill that gap.

Acknowledgements

The authors would also like to thank Yvonne Unruh for insightful discussions and general help in improving this manuscript. The authors would like to thank the anonymous referee for their careful review, which helped improve the presentation of results in this manuscript. TB is also grateful for access to the supercomputer Cobra at Max Planck Computing and Data Facility (MPCDF), on which all the simulations were carried out. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement Nos. 695075 and 101097844). VW and AIS acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant no. 101118581). TB and AIS also acknowledge support from DFG grant SH1489/1.

References

Appendix A Fluxes

The various fluxes mentioned in Sect. 4 are expressed as, $$F_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}=v_z(E+p)$$(A.1) $$F_{\mathrm{k}\mathrm{i}\mathrm{n}}=v_z(\rho v^2/2)$$(A.2) $$F_{\mathrm{p}\mathrm{o}\mathrm{y}}=v_z(B_x^2+B_y^2)-B_z(B_x{v}_x+B_y{v}_y)$$(A.3)

Here, E is the internal energy, p is the gas pressure, and the other symbols have their usual meanings (v for velocity and B for magnetic field). In the ideal scenario, the horizontally averaged vertical momentum flux in the simulation domain ρv_z should be zero. However, due to the presence of box modes, this term is, in general not zero at all time steps. This affects the calculation of the convective (Eq. A.1) and kinetic (Eq. A.2) fluxes. Even though the average over time averages out the box mode, the calculation of standard deviation over time is misleading in its magnitude. To account for the contribution of the box mode, the net convective flux is calculated as, $$⟨F_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}⟩=⟨v_z(E+p)⟩-⟨\rho v_z⟩⟨(E+p)/\rho ⟩$$(A.4) $$⟨F_{\mathrm{k}\mathrm{i}\mathrm{n}}⟩=⟨v_z(\rho v^2/2)⟩-⟨\rho v_z⟩⟨v^2/2⟩$$(A.5)

The Poynting flux is not strongly affected by this term since there is no ρ component in its expression.

Appendix B Opacity binning

The procedure leading to the selected opacity bins, briefly, is as follows: the formation height (that is, the value of τ_ref where τ_λ = 1 for each λ_i) for each wavelength point in the opacity distribution function (including continuum opacity) is calculated from a 1D atmosphere, which itself is computed from averaging a simulation timeseries horizontally and temporally. These opacity points are then averaged (either as Rosseland mean or Planck mean, depending on whether they form lower or higher in the atmosphere, respectively) over each bin to create a multibin opacity table. The simulation is run further with the new opacity table for enough time for most of the atmosphere to adjust to the new opacities (a couple of granulation timescales).

The distribution of wavelength points between different opacity bins depends on the mean stratification. Some points may fall in a different bin if the mean stratification changes, which can change the average opacity of that bin. The mean stratification couples back to the radiation via the radiative heating/cooling term in the energy equation. To make sure that the stratification that is used to compute the opacity bins is consistent with the structure on which the radiative scheme actually operates on, the procedure is iterated by using the evolved atmosphere (averaged horizontally and over 1–2 granulation timescales) at the end of each previous iteration. The changes in the effective temperature and thermodynamic structure near the surface remain small between the iterations and rapidly converged within 2–3 iterations per star.

Appendix C Additional figures

Fig. C.1 Geometric cuts of Bz (odd rows) and vz (even rows) for the 300 G simulations. First column shows vertical cuts and the columns to the right show horizontal cuts. The locations of the cuts are specified in dashed green lines. The yellow contour is the surface corresponding to τ = 1. The horizontal cuts in the second column correspond to z⟨τ⟩=1. The red contours enclose regions with $B_z>\sqrt{8\pi ⟨p_{\text{gas}}⟩}$. The blue contours enclose regions where Bz > ⟨ρv2/2⟩. Here, the averages ⟨pgas⟩ and ⟨ρv2/2⟩ are computed over the horizontal slice being plotted. The rows correspond, from top to bottom, to F-, G- and M-star, respectively.

Fig. C.2 Same as Fig. 1, but for the 100 G cases.

Fig. C.3 Same as Fig. 1, but for the 200 G cases.

Fig. C.4 Same as Fig. 1, but for the 500 G cases.

Fig. C.5 Same as Fig. 1, but for the SSD cases.

All Tables

All Figures

Fig. 1 Snapshot of the bolometric intensity I (in 1010 ergs per centimeter squared per second) for the 300 G case (left column) and the corresponding vertical magnetic field, Bz (in kilogauss) (middle column), and the vertical velocity, vz (in kilometers per second) (right column), for spectral types (from top to bottom) F, G, K, and M, respectively. The field and velocity plots correspond to the τ = 1 optical surface.

In the text

	Fig. 2 Teff as a function of field strength, normalized by the SSD Teff. The shaded regions represent changes of ±1%. The vertical bars at each point are 1σ average variation in Teff. The legend contains the average SSD Teff in Kelvin for each stellar type.
In the text

Fig. 3 Density, ρ (top), gas pressure, pgas (middle), and temperature, T (bottom), for all magnetic field strengths and stellar types. The horizontal axis is the number of pressure scale heights below the surface computed for the SSD case for each star (the left side of the plot is toward the bottom boundary and the right side is toward the top). The vertical dotted line marks the height at which ⟨τ⟩ = 1 in each case.

In the text

	Fig. 4 Changes in density, ρ (top), gas pressure, pgas (middle), and temperature, T (bottom), for the 300 G case relative to the SSD case for all four stellar types in percent (%). The error bands represent the 1σ variation in time for the horizontal means. All subsequent error bands have the same meaning.
In the text

	Fig. 5 Temperature, T, averaged over equal Rosseland tau τ surfaces instead of geometric horizontal slices for each stellar and magnetic field case. Absolute profiles (top) and difference in the imposed field temperature profiles relative to the SSD cases (bottom).
In the text

	Fig. 6 Changes in square of vertical velocity, $v_z^2$ (top), and in the ratio of horizontal to vertical velocities, vh,rms/vz,rms (bottom), for the 300 G case relative to the SSD case for all four stellar types.
In the text

	Fig. 7 PDF of logarithmic bolometric intensity for all cases (top) and the difference in PDFs between the imposed field cases and the SSD case (bottom).
In the text

Fig. 8 Zoom-in of the F-star pore-like structure in Fig. 1 showing umbral-dot-like features. Bolometric intensity of the zoomed-in region (left) with vertical cuts in convective velocity, vz (middle), and the vertical magnetic field, Bz (right), along the y and x axes. The dashed colored lines represent the locations of the cuts and the solid green line represents the τ = 1 surface. The bottom of the simulation box corresponds to z = 0 for reference.

In the text

Fig. 9 Top: PDFs of the strength of magnetic field B for the 300 G (dashed) and the SSD (solid) cases, calculated for the τ = 1 surface. Bottom: same as above, but calculated for the geometric surface z⟨τ⟩=1 corresponding to the height at which ⟨τ⟩ = 1. The vertical lines correspond to the pressure equipartition field, Beqp, for both field cases for all stars, calculated as described in Bhatia et al. (2024, Fig. 2 and Sect. 4.1).

In the text

	Fig. 10 Plasma β as a function of log(p/p0) for the 300 G (dashed) and SSD (solid) cases for all stars.
In the text

	Fig. 11 Fractional non-radiative flux, F; that is, sum of the convective Fconv, kinetic Fkin, and Poynting Fpoy flux normalized by the angleaveraged outgoing bolometric flux for the 300 G (dashed) and SSD (solid) cases.
In the text

	Fig. 12 Area fraction of regions with $B>\sqrt{8\pi p_{\text{gas}}}(top)$ and the corresponding entropy normalized by the upflow entropy at the bottom boundary, s/⟨sbot⟩up, in the strong-field regions calculated from single cubes for the 300 G (dashed) and SSD (solid) cases.
In the text

	Fig. 13 Horizontally averaged entropy normalized by the upflow entropy at the bottom boundary, s/⟨sbot⟩up, for the 300 G (dashed) and SSD (solid) cases.
In the text

	Fig. 14 Kinetic energy equipartition field strength as a function of pressure scale height for the 300 G (dashed) and SSD (solid) cases.
In the text

Fig. C.1 Geometric cuts of Bz (odd rows) and vz (even rows) for the 300 G simulations. First column shows vertical cuts and the columns to the right show horizontal cuts. The locations of the cuts are specified in dashed green lines. The yellow contour is the surface corresponding to τ = 1. The horizontal cuts in the second column correspond to z⟨τ⟩=1. The red contours enclose regions with $B_z>\sqrt{8\pi ⟨p_{\text{gas}}⟩}$. The blue contours enclose regions where Bz > ⟨ρv2/2⟩. Here, the averages ⟨pgas⟩ and ⟨ρv2/2⟩ are computed over the horizontal slice being plotted. The rows correspond, from top to bottom, to F-, G- and M-star, respectively.

In the text

	Fig. C.2 Same as Fig. 1, but for the 100 G cases.
In the text

	Fig. C.3 Same as Fig. 1, but for the 200 G cases.
In the text

	Fig. C.4 Same as Fig. 1, but for the 500 G cases.
In the text

	Fig. C.5 Same as Fig. 1, but for the SSD cases.
In the text