Stabilizing Effect of Magnetic Helicity on Magnetic Cavities in the Intergalactic Medium

To approach this problem and to provide a quantitative assessment we make use of direct numerical simulations, which we perform with the public code PencilCode.⁵ This code is particularly suitable for our study since it avoids using the magnetic field ${\boldsymbol{B}}$ as a primary variable and instead uses its vector potential ${\boldsymbol{A}}$, ensuring the fields stay solenoidal throughout the simulations. Moreover, it allows us to quantify the magnetic helicity $H=\int {\boldsymbol{A}}\cdot {\boldsymbol{B}}{\rm{d}}V$, where the integral is calculated over the whole computational domain. For our study we require an ideally conducting viscous medium. This is governed by the resistive magnetohydrodynamics equation together with the energy (temperature) equation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{A}}}{\partial t}={\boldsymbol{u}}\times {\boldsymbol{B}}+\eta {{\boldsymbol{\nabla }}}^{2}{\boldsymbol{A}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{D}}{\boldsymbol{u}}}{{\rm{D}}t} & = & -{c}_{{\rm{s}}}^{2}{\boldsymbol{\nabla }}\left(\displaystyle \frac{\mathrm{ln}T}{\gamma }+\mathrm{ln}\rho \right)+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{\rho }\\ & & -{\boldsymbol{g}}+{{\boldsymbol{F}}}_{\mathrm{visc}},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{Dln}\rho }{{\rm{D}}t}=-{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \mathrm{ln}T}{\partial t} & = & -{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}\mathrm{ln}T-(\gamma -1){\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}\\ & & +\displaystyle \frac{1}{\rho {c}_{v}T}\left({\boldsymbol{\nabla }}\cdot (K{\boldsymbol{\nabla }}T)+\eta {{\boldsymbol{J}}}^{2}\right.\\ & & \left.+2\rho \nu {\bf{S}}\otimes {\bf{S}}+\zeta \rho {\left({\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}\right)}^{2}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

with the magnetic vector potential ${\boldsymbol{A}}$, magnetic field ${\boldsymbol{B}}={\boldsymbol{\nabla }}\times {\boldsymbol{A}}$, fluid velocity ${\boldsymbol{u}}$, constant magnetic resistivity (diffusivity) η, advective derivative ${\rm{D}}/{\rm{D}}t=\partial /\partial t+{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}$, sound speed ${c}_{{\rm{s}}}\,=\gamma p/\rho$, adiabatic index $\gamma ={c}_{p}/{c}_{v}$, heat capacities ${c}_{p}$ and ${c}_{v}$ at constant pressure and volume, temperature T, density ρ, electric current density ${\boldsymbol{J}}={\boldsymbol{\nabla }}\times {\boldsymbol{B}}$, gravitational acceleration ${\boldsymbol{g}}$, viscous force ${{\boldsymbol{F}}}_{\mathrm{visc}}$, heat conductivity K, and bulk viscosity ζ. The viscous force is given as ${{\boldsymbol{F}}}_{\mathrm{visc}}={\rho }^{-1}{\boldsymbol{\nabla }}\cdot 2\nu \rho {\bf{S}}$, with the traceless rate of strain tensor ${S}_{{ij}}=\displaystyle \frac{1}{2}\left({u}_{i,j}+{u}_{j,i}\right)-\tfrac{1}{3}{\delta }_{{ij}}{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}$. The equation of state used here is for the ideal monoatomic gas and it appears implicitly in our equations, as we eliminated pressure p. Here the gas is monatomic with $\gamma =5/3$.

Our side boundaries (xy) are chosen to be periodic, while the bottom is closed (${\boldsymbol{u}}\cdot {\boldsymbol{n}}=0$) and the top open. This allows for outward fluxes. For the magnetic field the vertical boundaries are set to open, allowing magnetic flux.

As a simulation domain we choose a box of size $2.4\times 2.4$ in the horizontal (xy) plane and 9.6 in the vertical (z), using $480\times 480\times 1920$ meshpoints.

In order to reduce resistive magnetic helicity decay we choose a value of the magnetic resistivity as low as the resolution allows. Here we set it to $\eta =3\times {10}^{-4}$. Viscosity is either set to $\nu =1\times {10}^{-3}$ or to $\nu =2\times {10}^{-4}$ to check that the overall behavior of the simulation does not depend on this parameter, which prevents any accumulation of turbulent energy at small scales.

With the viscosity and magnetic resistivity we can then compute our Reynolds numbers:

$$\begin{eqnarray}&&\mathrm{Re}\,=\,\displaystyle \frac{{u}_{\max }d}{\nu },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\mathrm{Re}}_{{\rm{M}}}\,=\,\displaystyle \frac{{u}_{\max }d}{\eta },\end{eqnarray} \tag{ 6 }$$

