© The Authors 2026

1. Introduction

The art of magnetohydrodynamic (MHD) sunspot simulations has only recently been developed. In 2007, Heinemann et al. (2007) presented the first slab simulation of a sunspot, followed by similar simulations by Scharmer et al. (2008) and Rempel et al. (2009b). These simulations provided the first insights into the penumbral fine structure. Following a series of initial trials, particularly in (Rempel et al. 2009a), Rempel (2012) managed to run MHD simulations of a circular sunspot that qualitatively reproduced the observed properties of penumbral fine structure. In other words, the structure of the penumbral filaments is comparable to the observational analysis of (Tiwari et al. 2013).

However, the formation of an extended penumbra has been achieved in these MHD simulations by imposing artificial top boundary conditions that make the magnetic field more horizontal at the photospheric level. As shown by (Jurčák et al. 2020), the global magnetic properties of the resulting sunspot are not fully realistic. The magnetic field inclination on the umbral-penumbral boundary is around 60°, which is significantly more horizontal than in the observed sunspots, where the highest values for the inclination are around 45° and occur only in large sunspots with more flux than the simulated sunspot. Furthermore, the vertical magnetic field component, B_z, at the umbral boundary in these simulations is below the value expected from observations, a B_z^stable of around 1.8 kG. In Jurčák et al. (2020), it was also shown that the inclination distribution as a function of the fractional spot radius, as well as B_z at the umbral boundary of observed sunspots, is comparable to an MHD simulation of a sunspot with a potential field as an initial condition. The use of a potential field as an initial condition for sunspot simulations was first proposed by Nordlund (2015), using slab geometry. Based on this idea, we performed an MHD simulation of a circular sunspot, as described by Jurčák et al. (2020). To find the best agreement between the simulated and observed sunspots, we conducted a parameter study of sunspot simulations using potential field initial conditions.

The sunspot simulations of Panja et al. (2021) also adopted potential top boundary conditions. However, the penumbra of their sunspots was wide only between the nearby opposite polarity spots and very narrow in the perpendicular direction. Furthermore, their simulations are, by design, strongly affected by fluting instability, which makes them short-lived. The same authors also created simulations of starspots on other main sequence stars (Panja et al. 2020). Bhatia et al. (2025) simulated starspots, extending the work of Panja et al. (2020), but again forcing the magnetic field at the top to be inclined, resulting in too horizontal fields, at least when comparing the G2V starspot to sunspot observations.

Hideyuki Hotta investigated the effect of higher resolution on sunspot simulations with a potential-field top boundary condition. Matthias Rempel joined this investigation and also analysed what happens when the resolution is changed during the simulation runs. All simulations in their investigation use potential field top boundary conditions. While presentations of preliminary analysis (Rempel & Hotta 2025) are very encouraging, these results have not yet been published in a peer-reviewed journal. Furthermore, the high resolution used in these simulations results in prohibitively high costs.

In solar observations, the penumbra forms around existing flux concentrations (e.g. Schlichenmaier et al. 2010b). Various authors have done flux emergence simulations (e.g. Cheung et al. 2010; Chen et al. 2017; Hotta & Iijima 2020). These simulations show, consistently with solar surface observations, that deeply anchored bundles rise coherently up to around 2 Mm below the surface, to then decay in small parts that rise with individual convective cells. At the surface, they re-coalesce to larger structures that end up forming sunspots (see the ‘tethered-balloon’ model of Spruit 1981). However, they either manage to create a wide enough penumbra at most in some directions or have incorrect flow directions (inflows along the filaments and downflows at the umbral boundary) in most or all penumbral filaments.

In the following, Sect. 2 describes the simulation setups, including the initial conditions, their variables, and the values they take, as well as the data processing methods. Section 3 presents the results, focusing on intensities, surface flows, and magnetic field properties. Section 4 provides a summary and conclusion, including a discussion of how the results relate to observations.

2. Simulation set-ups and processing

Our simulations were run using the MPS/University of Chicago Radiative Magneto-hydrodynamics code (MURaM; Vögler et al. 2004; Rempel et al. 2009b). We used the updated numerical treatment from Rempel (2017), without using the coronal features introduced therein. The bottom boundary is open to mass flows and the magnetic fields are extrapolated symmetrically, as described by Rempel (2014, Sect. 2.2, item ‘OSb’). The top boundary is open to inflows and closed to outflows, with the magnetic field potential as in Cheung (2006, Appendix C.1) and Rempel (2012, Appendix B).

