Successful and Failed Flux Tube Emergence in the Solar Interior

TL;DR

This study uses 3D MHD simulations to explore how various parameters influence the emergence of magnetic flux tubes from the solar interior to the surface, revealing complex dependencies beyond initial magnetic strength.

Contribution

The paper introduces a detailed parametric analysis of flux tube emergence, highlighting the roles of geometry and twist, and identifies scaling laws relating magnetic field strength to local density.

Findings

Emergence efficiency depends on magnetic field strength and morphology.

Flux tube emergence is not solely determined by initial magnetic flux.

Scaling laws relate magnetic field strength to local density with a transition near the pressure scale height.

Abstract

We report on our three-dimensional (3D) magnetohydrodynamic (MHD) simulations of cylindrical weakly twisted flux tubes emerging from 18 Mm below the photosphere. We perform a parametric study, by varying the initial magnetic field strength ($B_0$), radius ($R$), twist ($\alpha)$ and length of the emerging part of the flux tube ($\lambda$) to investigate how these parameters affect the transfer of the magnetic field from the convection zone to the photosphere. We show that the efficiency of emergence at the photosphere (i.e. how strong the photospheric field will be in comparison to $B_0$) depends not only on the $B_0$ but also the morphology of the emerging field and the twist. We show that parameters such as $B_0$ and magnetic flux cannot alone determine whether a flux tube will emerge to the solar surface. For instance, high-$B_0$ (weak-$B_0$) fields may fail (succeed) to emerge at the photosphere, depending on their geometrical properties. We also show that the photospheric magnetic field strength can vary greatly for flux tubes with the same $B_0$ but different geometric properties. Moreover, in some cases we have found scaling laws, whereby the magnetic field strength scales with the local density as $B\propto \rho^\kappa$, where $\kappa \approx 1$ deeper in the convection zone and $\kappa <1$, close to the photosphere. The transition between the two values occurs approximately when the local pressure scale ($H_p$) becomes comparable to the diameter of the flux tube ($H_p\approx2R$). We derive forms to explain how and when these scaling laws appear and compare them with the numerical simulations.

Figures

Tables

Table 1. Table 1:

Case	B0d	$R$d	$\lambda$d	$\alpha$d	$\Phi$
	($\times B_c$)	($\times H_c$)	($\times H_c$)	($\times H_c^{-1}$)	(Mx)
1*	3.4	3.2	100	0.1	$1.1\times {10}^{19}$
2*	3.4	5	100	0.1	$2.6\times {10}^{19}$
3*	3.4	7.6	100	0.1	$6.0\times {10}^{19}$
4*	3.4	10.1	100	0.1	$1.1\times {10}^{20}$
5*	17	3.2	100	0.1	$5.3\times {10}^{19}$
6*	17	5	100	0.1	$1.3\times {10}^{20}$
7*	17	7.6	100	0.1	$3\times {10}^{20}$
8	17	10.1	100	0.1	$5.3\times {10}^{20}$
9*	34	3.2	100	0.1	$1.1\times {10}^{20}$
10	34	5	100	0.1	$2.6\times {10}^{20}$
11	34	7.6	100	0.1	$6.0\times {10}^{20}$
12	34	10.1	100	0.1	$1.1\times {10}^{21}$
13	340	3.2	100	0.1	$1.1\times {10}^{21}$
14	340	5	100	0.1	$2.6\times {10}^{21}$
15	340	7.6	100	0.1	$6.0\times {10}^{21}$
16	340	10.1	100	0.1	$1.1\times {10}^{22}$

Table 2. Table 2:

Case	B0d	$R$d	$\lambda$d	$\alpha$d	$\Phi$
	($\times B_c$)	($\times H_c$)	($\times H_c$)	($\times H_c^{-1}$)	(Mx)
1*	17	5	100	0.1	$1.3\times {10}^{20}$
2*	24	5	100	0.1	$1.8\times {10}^{20}$
3	34	5	100	0.1	$2.6\times {10}^{20}$
4	68	5	100	0.1	$5.1\times {10}^{20}$
5	34	3.2	100	0.1	$1.1\times {10}^{20}$
6	34	7.6	100	0.1	$6.0\times {10}^{20}$
7	34	10.1	100	0.1	$1.1\times {10}^{21}$
8*	34	5	20	0.1	$2.6\times {10}^{20}$
9	34	5	35	0.1	$2.6\times {10}^{20}$
10	34	5	50	0.1	$2.6\times {10}^{20}$
11	34	5	100	0.15	$2.6\times {10}^{20}$
12	34	5	100	0.25	$2.6\times {10}^{20}$
13	34	5	20	0.11	$2.6\times {10}^{20}$
14	34	5	20	0.15	$2.6\times {10}^{20}$
15	34	5	20	0.25	$2.6\times {10}^{20}$
16*	34	3.2	50	0.1	$1.1\times {10}^{20}$
17*	44	5	10	0.1	$3.4\times {10}^{20}$
18*	24	5	20	0.1	$1.8\times {10}^{20}$
19*	24	5	20	0.25	$1.8\times {10}^{20}$
20	24	7.6	100	0.1	$4.3\times {10}^{20}$
21*	17	7.6	100	0.1	$3.0\times {10}^{20}$

Equations

$$\frac{\partial v_x}{\partial x}=\alpha ,\frac{\partial v_y}{\partial y}=\alpha ,\frac{\partial v_z}{\partial z}=ϵ\alpha ,$$

$$\frac{D\mathbf{B}}{Dt}$$

$$\frac{D\rho }{Dt}$$

$$\kappa =\frac{1+ϵ}{2+ϵ}.$$

$$\frac{\partial \rho }{\partial t}=$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}=$$

$$-\rho g_0\hat{\mathbf{z}}+{\mathbf{S}}_{visc},$$

$$\frac{\partial \rho ϵ}{\partial t}=$$

$$\frac{\partial \mathbf{B}}{\partial t}=$$

$$ϵ=$$

$${\mathbf{S}}_{visc}=\frac{\partial }{\partial x_j}\left({\sigma }_{ij}+{\sigma }_{ij}^{shock}\right){\hat{\mathbf{e}}}_i,$$

$$T(z)=T_{ph}+\frac{T_{cor}-T_{ph}}{2}\left(\tanh \frac{z-z_{cor}}{w_{tr}}+1\right),$$

$$T(z)=T_{ph}-\frac{{\mu }_{m}g}{k_B}\frac{\gamma -1}{\gamma }z$$

$$B_y$$

$$B_{\varphi }$$

$$P_{exc}=\frac{1}{2\mu }\left[{\alpha }^2\left(\frac{R^2}{2}-r^2\right)-1\right]B_y^2.$$

$$\Delta \rho =-{\rho }_0\frac{P_{exc}}{P_0}{e}^{-y^2/{\lambda }^2},$$

$$B^2\kappa =B_x^2{\kappa }_x+B_y^2{\kappa }_y+B_z^2{\kappa }_z,$$

$$B\gtrsim \sqrt{\frac{H_p}{R}}{B}_{eq},$$

$$\frac{D{B}_x}{Dt}$$

$$\frac{D{B}_y}{Dt}$$

$$\frac{D{B}_z}{Dt}$$

$${\kappa }_x=1-\frac{1}{\nabla \cdot \mathbf{v}}\frac{\partial v_x}{\partial x}$$

$${\kappa }_y=1-\frac{1}{\nabla \cdot \mathbf{v}}\frac{\partial v_y}{\partial y}$$

$${\kappa }_z=1-\frac{1}{\nabla \cdot \mathbf{v}}\frac{\partial v_z}{\partial z}.$$

$$\frac{D{B}^2}{Dt}$$

$$B=B_0{\left(\frac{\rho }{{\rho }_0}\right)}^{\kappa },$$

$$\frac{D{B}^2}{Dt}$$

$$B^2\kappa =B_x^2{\kappa }_x+B_y^2{\kappa }_y+B_z^2{\kappa }_z.$$

$$\kappa ={\kappa }_x.$$

$$\frac{\partial v_x}{\partial x}=\chi ,\frac{\partial v_y}{\partial y}=\psi ,\frac{\partial v_z}{\partial z}=\zeta ,$$

$$\kappa =\frac{1}{2}({\kappa }_x+{\kappa }_y)=1-\frac{1}{2}{\kappa }_z.$$

$$\frac{\partial v_x}{\partial x}=\frac{\partial v_y}{\partial y}=\chi ,\frac{\partial v_z}{\partial z}=\zeta$$

$$\kappa =\frac{2}{3}.$$

$$\frac{\partial v_x}{\partial x}=\frac{\partial v_y}{\partial y}=\frac{\partial v_z}{\partial z}=\chi .$$

$$\frac{D{B}^2}{Dt}$$

Full text

Sun: activity – Sun: interior – Sun: Magnetic fields –Magnetohydrodynamics (MHD) –methods: numerical

††journal: ApJ

1 Introduction

It is believed that the origin of the magnetic field of the Sun is associated with the existence of a dynamo mechanism operating around the base of the deep convection zone (Parker, 1955a). The magnetic fields rise from the 200 Mm deep convection zone towards the photosphere due to buoyancy (Parker, 1955b), where they can emerge and form a variety of magnetic structures (from small scale pores to large scale active regions). The expansion of the flux tubes during their emergence within the solar interior, depends mostly on the local density of the convection zone. However, the density inside the convection zone drops by six orders of magnitude, and the main density decrease occurs mostly in the upper convection zone. For example, the density drops by approximately $10^{4}$ in the upper 20 Mm of the convection zone, of which a $10^{3}$ drop occurs only in the upper 10 Mm. So, the local pressure scale height ($H_{p}$) is large and decreases slowly at larger depths. Therefore, deeper in the solar interior, the motion of magnetic fields (e.g. a flux tube) is affected less by pressure variations than near the surface. This allows the flux emergence process there to be studied using either the thin-flux tube approximation (e.g. Spruit, 1981; Caligari et al., 1995; Fan et al., 1993; Weber et al., 2011) or the anelastic MHD approximation (e.g. Brun et al., 2004; Fan, 2008; Jouve & Brun, 2009; Fan & Fang, 2014). Closer to the photosphere, on the other hand, $H_{p}$ becomes small and decreases rapidly. Hence, the size of the emerging structures can become comparable to $H_{p}$ and full 3D compressible MHD is needed in order to study flux emergence in the upper convection zone.