with the maximum velocity ${u}_{\max }$ and bubble diameter d = 1.6.

As an initial condition we choose a stably stratified atmosphere in which we place an under-dense hot cavity of spherical shape. This can be found as part of the PencilCode under src/initial_condition/bubbles_init.f90. The stable atmosphere obeys the hydrostatic equilibrium

$$\begin{eqnarray}&&\rho g=-{\boldsymbol{\nabla }}p,\end{eqnarray} \tag{ 7 }$$

with the gravitational acceleration in the negative z-direction g. Our model atmosphere extends in the z-direction, which makes the gradient a derivative in z and we can rewrite Equation (7) as

$$\begin{eqnarray}&&\rho g=-\displaystyle \frac{{\rm{d}}p}{{\rm{d}}z}.\end{eqnarray} \tag{ 8 }$$

Density, pressure, and temperature are related through the ideal gas law

$$\begin{eqnarray}&&p=\displaystyle \frac{R}{\mu }\rho T,\end{eqnarray} \tag{ 9 }$$

where μ is the mass of one mol of gas and $R=8.31\,{\rm{J}}\,{{\rm{K}}}^{-1}\,{\mathrm{mol}}^{-1}$ is the ideal gas constant. With the ideal gas law we can express the hydrostatic equilibrium as

$$\begin{eqnarray}&&\displaystyle \frac{\mu g}{R}\displaystyle \frac{{\rm{d}}z}{T}=-\displaystyle \frac{{\rm{d}}p}{p}.\end{eqnarray} \tag{ 10 }$$

The gas is chosen to be adiabatic, i.e., 

$$\begin{eqnarray}&&{p}^{1-\gamma }{T}^{\gamma }=\mathrm{const}.,\end{eqnarray} \tag{ 11 }$$

which leads to the relation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}p}{p}=\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{{\rm{d}}T}{T}.\end{eqnarray} \tag{ 12 }$$

With Equation (10) we can write

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}T}{{\rm{d}}z}=-\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{\mu g}{R}.\end{eqnarray} \tag{ 13 }$$

We now integrate this equation and obtain

$$\begin{eqnarray}&&T=-\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{\mu g}{R}z+{T}_{0},\end{eqnarray} \tag{ 14 }$$

where T₀ is the temperature at an arbitrary height z₀.

We now choose this height to be the isothermal scale height

$$\begin{eqnarray}&&{z}_{0}=\displaystyle \frac{{{RT}}_{0}}{\mu g}.\end{eqnarray} \tag{ 15 }$$

This is justified if we assume that the gas is in isothermal equilibrium at z₀, which is a common assumption for the adiabatic atmosphere. Our temperature profile now obtains the form

$$\begin{eqnarray}&&T(z)={T}_{0}\left(1-\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{z}{{z}_{0}}\right).\end{eqnarray} \tag{ 16 }$$

With the temperature profile and Equation (12) we can compute the pressure profile to

$$\begin{eqnarray}&&p(z)={p}_{0}{\left(1-\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{z}{{z}_{0}}\right)}^{\gamma /(\gamma -1)}.\end{eqnarray} \tag{ 17 }$$

Taking Equation (8) we can also compute the density profile to

$$\begin{eqnarray}&&\rho (z)={\rho }_{0}{\left(1-\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{z}{{z}_{0}}\right)}^{1/(\gamma -1)}.\end{eqnarray} \tag{ 18 }$$