2.1. Initial condition: Potential with opposing flux

The initial vertical component of the magnetic field at the bottom of the box is given by

$$\begin{array}{cc}\hfill B_{\mathrm{z},0}(r)& =B_{0}exp(-\frac{r^2{B}_0\pi }{F_{\text{Gauss}}})-B_{\text{opp}},\hfill \end{array}$$(1)

where B₀ is the initial maximal magnetic field, r is the distance from the spot axis, and F_Gauss is the flux contributed by the Gaussian distribution. An alternative to a pure Gaussian is the subtraction of a uniform vertical magnetic field, B_opp. Performing a horizontal integration gives the total flux through the box,

$$\begin{array}{cc}\hfill F_{\text{tot}}& ={\displaystyle \int B_{\mathrm{z},0}}\left(\sqrt{x^2+y^2}\right)\mathrm{d}x\mathrm{d}y=F_{\text{Gauss}}-B_{\text{opp}}{w}^2,\hfill \end{array}$$(2)

where w is the width of the box in both the horizontal directions. This value is the same for all heights and all times (but not over the τ = 1 iso-surface). The initial field throughout the box is given by potential field extrapolation according to

$$\begin{array}{cc}\hfill B_{\mathrm{z}}& ={\mathcal{F}}^{-1}\left(\mathcal{F}(B_{z,0}){\mathrm{e}}^{-z|k|}\right)\text{and}{B}_{\mathrm{x}}={\mathcal{F}}^{-1}\left(\mathcal{F}(B_{z,0})\frac{-i{k}_{\mathrm{x}}}{|k|}{\mathrm{e}}^{-z|k|}\right),\hfill \end{array}$$(3)

where ℱ represents the horizontal Fourier transform. B_y can be obtained by replacing k_x by k_y, or in our axisymmetric case, via B_y(x, y) = B_x(y, x). This extrapolation is the same as the top boundary condition mentioned above, with the difference that there z = 0 stands for the top (extrapolation into the ghost cells during the simulation run) instead of the bottom of the box (creating the initial field).

The potential extrapolation results in the initial B_z falling much more slowly with increasing radius outside the spot (see Fig. A.1, top-left panel) than self-similar initial conditions as in Rempel (2012). For all initial conditions with a pure Gaussian in the smaller simulation box, B_z remained above 155 G everywhere in the box. Given that in observations, B_z outside the spot vanishes, we tested the effect of subtracting a uniform vertical field (B_opp). This subtraction can be performed before or after extrapolation and the result is the same.

2.2. Simulation dataset

The simulation runs described in this article are listed in Table 1, where the following parameters were varied:

• Initial maximal magnetic field strength B₀: 20 kG, 40 kG, 80 kG, or 160 kG (most runs). B₀ is the first number in run names.¹

• Box width, w: 49.152 Mm (most runs) or 98.304 Mm. Runs with w = 98.304 Mm have names ending in ‘w’.

• F_Gauss = 10²² Mx (most runs) or F_tot = 10²² Mx. Runs with F_Gauss > 10²² Mx have names ending in ‘c’.

• Opposing magnetic field B_opp: 0 G, 50 G, 100 G, 150 G, 200 G, and 300 G. This is the second number in the simulation names.

• Most runs use a 96 (32) km horizontal (vertical) resolution, except the high-resolution runs with 32 (16) km, which have a name ending containing ‘h’. All low-resolution simulations use a grey radiative transfer.

The simulations with names ending in ‘o’ were run with a previous version of the MURaM code and provided bolometric intensities instead of 500 nm continuum intensities. They also contain further subtle differences that we are discussed below, where necessary.

The box was 8.192 Mm deep for all our simulations, but we provided comparison simulations ‘alpha*’, which were only 6.114 Mm deep. These are continuations of those started in Rempel (2012), which were described in Jurčák et al. (2020) as ‘type II’. They have self-similar initial conditions and top boundary conditions, which force the field to be more horizontal than a potentially extrapolated field by a factor of α.

While the simulation ‘160_opp000o’ continued running at low resolution, after 4.8121 h, it was regridded to the higher 32 (16) km resolution; at the time 7.2163 h, it switched over to a non-grey radiative transfer. In that form, which we labelled as ‘160_opp000ho’, it was described as ‘type I’ in Jurčák et al. (2020). Simulation ‘160_opp000h’ was run at the higher resolution from the start.

