Issue |
A&A
Volume 700, August 2025
|
|
---|---|---|
Article Number | A28 | |
Number of page(s) | 9 | |
Section | The Sun and the Heliosphere | |
DOI | https://doi.org/10.1051/0004-6361/202554253 | |
Published online | 29 July 2025 |
Coriolis force acting on near-surface horizontal flows during simulations of flux emergence produces a tilt angle consistent with Joy’s law on the Sun
1
The University of Newcastle, University Dr, Callaghan NSW 2308, Australia
2
Max-Planck-Institut für Sonnensystemforschung, 37077 Göttingen, Germany
3
Institut für Astrophysik, Georg-August-Universität Göttingen, 37077 Göttingen, Germany
⋆ Corresponding author: hannah.schunker@newcastle.edu.au
Received:
25
February
2025
Accepted:
12
June
2025
Context. Joy’s law describes the tilt of bipolar active regions on the Sun away from an east-west orientation, where the flux of the polarity concentrated at the prograde side tends to be closer to the equator than the polarity on the retrograde side. Joy’s law is attributed to the Coriolis force because of the observed increase in the tilt angle at higher latitudes. This tilt plays a crucial role in some solar dynamo models.
Aims. Our goal is to model the effects of the Coriolis force on a flux tube as it rises through the near-surface convection zone.
Methods. We used a three-dimensional Cartesian magnetohydrodynamic simulation of an untwisted flux tube ascending from a depth of 11 Mm. We modelled the Coriolis effect using the f-plane approximation, which only considers and acts on horizontal flows. On the Sun, Joy’s law is weak and is only evident as an average over many active regions. To achieve a measurable effect in a single simulation, we considered a rotation rate 110 times faster than that of the Sun.
Results. The simulation shows that the flux tube emerges at the surface with a tilt angle consistent with Joy’s law when scaled to the Sun’s slower rotation, and the tilt angle does not substantially change after emergence.
Conclusions. This shows that the Coriolis force acting on flows horizontal to the surface within the near-surface convection zone is consistent with Joy’s law.
Key words: Sun: activity / Sun: granulation / Sun: magnetic fields / Sun: photosphere / sunspots
© The Authors 2025
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The orientation of magnetic active regions on the solar surface obeys Hale’s law and Joy’s law (Hale et al. 1919). Simple active regions consist of a bipole pair aligned in a predominantly east-west orientation. Hale’s law states that the sign of the leading polarity (relative to the Sun’s rotation) in each hemisphere is opposite, and the sign of the polarities switches from one solar cycle to the next. Joy’s law describes the statistical tendency of the leading polarity (in the direction of the Sun’s rotation) of an active region to be closer to the equator than the following polarity. Studies have shown that the tilt of the polarities increases with (absolute) solar latitude (Wang & Sheeley 1989). These two laws play a crucial role in the Babcock-Leighton dynamo model (Babcock 1961; Cameron & Schüssler 2015; Karak & Miesch 2017).
Joy’s law has conventionally been quantified using continuum intensity observations of sunspots (e.g. McClintock & Norton 2016), which limits the description to well-developed and stable active regions with established sunspots. To pinpoint the origin of active region tilt angles, however, observing the magnetic field is required to capture the earliest stages of flux emergence.
Monitoring missions, such as the Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory (SoHO; Scherrer et al. 1995) and the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) on board the Solar Dynamics Observatory (SDO; Scherrer et al. 2012), make routine observations of the Sun’s surface magnetic field at a high cadence and with a high duty cycle that are able to capture the early stages of active region emergence. From over 700 active regions, Kosovichev & Stenflo (2008) found that, statistically, active regions emerge roughly east-west aligned and with zero tilt angle. By analysing the motion of the individual polarities in SDO/HMI line-of-sight magnetic field observations, Schunker et al. (2020) found that the tilt angle develops during the emergence process (in 1–2 days; Weber et al. 2023) and that the latitudinal dependence of the tilt angle is due to the north-south separation of the polarities. The latitudinal dependence suggests that the Coriolis force is responsible for this north-south motion. This raises the question of which flows are influenced by the Coriolis force to generate the tilt angle.
One possibility involves the turbulent convective velocities (Parker 1955; Choudhuri & D’Silva 1990; Brandenburg 2005). Schmidt (1968) proposed that active region bipoles emerge with upwelling supergranulation cells, where the surface flows within the cell push the polarities outwards, further separating them. An alternative explanation arises from the buoyant ascent of magnetic flux tubes through the Sun’s convection zone (Wang & Sheeley 1991; D’Silva & Choudhuri 1993; Fisher et al. 1995; Weber et al. 2011). These types of studies usually employ the thin-flux-tube approximation, which is invalid in the top ∼20 Mm beneath the surface (see Fan 2009 and citations within).
The north-south component of the Coriolis force is associated with the creation of radial vorticity. This radial vorticity can introduce a north-south separation of the two magnetic polarities during flux emergence, giving rise to Joy’s law. In this study we investigated the effect of the radial vorticity generated by the Coriolis force on the orientation of an emerging flux tube. For this purpose, we considered the effect of the Coriolis force acting on horizontal flows (equivalently, the radial vorticity component).
We simulated a flux tube rising through the near-surface convective region of the Sun and then through the surface. In Sect. 2 we introduce our numerical approach, including a description of the implementation of the f-plane approximation of the Coriolis force and the initial condition for the flux emergence simulations. In Sect. 3 we describe the evolution of the tilt of the polarities in each simulation in relation to the flows. We discuss the implications of our results and our conclusions in Sect. 4.
2. Numerical approach
To simulate the interaction of a flux tube rising through the near-surface convection of the Sun, we used the Max-Planck-Institute for Aeronomy/ University of Chicago Radiation Magneto-hydrodynamics (MURaM) code (Vögler et al. 2005; and further developed by Rempel 2014, 2017). This code solves the magnetohydrodynamic (MHD) equations in a three-dimensional Cartesian box that includes a section of the surface of the Sun. The local plasma velocity, v, evolves with the density, ρ, pressure, p, and magnetic flux density, B, through the MHD equations:
Here: t is the time and T is the temperature; , the sum of the internal and kinetic energy; and FSR limits the Alfvén speed (for more detail about the definitions, see Rempel 2014, 2017). We added the Coriolis acceleration, acor, to the equations solved by MURaM, which we describe in Sect. 2.1. We note that since acor and v are perpendicular, acor has no effect on the energy.
The Lorentz force is given by
the Alfvén velocity is vA, and the identity tensor is I. The parameter c limits the speed of Alfvén waves, implemented specifically for extending the atmosphere of MURaM into the chromosphere and corona. The case when c is set to the speed of light corresponds to relativistic Alfvén waves. We limited the Alfvén wave speed to 100 km s−1, which should be inconsequential in our simulations of the convection zone and photosphere. We used a constant value for the gravitational acceleration g = −2.74 × 104 cm2 s−1 in the vertical direction.
In the energy equation (Eq. (3)), Qrad is the radiative heating (and cooling) term. The radiation transport equations are solved using the short characteristics method (Kunasz & Auer 1988; Vögler et al. 2005). We treated the radiation as monochromatic (grey) with opacities from the MPS ATLAS package (Witzke et al. 2021). We defined the solar surface in the simulation to be the τ500 = 1 layer.
An equation of state is required to close the system. We assumed the plasma is in local thermodynamic equilibrium. The equation of state was determined numerically by relations specifying T(ρ, Eint) and p(ρ, Eint), where Eint is the internal energy density, which are read from tables generated using the FreeEOS package (Irwin 2012).
MURaM has a free parameter, h, that controls the numerical diffusion (Rempel 2014). Following Rempel (2017), we used a value of h = 2 for the lower boundary and regions with density ρ > 10−11 g/cm3 (mostly below the surface), h = 1 for regions with density below ρ ≤ 10−11 g/cm3 (above the surface), and h = 0 for the upper boundary. A discussion of the implications of the h parameter and diffusion scheme is given in Rempel (2017).
The top boundary condition of the domain is semi-transparent, permitting outflows but not inflows, and the magnetic field is forced to be purely vertical. We used the ‘Open boundary/symmetric field’ condition (OSb. in Rempel 2014) at the bottom boundary, which is open to flows and is symmetric for the magnetic field. We set a constant entropy in the inflows at the bottom boundary and extrapolated the pressure across the boundary cells. The horizontal boundaries are periodic. With MURaM, we simulated the near-surface of the Sun from 15.128 Mm below the surface to 1 Mm above (a total of 504 grid points with a resolution of 0.032 Mm in the z-direction) and extending 96.768 Mm (over 1008 grid points with a resolution of 0.096 Mm) in each horizontal direction (x, y).
To initiate the simulations, we prescribed the density, temperature, sound-speed and pressure of a standard solar model (Christensen-Dalsgaard et al. 1996, based on Model S) in the box that varies only in height, and initiated the convection with small amplitude random velocity perturbations uniformly distributed between ±100 cm s−1. Importantly, we initiated the convection without rotation, i.e. acor = 0. Starting from our initial condition the simulation required about 8 solar hours until convection reached a steady state. After reaching this steady state, we selected a single snapshot as the initial condition for all following simulations.
2.1. Implementing the Coriolis force in the f-plane approximation
Our goal is to understand the effects of the radial vorticity introduced by the Coriolis force. The Coriolis force is a pseudo-force felt by objects in a rotating frame of reference. The Coriolis acceleration is given by 2Ω × v where Ω is the rotation vector and v = (vϕ, vθ, vr) is the velocity within the reference frame in spherical coordinates where θ is latitude, ϕ is longitude, and r is the radial direction. The Coriolis acceleration is the cross product of the flow velocity, and the rotation rate,
where the first term is in the longitudinal direction, , the second term is in the latitudinal direction,
, and the third is in the radial direction,
. In our case, we were simulating a localised box covering a small section of the Sun’s surface that we considered to be located at a single latitude. We approximated the plasma velocity in spherical coordinates to Cartesian coordinates in our box as (vϕ, vθ, vr)→(vx, vy, vz).
Because we are interested in the radial vorticity introduced by the Coriolis force, we considered the terms generating radial vorticity from horizontal flows via Coriolis acceleration in the f-plane approximation
We plan to include the effects of the other terms in the f-plane approximation in a future paper. In our simulations we chose a mid-latitude of solar activity in the northern hemisphere, θ = 15°, and a rapid rotation rate of Ω = 47.75 μ Hz ≈110 Ω⊙ (where the Sun has a surface rotation rate of about 2 km s−1, or 430 nHz). We selected a fast rotation rate to ensure we would have a measurable effect from the Coriolis force on the typical convective velocities in a practical time frame of the order of about 10 hours of solar time. We ignored any perturbations to the asphericity or local gravity since we are only interested in the effect of the Coriolis force on the horizontal flows, as in the case for the Sun. We added the Coriolis term to the MHD equation of motion (Eq. (2)) and energy equation in the MURaM code (Eq. (3)). The local plasma velocity then evolves with the MHD equations, with the added effect from the rotation in the f-plane approximation with the acor term. To quantify the effect of the Coriolis force on the supergranulation scale convective flows, we simulated the behaviour of the convection with, and without, the Coriolis acceleration for ≈24 hours starting from identical initial conditions. We refer to the simulations with Ω = 110 Ω⊙ as the ‘rotational simulation’, and Ω = 0 Ω⊙ as the ‘non-rotational simulation’.
We validated the implementation of the f-plane approximation by comparing the divergence and vorticity of the flows in the hydrodynamic simulation to an analytic estimate (see Appendix A). A visual representation of the effects of the Coriolis force on the flows can be seen in surface supergranulation flows. We averaged the vx and vy surface velocity maps in 6-hour-long segments and filtered out spatial scales less than 15 Mm and greater than 25 Mm (the approximate size of supergranules that could form in the limited size of the simulation box). From each of these four segments we computed maps of the horizontal flow divergence, ∇h ⋅ vh (see Fig. 1).
![]() |
Fig. 1. Example of a surface horizontal flow divergence map averaged over 6 hours and filtered between 15 Mm and 25 Mm from the non-rotational hydrodynamic simulation. The cyan contours indicate the identified supergranules, and the blue points indicate the location of the peak divergence within the boundary. |
To identify supergranules, in each temporally averaged map of flow divergence we set a threshold of 0.8 standard deviations from the mean of all four flow divergence maps, following the definition used by Hirzberger et al. (2008) and Langfellner et al. (2015). We contoured the boundaries of the features at the threshold value, and defined the centre of each feature as the pixel with peak horizontal divergence within the boundary. If the peak was less than 1 Mm distance from any part of the contour, we excluded it from our selection of supergranules (see the example in Fig. 1).
We averaged the horizontal flow divergence maps centred at each supergranule peak to create the average supergranule (see Fig. 2, top). We repeated this for the horizontal velocity, vx and vy, maps. The average supergranule in the simulations with the f-plane approximation is shown in the middle panel of Fig. 2. As expected, the average supergranule in the non-rotational simulation has purely diverging, radial velocities, and the average supergranule in the rotational simulation has an azimuthal component with a similar amplitude to the radial component (bottom panel, Fig. 2). The average supergranule in the hydrodynamic simulations is qualitatively consistent with the maps of average supergranules on the Sun. The radial velocities have similar amplitudes (200–300 ms−1), and the radial and polar components have a similar sine-like profile with a peak about half a supergranule radius from the centre (e.g. see Figs. 12 and 13 in Langfellner et al. 2015, and recall that we are simulating the effect from a rotation rate more than 100 times that of the Sun).
![]() |
Fig. 2. Top panel: horizontal flow divergence of the average supergranule in the non-rotational hydrodynamic simulation. The cyan arrows represent the horizontal flows at the surface. Middle panel: same as the top panel but for the rotational hydrodynamic simulation. Bottom panel: azimuthal averages of the radial velocity ( |
2.2. Flux tube initial condition
To simulate flux emergence, we placed a magnetic flux tube in the simulation box with an axis aligned in the -direction at yc = 0 and zc = −11.096 Mm relative to z = 0 Mm at the surface. The magnetic field is given by
where d2 = (y − yc)2 + (z − zc)2 and yc and zc are the y-coordinate and z-coordinate of the tube axis, respectively, the maximum field strength was B0 = 50 kG, and the total flux was ΦT = 1021 Mx, typical of the observed flux in a small active region at the surface of the Sun. The flux tube had no radial or azimuthal component (twist) and satisfies flux conservation and the solenoidal constraint. The full-width-half-maximum of the radial Gaussian profile was 1.33 Mm, and we set the magnetic field to zero at a radius of d = 1.71 Mm, which enclosed 99% of the total flux. The width of the flux tube was less than the pressure scale height at this depth. To maintain pressure balance inside and outside the flux tube, we subtracted the magnetic pressure to make a new initial condition from the non-magnetic snapshot by modifying the pressure to be, pnew = p0 − B2/2μ0 where p0(x, y, z) is the pressure before the tube was added.
Because we only want one section of the tube to emerge, we modified the entropy, s, in the simulation according to
where s0(x, y, z) is the entropy before we insert the flux tube,
scool = 679444350 K−1 (corresponding to very cool material), N1 = 0.009, and σ = 7.4 Mm. The parameters N1, scool and σ were chosen so that only the part of the tube at the horizontal centre of the box (x = 0, y = 0) rises (see Rempel 2017, for more details on setting the entropy). We note that the horizontal centre of the domain is located in an upflow (see Figs. 4 and 5). Because we are assuming local thermodynamic equilibrium, specifying p and s specifies the plasma properties. The velocity field v is left unchanged. Together these specify the flux tube initial condition. We began with this initial condition and simulated the evolution of the flux tube, first without rotation (Fig. 4) and then with the f-plane approximation (Fig. 5).
![]() |
Fig. 3. Left: total magnetic field strength across the centre of the tube as a function of height. Right: Vertical slice through the centre of the flux tube showing the difference in density between the initial condition with the flux tube and the purely hydrodynamic initial condition. A reduced density (negative) will cause the tube to rise, and an increased density (positive) will cause it to sink. The contour of the magnetic field value of 25 kG is shown in green. The relative density perturbation is constant as a function of height in the background Model S atmosphere. |
![]() |
Fig. 4. Slices through y = 0 of the magnetic field strength and vertical velocity of the initial condition (top two panels), and then at two subsequent time steps of the non-rotational (Ω = 0 Ω⊙) MHD simulation (bottom four panels). The flux tube rises in an arch structure before forming two dominant legs. In the velocity slices, blue indicates downward flow and red upward flow. See the movies of the simulation online and here. |
![]() |
Fig. 5. Slices through y = 0 of the magnetic field strength and vertical velocity of the initial condition (top two panels), and then at two subsequent time steps of the rotational (Ω = 110 Ω⊙) MHD simulation (bottom four panels). In the velocity slices, blue indicates downward flow and red upward flow. See the movies of the simulation online and here. |
2.3. Emergence of the flux tube at the surface
In both the non-rotational and rotational MHD simulations, the emergence process occurs between t = 4 and t = 6 hours (Fig. 6), which is much faster than that observed on the Sun (about 2 days Weber et al. 2023). A general bipole magnetic field structure appeared at the surface at about t = 4 hours, and coherent polarities with intensity darkening formed at about t = 7 hours in the non-rotational MHD simulation (Fig. 7) roughly aligned along the flux tube axis (along y = 0). The unsigned flux at the surface is about twice the flux in the initial condition (as expected), and comparable to observed active regions at the very beginning of their emergence at the surface (e.g. Schunker et al. 2016, 2019).
![]() |
Fig. 6. Total unsigned magnetic flux at the surface as a function of time for both non-rotational and rotational MHD simulations. The vertical line indicates the approximate end time of the flux emergence (at 6 hours). |
![]() |
Fig. 7. Vertical magnetic field at the surface (grey-scale maps) and slices through y = 0 Mm showing the unsigned magnetic field strength and vertical flows at three different times of the Ω = 0Ω⊙ non-rotational MHD simulation. The red and yellow symbols show the centres of the positive and negative magnetic poles, respectively. The crosses indicate the flux-weighted centres of the polarities over a threshold of |Bz| = 800 G, and the dots indicate the centres without a threshold. See the movies of the simulation online and here. |
The apex of the rising flux tube is accompanied by a significant horizontal outflow (a positive divergence) at the surface, of the order of 2 km s−1 (Fig. 9), which we note is inconsistent with observations (e.g. Birch et al. 2016; Schunker et al. 2024). The outflow is due to the fast rise speed of the flux tube (see the middle panel of Figs. 4 and 5). A future goal is to develop an initial condition so that the simulation produces a realistic active region that rises more slowly and does not produce these strong surface outflows.
We note that there are two stages of emergence in the rotational simulation (see Fig. 6). The magnetic flux that reaches the surface in the rotational simulation is initially about 2 × 1021 Mx, not substantially less than the non-rotational simulation at about 2.5 × 1021 Mx. In the rotational simulation the flux grows again from about 10 hours, peaking at about 2.7 × 1021 Mx. This difference in flux between the simulations may be due to the effects of rotation or the non-linear nature of the simulations; however, a number of additional simulations with different convective states would be required to determine that.
3. Relationship between the motion of the polarities and surface flows
The positions of the surface polarities begin close together and separate as more flux emerges, as observed for active regions emerging on the Sun. The fast rise-speed of the flux tube generates outflows (Fig. 9) that are an order of magnitude greater (∼2 km s−1) than the convective flows (∼100 ms−1) during the emergence phase. Such outflows are not observed on the Sun (Birch et al. 2019). However, the fast rise speed coupled with the fast rotation rate produces a clear effect of the Coriolis force on the horizontal flows associated with a rising flux tube. In the non-rotational simulation, the polarities emerge at the surface aligned with the flux tube axis (Fig. 7). Towards the end of the simulation a small counter-clockwise tilt develops. In the rotational simulation, the flux tube emerges with a larger tilt angle in the clockwise direction (Fig. 8).
![]() |
Fig. 8. Same as Fig. 7 but for the110 Ω⊙ rotational MHD simulation. The red and yellow symbols show the centres of the positive and negative magnetic poles, respectively, used to calculate the geometric tilt angle, γg. See movies of the simulation online and here. |
In this section we measure the geometric tilt angle and compute the component of the tilt angle expected from the Coriolis force acting on the surface horizontal flows.
3.1. Geometric tilt angle
To track the location of the polarities, we computed the flux-weighted centre of each polarity, at the surface as a function of time (see Figs. 7 and 8). We defined the geometric tilt angle as the angle between the x-axis aligned with the flux tube initial condition and the line connecting the two polarities at sometime,
where Δxg and Δyg are the separation between the polarities (i.e. Δxg = x− − x+ and Δyg = y− − y+) and we defined a positive tilt in the clockwise direction, consistent with the direction of Joy’s law in the northern hemisphere on the Sun.
By design, in the rotational simulations we imposed a rotation rate more than a hundred times faster than that of the real Sun, and the polarities emerge with a significant tilt in the clockwise direction of about 30° (see Figs. 8 and 9). In both the non-rotational and rotational cases, most of the flux has emerged by t = 6 hours (Fig. 6) and, by inspection, any change in tilt angle after that time is not due to significant additional flux emerging from below, but from the motion of existing surface flux. Figure 10 shows that the motion is predominantly in the y-direction (consistent with Schunker et al. 2020), since Δxg remains roughly constant after emergence. The question is then how much of the motion of the polarities can be explained by the near-surface horizontal flows.
![]() |
Fig. 9. Surface flow velocities (cyan arrows) and vertical magnetic field averaged over time from t = 4.25 h to t = 5.44 h for the non-rotational (0 Ω⊙) hydrodynamic simulation (top) and the rotational (110 Ω⊙) hydrodynamic simulation (bottom). The velocity has also been spatially smoothed for clarity by a Gaussian with standard deviation of 2.4 Mm. Note that this has been computed at a constant geometrical height, at z = 480 pixels = 15.360 Mm from the bottom of the box. |
![]() |
Fig. 10. Separation of polarities for the non-rotational simulation (blue) and the rotational simulation (red). We show the geometric separation without a Bz threshold (solid curves, Δxg) and the separation expected from the flows (Eqs. (11) and (12); dashed curves, Δxf). The vertical line at 6 hours is close to the time when most of the magnetic flux has emerged (see Fig. 6). |
We note that when measuring the tilt angle in observations it is usually the tilt angle between coherent magnetic structures such as sunspots, and thresholds of magnetic field strength have been used to isolate dominant coherent structures in magnetic field maps (e.g. Schunker et al. 2016; Sreedevi et al. 2023). To demonstrate the equivalence, we include an example of the polarity locations and tilt angle computed only for |Bz|> 800 G (see Figs. 7, 8, and 11) that isolates the small pores that form. The polarity separation is larger in both directions; but, the tilt angle is not significantly different.
3.2. Tilt angle due to surface horizontal flows
The geometric positions of the polarities are affected by the motion of the vertical magnetic field due to the horizontal flows, the appearance of vertical field from flux emerging through the surface, and a diffusive component. In this section we measure the expected contribution to the motion of the polarities from the surface horizontal flows, Δxf and Δyf.
We computed the separation of the polarities in the x-direction due to the flows, Δxf, by integrating the horizontal flow velocities weighted by the magnetic field:
where the angle brackets indicate the mean over the surface area of the simulation and the integral is cumulative in time from t0 = 3 hours, after which there is consistently non-zero magnetic field at the surface (Fig. 10). Our flux tube lies along y = 0 and the magnetic field is pointed in the positive x-direction. Our sign convention is that a positive Δx corresponds to the negative polarity lying in the positive x direction relative to the positive polarity, and a positive Δy corresponds to the negative polarity lying in the negative y direction relative to the positive polarity, which explains the negative sign in Eqs. (11) and (12).
The tilt angle imparted by the flows is then
defined as positive in the clockwise direction. In both the rotational and non-rotational simulations, the tilt angle is not significantly different to the geometric tilt angle (Fig. 11), despite the differences in separation components. Figure 10 shows that the magnitude of the separation due to the horizontal flows is about one-third less than the geometric separation, and so there is some remaining component of motion not due to the horizontal flows at the surface.
![]() |
Fig. 11. Tilt angle of the polarities for the non-rotational simulation (blue) and the rotational simulation (red). The geometric tilt angle is shown by the solid curves, and the tilt angle expected from the flows is shown by the dashed curves. The tilt angle resulting from the location of the polarities with a threshold of |Bz|> 800 G applied is shown by the thick shaded curves. |
3.3. Remaining tilt angle component
We defined the difference between the geometric and horizontal flow components as the remaining component, Δxr = Δxg − Δxf and Δyr = Δyg − Δyf, and the corresponding tilt angle as γr = γg − γf. This remaining component of the separation is then due to the emergence of magnetic field through the surface and diffusion. On the Sun, the diffusive component is expected to be small, but not necessarily in the simulations.
In the non-rotational simulation, the tilt angle due to the flow component, γf, can fully explain the geometric tilt angle (Fig. 10). This is true for the rotational simulation only during the emergence phase (from 4 to 6 hours). After the emergence phase the tilt angle from the flows is systematically smaller than the geometric tilt angle by about γr ≈ 5°, but within the estimated uncertainty. The tilt angle in the non-rotational simulation is an estimate of the uncertainty in the tilt angle measurement, δγ ∼ 10°, by buffeting from convective flows (e.g. Wang & Sheeley 1989; Will et al. 2024).
4. Discussion
By implementing the f-plane approximation in the MURaM code, we have shown that the Coriolis force acts on the horizontal flows at the surface of the Sun, inducing radial vorticity in supergranules. We simulated a magnetic flux tube rising from 11 Mm beneath the solar surface that formed a magnetic bipole at the surface and generated significant outflows of about 2 km s−1 out to a radius of about 20 Mm from the initial emergence location. Our simulations with and without rotation formed small, low-flux magnetic bipoles that are consistent with the early stages of active region emergence.
By tracking the motion of the two polarities at the surface, we measured the separation and tilt angle as a function of time. Without rotation, the polarities emerged aligned with the initial flux tube axis. With a rotation rate about 100 times the Sun’s rotation rate, the polarities emerged tilted close to 30° in a clockwise direction away from the initial flux tube axis. In both cases, the tilt angle remained constant during the emergence process.
The tilt angle of the bipole formed in the rotational simulation, while consistent with Joy’s law, is inconsistent with observations in two respects. Firstly, the mean tilt angle of active regions on the Sun is about ∼6° with an uncertainty a factor of two larger (e.g. Hale et al. 1919; Dasi-Espuig et al. 2010). If we apply a scaling argument for the Sun with a 110 times slower rotation rate, surface velocities of about 100 m s−1 (compared to 2000 m s−1), and an emergence time of 2 days (compared to 2 hrs), the y separation (Δy ∼ ΩvxΔt2) becomes
Thus, γ⊙ is ∼7°, which is of the order of Joy’s law on the Sun (e.g. Hale et al. 1919; Wang & Sheeley 1989; Will et al. 2024). Secondly, although we have demonstrated that a tilt angle of the order of Joy’s law can be produced in the top 10 Mm of the convection zone, our simulations have not captured the observed time evolution through emergence. In the observations, active regions on the Sun emerge east-west aligned on average, and the tilt angle develops throughout the emergence process (Schunker et al. 2020).
Presumably, active regions on the Sun originate from the global toroidal magnetic field oriented in the east-west direction. In our simulations, the tilt away from the flux tube axis (along the east-west direction) generates a perpendicular component of the magnetic field (north-south direction), and this is consistent with being partly due to the effect of rotation on the near-surface horizontal convective flows.
We have shown that a bipole can emerge tilted away from the sub-surface orientation of the flux tube axis due to rotational forces acting on a flux tube as it rises through the near-surface convective flows. The change in tilt angle once the flux has reached the surface is consistent with a component due to the flows and a remaining component coming from a combination of the vertical velocity acting on the horizontal magnetic field and the effects of numerical diffusion in the simulation. Future simulations of flux emergence should improve upon our experiment by simulating a flux tube with a lower rise speed that ideally forms a stable, higher flux active region and by using a smaller grid size. To isolate the effect of the Coriolis force, we chose not to impose a twisted magnetic field in our flux tube; however, it may be an important component for the coherency of the flux tube and even for Joy’s law (e.g. Schuessler 1979; Longcope & Fisher 1996). These simulations are required to interpret the observational flow signatures as active regions emerge and the tilt angle develops.
Data availability
The simulations can be fully reproduced following the description in this paper. The version of the code including the f-plane approximation is available in the MURaM repository. The initial conditions for the simulations can be provided on request. The digital data presented in this paper are available through private communication with the authors. Movies associated to Figs. 4, 5, 7, and 8 are available at https://www.aanda.org
Acknowledgments
We acknowledge the sincere and helpful comments from the anonymous referee towards improving this paper. HS, WRB and LG acknowledge support from the project “Preparations for PLATO asteroseismology” DAAD project 57600926. This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. HS is the recipient of an Australian Research Council Future Fellowship Award (project number FT220100330) and this research was funded by this grant and an Australian Government Research Training Program (RTP) Scholarship from the Australian Government. RHC and LG acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 810218 WHOLESUN). HS, DIP and WRB acknowledge the Awabakal people, the traditional custodians of the unceded land on which their research was undertaken at the University of Newcastle.
References
- Babcock, H. W. 1961, ApJ, 133, 572 [Google Scholar]
- Birch, A. C., Schunker, H. S., Braun, D. C., et al. 2016, Sci. Adv., 2 [Google Scholar]
- Birch, A. C., Schunker, H., Braun, D. C., & Gizon, L. 2019, A&A, 628, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Brandenburg, A. 2005, ApJ, 625, 539 [NASA ADS] [CrossRef] [Google Scholar]
- Cameron, R., & Schüssler, M. 2015, Science, 347, 1333 [Google Scholar]
- Choudhuri, A. R., & D’Silva, S. 1990, A&A, 239, 326 [NASA ADS] [Google Scholar]
- Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [Google Scholar]
- Dasi-Espuig, M., Solanki, S. K., Krivova, N. A., Cameron, R., & Peñuela, T. 2010, A&A, 518, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- D’Silva, S., & Choudhuri, A. R. 1993, A&A, 272, 621 [NASA ADS] [Google Scholar]
- Fan, Y. 2009, Liv. Rev. in Sol. Phys., 6, 4 [Google Scholar]
- Fisher, G. H., Fan, Y., & Howard, R. F. 1995, ApJ, 438, 463 [NASA ADS] [CrossRef] [Google Scholar]
- Hale, G. E., Ellerman, F., Nicholson, S. B., & Joy, A. H. 1919, ApJ, 49, 153 [NASA ADS] [CrossRef] [Google Scholar]
- Hirzberger, J., Gizon, L., Solanki, S. K., & Duvall, T. L. 2008, Sol. Phys., 251, 417 [NASA ADS] [CrossRef] [Google Scholar]
- Irwin, A. W. 2012, Astrophysics Source Code Library [record ascl:1211.002] [Google Scholar]
- Karak, B. B., & Miesch, M. 2017, ApJ, 847, 69 [NASA ADS] [CrossRef] [Google Scholar]
- Kosovichev, A. G., & Stenflo, J. O. 2008, ApJ, 688, L115 [NASA ADS] [CrossRef] [Google Scholar]
- Kunasz, P., & Auer, L. H. 1988, J. Quant. Spectr. Rad. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
- Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Longcope, D. W., & Fisher, G. H. 1996, ApJ, 458, 380 [Google Scholar]
- McClintock, B. H., & Norton, A. A. 2016, ApJ, 818, 7 [NASA ADS] [CrossRef] [Google Scholar]
- Parker, E. N. 1955, ApJ, 121, 491 [Google Scholar]
- Rempel, M. 2014, ApJ, 789, 132 [Google Scholar]
- Rempel, M. 2017, ApJ, 834, 10 [Google Scholar]
- Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129 [Google Scholar]
- Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207 [Google Scholar]
- Schmidt, H. U. 1968, in Structure and Development of Solar Active Regions, ed. K. O. Kiepenheuer, IAU Symp., 35, 95 [Google Scholar]
- Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229 [Google Scholar]
- Schuessler, M. 1979, A&A, 71, 79 [Google Scholar]
- Schunker, H., Braun, D. C., Birch, A. C., Burston, R. B., & Gizon, L. 2016, A&A, 595, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Schunker, H., Birch, A. C., Cameron, R. H., et al. 2019, A&A, 625, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Schunker, H., Baumgartner, C., Birch, A. C., et al. 2020, A&A, 640, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Schunker, H., Roland-Batty, W., Birch, A., et al. 2024, MNRAS, 533, 225 [Google Scholar]
- Sreedevi, A., Jha, B. K., Karak, B. B., & Banerjee, D. 2023, ApJS, 268, 58 [Google Scholar]
- Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [Google Scholar]
- Wang, Y.-M., & Sheeley, N. R., Jr 1989, Sol. Phys., 124, 81 [NASA ADS] [CrossRef] [Google Scholar]
- Wang, Y.-M., & Sheeley, N. R., Jr 1991, ApJ, 375, 761 [NASA ADS] [CrossRef] [Google Scholar]
- Weber, M. A., Fan, Y., & Miesch, M. S. 2011, ApJ, 741, 11 [NASA ADS] [CrossRef] [Google Scholar]
- Weber, M. A., Schunker, H., Jouve, L., & Işık, E. 2023, Space Sci. Rev., 219, 63 [NASA ADS] [CrossRef] [Google Scholar]
- Will, L. W., Norton, A. A., & Hoeksema, J. T. 2024, ApJ, 976, 20 [Google Scholar]
- Witzke, V., Shapiro, A. I., Cernetic, M., et al. 2021, A&A, 653, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Validation of the f-plane approximation in MURaM
To test the implementation of the f-plane approximation in the MURaM code, we computed the expected effect on the horizontal flows of a supergranule at the surface. We described the radial profile of velocity of a supergranule at the surface as a sine function in time and space with a lifetime of T = π/ω and spatial scale of L = π/k, where k is the spatial scale and ω is the temporal scale. The velocity in the radial direction vr, is given by
for 0 ≤ t ≤ T, where r is the radius from the centre of the supergranule and V0 is the peak velocity. The Coriolis force acts perpendicular to the radial velocity in the plane, so we can write the acceleration due to the Coriolis force as
Integrating with respect to t gives
If we define vθ = 0 at t = 0 for all r, then
which leads to
The horizontal flow divergence (in polar coordinates centred on the supergranule) is given by
and the vertical component of the vorticity of the horizontal flows is
We weighted the vorticity by the horizontal flow divergence to favour the approximation in the supergranules, by multiplying Eq. A.7 by Eq. A.6 to give
The average product of divergence and vorticity over the lifetime of a supergranule is
Averaging over all surface area of the simulation gives
where we solved the integral was numerically. Assuming the lifetime of a supergranule is the time it takes for fluid to cross the diameter once, V0 = L/T = ω/k, and Eq. A.10 becomes
To obtain an estimate for the average product independent of the size or lifetime of a supergranule, we expressed the analytic product in terms of the root-mean-square horizontal flow divergence. We derived the mean-square divergence (from Eq. A.6) to be
leading to the root-mean-square divergence given by
We measured the root-mean-square divergence of the simulations without the f-plane approximation, and then substituted that value into Eq. A.11 to get
Figure A.1 shows that this analytic estimate for the spatially averaged product Eq. A.13 agrees within the uncertainties with the product computed for the simulations with the f-plane approximation. Note that the average product is a negative quantity since the value of f is negative due to the direction of rotation.
![]() |
Fig. A.1. Mean product of the vorticity and divergence at the surface for the MURaM hydrodynamic simulations with (red) and without (blue) the f-plane approximation. The expected analytic value of the product (green; Eq. A.13) is systematically offset but within a standard deviation of the mean of the product for the rotational hydrodynamic simulation with the f-plane approximation. |
All Figures
![]() |
Fig. 1. Example of a surface horizontal flow divergence map averaged over 6 hours and filtered between 15 Mm and 25 Mm from the non-rotational hydrodynamic simulation. The cyan contours indicate the identified supergranules, and the blue points indicate the location of the peak divergence within the boundary. |
In the text |
![]() |
Fig. 2. Top panel: horizontal flow divergence of the average supergranule in the non-rotational hydrodynamic simulation. The cyan arrows represent the horizontal flows at the surface. Middle panel: same as the top panel but for the rotational hydrodynamic simulation. Bottom panel: azimuthal averages of the radial velocity ( |
In the text |
![]() |
Fig. 3. Left: total magnetic field strength across the centre of the tube as a function of height. Right: Vertical slice through the centre of the flux tube showing the difference in density between the initial condition with the flux tube and the purely hydrodynamic initial condition. A reduced density (negative) will cause the tube to rise, and an increased density (positive) will cause it to sink. The contour of the magnetic field value of 25 kG is shown in green. The relative density perturbation is constant as a function of height in the background Model S atmosphere. |
In the text |
![]() |
Fig. 4. Slices through y = 0 of the magnetic field strength and vertical velocity of the initial condition (top two panels), and then at two subsequent time steps of the non-rotational (Ω = 0 Ω⊙) MHD simulation (bottom four panels). The flux tube rises in an arch structure before forming two dominant legs. In the velocity slices, blue indicates downward flow and red upward flow. See the movies of the simulation online and here. |
In the text |
![]() |
Fig. 5. Slices through y = 0 of the magnetic field strength and vertical velocity of the initial condition (top two panels), and then at two subsequent time steps of the rotational (Ω = 110 Ω⊙) MHD simulation (bottom four panels). In the velocity slices, blue indicates downward flow and red upward flow. See the movies of the simulation online and here. |
In the text |
![]() |
Fig. 6. Total unsigned magnetic flux at the surface as a function of time for both non-rotational and rotational MHD simulations. The vertical line indicates the approximate end time of the flux emergence (at 6 hours). |
In the text |
![]() |
Fig. 7. Vertical magnetic field at the surface (grey-scale maps) and slices through y = 0 Mm showing the unsigned magnetic field strength and vertical flows at three different times of the Ω = 0Ω⊙ non-rotational MHD simulation. The red and yellow symbols show the centres of the positive and negative magnetic poles, respectively. The crosses indicate the flux-weighted centres of the polarities over a threshold of |Bz| = 800 G, and the dots indicate the centres without a threshold. See the movies of the simulation online and here. |
In the text |
![]() |
Fig. 8. Same as Fig. 7 but for the110 Ω⊙ rotational MHD simulation. The red and yellow symbols show the centres of the positive and negative magnetic poles, respectively, used to calculate the geometric tilt angle, γg. See movies of the simulation online and here. |
In the text |
![]() |
Fig. 9. Surface flow velocities (cyan arrows) and vertical magnetic field averaged over time from t = 4.25 h to t = 5.44 h for the non-rotational (0 Ω⊙) hydrodynamic simulation (top) and the rotational (110 Ω⊙) hydrodynamic simulation (bottom). The velocity has also been spatially smoothed for clarity by a Gaussian with standard deviation of 2.4 Mm. Note that this has been computed at a constant geometrical height, at z = 480 pixels = 15.360 Mm from the bottom of the box. |
In the text |
![]() |
Fig. 10. Separation of polarities for the non-rotational simulation (blue) and the rotational simulation (red). We show the geometric separation without a Bz threshold (solid curves, Δxg) and the separation expected from the flows (Eqs. (11) and (12); dashed curves, Δxf). The vertical line at 6 hours is close to the time when most of the magnetic flux has emerged (see Fig. 6). |
In the text |
![]() |
Fig. 11. Tilt angle of the polarities for the non-rotational simulation (blue) and the rotational simulation (red). The geometric tilt angle is shown by the solid curves, and the tilt angle expected from the flows is shown by the dashed curves. The tilt angle resulting from the location of the polarities with a threshold of |Bz|> 800 G applied is shown by the thick shaded curves. |
In the text |
![]() |
Fig. A.1. Mean product of the vorticity and divergence at the surface for the MURaM hydrodynamic simulations with (red) and without (blue) the f-plane approximation. The expected analytic value of the product (green; Eq. A.13) is systematically offset but within a standard deviation of the mean of the product for the rotational hydrodynamic simulation with the f-plane approximation. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.