For the cavity we choose a radius of 0.8 and place it centrally in x and y and at 0.8 in z. Its initial peak temperature is ${T}_{\mathrm{cavity}}=4$ and density ${\rho }_{\mathrm{cavity}}=0.25$. In order to avoid sudden wave formation we apply (sharp) smoothing for ${T}_{\mathrm{cavity}}$ and ${\rho }_{\mathrm{cavity}}$ near the edges using a tanh profile that matches the surrounding values. With that the temperature and density are constant up to ca. 72% of the radius of the cavity.

This should be contrasted with the surrounding medium, which is stably stratified with gravitational acceleration of g = 0.1. For that we choose ${z}_{0}=4$, ${\rho }_{0}=1$, and ${T}_{0}=1$ such that the cavity is in approximate pressure balance with its surrounding medium and its expansion or compression is insignificant.

For the simulations that include a magnetic field inside the bubble we consider two kinds of initial magnetic conditions. In the first case we insert a helical magnetic field of the Arnold–Beltrami–Childress (ABC) flow type:

$$\begin{eqnarray}&&{\boldsymbol{A}}=f(r){A}_{0}\left(\begin{array}{c}\cos ((y-{y}_{\mathrm{cavity}})k)+\sin ((z-{z}_{\mathrm{cavity}})k)\\ \cos ((z-{z}_{\mathrm{cavity}})k)+\sin ((x-{x}_{\mathrm{cavity}})k)\\ \cos ((x-{x}_{\mathrm{cavity}})k)+\sin ((y-{y}_{\mathrm{cavity}})k)\end{array}\right),\end{eqnarray} \tag{ 19 }$$

with the function

$$\begin{eqnarray}&&f(r)=1-{\left(r/{r}_{{\rm{b}}}\right)}^{{n}_{\mathrm{smooth}}}\end{eqnarray} \tag{ 20 }$$

that makes sure that the magnetic field vanishes outside the bubble without any strong current sheets. Here $r\,=\sqrt{{(x-{x}_{\mathrm{cavity}})}^{2}+{(y-{y}_{\mathrm{cavity}})}^{2}+{(z-{z}_{\mathrm{cavity}})}^{2}}$ and ${r}_{{\rm{b}}}$ is the radius of the bubble. We use as a smoothing coefficient ${n}_{\mathrm{smooth}}=2$. For values of $r\gt {r}_{{\rm{b}}}$ we set the field to 0.

Taking the curl we obtain ${\boldsymbol{B}}=k{\boldsymbol{A}}$. The magnetic energy scales like ${A}_{0}^{2}{k}^{2}$, while the magnetic helicity scales like ${A}_{0}^{2}k$. By inversely scaling the amplitude A₀ and the inverse scale k we can change the magnetic helicity while keeping the magnetic energy fixed. For our low-helicity setup we choose ${A}_{0}=2.5\times {10}^{-2}$ and k = 20. For the high-helicity case we choose ${A}_{0}=0.1$ and k = 5. For simplicity we will occasionally write H = 1 and H = 4 for the two cases.

To make sure our results are not dependent on a specific initial geometry of the internal magnetic field, but only the topology described by the amount of magnetic helicity, we perform a second set of simulations that makes use of a spheromak magnetic field. For that we use the form given by Aly & Amari (2012). There the magnetic field is given in spherical coordinates $(r,\theta ,\phi )$, with the origin at the bubble's center, as

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{B}} & = & 2{A}_{0}\displaystyle \frac{g(\alpha r)}{{\left(\alpha r\right)}^{2}}\cos (\theta ){\hat{e}}_{r}\\ & & -{A}_{0}\displaystyle \frac{g^{\prime} (\alpha r)}{\alpha r}\sin (\theta ){\hat{e}}_{\theta }\\ & & +{A}_{0}\displaystyle \frac{g(\alpha r)}{\alpha r}\sin (\theta ){\hat{e}}_{\phi },\end{array}\end{eqnarray} \tag{ 21 }$$

with $\alpha =\tau /{r}_{{\rm{b}}}$. The function g is given as