All high-resolution simulations presented here (‘160_opp000{h,ho}’, ‘alpha*’) use non-grey radiative transfer (Nordlund 1982; Ludwig 1992; Vögler et al. 2004) with four opacity bins. The main benefit of this is that it allows forward radiative transfer, degradation to instrumental resolution, adding noise, and subsequent inversions, as done in Jurčák et al. (2020). The simulation ‘160_opp000b’ suppresses mass flows through the bottom in regions where |B_z|≳10 kG. The azimuthal averages of the initial conditions at photospheric heights are shown in Fig. A.1. In addition, we used the Helioseismic and Magnetic Imager (HMI/SDO) (Scherrer et al. 2012) observations of NOAA AR 11591 for comparison. For more details on this spot, in particular, its umbral and spot boundary, see Schmassmann et al. (2018, 2026).

2.3. Data processing

During the simulation run, every 50 time steps, the intensity, I_c, as well as velocities, v_xyz, and magnetic field vectors, B_xyz, at the location of the τ = 1 layer were saved. This corresponds to one slice every ≈36 s or ≈43 s. The larger time steps are for the older simulations (with names ending in ‘o’), where the Alfvén speed limiter was set to a lower value. For every simulation, we evaluated 101 such slices, thus covering ≈1 h or ≈1.2 h. The time ranges covered are listed as t_ref in Table 1. The different Alfvén speed limiter values only change the length of the computation time steps and affect the low-density regions above the umbra, which are not relevant to this study.

The quiet Sun intensity, I_qs, as well as the umbral and spot boundaries, were calculated using the same methods as in Schmassmann et al. (2021), where the umbral boundary corresponds to I_c = 0.5I_qs and the spot boundary to I_c = 0.9I_qs. For consistency, we calculated a single I_qs for all new simulations (not ending in ‘o’) by averaging over the quiet sun data of all simulations and all time steps and we did the same, separately, for the old simulations (with names ending in ‘o’). These two sets of simulations were treated separately because the changes in MURaM between these runs resulted in incompatible intensity values (500 nm continuum vs bolometric). Maps of the intensity and velocities indicate that the umbral boundary is sharp and well retrieved by this procedure, whereas an accurate determination of the spot boundary would require identifying individual convective cells and classifying them as belonging to the spot or quiet Sun. Using a contour on a degraded intensity map (as we do here) will only result in an approximate result, but it is good enough to retrieve the spot areas and radii.

In the next step, the vectors were transformed to obtain the radial field², B_r, total field, |B|, inclination to the surface normal $\gamma =arctan(\sqrt{B_{\mathrm{x}}^2+B_{\mathrm{y}}^2}/B_{\mathrm{z}})$ and corresponding results for the velocities. The values were then remapped from Cartesian to cylindrical coordinates, after which an azimuthal average was performed. In addition to the averages of the values themselves (e.g. mean vertical velocity, v_z), the fraction of positive terms (upflow filling factor n_vz > 0) and signed averages (mean upflow velocity, v_z, +, mean downflow velocity, v_z, −) were calculated.

In the final step, we performed a least-squares linear fit of all time-dependent parameters, including the azimuthal averages. Then we evaluate the value of this linear function at 4 h (‘160_opp000ho’: 7.5 h, ‘alpha*’: $6.\overline{3}$h). For simplicity, hereafter in this article, these temporal fits are implied unless mentioned otherwise. For instance, when we refer to azimuthal averages, we mean these values of the linear functions at 4 h fit to the azimuthal averages. When we refer to maps in this work, we mean those of the last time in the reference range (last column of Table 1).

3. Results

3.1. Intensity

Figure 1 shows the bolometric intensity maps of sunspot simulations with four values of the initial maximal magnetic field strength B₀: 20 kG, 40 kG, 80 kG, and 160 kG, respectively. All runs contain the same flux (10²² Mx), with no opposing field B_opp, and situated in the small (49 Mm wide) box. The green line shows the umbral boundary as the longest contour of I = 0.5I_qs, in yellow other contours at the same intensity level, and in red the spot boundary as the longest contour at I_c = 0.9I_qs. The other contours at this level are not shown.