Toriumi & Yokoyama (2010) performed 2D MHD simulations of a magnetic flux sheet positioned at $z=-20$ Mm below the photosphere. They reported that in order for the magnetic field to emerge into the photosphere and above, its flux needs to be $10^{21}-10^{22}$ Mx. However, these fluxes were derived by assuming the length of the magnetic flux sheet along the third dimension. Toriumi & Yokoyama (2013) performed 3D MHD simulations of a magnetic flux tube originating from the same depth. They varied the initial magnetic field strength, twist and length of the buoyant part of the flux tube. They found that for higher values of the magnetic field strength and twist, the flux tube emerges faster inside the solar interior, and expands more dynamically above the photosphere. On the other hand, the flux tube emerges faster inside the solar interior, but expands less dynamically above the photosphere, when the buoyant part of the flux tube is longer. The above results are important to understand the emergence process of flux tubes in the solar interior. However, many questions remain open. For example, how the parameters of the initial sub-photospheric field affect the amount of flux emerging at the photosphere is still unknown.

Another interesting question is how the magnetic field strength ($B$) scales with the local plasma density ($\rho$) during the emergence of the flux tube within the convection zone. A simple scaling law can be derived if we assume a flux tube with a uniform axial magnetic field of strength $B$ and then vary its cross section ($A$), while keeping its length ($L$) constant (Fig. 1a). From conservation of mass ($M=AL\rho$) and magnetic flux ($\Phi=AB$), it is easy to show that $B\propto\rho$ (or $B\propto\rho^{\kappa}$ with $\kappa=1$).

Another scaling law can be derived if we vary the length of the flux tube while keeping its cross section constant (Fig. 1b). Conservation of mass and flux suggests that $B$ and $\rho$ will depend on the length of the flux tube. Useful information about the scaling can be derived by assuming that $B\propto\rho^{\kappa}$. Then, $\kappa$ becomes constrained to $0<\kappa<1$ (Pinto & Brun, 2013).

Finally, the scaling of the magnetic field strength with the local plasma density can be affected by the action of velocity gradients on the magnetic field. To show that, Cheung et al. (2010) assumed a horizontal magnetic field, $\mathbf{B}=B\;\mathbf{\hat{x}}$, (Fig. 1c). This field was then distorted by an asymmetric velocity gradient, defined by:

[TABLE]

where $\alpha$ is the horizontal expansion rate and $\epsilon$ is a parameter describing the asymmetry of the flow in the vertical direction. Combining the ideal induction and the continuity equations,

[TABLE]

they found that the scaling of $B$ with $\rho$ is indeed affected by the velocity gradients and that the power $\kappa$ depends on the degree of the asymmetry of the velocity gradients as

[TABLE]

For a purely horizontal expansion ($\epsilon=0$), $\kappa=0.5$ and for expansion transverse to the field ($\epsilon>>1$), $\kappa=1$, as expected from the conservation of flux and mass (as in Fig. 1a). For isotropic expansion ($\epsilon=1$), they found $\kappa=2/3$.

But which scaling is more suitable at different depths inside the convection zone? In the deeper parts of the convection zone, the local pressure scale height is large. The characteristic length of the flux tube’s cross section (e.g. radius, $r$) is $r<<H_{p}$ and the characteristic length of the emerging part (e.g. an axial perturbation, $l$) is $l>>H_{p}$ (i.e. a thin flux tube). So, the flux tube’s axial expansion would be small and its cross sectional expansion would be gradual and symmetric. As a result, velocity gradients along the axis would be small and the scaling of the magnetic field with the local density should follow $\kappa=1$.

Pinto & Brun (2013) studied the emergence of twisted flux tubes in a global dynamo model, using 3D anelastic MHD simulations. They found that $B\propto\rho^{\kappa}$, $\kappa=0.998\pm 0.001$. Similar behaviour was found in cases without a dynamo. When the emerging field was less buoyant, they found steeper, but similar slopes. Overall, they found values of $0.998<\kappa<1.002$ during the emergence of the flux tubes from 0.8 $R_{\odot}$ to 0.92 $R_{\odot}$. These results suggested that the poloidal component of the magnetic field dominated over the toroidal component and that the perturbations along the axis of the tube where indeed small.

In the upper parts of the convection zone, the length of the flux tube can increase significantly as the $\Omega$-loop shaped flux tube rises towards the photosphere. Moreover, the flux tube expands radially, as its cross section becomes comparable to $H_{p}$. Close to the photosphere, the flux tube experiences a significant horizontal expansion (Spruit et al., 1987) due to the rapid decrease of $H_{p}$. The flux tube cannot emerge above the photosphere until its magnetic forces dominate over the gas pressure forces and trigger a magnetic buoyancy instability (Acheson, 1979; Archontis et al., 2004). So, when the flux tube is underneath the photosphere, it becomes compressed and further expands horizontally, increasing its magnetic field strength until it is large enough to trigger the buoyancy instability. In addition to the above, the velocity gradients of the convective flows are larger closer to the photosphere than deeper in the solar interior. All the above should lead to a decrease of $\kappa$. Thus, the scaling of the magnetic field with the local density should follow $\kappa<1$ in the upper convection zone.

Cheung et al. (2010) compared their analytical result on $\kappa$, with a 3D radiative MHD simulation of the emergence of a toroidal flux tube, positioned $7.5$ Mm below the photosphere, inside a convective layer. They found that a value of $\kappa=0.5$ for the scaling of the magnetic field strength with the local density.

Cheung & Isobe (2014) suggested that the transition from $\kappa\approx 1$ (similar to Pinto & Brun, 2013) to $\kappa<1$ occurs somewhere during the rise of a flux tube from the deeper parts of the convection zone to the surface.

In this paper, we address a series of questions on the emergence of flux tubes from the convection zone to the photosphere. For this, we use 3D resistive and compressible MHD and assume a horizontal flux tube positioned at 18 Mm below the photosphere. The free parameters of our model are a) the initial magnetic field strength at the center of the flux tube, b) the twist, c) the radius and d) the length of the buoyant part of the flux tube. We perform a detailed parametric study to identify i) how $\kappa$ behaves when each of these parameters are varied, ii) where does the transition from $\kappa\approx 1$ to $\kappa<1$ occur, iii) what is the efficiency of emergence, namely what is the ratio of the photospheric field strength to the initial field strength, and iv) how to use the above in order to understand the initial conditions leading to “successful” flux emergence above the photosphere. Furthermore, we derive analytically the conditions under which $\kappa$ is constant.

In Sec. 2, we present the model and the initial conditions. In Sec. 3.1, we vary only the magnetic field strength and radius of the flux tube, in order to explore the parameter space and identify combinations of those two parameters leading to “successful” emergence of magnetic field above the photosphere. In Sec. 3.2, we describe analytically conditions under which $B$ scales with $\rho$, and compare our analysis with a numerical simulation of a “successful” emergence. In Sec. 3.3, we focus on one of the “successfull” emergence cases of Sec. 3.1 and perform a large parametric study to identify how each parameter affects the emergence of the field (magnetic field strength (Sec. 3.3.1), radius (Sec. 3.3.2), length of the buoyant part (Sec. 3.3.3) and twist (Sec. 3.3.4) ). In Sec. 3.3.5 we discuss all the results together, and present a “border-line” case that separates the “successful” and the “failed” emergence cases. In Sec. 4 we summarize and discuss our results.

2 Numerical Setup

To perform the simulations, we numerically solve the 3D time-dependent, resistive, compressible MHD equations in Cartesian geometry using the Lare3D code of Arber et al. (2001). The equations in dimensionless form are:

[TABLE]

where $\rho$, $\mathbf{v}$, $\mathbf{B}$ and P are density, velocity vector, magnetic field vector and gas pressure. Gravity is $g_{0}=274$ m s*-2*. We assume a perfect gas with ratios of specific heat $\gamma=5/3$. Viscosity is included through:

[TABLE]

where $\sigma_{ij}=2\nu\left(\varepsilon_{ij}-\frac{1}{3}\delta_{ij}\nabla\cdot\mathbf{v}\right)$ is the viscous stress tensor and $\sigma_{ij}^{shock}=\rho l(\nu_{1}c_{ms}+\nu_{2}l|s|)\left(\varepsilon_{ij}-\frac{1}{3}\delta_{ij}\nabla\cdot\mathbf{v}\right)$ is the shock tensor. In these tensors, $\varepsilon_{ij}=\frac{1}{2}\left(\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}\right)$ is the strain rate tensor, $\delta_{ij}$ is the Kronecker delta, $l$ is the distance across a grid cell in the direction normal to the shock front, $|s|$ is the rate of the strain tensor in the direction normal to the shock front and $c_{ms}=\sqrt{c_{s}^{2}+v_{A}^{2}}$ is the magnetosonic speed, with $c_{s}$ being the sound speed and $v_{A}$ being the Alfvén speed (more details in e.g. Arber et al., 2001; Bareford & Hood, 2015). The viscosity coefficients take the values $\nu=622$ kg m*-1* s*-1* ($0.01$ in non-dimensional units), and $\nu_{1}=0.1$ and $\nu_{2}=0.5$ (in non-dimensional units). Viscous heating is added through $Q_{visc}=\varepsilon_{ij}(\sigma_{ij}+\sigma_{ij}^{shock})$.