$$\begin{eqnarray}&&g(t)=\displaystyle \frac{{t}^{2}}{{\tau }^{2}}-\displaystyle \frac{3}{\tau \sin (t)}\left(\displaystyle \frac{\sin (t)}{t}-\cos (t)\right)\end{eqnarray} \tag{ 22 }$$

and its derivative is

$$\begin{eqnarray}&&g^{\prime} (t)=\displaystyle \frac{2t}{{\tau }^{2}}-\displaystyle \frac{3}{\tau \sin (t)}\left(-\displaystyle \frac{\sin (t)}{{t}^{2}}+\displaystyle \frac{\cos (t)}{t}+\sin (t)\right),\end{eqnarray} \tag{ 23 }$$

with τ being the zeros of the equation $g(\tau )=0$ , which can be calculated numerically to ${\tau }_{1}=5.763$, ${\tau }_{2}=9.095$, ${\tau }_{3}=12.229$, ${\tau }_{4}=15.515$, and ${\tau }_{5}=18.689$ for the first five roots. We then perform the transformation into Cartesian coordinates using the relations

$$\begin{eqnarray}\begin{array}{rcl}{B}_{x} & = & {B}_{r}(x-{x}_{\mathrm{cavity}})/r\\ & & +{B}_{\theta }(x-{x}_{\mathrm{cavity}})(z-{z}_{\mathrm{cavity}})/({r}_{{xy}}r)\\ & & -{B}_{\phi }(y-{y}_{\mathrm{cavity}})/{r}_{{xy}},\end{array}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}\begin{array}{rcl}{B}_{y} & = & {B}_{r}(y-{y}_{\mathrm{cavity}})/{r}_{{\rm{b}}}\\ & & +{B}_{\theta }(y-{y}_{\mathrm{cavity}})(z-{z}_{\mathrm{cavity}})/({r}_{{xy}}r)\\ & & +{B}_{\phi }(x-{x}_{\mathrm{cavity}})/{r}_{{xy}},\end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{B}_{z}={B}_{r}(z-{z}_{\mathrm{cavity}})/{r}_{{\rm{b}}}-{B}_{\theta }{r}_{{xy}}/r\end{eqnarray} \tag{ 26 }$$

with ${r}_{{xy}}=\sqrt{{(x-{x}_{\mathrm{cavity}})}^{2}+{(y-{y}_{\mathrm{cavity}})}^{2}}$.

Aly & Amari (2012) derived this field for a cylindrically symmetric force balance between magnetic and pressure forces solving the Grad–Shafranov equation. However, we note that in our case the hydrostatic pressure is not constant in space. Nevertheless, since the pressure deviation from p₀ is of the order of 0.5% we assume it is constant in our bubble.

We restrict the field within the bubble by setting it to 0 outside. We can do this, as the field naturally and smoothly vanishes at $r={r}_{{\rm{b}}}$ without generating any strong surface currents. To obtain the magnetic vector potential to this field we solve the equation ${{\boldsymbol{\nabla }}}^{2}\,{\boldsymbol{A}}=-{\boldsymbol{\nabla }}\times {\boldsymbol{B}}$ in Fourier space, where we use the Coulomb gauge (${\boldsymbol{\nabla }}\cdot {\boldsymbol{A}}=0$) and the fact that the field is periodic, since it vanishes at the boundaries. A rendering of the magnetic field can be seen in Figure 1.

Zoom In Zoom Out Reset image size

The third magnetic case we take into consideration is that of an initial vertical magnetic field

$$\begin{eqnarray}&&{\boldsymbol{B}}={B}_{0}{{\boldsymbol{e}}}_{z},\end{eqnarray} \tag{ 27 }$$

where B₀ is the initial amplitude of the field in the entire computational domain.