Fig. 1. Quadrants of bolometric intensity maps for four selected simulation runs with different initial magnetic field strengths. From top-left clock-wise: B0 = 20 kG, 40 kG, 80 kG, and 160 kG. The magnetic flux is the same for all four runs: F = 1022 Mx.

As a global property, the size of the penumbrae increased with increasing initial magnetic field strengths, B₀, up to 160 kG. Simulation runs with B₀ = 20 kG and 40 kG show elongated convective cells resembling penumbral filaments. The B₀ = 160 kG run shows the slenderest and longest penumbral filaments, and B₀ = 80 kG shows an intermediate state.

Figure 2 and the bottom-right panel of Fig. 4 show the azimuthal averages of the intensity for selected simulations. The vertical dashed lines show the average spot and umbral radii, r, calculated from the contour areas using $r=\sqrt{A/\pi }$. For all simulations, the umbral and spot sizes, along with their penumbral fractions are given in Table 2 (columns r_u, r_s, and r_pu/r_s = 1 − r_u/r_s). Increasing F_Gauss beyond 10²² Mx increases the sizes of the umbrae and spots. Based on sunspot statistics, Kiess et al. (2014) showed that the umbral radii increase linearly with the spot flux. In our simulations, the umbral radii increased with spot flux, consistent with their linear fit. However, for B₀ < 160 kG, the penumbral fractions are smaller than for those with B₀ = 160 kG, which, at 47%–52%, are already below the typical value observed in sunspots, where the penumbra covers about 60% of the spot size (e.g. Keppens & Martinez Pillet 1996; Westendorp Plaza et al. 2001; Mathew et al. 2003; Borrero et al. 2004; Bellot Rubio et al. 2004; Sánchez Cuberes et al. 2005; Beck 2008; Borrero & Ichimoto 2011; Chowdhury et al. 2024, and references therein). In summary, our simulations provide sunspots with a realistic umbral size but diminished penumbral length; that is, our simulations provide smaller sunspots than the observed stable ones. However, the simulated spot sizes resemble those of developing sunspots during penumbra formation (see e.g. Jurčák et al. 2014).

Fig. 2. Azimuthally averaged intensity, I/Iqs, (solid lines) as a function of radius in Mm and fraction of the spot radius. The colours indicate the simulations listed in the inset legend. For the older runs (with names ending in ‘o’ in red), the bolometric intensity is shown, and for the other runs (blue), the continuum intensity is shown. Vertical dashed lines show the position of the umbral (left panel with R < 10 Mm, right panel R < 1) and spot boundaries (left panel R > 10 Mm, right panel R = 1 black). The horizontal dashed lines in the right panel show the average umbral intensity. Values corresponding to the location of the vertical dashed lines are given in Table 2, columns ru, rs, ru/s.

The difference in brightness in the umbrae between older runs (names ending in ‘o’, red lines in Figs. 2) and the newer ones (blue lines) are caused by the former giving bolometric intensities, I = I_bol; whereas the latter are giving the continuum intensity at 500 nm, I = I_c500. The differences between these intensities are consistent with the expectations of blackbody radiation.

3.2. Surface magnetic fields

In Fig. 3, we plot the radial dependence of the azimuthally averaged B_z for the runs given in the figure’s legend. Simulations with F_Gauss > 10²² Mx (names ending in ‘c’) have stronger vertical magnetic fields in the spot axis B_z(r = 0) at the average τ = 1 height in the initial condition (see top-left panel of Fig. A.1, three light blue lines) than those with F_Gauss = 10²² Mx (dark blue line). This results in unrealistically strong fields in the later stages of the simulation (see Fig. 3, blue lines, and the ‘maxB_z’ column in Table 2). By unrealistic, we mean that they are stronger than observed on the sun in spots of comparable fluxes (Rezaei et al. 2012; Kiess et al. 2014).

Fig. 3. Vertical magnetic field component, Bz, of selected simulations. The Bz value from an observation (green) is also displayed for comparison.