We use constant explicit resistivity of $\eta=0.01$ (non-dimensional units). Joule dissipation is added through $Q_{\mathrm{joule}}=\eta j^{2}$. The normalization is based on the photospheric values of density $\rho_{\mathrm{c}}=1.67\times 10^{-7}\ \mathrm{g}\ \mathrm{cm}^{-3}$, length $H_{\mathrm{c}}=180\ \mathrm{km}$ and magnetic field strength $B_{\mathrm{c}}=300\ \mathrm{G}$. From these we get temperature $T_{\mathrm{c}}=6230~{}\mathrm{K}$, pressure $P_{\mathrm{c}}=7.16\times 10^{3}\ \mathrm{erg}\ \mathrm{cm}^{-3}$, velocity $v_{\mathrm{0}}=2.1\ \mathrm{km}\ \mathrm{s}^{-1}$ and time $t_{\mathrm{0}}=85.7\ \mathrm{s}$.

The computational domain has a physical size of $72^{3}\mathrm{Mm}$ on a $600^{3}$ uniform grid. We assume periodic boundary conditions in the $y$ direction. Open boundary conditions are used in the $x$ direction. Open (closed) boundary conditions are assumed and at the top (bottom) of the numerical domain.

The temperature of the atmosphere ($z>0$) follows a tangential profile,

[TABLE]

where $T_{ph}=6100$ K, $T_{cor}=0.92$ MK, $z_{cor}=2.38$ Mm and $w_{tr}=0.18$ Mm. This results in an isothermal photospheric-chromospheric layer at $0\ \mathrm{Mm}\leq z<1.8\ \mathrm{Mm}$, a transition region at $1.8\ \mathrm{Mm}\leq z<3.2\ \mathrm{Mm}$ and an isothermal coronal at $3.2\ \mathrm{Mm}\leq z<45\ \mathrm{Mm}$. The atmospheric density is derived by numerically solving the hydrostatic equation $dP/dz=-g\rho$, having as boundary condition $\rho_{ph}=1.67\times 10^{-7}$ g cm*-3*. The atmosphere is field-free.

The solar interior ($-27\ \mathrm{Mm}\leq z<0\ \mathrm{Mm}$) is convectively stable and in hydrostatic equilibrium. The temperature profile of the interior increases linearly, with depth with the constant temperature gradient given by:

[TABLE]

where $\mu_{m}=m_{f}m_{p}$ is the reduced mass, $m_{p}$ is proton mass, $m_{f}=1.2$, and $k_{B}$ is Boltzmann’s constant. The density in the interior is calculated by solving the hydrostatic equation with boundary condition $\rho_{ph}$. This stratification (sometimes with different $m_{f}$) is commonly used in flux emergence simulations of a fully ionized convectively stable solar interior (e.g. Fan, 2001; Manchester et al., 2004; Archontis et al., 2004; Moreno-Insertis et al., 2008; Toriumi & Yokoyama, 2011; Leake et al., 2013; Syntelis et al., 2015, 2017). The initial distribution of temperature ($T$), density ($\rho$), gas pressure ($P_{\mathrm{g}}$) of the interior and the atmosphere is shown in Fig. 2. The gas pressure of the interior at $-20$ Mm ($-10$ Mm) is 1.1$\times 10^{4}$ (3$\times 10^{3}$) larger than the photospheric one.

We place a cylindrical flux tube at $z=-18$ Mm, oriented along the $y$-axis. The magnetic field of the flux tube is defined as:

[TABLE]

where $R$ is a measure of the flux tube’s radius, $r$ the radial distance from the flux tube’s axis and $\alpha/2\pi$ is the twist per unit of length. The magnetic pressure ($P_{m}$) of a flux tube with $B_{0}=34$, $R=5$ is over-plotted in Fig. 2 (black line).

The background solar interior has a pressure, temperature and density profile of $P_{0}$, $T_{0}$ and $\rho_{0}$. When adding the flux tube, we introduce a pressure excess due to the magnetic field. By requiring the flux tube to be in radial force balance (i.e. ($\nabla P)\cdot\mathbf{\hat{e}_{r}}=(\mathbf{j}\times\mathbf{B})\cdot\mathbf{\hat{e}_{r}}$), we find the excess pressure $P_{exc}$ to be (see details in Murray et al., 2006):

[TABLE]

So, to set the flux tube in force balance with the background, we set the gas pressure inside the flux tube ($P_{i}$) to be $P_{i}=P_{0}-P_{exc}$. To initiate the flux tube emergence, we assume the flux tube is in thermal equilibrium with the background ($T_{i}=T_{0}$) and this leads to a difference of density of $\Delta\rho=\rho_{i}-\rho_{0}=-\rho_{0}P_{exc}/P_{0}$ between the flux tube interior and the non-magnetized background plasma (density deficit), which makes the flux tube buoyant. To avoid emerging the whole length of the flux tube, we reduce the density deficit towards the flanks of the flux tube by (Fan, 2001):

[TABLE]

where $\lambda$ is thus a measure of the length of the buoyant part of the flux tube. The above ensures that the middle part of the flux tube will be buoyant, while the flanks will not. Therefore, during the emergence, the flux tube will adopt an $\Omega$-loop shape during its emergence.

The values of the parameters used for our parametric numerical study are show in Tables 1, 2. From now on, we will refer to the dimensionless values of a variable using the subscript “d” (e.g. $B_{0}$d will be the dimensionless initial magnetic field strength).

3 Results

3.1 Magnetic flux and emergence

First we study the emergence of flux tubes by varying their initial magnetic flux from 1019 Mx up to 1021 Mx. To change the initial magnetic flux we vary both the magnetic field strength and the radius of the flux tube. We select $B_{0}=1,\,5,\,10,\,100$ kG ($B_{0}$d$=3.4,\,17,\,34,\,340$) and $R=0.6,\,0.9,\,1.4,\,1.8$ Mm ($R$d$=3.2,\,5,\,7.6,\,10.1)$. The combination of these values produce 16 cases, shown in Table 1. Fig. 3 shows the resulting values of the initial flux as a function of $B_{0}$. In all cases, the initial twist is low ($\alpha$d=0.1). In general, increasing the twist keeps the flux tube more coherent and assists the emergence process (e.g. Moreno-Insertis & Emonet, 1996; Murray et al., 2006; Toriumi & Yokoyama, 2011). We choose the length scale of the buoyant part of the flux tube to be relatively large ($\lambda$d$=100$). As a result, the apex of the flux tube will adopt a horizontal-like shape during its emergence.

Archontis et al. (2004) found that a flux tube will emerge from the solar interior into the solar atmosphere when a magnetic buoyancy instability (Acheson, 1979) is triggered. Before the instability is triggered, the emerging field cannot penetrate the solar surface and, instead compresses significantly just below the photosphere. The increase in magnetic field strength at this location reduces the plasma $\beta$ inside the flux tube and the instability is triggered when the plasma $\beta$ drops below unity. In our simulations, to classify a case as “successful” or “failed” emergence at and above the photosphere, we use the following criteria. If the buoyancy instability criterion (see e.g. Acheson, 1979; Archontis et al., 2004) measured at the photosphere is satisfied and the photospheric magnetic field is at least 100 G, then we consider the emergence as “successful”. If the rising flux tube reaches the photosphere, but its plasma $\beta$ remains very high (e.g. $\geq 100$) and does not decrease considerably over several (at least 100) Alfvén times, we consider the emergence as “failed”.

Not surprisingly, all the cases with $B_{0}=100$ kG emerge “successfully” (Fig. 3). Most of the flux tubes with $B_{0}=10$ kG, manage to emerge “successfully”, with the exception exception being the thinnest of these flux tubes ($R$d$=3.2$). The magnetic field strength at its center drops significantly, resulting in an internal $\beta\approx 4000$ plasma just below the photosphere. For two flux tubes with same $B_{0}$ but different radius, at $r=R_{large}$ and $r=R_{small}$, the pressure difference between the interior and the exterior of the tube will be the same. Thus, the expansion rate of the flux tubes will be, at least initially, the same. However, as the flux tubes expand, the cross-sectional area of the smaller flux tube grows more (as a percentage of the cross-sectional area at $t=0$). Due to conservation of flux, the magnetic field strength of a smaller radius flux tube will decrease more in comparison to a larger radius flux tube. Therefore, its magnetic pressure will decrease faster and it will bring higher plasma $\beta$ material close to the photosphere.

For $B_{0}=5$ kG, only the largest flux tube radius manages to emerge above the photosphere. All the $B_{0}=1$ kG cases fail to emerge. They rise very slowly and end up reaching force balance inside the solar interior, with a very large $\beta$. In these cases, the magnetic field brought below the photosphere is very low and the buoyancy instability is never triggered.

Notice that some cases “successfully” emerge (e.g. $B_{0}=10$ kG and $R$d$=5$) while other cases with a similar flux but different $B_{0}$ and $R$ fail to emerge ($B_{0}=1$ kG and $R$d$=10.1$, $B_{0}=5$ kG and $R$d$=7.6$). Despite the substantial initial flux (greater than $10^{20}$ Mx), these two flux tubes are not buoyant enough to emerge “successfully”. Therefore, we conclude that the initial magnetic flux within the rising magnetic structure cannot indicate directly whether the magnetic structure will emerge.

From Fig. 3 we find that in some cases an increase of $R$ (for constant $B_{0}$) leads to “successful” emergence. Such cases are the $B_{0}=5$ kG and $R$d$=10.1$ $B_{0}=10$ kG and $R$d$=5$. In Sec. 3.3 we will present the results of a parametric study on $B$, $R$, $\alpha$, and $\lambda$ in the latter case, to identify how these parameters affect the emergence. However, it is important first to show how the magnetic field strength varies with the local density, during the emergence of the flux tubes in the solar interior. This is discussed in the next section.