For a fair comparison, and a quantitative evaluation of the role of magnetic fields in the three cases described in Sections 2.3, 2.4, and 2.5, we first perform calculations of a purely hydrodynamical case with no magnetic fields. We compare this with the two cases with an internal magnetic field within the cavities. The two ABC cases have almost the same magnetic energy of ${E}_{{\rm{m}}}=0.2182$, while the two spheromak cases have energies of ${E}_{{\rm{m}}}=0.264$ for the high-helicity and ${E}_{{\rm{m}}}=0.189$ for the low-helicity case. This difference is a result of adjustments of the parameters that are chosen to match the magnetic helicity content of the ABC case as closely as possible, while keeping the difference with the magnetic energy minimal. Here the high-helicity case has a magnetic helicity content that is four times larger. An overview of the parameters used in the simulations is shown in Table 1. With these initial conditions we can compute the plasma β where the magnetic field is strongest using

$$\begin{eqnarray}&&\beta =\min \left(\displaystyle \frac{2(R\rho T/\mu )}{{B}^{2}}\right).\end{eqnarray} \tag{ 28 }$$

We compute this to be $\beta \approx 0.6$ for the helical ABC cases, $\beta =0.44$ for the high-helicity spheromak, and $\beta =0.038$ for the low-helicity spheromak case. However, we have to note that this is the minimum value β attains in our model, since with the presence of magnetic null points within the domain, β is not constant and even diverges at some points. In the external-field cases we have instead $\beta =20$ (weak field) and $\beta =1.25$ (strong field).

Table 1.  Simulation Parameters

Model	${\boldsymbol{B}}({A}_{0})$	${H}_{{\rm{m}}}$	ν	η	k/τ	Re	Re${}_{{\rm{M}}}$
hydro	⋯	⋯	10−3	⋯	⋯	960	⋯
hydro2	⋯	⋯	2 × 10−4	⋯	⋯	4800	⋯
hel_l	0.025	1	10−3	$3\times {10}^{-4}$	20	1280	4200
hel_h	0.1	4	10−3	$3\times {10}^{-4}$	5	1280	4200
hel_l2	0.025	1	2 × 10−4	$3\times {10}^{-4}$	20	5600	3700
hel_h2	0.1	4	2 × 10−4	$3\times {10}^{-4}$	5	6400	4200
sph_l	6.39	1	2 × 10−4	$3\times {10}^{-4}$	21.85	7200	4800
sph_h	1.7	4	2 × 10−4	$3\times {10}^{-4}$	5.76	11000	7500
ex_low	0.2	0	10−3	$3\times {10}^{-4}$	⋯	320	1000
ex_high	0.8	0	10−3	$3\times {10}^{-4}$	⋯	320	1000

Note. The magnetic field intensity, measured by the parameter A₀ in Equation (19), the magnetic helicity H_m, the viscosity ν, the magnetic diffusivity η, the parameter k of the ABC flow in Equation (19) (parameter τ of the spheromak configuration in Equation (21)), the Reynolds number Re and the magnetic Reynolds number Re_M, as defined in Equations (5) and (6).

Download table as:  ASCIITypeset image

We choose the conversion rate from code units to physical units such that the dimensions of the setups correspond to physically observed numbers in the intergalactic medium (see Table 2).

Table 2.  Code Unit Conversion Table

One Code Unit of	Physical Unit
length	10 kpc
time	10 Myr
density	$\times {10}^{-25}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$
temperature	$\times {10}^{6}\,{\rm{K}}$
magnetic field	$\times {10}^{-4}\,{\rm{G}}$

Download table as:  ASCIITypeset image

With these code unit conversions our setup has a size of ${L}_{{xy}}\,=24\,\mathrm{kpc}$ horizontally and ${L}_{z}=96\,\mathrm{kpc}$ vertically. The cavity has a radius of ${r}_{{\rm{b}}}=8\,\mathrm{kpc}$ with a density of ${\rho }_{{\rm{b}}}=2.5\times {10}^{-26}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and temperature of ${T}_{{\rm{b}}}=4\times {10}^{6}\,{\rm{K}}$. The surrounding medium at z = 0 has a density of ${\rho }_{0}=1\times {10}^{-25}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and temperature of ${T}_{0}=1\times {10}^{6}\,{\rm{K}}$. The system experiences a gravitational acceleration of $g=3.0985\times {10}^{-7}\,\mathrm{cm}\,{{\rm{s}}}^{-2}$. Our simulations then run between $200\,\mathrm{Myr}$ and $250\,\mathrm{Myr}$. The amplitude of the magnetic field is between ${B}_{0}=2.5\times {10}^{-6}\,{\rm{G}}$ and ${B}_{0}\,=6.39\times {10}^{-4}\,{\rm{G}}$. For the viscosity we obtain $\nu =3.0172\,\times {10}^{27}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ and resistivity $\eta =9.0516\times {10}^{26}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$.