The initial phase of the simulations had both concentrating and dispersing effects in the magnetic fields at the spot centre. For the B₀ = 160 kG simulations starting with the lower resolution (96 / 32 km), the concentrating effects are slightly stronger than the dispersing effects, resulting in a moderate increase in B_z. For the simulation starting with a higher resolution (‘160_opp000h’), the concentrating effects were stronger, resulting in an excessively high B_z. For B₀ = 40 kG and 80 kG, the dispersing effects are weaker than those for the B₀ = 160 kG simulations, resulting in an overly strong maxB_z. This is consistent with a simulation with B₀ = 160 kG, but with a closed bottom boundary at the spot foot-point (‘160_opp000b’), which also lacks the very dynamic initial phase with its dispersion effects and has a higher maxB_z. The simulation with a closed bottom boundary also lacked a properly formed penumbra and the few elongated granules had inflows along the full length and downflows outside the umbral boundary.

Low-resolution simulations with F_Gauss = 10²² Mx and B₀ = 160 kG have maxB_z = 4240 − 4478 G, while B_z drops with increasing radius to values of 200 G similar for all such simulations. Outside the spot, the values for B_z fluctuate, whereas simulations in larger boxes or higher B_opp fluctuate around lower values, as shown in the top left panel of Fig. 4.

Fig. 4. Azimuthal averages of the vertical magnetic field component, Bz (top left), the radial magnetic field component, Br (top right), the inclination (bottom left), and the intensity, Iout (bottom right), of our more realistic simulations (blue, red, black and grey), and of an observation (green), and an old reference simulation (gold). Horizontal dashed lines refer to Bz = 1867 G, the observed critical value for penumbral-type convection to operate (top left), and to an inclination of 90°, meaning a horizontal field (bottom-left).

Simulations with F_Gauss > 10²² Mx have increased B_r compared to those with F_Gauss = 10²² Mx. Beyond that, within the spot, all simulations with B₀ = 160 kG have a similar dependence of B_r on the fractional spot radius, as depicted in the upper right panel of Fig. 4. The B₀ = 20 kG simulation has the maximum B_r at a larger radius, consistent with the higher umbral fraction. All B₀ < 160 kG simulations have a faster drop of B_r outside the spot than the B₀ = 160 kG simulations. Of the B₀ = 160 kG simulations, B_r decreases more slowly outside the spot for those in larger boxes, indicating that the drop-off of B_r in smaller boxes may be caused by the horizontal periodic boundary conditions.

The increase in field inclination, γ, with fractional spot radius is similar for all B₀ = 160 kG simulations up to γ ≈ 80° (Fig. 4, bottom left panel), whereas γ rises more slowly with the fractional spot radius for the B₀ < 160 kG simulations, which is consistent with their larger umbral fraction. The other differences in the initial conditions are shown in the bottom right panel of Fig. A.1 within the spot were removed during the dynamic initial phase. Outside the spot, the average inclination fluctuated around 75° −105°, depending on the average B_z there.

3.3. Surface flows

The simulations with B₀ = 20 kG and 40 kG have inflows directed towards the umbra along the full length of all elongated convective cells, with downflows at umbral boundaries. Since this is clearly inconsistent with what is observed on the sun, we did not emark on a more detailed discussion of the surface flows in these runs. A comparison of the other runs with observations is given in Sect. 4.

Simulations with B₀ = 80 kG and 160 kG and low resolution (96/32 km) have strong in- and downflows in most umbral ends of penumbral filaments (filament heads), with out- and upflows in the middle and outer parts of the filaments. Downflows were also observed between the filaments and at the outer ends. An example map of the radial flows (in- and outflows) and vertical flows (up and down) is presented in Fig. 5.

Fig. 5. Surface velocity map of sunspot run ‘160_opp000’. Top-left: Vertical flows, vz, downflows in red, upflows in blue. Bottom-right: Radial velocities, vr, inflows in red (towards the spot centre), outflows in blue. The black contours are the same as shown in Fig. 1.

In the umbra, the dominant contribution to the vertical flows, v_z, is an overall oscillation, which at the time shown here is near a zero-crossing and, hence, it is not visible in the map. Our azimuthally averaged and temporally fitted values of v_z have these oscillations averaged out over time. This was verified in some simulations by subtracting the τ = 1 surface oscillations from the vertical velocities. Outside the umbra, the averages of the upflow filling factor and separate averages over up- and downflow speeds change quantitatively, but the radial distance of the minima and maxima is not affected. Inside the umbra, these values are not usable because vertical oscillations are the dominant contributors to the velocity.

The following discussion focuses on the common flow behaviour of the B₀ = 160 kG simulations. We find that the simulation run with B₀ = 80 kG is similar in many aspects, but the shorter penumbral filaments result in slower outflows.