3.2 Scaling of magnetic field strength with the local density

To study how the magnetic field strength scales with the local plasma density, we use the following approach. We examine only the field at the $xz$-midplane, which is the cross section of the middle part of the flux tube. Notice how the histogram of $B$ and $\rho$, at that plane and shown in Fig. 4, evolves in time. As the flux tube emerges (panels a to c), this histogram is shifted towards lower values of density and field strength. To track the overall change of the field strength with the local density, we make the histogram of all the values of $B$ and $\rho$, from $t=0$ until the end of the simulation (an example is shown in Fig. 5a). We then plot the line that outlines the uppermost part of the histogram (black line). This line highlights how the maximum magnetic field strength scales with the local density (the undulations of this line are due to the snapshots frequency of the simulation: the higher the frequency, the smoother the line). We will refer to such lines as scaling curves.

We follow this process for the Table 2 cases 1-15, which explore the $B_{0}$, $R$, $\alpha$, and $\lambda$ parameter space around the $B$d$=34$ ($B_{0}=10$ kG), $R$d$=5$ case of Fig. 3. We show their scaling curves in Fig. 5b-f. Notice that in most of the cases shown in Fig. 5, the scaling curves consist of two major parts: a less steep part ($\log\rho\in[-7,-6]$, i.e. $-7<\log\rho<-6$, where $\rho$ is in g cm*-3*) and a more steep part ($\log\rho\in[-5,-4.3]$). We identify the mean inclination ($\kappa$) of these parts by performing linear fits ($\log B=\kappa\log\rho+c$). The value of $\kappa$ is shown in each panel, inside the parenthesis next to the value of the varied parameter.

As discussed in the Introduction, Cheung & Isobe (2014) suggested that the scaling curve will change from a steeper ($\kappa=1$) to a less steep ($\kappa<1$) power law during the emergence within the solar convection zone. Fig. 5b-e shows a similar transition in our numerical experiments. However, Fig. 5f shows a number of simulations where the scaling curves are not linear at all. Why do some of the scaling curves develop power-laws while others behave non-linearly? How does the steep and less steep part of the scaling curves develop? We first address these questions and then discuss how the variation of each flux tube parameter modifies the scaling curves and affects the emergence.

3.2.1 Derivation of scaling laws and comparison with simulation

In Appendix A, we derive forms to explain under which conditions the magnetic field strength scales with the local density as $B\propto\rho^{\kappa}$. Here we summarize these results.

First, we assume a velocity field with no shearing terms. We demand the magnetic field strength to be written as $B\propto\rho^{\kappa}$, where $\kappa$ is a constant. Combining the induction and continuity equation, we get Eq. A6, which we write again here:

[TABLE]

where $\kappa_{x}$, $\kappa_{y}$, $\kappa_{z}$ (Eq. A2a, A2b, A2c respectively) are functions of the non-shear velocity gradients. We then find solutions for $\kappa$ that satisfy the above equation. To do so, we focus on the following three cases.

The first case is when the magnetic field can be described with one component (e.g. $B_{i}\gg B_{j},B_{k}$, where $i,\,j,\,k$ are the components of the field). Then, the magnetic field strength scales as $B\propto\rho^{\kappa_{i}}$. For instance, if $B_{x}\gg B_{y},B_{z}$, then $\kappa=\kappa_{x}$, whereas if $B_{y}\gg B_{x},B_{z}$, then $\kappa=\kappa_{y}$. The assumption of $\kappa$ being constant constrains $\kappa_{i}$ to be constant too. Therefore, the velocity field is constrained to be $\frac{\partial v_{x}}{\partial x}=\chi,\,\frac{\partial v_{y}}{\partial y}=\psi,\,\frac{\partial v_{z}}{\partial z}=\zeta,$ where $\chi,\,\psi,\,\zeta$ are constants. Notice that the expression of $\kappa$ from Cheung et al. (2010) (see Introduction) is a special case of the above. For their assumption that $\mathbf{B}=B\,\mathbf{\hat{x}}$, the field will scale with $\kappa_{x}$. Assuming the velocity gradients to be in the form of Eq. 1, then Eq. A2a will give Eq. 4.

The second case is when the magnetic field can be described with two components of the full magnetic field vector (e.g. $B_{i},\,B_{j}\gg B_{k}$). Then, the magnetic field strength scales as $B\propto\rho^{\kappa},\,\kappa=1-\frac{1}{2}\kappa_{k}$. For instance, if $B_{x},\,B_{y}\gg B_{z}$, then the magnetic field strength will scale with $\kappa=1-\frac{1}{2}\kappa_{z}$ (Eq. A9). In this case, the constraint imposed on the velocity field will be stricter ($\frac{\partial v_{x}}{\partial x}=\frac{\partial v_{y}}{\partial y}=\chi,\quad\frac{\partial v_{z}}{\partial z}=\zeta$, where $\chi,\,\zeta$ are constants). This constraint guarantees that the magnetic field vector described by $B_{x}$ and $B_{y}$ does not change direction. As a result, in a high $\beta$ plasma (like a flux tube in the solar interior), such a velocity field will force the two-component field to behave as a one-component field. Therefore, in the field aligned coordinate system, this case is a special case of the first one.

The third case is when all three components of the magnetic field vector are important to describe the field. Then all $\kappa_{x}$, $\kappa_{y}$, $\kappa_{z}$ are equal. The velocity field constraint then becomes $\frac{\partial v_{x}}{\partial x}=\frac{\partial v_{y}}{\partial y}=\frac{\partial v_{z}}{\partial z}=\chi$, where $\chi$ is constant. So, the strength of the field will scale with the local density only if the field expands isotropically. In that case, $\kappa=\frac{2}{3}$ (Eq. A11). In the field aligned system the three-component field behaves as a one-component field. Therefore this is also a special case of the first one.

From the above, we infer that in order to express the magnetic field strength as $B\propto\rho^{\kappa}$, the magnetic field needs to have one dominant direction. Otherwise, a power law cannot be derived in general.

In the above, we assumed that the shearing terms of the velocity field equal to zero, i.e. $\frac{\partial v_{i}}{\partial x_{j}}=0$, $i\neq j$. Assuming constant $\kappa$ and a velocity field with non-zero shearing terms, we deduce that the magnetic field strength will scale with the local density with $\kappa=\frac{1}{B^{2}}\kappa_{ij}B_{i}B_{j}$ (Eq. A14), where $\kappa_{ij}$ is a tensor describing the deformation of the velocity field, given by Eq. A15. Because we assumed constant $\kappa$, the components of the tensor are required to be constant too. Note that the previous expressions derived for zero shear velocities are special cases of this general expression.

So far, we have assumed a strict power law between $B$ and $\rho$ (i.e. constant $\kappa$), which led to the constraint that the gradients of the velocity field components are constants. In general, the velocity gradients would be expected change during the emergence of a field. Assuming a non constant $\kappa$, we derived that $\kappa$ can be described by Eq. A17. However, if $\kappa$ changes slowly both in space and time ($\frac{D\kappa}{Dt}\approx 0$), we get that $\kappa\approx\frac{1}{B^{2}}\kappa_{ij}B_{i}B_{j}$. Therefore, the latter expression for $\kappa$ can describe the scaling of the $B$ with $\rho$, not only when $\kappa$ is constant, but also when $\kappa$ is changing slowly.

Our analysis suggests that when the velocity gradients change rapidly in space and/or time, or when the magnetic field cannot be adequately described by one component of the full magnetic field vector, $\kappa$ will not be constant. In that case, a power law between $B$ and $\rho$ should not be expected to occur.

We will now discuss where in our simulations we find conditions favouring the formation of power laws. We focus on where the magnetic field has a dominant component during its emergence within the solar interior. Notice that closer to the center of the flux tube, the poloidal component of the field becomes less significant than the axial one ($B_{\phi}/B_{y}$ decreases for smaller $r$, Eq. 14, 15). During the emergence of the flux tube, the shape of the field is crucial for the development of dominant field directions. In our numerical experiments, the length scale of the buoyant part of the flux tube ($\lambda$) is the parameter that affects the shape of the apex of the emerging flux tube the most. In Fig. 6a, top, we show a cartoon-like illustration of the upper part of a flux tube (oriented along the $y$-axis) with large $\lambda$. In this case, because of the large value of $\lambda$ the apex of the emerging tube is almost horizontal, oriented along the $y$-axis. Close to the axis of the flux tube (grey line in Fig. 6a top, grey shaded cross sectional region in Fig. 6a bottom left) the axial ($B_{y}$) component of the field will be dominant across a length $h_{large}$. The magnetic field in this region is, therefore, expected to scale with $\kappa=\kappa_{y}$, if the velocity field is changing slowly. It is important to note that the magnitude of the magnetic field is stronger close to the flux tube axis (Eq. 14). Consequently, the region around the axis contains the bulk of the magnetic energy of the flux tube and, therefore, has the most important role in the transfer of that energy to the photosphere. Away from the center, (thick black line in Fig. 6a top, black shaded region Fig. 6a bottom right) the poloidal component of the field becomes important. There, the field strength should not be expected to scale with the local density in general.

During the emergence process, the tube expands and, hence, its radius increases. Parker (1974) showed that the radial expansion of a flux tube causes the poloidal component of the field to increase when the twist remains constant (i.e. for a tube oriented along the $y$-axis, $B_{\phi}/B_{y}$ increases). We do find $B_{\phi}/B_{y}$ increases, in agreement with Parker. Eventually, this effect would result in a decrease of the size of the region close to the axis that scales with $\kappa_{y}$ (shaded region in Fig. 6a bottom left).