The first case we take into account is the purely hydrodynamical one, to which we will compare the simulations including magnetic fields. As the cavity rises the Kelvin–Helmholtz instability acts on its interface with the surrounding medium and eventually affects the entire cavity by non-linear Kelvin–Helmholtz instability-induced turbulent mixing, as illustrated in the left panels of Figures 2, 3, and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We then perform two different simulations including an ABC helical magnetic field (see Section 2.3) inside the cavity. There it becomes evident that magnetic helicity contributes to keeping the cavity in a significantly more coherent state and prevents its disruption, provided a minimum amount of helicity is present in the field. This is evident from the central and left panels of Figure 2, showing two helical cases: the amount of magnetic helicity in the simulation shown in the right panel is four times larger than that in the central panel.

In order to make it easier to compare our results with observations we create artificial emission measures, assuming an optically thin medium. We place the observer along the y axis at an infinite distance. The emission measure E is then simply the line integral of the temperature to the fourth power,

$$\begin{eqnarray}&&E(x,z)=\int {T}^{4}\ {\rm{d}}y.\end{eqnarray} \tag{ 29 }$$

Similar to the slices plotted in Figure 2, we observe from the emission measures (Figure 3) an increased stability when a helical magnetic field is present.

We can also observe this behavior in the volume rendering of the temperature (Figure 4). While the purely hydrodynamical case (left panel) and the low-helicity case (center) disintegrate after two bubble diameter crossings, the strongly helical case (right) remains largely intact. Furthermore, in the purely hydrodynamical case, the temperature remains symmetric about the central axis of the domain, as expected from symmetry properties of our configuration. This is in agreement with simulations of Ruszkowski et al. (2007), as well as Dong & Stone (2009), who both produced perfectly symmetric configurations in the hydrodynamical case. Conversely, the disruption of the cavity in the low-helicity case develops a very asymmetric, chaotic structure. Also in the high-helicity case, although the disruption is only marginal, we can still observe an asymmetric evolution of the cavity.

To test the stability of the cavities we measure their coherence. In order to do so, we start by defining the space filled by the cavity as the loci for which ${\mathrm{log}}_{10}(T)\gt 1.5$. We choose this threshold because the surrounding cold medium has a significantly lower temperature and in the simulated times we do not observe a high enough temperature diffusion or conduction that would reduce the cavity temperature below this value. We then measure the mean distance ${d}_{\mathrm{mean}}$ of all the points in the cavity

$$\begin{eqnarray}&&{d}_{\mathrm{mean}}=\langle \left|{{\boldsymbol{r}}}_{\mathrm{cavity}}-{{\boldsymbol{r}}}_{\mathrm{CM}}\right|\rangle ,\end{eqnarray} \tag{ 30 }$$

where ${{\boldsymbol{r}}}_{\mathrm{CM}}$ is the position vector to the center of mass defined as

$$\begin{eqnarray}&&{{\boldsymbol{r}}}_{\mathrm{CM}}=\displaystyle \frac{{\int }_{{\mathrm{log}}_{10}(T)\gt 1.5}T{\boldsymbol{r}}{\rm{d}}V}{{\int }_{{\mathrm{log}}_{10}(T)\gt 1.5}T{\rm{d}}V},\end{eqnarray} \tag{ 31 }$$

where ${{\boldsymbol{r}}}_{\mathrm{cavity}}$ is the position vector to a point within the cavity, i.e., ${\mathrm{log}}_{10}(T)\gt 1.5$, and we take the average over the entire domain. We can therefore study the evolution of ${d}_{\mathrm{mean}}$ with time. Since some of the bubbles rise at different speeds, due to the different magnetohydrodynamic parameters of different models, we can also study the behavior of ${d}_{\mathrm{mean}}$ as a function of the height reached by the bubble, i.e., by its center of mass. For that we also compute the mean height of the bubbles as