In terms of azimuthal averages, the outflow speed v_r, + increases with increasing radius until the spot boundary and then slowly decreases. Inflow speeds of −v_r, − increased until the umbral boundary, then converged to the lower quiet sun value. The outflow filling factor, n_vr > 0, has values of 39–48% at a local minimum near the umbral boundary, then rapidly increases and tapers off to a maximum around the spot boundary of 88–93% for the B₀ = 160 kG simulations. This results in the radial velocity, v_r, having a minimum near the umbral boundary and a maximum near the spot boundary. The minima, maxima, and their locations in our simulation runs are listed in Table 2, with v_r and v_z shown in Fig. 7, where the extrema are marked with + symbols.

Upflow speeds, v_z, +, increased steadily with radius until they reached the quiet-sun average at around 1.2 r_s. Downflow speeds, −v_z, −, have a local speed maximum near the umbral boundary, followed by a decrease in speed over the inner penumbra to increase to the quiet Sun average. Outside the umbral boundary, the upflow filling factor, n_vz > 0, increased over the same range as the downflows slowed with increasing radius. Outside this narrow range, any trend in the radial dependence of the filling factor is sufficiently weak to be hidden by any fluctuations that might appear random. The lower panel of (‘160_opp000h’) shows that inside the central umbra, vertical velocities, v_z, are, on average, near zero, then drop to a local minimum near the umbral boundary. This is followed by a quick increase over the inner part of the penumbra and a decrease to the quiet sun average over the width of the rest of the penumbra, accompanied by an increase in fluctuations.

Having a negative average quiet-sun vertical velocity does not stand in contradiction to the convective blue shift, since the convective blue shift results from weighting the average of the light, meaning that the bright granular upflows outshine the darker inter-granular downflows. Our average is negative because the τ = 1 surface in the inter-granular lanes is lower; hence, the downflows are faster than the upflows.

The maximal total velocity, v, was found to be outside, but near the spot boundary, beyond which it dropped to the quiet-sun average. Most simulation runs have a minor local maximum near the umbral boundary.

In summary, as shown in Fig. 5, a pattern that appears characteristic of the potential field simulations with B₀ ≥ 80 kG presented here is manifested by the bi-directional flows exhibited by the penumbral structure: radial outflows (blue) resembling the observed Evershed flows are present in the mid and outer penumbra, yet the inner penumbra is characterised by strong inflows (red) combined with downflows. Bi-directional flows were first reported in high-resolution observations by Schlichenmaier et al. (2011) as counter-Evershed flows. Recently, García-Rivas et al. (2024) showed that some of these counter-Evershed flows correspond to bi-directional flows that act as a precursor to the formation of a stable penumbra. The simulation with a closed bottom boundary inside the magnetic trunk (‘160_opp000b’) has elongated granules instead of penumbral filaments and, hence, no outflows.

Figure 6 shows that the high-resolution (32/16 km) simulations have, however, both filaments with bi-directional flows and those with normal Evershed flows, meaning upflows in the filament heads, outflows along the whole lengths of the filaments, and downflows at the outer ends and between the filaments. The simulation that started with high-resolution (‘160_opp000h’) has faster outflows (Fig. 7, upper panel), consistent with the resolution dependence of the Evershed flow speed reported by Rempel (2012). The simulation that started at low resolution and was later regridded to high resolution (‘160_opp000ho’) has slower outflows, because the flow speeds decrease over time and this simulation is evaluated at a later time.

Fig. 6. Same as Fig. 5, but for sunspot run ‘160_opp000h’ (high res.), with no spot boundary contour.

Fig. 7. Azimuthally averaged radial, vr (top panel), and vertical, vz (bottom panel) velocities of our more realistic simulations.

3.4. Flux emergence

Our low-resolution (96/32 km) simulations with B₀ = 160 kG and open bottom boundary (names starting with ‘160’ and an ending that does not contain ‘h’ or ‘b’) keep increasing their flux during our observation time (Table 2, column d_tF/F, 8.5–13.0%/h). This corresponds to flux emergence. In the only simulation run for substantially longer than 4 h (‘160_opp000’), the stalling of flux emergence at 8–10 h coincides with a narrowing of the penumbra and the disappearance of most outflows.