Fig. 6b shows an illustration of the flux tube when its apex reaches the photosphere. The upper part of the tube (shaded region) undergoes compression and horizontal expansion. If compressed enough, this region will develop locally a strong $B_{x}$ component. The $B_{x}$ component can eventually become much stronger than the local $B_{y}$ component of the field ( $B_{x}\gg B_{y},\,B_{z}$). Then, the magnetic field strength inside the compressed region will scale with $\kappa=\kappa_{x}$ (and not with $\kappa=\kappa_{y}$, which was the scaling exponent during the rise of the flux tube deeper in the solar interior, where $B_{y}\gg B_{x},\,B_{z}$ ). Note that the axis of the flux tube might not be inside the compressed region. In fact, in our simulations the center of the tubes is found at lower heights.

Fig. 6c (top) shows a flux tube that develops a highly bent apex when $\lambda$ is small. In this case, only a small segment of the apex (with a horizontal size $h_{small}$, Fig 6c bottom) could be oriented parallel to the photosphere, adopting a horizontal-like configuration. Also, due to the highly bent apex, plasma draining is expected to be more profound in this case, which could develop strong variations in the velocity gradients. Thus, in this case, we should not expect $B$ and $\rho$ to scale with a power law, except if the tube undergoes significant compression at the photosphere and adopts a similar configuration to the case described in Fig. 6b.

To use the above analysis towards studying the results of the simulations, we select an experiment with strong $B_{0}$ and large $\lambda$ (case with $B_{0}=20$ kG ($B_{0}$d$=68$), $\lambda$d$=100$, $\alpha$d$=0.1$ and $R$d$=5$) and find its scaling curve (Fig. 7a, solid line). We shift the scaling curve by $\Delta\log B=0.1$ (dashed line) and we take into account all the points between the two curves. These are the points with very high $B$. For these points, we plot the distributions of the absolute value of each individual component of the field ($|B_{x}|$ is the blue, $|B_{y}|$ is the magenta and $|B_{z}|$ is the yellow distribution, Fig. 7b). Then, we plot the distribution of $\kappa_{x}$, $\kappa_{y}$, $\kappa=1-\frac{1}{2}\kappa_{z}$ and $\kappa=\frac{1}{B^{2}}\kappa_{ij}B_{i}B_{j}$ (Eq. A2a, A2b, A9, A14), calculated directly from the velocity field (oragne color, Fig. 7c, d, e, f). The diamond symbols show the mean value of the distributions (i.e. mean value of $\kappa$) at each density bin, and the error bars show the standard deviation. The black line in panels (c)-(f) show the derivative of the scaling curve (i.e. the $\kappa$ measured from the histogram).

Fig. 7b shows that, for $\log\rho\in[-5.4,-4.35]$ (meaning $-5.4<\log\rho<-4.35$, where $\rho$ is in  g cm*-3*) or in terms of height $z\in[-18,-3]$ Mm), the strongest component of the magnetic field is $B_{y}$ (purple between second and third vertical line, in comparison to cyan and yellow). The large $\lambda$ ensures that the apex will be locally horizontal along a relatively large region, similar to Fig. 6a Therefore, when $\log\rho\in[-5.4,-4.35]$, the steep power law segment of the scaling law should be described by $\kappa=\kappa_{y}$. Indeed, in Fig. 7d, the values of $\kappa_{y}$ measured from the velocity field (orange) and $\kappa$ measured from the gradient of the scaling curve (black line) are in a relatively good agreement.

The small buildup of $B_{x}$ and $B_{z}$ when $\log\rho\in[-5.4,~{}-4.9]$ (cyan and yellow, Fig. 7b) is due to the expansion of the flux tube (which increases $B_{\phi}/B_{y}$). However, not many points have comparable $B_{x}$ and $B_{y}$. We should highlight that Fig. 7b is a true color image and the colors blend proportional to the value on the histogram. When comparable number of points have similar $B_{x}$ and $B_{y}$, cyan becomes purple. Therefore, when $\log\rho\in[-5.4,-4.35]$, the $B_{\phi}/B_{y}$ increase during the expansion of the flux tube is not significant, and it does not affect the steep power law.

For values in the range $-6<\log\rho<-5.4$ (or $z\in[-3,-1]$ Mm) the steepness of the scaling curve changes, revealing a transition to another regime with a different power law dependence between $B$ and $\rho$ (Fig. 7a, between first and second vertical line). During that transition, $B_{x}$ increases and becomes comparable to $B_{y}$ (purple color). The comparison between $\kappa$, deduced from the scaling curves, and the expression $\kappa=1-\frac{1}{2}\kappa_{z}$ (which is derived when both $B_{x}$ and $B_{y}$ are important) is in agreement at these depths (Fig. 7e, orange histogram and black line between first and second vertical lines).

For $\log\rho\lesssim-6$ or $z\gtrsim-1$ Mm, the apex of the flux tube is compressed significantly and $B_{x}$ becomes the strongest component of the magnetic field (Fig. 7b, cyan before first vertical line, in comparison to purple and yellow), as it is schematically illustrated Fig. 6b. There, we find the less steep power law of the scaling curve. Since $B_{x}$ is significantly stronger than the other components of the field, the magnetic field strength is expected to scale with the local density raised to the power $\kappa=\kappa_{x}$. Indeed, in Fig. 7c, the values of $\kappa_{x}$ measured from the velocity field and the $\kappa$ measured from the scaling curve are in good agreement.

Finally, we compare $\kappa$ measured from the scaling curve with the more general expression which includes velocity shear, $\kappa=\frac{1}{B^{2}}\kappa_{ij}B_{i}B_{j}$, calculated directly from the velocity field (Fig. 7f). We find that they are overall in agreement.

It is important to note that, for the derivation of the expressions of $\kappa$, we assumed that the velocity gradients are either constant or change slowly. For $\kappa_{x}$, $\kappa_{y}$ and $\kappa_{z}$, we also assumed zero shearing velocities. In the simulation, the velocity gradients are not changing slowly close to the photosphere. Also, the $\frac{\partial v_{y}}{\partial x}$ shear is significant when $\log\rho\lesssim-5.4$. However, the expressions of $\kappa$ shown in Fig. 7c, d, e, f, are in agreement with the values measured from the scaling curve. Therefore, we conclude that the most important parameter for the development of the power laws is a strong locally horizontal field across a large region, and not the strict velocity field constraints. However, we expect that for significant variations of the velocity gradients, which can perturb the direction of the magnetic field, it is not possible to form at power law.

The effects that the resolution, resistivity and viscosity have on the scaling curve of the studied case is discussed in Appendix B.

3.3 Height-time profiles and scaling curves of parametric study

In the following, we will study how the initial parameters of the flux tube (e.g. $B_{0}$, $R$, $\lambda$, $\alpha$) affect the emergence to the photosphere and above. For that, we focus on the cases 1-15 of Table 2, which explore the parameter space around the $B_{0}=10$ kG and $R$d$=5$ point of Fig. 3. We will study the emergence, focusing on the height-time profiles and the scaling curves of the emerging fields.

To plot the height-time profiles, we follow the rising motion of two points of the emerging flux tubes. The first one is the center of the flux tube, which is the point where $B_{y}$ is maximum and $B_{x}$ changes sign, along the $z$-axis at the center of the numerical box. The second one is the apex of the rising flux tube, which we consider to be the uppermost point along the $z$-axis at the center of the numerical domain, where $B>0.001B_{0}$. The profiles are plotted in Fig. 8 with solid (apex) and dashed (center) lines.

3.3.1 Variation of magnetic field strength

Firstly, we focus on the dynamics of the emerging flux tube when the magnetic field strength is varied and the other parameters are kept constant. We select $B_{0}=5,\,7.2,\,10,\,20$ kG ($B_{0}$d$=17,\,24,\,34,\,68$) and $\alpha$d$=0.1$, $\lambda$d$=100$, $R$d$=5$ (Table 2, cases 1-4).

The height-time profiles are shown in Fig. 8a. It is clear, that the stronger the field strength the faster the flux tube will rise inside the solar interior. Notice that the $B_{0}$d$=68$ case emerges above the photosphere almost immediately. In comparison, the $B_{0}$d$=34$ case exhibits a phase of deceleration before it emerges above the photosphere (during which the magnetic field at the apex is locally compressed). This is consistent with the results reported in previous studies (e.g. Fan, 2001; Archontis et al., 2004; Toriumi & Yokoyama, 2013). For lower $B_{0}$, the buoyancy of the flux tubes decreases and, thus, their center reach lower heights in the convection zone. The cases where $B_{0}$d$=17,\,24$ “failed” to emerge.

If we scale the time as $t\times B_{0}$ (Fig. 8b), we find that the height-time profiles are “clustering” closer together indicating self-similar behaviour (Murray et al., 2006; Sturrock & Hood, 2016). Still, the “clustering” is not as “compact” as in the previous studies. In our simulations, the flux tubes emerge from much deeper down in the solar interior. Thus, the downwards tension force becomes higher, reducing the upwards buoyancy force.

Fig. 5b shows the scaling of $B$ with $\rho$. We focus on the steeper part of the scaling curves. Increasing $B_{0}$ decreases the value of $\kappa$ (from $\kappa$=2.2 in the $B_{0}$d$=17$ case to $\kappa$=0.76 in the $B_{0}$d$=68$ case). Therefore, flux tubes with higher $B_{0}$ emerge more efficiently. In the “failed” emergence cases, the central part of the emerging fields reach moderate heights within the convection zone (around -9 Mm and -5 Mm for $B_{0}$d$=17$ and $B_{0}$d$=24$ respectively, Fig. 8a). The apexes move slowly upwards, but never emerge through the photosphere. Because of the lower $B_{0}$, these flux tubes do not undergo a 3D full expansion but mainly experience a vertical stretching in the following manner. The lower segments of the buoyant part of the flux tubes remain almost anchored at the initial depth. The rest of the tube emerges slowly, causing the vertical stretching. This stretching leads to faster decrease of the axial field strength and as a result a higher $\kappa$ ($\kappa>1)$.