$$\begin{eqnarray}&&{z}_{\mathrm{mean}}=\langle \left|{z}_{\mathrm{cavity}}-{z}_{\mathrm{CM}}\right|\rangle ,\end{eqnarray} \tag{ 32 }$$

where we use the loci for which ${\mathrm{log}}_{10}(T)\gt 1.5$, similar to our process for ${d}_{\mathrm{mean}}$.

Figure 5 (upper panel) depicts ${d}_{\mathrm{mean}}$ as function of ${z}_{\mathrm{mean}}$ in our simulations. We observe that without a magnetic field the cavity is dispersed without extending far from its starting position. An internal magnetic field with low helicity content (model hel_l in Table 1) does not improve significantly the stability for almost the whole simulation, but it has a stabilizing effect toward the end. However, an internal field with higher magnetic helicity (model hel_h in Table 1), and the same magnetic energy has some clear positive effect on the bubble's stability. In this latter case we observe a relatively stable rise throughout. Figure 5 (lower panel) depicts instead the time evolution of ${d}_{\mathrm{mean}}$. Here too it is evident how the strong-helical case is the most stable, while the low helical one only marginally stabilizes the cavity. From Figure 5 we see that high helicity stabilizes the cavity for at least 250 Myr, while the hydrodynamical case and the low-helicity case disrupt after about 100 Myr.

Zoom In Zoom Out Reset image size

We test the effect of a higher fluid Reynolds number on our results by performing simulations with a lower viscosity, $\nu =2\times {10}^{-4}$, for the hydrodynamical and the helical cases (models hydro2, hel_l2, and hel_h2 in table Table 1). There we observe a clearer onset of the Kelvin–Helmholtz instability. However, the coherence measure does not change significantly for any of the models compared to the low Reynolds number case, and the relative behavior between the hydrodynamical case and the two cases with internal helical magnetic fields resembles that obtained with higher viscosity (Figure 6). Therefore, this confirms our results in the lower Reynolds number regime.

Zoom In Zoom Out Reset image size

In order to eliminate effects from the particular geometry of the magnetic field, we test the stabilizing properties of the spheromak magnetic field configuration from Equation (21). To obtain a fair comparison with the ABC field case, we want to analyze a scenario with the same amount of magnetic energy and magnetic helicity. To do so, we choose appropriate values for the field amplitude and parameter τ, which in this case plays the same role as the parameter k for the ABC field. Since we have two degrees of freedom we can match both the magnetic energy and the helicity contents, such that our low (high) helicity spheromak configuration has almost the same magnetic energy and helicity as the low (high) helicity ABC configuration (see sph_l and sph_h in Table 1).

From the parameters in Table 1 we can see that our spheromak simulations belong to a high $\mathrm{Re}$ and ${\mathrm{Re}}_{{\rm{M}}}$ scenario. Since the parameters ν and η are the same used in the high Reynolds number cases for the ABC field, we can deduce that higher values of ${u}_{\max }$ are reached. Due to these localized high velocities of the gas that we attribute to local currents, we are able to run the high-helicity simulation only to a simulation time of ca. 9. Nevertheless, this gives us enough data to confirm the results obtained through the ABC configuration. That is, a high magnetic helicity internal magnetic field can stabilize the bubbles.

We observe this stabilizing effect from the slices of the temperature (Figure 7) and the emission measure (Figure 8). The high-helicity spheromak case clearly stays more stable. From the coherence calculations depicted in Figure 9, where we compare the hydrodynamical case with the low- and high-helicity spheromak cases, this stabilizing effect is evident as well, which is clearly visible in the high-helicity configuration. Conversely, the low-helicity configuration is characterized by a coherence similar to that of the hydrodynamic case. However, from the lower panel of Figure 9 we notice how in the low-helicity case the bubble can reach higher parts of the domain.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size