In simulation ‘160_opp000ho’, for which the data were evaluated 3.5h later than in most other simulations, the flux emergence appears to be tailing off. In the simulation with a closed bottom boundary inside the magnetic trunk (‘160_opp000b’), only slight flux emergence is observed.

4. Summary and conclusion

As the aim of this study is to explore more realistic sunspot simulations by extending the ideas introduced by Nordlund (2015) and examined in Jurčák et al. (2020), we carried out a parameter study in MURaM sunspot simulations under the potential-field initial condition. Here, we varied the following parameters:

• Initial magnetic field strength, B₀: 20, 40, 80, and 160 kG;

• Magnetic flux content F_Gauss and opposing flux, B_opp;

• Box width;

• Resolution (96 / 32 km versus 32 / 16 km).

Our main results are:

1. The penumbral extent increases systematically with B₀, from no penumbra at 20 kG to long, slender filaments at 160 kG.

2. Runs with F_Gauss > 10²² Mx develop unrealistically large umbral vertical fields, B_z.

3. The most realistic sunspots come from B₀ = 160 kG, F_Gauss = 10²² Mx, which produce (1) maximum B_z = 4.2 − 4.5 kG (still larger than found in observations) and (2) a realistic profile of B_z and B_r with the radius.

4. Closing the bottom boundary inside the spot trunk does not create a realistic sunspot.

5. Varying B_opp or the box width changes the background B_z outside the spot but does not significantly affect the internal spot structure or dynamics.

6. For B₀ ≤ 40 kG, no Evershed flows develop; there are only in- and downflows, with only slight flux emergence.

7. For B₀ = 80 − 160 kG and low resolution (96/32 km), bi-directional flows, resembling those found in observations of the early stages of penumbra formation, are found during ongoing flux emergence.

8. High-resolution (32/16 km) runs produce both filaments with bi-directional flows and filaments with Evershed flows.

From this parameter study, we conclude that simulations initialised with potential-field conditions, sunspot fluxes on the order of 10²² Mx, and strong (160 kG) magnetic fields and open bottom boundary reproduce key aspects of the observed morphology and dynamics of sunspots during the early phases of penumbra formation (see e.g. Schlichenmaier et al. 2010a; Jurčák et al. 2014). Specifically, our simulations have early penumbral-like regions characterised by slender granulation of reduced intensity; rather than exhibiting the classical radial velocities of the Evershed flow (Schlichenmaier & Schmidt 2000; Löhner-Böttcher & Schlichenmaier 2013), they show bi-directional motions that are likely driven by flux emergence at the outskirts of the spot (García-Rivas et al. 2024). These results highlight the value of these simulations in investigating the initial stages of sunspot formation.

Comparisons with failed simulations show that – at least for low-resolution simulations with potential field top boundary conditions – we only get a realistic-looking sunspot with outwards flow in a sufficiently wide penumbra, while flux keeps emerging in the outer parts of the spot.³ We speculate that any low-density plasma in ∩-shaped magnetic field lines below the surface might be sufficient to drive flux emergence and create penumbrae as in our simulations.

However, observations suggest that an overlying canopy is required to form a stable penumbra with Evershed flows (Lindner et al. 2023; Chifu et al. 2025). Older simulations have achieved this by forcing the field to be more inclined at the top boundary, so that the inclined field drives the Evershed flows. The presence of Evershed flows in some of the filaments of our high-resolution simulations, coupled with their absence in low-resolution simulations, indicates that insufficient resolution is an inhibiting factor for Evershed flows, in line with preliminary reports from Rempel & Hotta (2025). We speculate that embedding the spot in a much larger simulation box (across all dimensions above and below the photosphere), while including a small-scale dynamo run that allows a canopy to form naturally, might also force the field to be sufficiently inclined. As a result, Evershed flows might indeed form even in low-resolution simulations if the spot is placed in the right location.

Acknowledgments

We thank the anonymous reviewer for their valuable and insightful comments, which helped improve this manuscript. This work was supported by the Czech-German common grant, funded by the Czech Science Foundation under project 23-07633K and by the Deutsche Forschungsgemeinschaft under project BE 5771/3-1 (eBer-23-13412), and the institutional support ASU:67985815 of the Czech Academy of Sciences. Most of the simulations leading to the results obtained were run on Piz Daint at the Swiss National Supercomputing Centre, Switzerland, financed through the ACCESS programme of the SOLARNET project, which has received funding from the European Union Horizon 2020 research and innovation programme under grant agreement no 824135. For the older simulations, we would like to thank Matthias Rempel for providing us with their data and acknowledge the high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by the NFS NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Initial magnetic field configurations