The transition to the less steep part of the scaling curves occurs when the flux tubes are close to the photosphere (around -4 Mm or $\log(\rho)=-5.3$). There, the scaling curve transitions from scaling with $\kappa=\kappa_{y}$ to scaling with $\kappa=\kappa_{x}$. The “failed” emergence cases with $B_{0}$d$=17$ and $B_{0}$d$=24$ do not experience significant compression and, therefore, do not develop the less steep slope.

3.3.2 Variation of radius

Next, we focus on the dynamics of the emerging flux tubes when their radius is varied. We select $R$d$=3.2,\,5,\,7.6,\,10.1$, and $\alpha$d$=0.1$, $\lambda$d$=100$, $B_{0}=10$ kG ($B_{0}$d$=34$) (Table 2, cases 3, 5, 6, 7). The height-time profiles of these cases are shown in Fig. 8c. Notice that the larger the radius the faster and higher the flux tube’s apex and center will rise. At $t=0$, all the $B_{0}=10$ kG flux tubes are equally buoyant (buoyancy $\propto B^{2}$) at their centers (where $B=B_{0}$). However, they are not equally buoyant away from their centers, as $B_{R_{large}}>B_{R_{small}}$ when $r>0$ (see Eq. 14, 15). Thus, a larger radius tube will be more buoyant across its whole cross section.

As discussed in Sec. 3.1, the magnetic pressure of smaller radius flux tubes will decrease faster. This can be seen in the steeper part of the scaling curves in Fig. 5c. Notice that $\kappa$ decreases as $R$d increases, both at the steeper ($\kappa=2.3,\,0.98,\,0.86,\,0.77$ for the $R$d$=3.2,\,5,\,7.6,\,10.1$ cases respectively) and at the less steep part of the scaling curve ($\kappa=0.23,\,0.16,\,0.15$ for the $R$d$=5,\,7.6,\,10.1$ cases). Overall, higher $R$ leads to more efficient emergence. This is also reflected in the time needed for the flux tube to emergence above the photosphere. For instance, in Fig. 8, the $R$d$=10.1$ (green) flux tube emerges almost directly in comparison to the $R$d$=5$ (blue) case. Therefore, the radius of the tube is an important parameter affecting the dynamics of the emergence.

Notice that the point where the scaling curves transitions from the steep to the less steep power law is different for each case. It can be traced approximately at $\log(\rho)=-5.3,\,-5,\,-4.7$, where $H_{p}\approx 2R$ (local pressure scale of $H_{p}$d$=9,\,14,\,22$) for the $R$d$=5,\,7.6,\,10.1$ cases respectively).

3.3.3 Variation of $\lambda$

We now focus on the dynamics of the emerging flux tubes when the length of their buoyant part is varied. We select $\lambda$d$=20,\,35,\,50,\,100$ and $\alpha$d$=0.1$, $B_{0}$d$=34$ and $R$d$=5$ (Table 2, cases 3, 8, 9, 10).

The buoyant part of the flux tube becomes more bent for smaller $\lambda$, resulting to higher downward magnetic tension at its apex. Due to that, smaller $\lambda$ flux tubes emerge slower (Fig. 8d) (e.g. Schuessler, 1979; Longcope et al., 1996; Moreno-Insertis & Emonet, 1996). The $\lambda$d$=35,\,50,\,100$ results are consistent with the results of previous studies (e.g. Fan, 2001; Syntelis et al., 2015).

However, the $\lambda$d$=20$ case (Fig. 8d , black line) behaves differently. This is a case of a “failed” emergence. Initially, the flux tube rises for a time period of about $t=1000$ minutes (black solid and dashed line). Then, the emerging flux system enters a short phase of deceleration (i.e. from $t=1000$ minutes until $t=1400$ minutes), during which the downward tension force of the envelope field lines becomes comparable to the magnetic pressure force. At the same time, plasma draining from the apex of the tube towards its flanks, becomes more efficient due to the highly curved shape of the flux tube. The draining makes the flanks significantly heavier than the surrounding material. Thus, while the apex continues to emerge, the flanks start to submerge. The submergence modifies the geometrical shape of the emerging field further, making the apex of the flux tube even more curved and, therefore, further enhancing the plasma draining. Eventually (after $t=1400$ minutes), the middle part of the flux tube loses enough mass to become buoyant again, and, hence, continues to rise and to expand. This complicated process affects the overall horizontal and vertical expansion the flux tube, resulting in a reduced magnetic field strength. Thus, when the flux tube reaches the photosphere, it carries very high $\beta$ plasma. Furthermore, the compression rate of the field below the photosphere is very low. As a result, the field fails to emerge.

The biggest difference between $\lambda=20$ and higher $\lambda$ cases is found at the scaling curves (Fig. 5d). The $\lambda=20$ case does not scale with a power law. Notice also that the magnetic field strength, $B$ is significantly reduced as the field rises. Interestingly, the variation of $\lambda$ from 35 to 100 does not effect dramatically the scaling curves.

3.3.4 Variation of $\alpha$

We now focus on the dynamics of the emerging flux tubes when the twist is varied. We showed that small $\lambda$ affects significantly the plasma draining along the field lines. The twist is a parameter that affects the efficiency of the draining, as higher twisted field lines has ‘dips’that can trap dense plasma. We will study the effects of the variation of twist by using both large and small $\lambda$ to capture the effect of the twist on the draining along the field lines.

First, we select values of $\alpha$d$=0.1,\,0.15,\,0.25$ and $\lambda$d$=100$, $B_{0}$d$=34$, $R$d$=5$ (Table 2, cases 3, 11, 12). The larger, twist flux tubes emerge slightly faster (Fig. 8e). This is consistent with previous studies (e.g. Murray et al., 2006). Notice also that their scaling curves behave similarly, deeper in the convection zone (steeper slopes in Fig. 2e). Closer to the photosphere, the higher the twist, the smaller the value of $\kappa$ ($\kappa=0.23,\,0.13,\,0.11$ for $\alpha$d$=0.1,\,0.15,\,0.25$). This is expected as i) the radial magnetic tension from the twist keeps the flux tube more coherent, bringing stronger field below the photosphere and ii) higher twist flux tubes have stronger poloidal field component, which is further enhanced during the compression below the photosphere. So, overall, the flux tubes with higher twist emerge more efficiently.

Now, we select cases with a smaller $\lambda$d$=20$ and $\alpha$d$=0.1,\,0.11,\,0.15,\,0.25$, $B_{0}$d$=34$, $R$d$=5$ (Table 2, cases 8, 13, 14, 15). In the low $\lambda$ cases we find some unexpected results.

The cases with $\alpha$d$=0.1,\,0.11,\,0.15$, at $t\approx 1000$ minutes, stop rising for a small time period (Fig. 8f). Then they start rising again until they become decelerated by the photosphere. This is similar to the $\lambda=20$ case of Sec. 3.3.3. The net effect of this motion is enhanced plasma draining, leading to a complicated horizontal and vertical expansion. However, the $\alpha$d$=0.25$ flux tube behaves differently (green lines). There, the higher twist prohibits the enhanced draining that occurs in the lower $\alpha$d cases. This flux tube emerges without the complicated horizontal and vertical expansion that is present in the lower twist cases. As a result, its internal magnetic pressure is reduced less during its emergence. However, the high downwards magnetic tension and the lack of draining eventually reduces the rate of emergence of the $\alpha$d$=0.25$ case ($t=1200-2500$ min, green line).

For the cases with the enhanced draining ($\alpha$d$=0.1,\,0.11,\,0.15$), increasing $\alpha$ led to more efficient emergence (Fig. 5f). We do not find power laws deep in the convection zone for these cases. A less steep linear power law appears only for $\alpha$d$=0.11,\,0.15$, when they compress below the photosphere. However, for $\alpha$d$=0.25$, due to the deceleration of the flux tube, the compression below the photosphere is less. Thus, less steep part of the scaling curve has higher $\kappa$ value than the values for the less twisted cases. Therefore, we find that for $\lambda$d$=20$, the higher twisted flux tube emerges less efficiently than the less twisted cases.

3.3.5 All cases

We now plot all the scaling curves (Table 2, cases 1-15) in Fig. 9a. We also plot some additional cases that mostly describe “failed” emergence (cases 16-21). The blue lines are the cases that “successfully” emerge above the photosphere (non-asterisk cases in Table 2) and the red lines are the cases that “fail” to emerge above the photosphere (asterisk cases). Notice that there is a clear separation and clustering of the blue and the red lines. The green line is the scaling curve of the “failed” emergence of case 8 (discussed in Sec. 3.3.3), which acts as a “border-line” between the bulk of the “successfully” emerged cases and the ones that “failed” to emerge.

In Fig. 9b, we plot again the “border-line” case (green line). We color the region above that line with blue and below with red. We noticed that if the left-most part of flux tube’s scaling curve is located inside the blue region, then parts of the flux tube will eventually emerge above the photosphere (case 20, blue line, Fig. 9b). If it ends inside the red region, it will eventually fail to emerge (case 17, red line). Using the above comments and how $\kappa$ behaves when varying $B_{0}$, $R$, $\alpha$, we were able to estimate flux tube parameters needed for a “successful” or “failed” emergence, and roughly estimate the value of the magnetic field below the photosphere.

An interesting result is the following. In Fig. 9a, most of the blue lines originate from the same point, as they initially have a $B_{0}=10$ kG field. However, the photospheric field is much different. Therefore, the efficiency of the emergence (ratio of maximum photospheric field strength over $B_{0}$) is different. For instance, case 15 has an efficiency of 0.02, while case 7 has an efficiency of 0.1. The rest of the “successful’ emergence cases starting with $B_{0}=10$ kG have intermediate values of efficiency. This difference in the efficiency, is due to effects of the geometry of the field (twist, radius, curvature) on the emergence. Fig. 9b is another example of the effect of the geometry on the emergence. Case 17 (red line) is a case with $B_{0}$d$=44$ ($B_{0}=13200$ G), that fails to emerge because of its very small $\lambda$d$=10$. However, case 20 (blue line), which has similar flux to case 17, but almost half the magnetic field strength ($B_{0}$d$=24$ ($B_{0}=7200$ G)), emerges “successfully” due to the larger radius and the larger $\lambda$.