Figure A.1 shows the azimuthal averages of the initial magnetic field at the average height of the τ = 1 iso-surface in the quiet Sun. The top two panels (B_z + B_opp and B_r), show the initial conditions for all our simulations (except ‘80_opp000o’ and ‘alpha*’), because introducing a B_opp ≠ 0 only offsets B_z uniformly, and does not affect B_r. In the bottom right panel, all simulations (except ‘alpha*’) are plotted. Colours and line-styles are described in Table A.1, where the ‘?’ in ‘160_opp000?’ is a wildcard, standing for any of ‘o’, ‘ho’, ‘h’, ‘’, or ‘b’, since these simulations share the same initial conditions.

Fig. A.1. Initial magnetic field. Vertical field Bz + Bopp (top left), radial field Br (top right), total field |B| (bottom left), and inclination to the surface normal γ. The legend in the top left panel shows the simulations plotted in the first three panels. The legend in the bottom right panel gives a selection of the simulations shown there, the rest is given in Table A.1; for each line in that table (colour), the further right the column is (the higher B0 or Bopp, the changing numbers in the run names), the higher the inclination.

While this figure shows many details, primarily the following aspects are pertinent:

• 1)

  F_Gauss > 10²² Mx leads to strong B_z on the spot axis.

• 2)

  w = 98.304 Mm leads to the field being distributed over a wider box, resulting in a weaker B_z far from the spot axis and more extended B_r.

• 3)

  The inclinations in the initial conditions show a wide range of behaviour, in particular outside the spot.

All Tables

All Figures

	Fig. 1. Quadrants of bolometric intensity maps for four selected simulation runs with different initial magnetic field strengths. From top-left clock-wise: B0 = 20 kG, 40 kG, 80 kG, and 160 kG. The magnetic flux is the same for all four runs: F = 1022 Mx.
In the text

Fig. 2. Azimuthally averaged intensity, I/Iqs, (solid lines) as a function of radius in Mm and fraction of the spot radius. The colours indicate the simulations listed in the inset legend. For the older runs (with names ending in ‘o’ in red), the bolometric intensity is shown, and for the other runs (blue), the continuum intensity is shown. Vertical dashed lines show the position of the umbral (left panel with R < 10 Mm, right panel R < 1) and spot boundaries (left panel R > 10 Mm, right panel R = 1 black). The horizontal dashed lines in the right panel show the average umbral intensity. Values corresponding to the location of the vertical dashed lines are given in Table 2, columns ru, rs, ru/s.

In the text

	Fig. 3. Vertical magnetic field component, Bz, of selected simulations. The Bz value from an observation (green) is also displayed for comparison.
In the text

Fig. 4. Azimuthal averages of the vertical magnetic field component, Bz (top left), the radial magnetic field component, Br (top right), the inclination (bottom left), and the intensity, Iout (bottom right), of our more realistic simulations (blue, red, black and grey), and of an observation (green), and an old reference simulation (gold). Horizontal dashed lines refer to Bz = 1867 G, the observed critical value for penumbral-type convection to operate (top left), and to an inclination of 90°, meaning a horizontal field (bottom-left).

In the text

	Fig. 5. Surface velocity map of sunspot run ‘160_opp000’. Top-left: Vertical flows, vz, downflows in red, upflows in blue. Bottom-right: Radial velocities, vr, inflows in red (towards the spot centre), outflows in blue. The black contours are the same as shown in Fig. 1.
In the text

	Fig. 6. Same as Fig. 5, but for sunspot run ‘160_opp000h’ (high res.), with no spot boundary contour.
In the text

	Fig. 7. Azimuthally averaged radial, vr (top panel), and vertical, vz (bottom panel) velocities of our more realistic simulations.
In the text

Fig. A.1. Initial magnetic field. Vertical field Bz + Bopp (top left), radial field Br (top right), total field |B| (bottom left), and inclination to the surface normal γ. The legend in the top left panel shows the simulations plotted in the first three panels. The legend in the bottom right panel gives a selection of the simulations shown there, the rest is given in Table A.1; for each line in that table (colour), the further right the column is (the higher B0 or Bopp, the changing numbers in the run names), the higher the inclination.

In the text