The physical meaning of the above is that, in order for a flux tube to emerge above the photosphere, it must bring with it the necessary amount of magnetic field strength and flux. If its geometry and twist are not favouring the efficient emergence of this field, then even an initially strong field will fail to emerge. On the other hand, weaker fields can emerge above the photosphere if their geometry results to a more efficient emergence. Our “border-line” case is a numerically derived limit that separates the two states.

4 Summary and Discussion

In this work we studied the emergence of flux tubes from 18 Mm below the photosphere, using 3D MHD numerical simulations. We performed a detailed parametric study on: (i) the magnetic field strength; (ii) the twist; (iii) the radius and (iv) the length of the buoyant part of a flux tube. Initially, we varied the radius and the magnetic field strength (while keeping the twist and the length of the buoyant part constant), to study whether the initial amount of sub-photospheric magnetic flux is a good indicator for “successful” emergence (Fig. 3). Then, we focused on the scaling of the maximum magnetic field strength with the local density. We identified the curve that describes the maximum $B$ as a function of $\rho^{\kappa}$ (scaling curve). The scaling curve had a part with steeper slope (larger $\kappa$, where $\kappa$ is the power of the density such that $B\propto\rho^{\kappa}$) and this developed in the deeper part of the solar convection zone. Close to the photosphere, $B$ scales with $\rho^{\kappa}$ with a smaller $\kappa$. However, in a few cases, the curves did not follow such a power laws. We identified under which conditions the scaling curve can form a power law, and derived expressions for $\kappa$ that describe approximately the scaling. Finally, we studied the scaling curves and the height-time profiles for a number of different initial conditions (Table 2) by keeping constant three of the $B_{0},\,R,\,\alpha,\,\lambda$ and varying the remaining variable (Fig. 5, 8).

Our results are summarized as follows:

Magnetic flux alone is not sufficient to estimate whether the magnetic field will emerge, especially below $10^{21}$ Mx. 2. 2.

$B$ scales as $\rho^{\kappa}$ when the magnetic field has one dominant direction (the apex of the emerging flux tube is locally horizontal along a large enough segment) and the spatial/temporal changes of the velocity gradients and shear are not significant. In its most general form, a constant $\kappa$ can be described by Eq. A14. 3. 3.

The steeper part of the scaling curves develops when the flux tube apex is horizontal-like and is located deeper in the solar interior (similar to Fig. 6a). The less steep part of the scaling curves develops due to the compression of the flux tube just below the photosphere (similar to Fig. 6b). The transition from the less to the more steep part of the scaling curve occurs approximately when the characteristic radial size of the emerging tube is similar to the local pressure scale height ($2R\approx H_{p}$ in our case). Some parameters (like twist) can affect this characteristic length as they affect the rate of the flux tube expansion. For flux tubes whose apex is not horizontal-like, the field strength does not scale with the local density deeper in the solar interior (similar to Fig. 6c). However, a power law can be developed below the photosphere if such a flux tube compresses significantly. 4. 4.

The magnetic field is transferred upwards more efficiently when $B_{0}$ or $R$ is increased. In most cases, this applies also for the twist. 5. 5.

A highly curved flux tube (small $\lambda$) with low twist emerges less efficiently in comparison to a lower curvature flux tube (large $\lambda$) with similar twist. 6. 6.

In a highly curved flux tube, increasing the twist increases the efficiency of the emergence to a certain extent. Eventually, the higher twist obstructs the plasma draining by maintaining a local dip in the magnetic field, the flux tube remains heavy and the efficiency of emergence is reduced. As a result, a higher twisted flux tube can eventually bring less magnetic field closer to the photosphere in comparison to a less twisted one. 7. 7.

The combined effect of all the above (Fig. 9a) shows that the efficiency with which the magnetic field is brought upwards is a significant aspect of the emergence of buoyant magnetic fields in the solar interior For instance, high-$B_{0}$ (weak-$B_{0}$) fields may fail (succeed) to emerge to the photosphere, depending on their geometrical properties.

Based on our results, it is clear that there is neither a specific $\kappa$ for which $B\propto\rho^{\kappa}$ everywhere in the solar interior nor a specific $\kappa$ that describes the field close to the photosphere.

Deep in the solar interior, Pinto & Brun (2013) found in their dynamo simulation that $\kappa\approx 1$. They showed that the poloidal expansion dominated over the axial expansion. This is in agreement with our analysis. Assuming a strong axial field oriented along the $x$-axis, the field would scale with $\kappa_{x}$. Then, for negligible axial expansion ($\frac{\partial v_{x}}{\partial x}\approx 0$), from Eq. A2a we derive that $\kappa\approx 1$. If the axial expansion is not negligible in comparison to the poloidal one, the value of $\kappa$ can be different.

Cheung et al. (2010) studied $\kappa$ in the case of the emergence of a highly twisted toroidal flux tube inside a convective layer and found $\kappa=0.5$. We consistently find lower values than that, meaning that in our simulations the magnetic field is transferred more efficiently upwards. It is possible that this is due to the lack of a fully developed convective envelope in our or their simulations. Convective motions should deform to an extent the flux tubes, and reduce the efficiency of emergence. Thus, the effect that convective motions have on emerging flux tubes is very important for the study of the scaling of $B$ with $\rho$.

Such effects cannot be easily estimated. However, the comparison between the buoyancy force and the drag force has been proposed as a measure for identifying whether convective motions will have a destructive effect on a flux tube or not. Moreno-Insertis (1983); Fan et al. (2003); Cheung et al. (2007) showed that the flux tube will not be fragmented by the convective motions if its magnetic field strength is:

[TABLE]

where $B_{eq}$ is the equipartition value of the magnetic field strength with the local kinetic energy density ( $B_{eq}=\sqrt{\mu\rho u_{downflow}}$, where $u_{downflow}$ is the local velocity of downdrafts). To estimate $B_{eq}$, we need the velocities of the local vertical flows.

In helioseismology, vertical velocities are calculated by averaging data across large regions (e.g. Komm et al., 2004, 2011). Therefore, the local fast upflows and downflows are smoothed and such vertical velocities estimate the mean value across these large regions. Also, comparisons between models and helioseismology methods have posed questions about the accuracy of vertical velocity measurements below certain depths (Zhao et al., 2010). As a result, we cannot use vertical velocities from helioseismology to estimate $B_{eq}$ in Eq. 18. To estimate $B_{eq}$, we could assume some values for the vertical velocities. For instance, we can assume that the vertical velocities are of the order of the horizontal velocities derived from helioseismology (e.g. Greer et al., 2015). Another approach could be to use the vertical velocities at different depths given from models (Stein et al., 2011). Using either the vertical velocity root mean square from Stein et al. (e.g. 2011), or the horizontal velocity root mean square from Greer et al. (2015) we find that our selected $B_{0}$ values satisfy Eq. 18. So, our flux tubes would not be fragmented from the downdrafts, at least deeper in the interior. However, as the flux tubes expand closer to the photosphere, we expect that the convective motions will deform these flux tubes, reducing the magnetic field strength. We would expect also to find the “opposite” of the deformation. In a fully developed convective layer weaker fields can intensify locally due to “convective intensification” (e.g. Parker, 1978; Spruit, 1979). The actual degree of deformation and intesification, their effect on $\kappa$ and whether they could significantly impact the emergence of magnetic elements above the photosphere, is unknown. To estimate these effects requires 3D compressive simulations with fully developed convection zones.

Note that, in this work, we do not aim to identify conditions where flux tubes will form an active region (of any size). Our aim is to study the scaling, and identify cases where the field emerges above the photosphere, even if the photosheric magnetic field strength is small.

In most 3D flux emergence simulations, the flux tube is initially located close to the photosphere, around $-5$ Mm to $-1$ Mm (e.g. Fan, 2001; Magara & Longcope, 2001; Manchester et al., 2004; Archontis et al., 2004; Murray et al., 2006; Hood et al., 2009; MacTaggart & Hood, 2009; Fan, 2009; Moreno-Insertis & Galsgaard, 2013; Toriumi & Yokoyama, 2013; Leake et al., 2013; Fang et al., 2014; Lee et al., 2015; Takasao et al., 2015; Syntelis et al., 2015, 2017). For instance, in the parametric study of Murray et al. (2006), the flux tube is placed at $-1.7$ Mm, whereas in our simulation the flux tube is placed at $-18$ Mm. We found that the previous results are consistent with the results of flux tubes placed deeper in the interior. However, our work shows that additional effects are also important during the emergence of flux tubes from deeper in the interior, associated mostly with the plasma draining along the field lines. Toriumi & Yokoyama (2013) performed 3D simulations of flux tubes placed at -20 Mm. They did not find the effects on the plasma draining that we identified in our simulations when varying $\lambda$ and the twist. This is probably because they did not explore the same parameter space of low $\lambda$ and twist. For flux tubes similar to theirs, our results are in agreement. On the other hand, they showed that for higher values of $\lambda$ (e.g. $\lambda$d$=400$) than the ones we used, the flux emerges slightly slower compared to lower $\lambda$ cases (e.g. $\lambda$d$=100$). They attributed this behaviour to very slow plasma draining. We do not find a similar behaviour in our simulations, but it is possible that the further increase of $\lambda$ could lead to similar results.

Notice that in Fig. 9a, the vast majority of the “successfully” emerged cases (blue lines) start with the same $B_{0}$, and differ in flux, twist and the $\lambda$. Just below the photosphere, though, the magnetic field strength ranges from $200-1000$ G. Therefore, the magnitude of the photospheric magnetic field does not contain sufficient information to infer the magnetic field strength of the initial flux tube. To estimate $B_{0}$, information about the radius and the shape of the flux tube are needed. This information, along with some estimate of the subphotospheric velocity vector, can assist in estimating the value of $\kappa$ close to the photosphere and the depth where the scaling curve changes behaviour. Hence, this could be used to estimate the magnetic field strength deeper in the interior. For such a calculation, further work is needed in many aspects. For instance, using 3D flux emergence models, it is important to identify whether the photospheric values of twist and the length scale of the emerged field can be correlated with the corresponding subphotospheric values. If no such relation exists (similar to our result for $B_{0}$ and photospheric $B$), then using the photospheric values of twist, $B$, and the size of an active region would provide little information about the conditions below the photosphere. To understand the nature of the magnetic fields below the photosphere, information about the sub-photospheric magnetic field strength and the sizes of the typical emerging structures is required. Such parameters are essential to further develop our understanding of solar flux emergence and to pose constrains on numerical models.

This project has received funding from the Science and Technology Facilities Council (UK) through the consolidated grant ST/N000609/1. The authors acknowledge support by the Royal Society. This work was supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facility - ARIS.

Appendix A Derivation of scaling laws

A.1 Velocity field without shearing terms

Assuming a velocity field with no shearing $(\frac{\partial v_{i}}{\partial x_{j}}=0,\quad i\neq j)$, the components of the ideal induction equation (Eq. 2) can be written as:

[TABLE]

where we define

[TABLE]

Combining Eq. A1a, A1b, A1c we get:

[TABLE]

To study the conditions under which the magnetic field strength will scale with a power of the local density, we assume that the magnetic field strength can be written as

[TABLE]

where $\kappa$ is constant and $B_{0}$, $\rho_{0}$ are the values of $B$ and $\rho$ at $t=0$. By solving the above for $\rho$, substituting that expression into the continuity equation (Eq. 3), and then multiplying by $2B$, Eq. 3 becomes:

[TABLE]

Eq. A3 and Eq. A5 are consistent only if:

[TABLE]

We will now identify the possible solutions of the above equation.

A.1.1 Case 1: $B_{x}\gg B_{y},B_{z}$

In this case, we assume that the magnetic field has one dominant direction, say along the $x$-axis. Then, the magnetic field can be described locally only by the $B_{x}$ component of the full magnetic field vector. Hence, Eq. A6 suggests that the magnetic field strength will indeed scale with the local density and that:

[TABLE]

Since we have assumed $\kappa$ to be constant, $\kappa_{x}$ needs to be constant too and the velocity field is constrained such that (see Eq. A2a):

[TABLE]

where $\chi,\,\psi,\,\zeta$ are constants. The constant velocity gradients guarantee that the magnetic field will have the same direction as the initial field.

If $B_{x}$ is the dominant magnetic field component, then the field will scale with $\kappa_{x}$. If $B_{y}$ or $B_{z}$ is the dominant magnetic field component, the field will scale with $\kappa_{y}$ or $\kappa_{z}$ respectively.

A.1.2 Case 2: $B_{x}\sim B_{y}\gg B_{z}$

Now we assume that the magnetic field has two dominant directions (e.g. along the $x$-axis and $y$-axis), then the magnetic field can be described locally by two components of the full magnetic field vector. Then, from Eq. A6, we get that $\kappa=\kappa_{x}=\kappa_{y}$, since the magnitude of the magnetic field to be $B^{2}=B_{x}^{2}+B_{y}^{2}$. Adding them together, we get:

[TABLE]

In order for $\kappa=\kappa_{x}=\kappa_{y}$ and $\kappa$ to be constant the velocity gradients are required to be:

[TABLE]

where $\chi,\,\zeta$ are constants.

Therefore, if the magnetic field can be described using two components of the full magnetic field vector, its strength will scale with the local density only when the field expands at the same rate along these two directions. If the important components are $B_{x}$ and $B_{y}$, then the field will scale as $\kappa=1-\frac{1}{2}\kappa_{z}$. For $B_{x}$ and $B_{z}$ ($B_{y}$ and $B_{z}$) the field will scale with $\kappa=1-\frac{1}{2}\kappa_{y}$ ($\kappa=1-\frac{1}{2}\kappa_{x}$).

In general, such a “2D” field will not scale with the local density. The velocity field restriction is such that it forces the magnetic field to maintain the direction of the total magnetic field vector. As a result, the restriction forces the “2D” magnetic field to behave as “1D” in the field aligned coordinate system. Thus, in the field aligned system, the magnetic field behaves according to Case 1.

A.1.3 Case 3: $B_{x}\sim B_{y}\sim B_{z}$

Now we assume that all the magnetic field components are needed to describe the magnetic field. Then, from Eq. A6 we get that $\kappa=\kappa_{x}=\kappa_{y}=\kappa_{z}$, in order for the magnitude of the magnetic field to be $B^{2}=B_{x}^{2}+B_{y}^{2}+B_{z}^{2}$. Adding these terms together gives that

[TABLE]

In order for $\kappa=\kappa_{x}=\kappa_{y}=\kappa_{z}$ and $\kappa$ to be constant the velocity gradients are required to be:

[TABLE]

where $\chi$ is constant.

Therefore, a general magnetic field can scale with the local density with a constant $\kappa$ only if it expands isotropically. In general, a “3D” field will not scale with the local density. As in Case 2, the velocity field restriction is such that the magnetic field maintains the direction of the total magnetic field vector. Therefore, this restriction makes the magnetic field in the field aligned coordinate system to behave as “1D”.

A.2 Velocity field with shearing terms

We now include the shearing terms of the velocity field. Following the same steps as before, we get from the induction equation that:

[TABLE]

Combining with Eq. A5, we get the generalized expression of Eq. A6:

[TABLE]

The above can be written compactly as:

[TABLE]

where $i,j$ are indices corresponding to the $x,y,z$ coordinates, and $\kappa_{ij}$ is

[TABLE]

where $E_{ij}=\frac{1}{2}\left(\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}\right)$ is the strain rate tensor. As $E_{ij}$, $\kappa_{ij}$ is symmetric. Its diagonal elements can also be expressed in terms of the strain rate tensor, so that $\kappa_{ij}=I-\frac{1}{\mathbf{\nabla}\cdot\mathbf{v}}\,E_{ij}$, where $I$ is the identity matrix. Therefore, $\kappa_{ij}$ is, as $E_{ij}$, a metric of the deformation of the velocity field. Eq. A7, A9, A11 are special cases of Eq. A14.

The only assumption made to derive Eq. A14 was that $\kappa$ is constant (in Eq. A5). In order for $\kappa$ to be constant, all terms of $\kappa_{ij}$ and $B_{i}$ have to be independent of position and time.

A.3 Non-constant $\kappa$

We now assume that $\kappa$ is not constant, but a general function of $x$, $y$, $z$ and $t$ (i.e. $\kappa(x,y,z,t)$). Substituting Eq. A4 in Eq. 3 and multiplying by $2B$ we get:

[TABLE]

This is the generalization of Eq. A5. We now write Eq. A12 compactly as $\frac{DB^{2}}{Dt}=-2\kappa_{ij}B_{i}B_{j}\mathbf{\nabla}\cdot\mathbf{v}$ and substitute it in the previous equation. We get that:

[TABLE]

Notice that in order for $\kappa$ to be constant, the term on the right hand side needs to be zero, which gives the previous result for constant $\kappa$ (Eq. A14). If $\kappa$ is changing slowly, so that $\frac{D\kappa}{Dt}\approx 0$, then $\kappa\approx\frac{1}{B^{2}}\kappa_{ij}B_{i}B_{j}$. Therefore, Eq. A14 can describe both the scaling of the magnetic field strength with the local density when $\kappa$ is constant, but also when $\kappa$ is changing slowly in time and deviating slightly from a power law.

Appendix B Resolution, resistivity and viscosity effects

We examine the effect of the resolution on the scaling curves. To do so, we choose case 4 Table 2 to be our reference simulation, as this case is examined in detail in Sec. 3.2.1 when explaining the formation the scaling laws. We perform simulations with the same initial conditions and physical domain, and vary the number of grid points. In Fig. 10a, we plot the scaling curve of the reference simulation (solid line, $600^{3}$ grid points) and compare it with simulations of lower (dotted-dashed line is $400^{3}$ grid pints and dotted line $500^{3}$ grid points) and higher (dashed line, $700^{3}$ grid points) resolution. We find that the scaling curve is minorly affected by the change of the resolution.

We also examine the effect of the viscosity on the scaling curves, by performing a simulation with $\nu_{3}$=0. We do not set the shock viscosity coefficients ($\nu_{1}$ and $\nu_{2}$) to zero to ensure the numerical stability of the simulation. In Fig. 10b, we plot the reference simulation (solid line) and the $\nu_{3}$=0 simulation (dashed line). The two curved mostly overlap.

Finally, we examine effect of the resistivity on the scaling curves by performing a simulation with no explicit resistivity ($\eta=0$). We plot its scaling curve in Fig. 10b (dotted line). We find that the scaling curve of this simulation is different from the reference simulation by some degree.

The effect of the resistivity on the scaling curves can be estimated analytically by extending the analysis of Appendix A. Instead of the ideal induction equation, we use the non-ideal induction equation with uniform resistivity:

[TABLE]

Using the above in our analysis, Eq. A12 becomes:

[TABLE]

Combining the above with Eq. A5 gives:

[TABLE]

So, resistivity will have an effect on the value of $\kappa$. The second term of the above equation for the reference simulation is of the order of $0.01-0.1$ during the simulation and this is approximately the difference we find between the solid and the dotted curves in Fig. 10b.