Relativistic magnetic explosions

Many explosive astrophysical events, like magnetars' bursts and flares, are magnetically driven. We consider dynamics of such magnetic explosions - relativistic expansion of highly magnetized and highly magnetically over-pressurized clouds. The corresponding dynamics is qualitatively different from fluid explosions due to the topological constraint of the conservation of the magnetic flux. Using analytical, relativistic MHD as well as force-free calculations, we find that the creation of a relativistically expanding, causally disconnected flow obeys a threshold condition: it requires sufficiently high initial over-pressure and sufficiently quick decrease of the pressure in the external medium (the pre-explosion wind). In the subcritical case the magnetic cloud just"puffs-up"and quietly expands with the pre-flare wind. We also find a compact analytical solution to the Prendergast's problem - expansion of force-free plasma into vacuum.


Introduction
Magnetically powered outflows form the basis of our understanding of pulsars (Goldreich & Julian 1969;Michel 1973;Masada et al. 2010) and Active Galactic Nuclei (AGNe, Blandford & Znajek 1977;Komissarov et al. 2007). In these cases the flows are typically approximated as (quasi)-stationary. Conservation of the magnetic flux then requires that the toroidal magnetic field scales as B φ ∝ 1/r while the poloidal B p ∝ 1/r 2 .
Dynamics of magnetically-driven explosions is qualitatively different from fluid explosions, quasi-stationary MHD flows, and from the 1D magnetized models. The two key ingredients are: (i) non-stationarity; (ii) the requirement of conservation of the magnetic flux. They change completely the overall dynamics of magnetic explosions, if compared with magnetized steady flows, fluid explosions and/or 1D models. Describing dynamics of multi-dimensional magnetic explosions is the main goal of the present paper.
First, the acceleration properties are different between steady-state and explosive events. For example, in steady-state pulsar winds the acceleration stops when the wind four-velocity equals approximately the fast magneto-sonic four-velocity (Michel 1969;Goldreich & Julian 1970;Heyvaerts & Norman 2003;Lyubarsky 2009). At this point the flow remains magnetically dominated, hence most of the initial magnetic power is not spent on acceleration. If the initial magnetization was σ 0 ≡ B 2 /(4πρc 2 ) 1, the flow reaches terminal Lorentz factor Γ w ∼ σ 1/3 0 , while the final magnetization of the flow has σ w ∼ σ 2/3 0 -magnetization decreases during acceleration, but remains high in the coasting stage.
Non-steady flows behave differently. Lyutikov (2010) found a fully analytic one-dimensional solution, a simple wave, for expansion of a highly magnetized plasma into vacuum and/or low density medium. The result is somewhat surprising: initially the plasma accelerates as Γ ∝ t 1/3 and can reach (locally, near the head part of the flow) terminal Lorentz factor Γ f = 1 + 2σ 0 . (see also Levinson 2010). Thus, it was shown, that time-dependent explosions can achieve much larger Lorentz factors that the steady state flows: 1 + 2σ 0 for non-steady versus σ 1/3 0 for steady-state. (The cold MHD solution were generalized to include pressure effects using a mathematically tricky hodograph transformation by Lyutikov & Hadden 2012).
The second important point, missed by 1D magnetized models, is the requirement of conservation of the magnetic flux. Consider a magnetized ball ejected during a burst, Fig. 1. Fields are tangled, and there is no clear separation between poloidal and toroidal components: they are linked and must evolve similarly. The flux conservations then requires that field scales as B ∝ 1/R(t) 2 , where R(t) is the overall size of an expanding cloud.
The scaling B ∝ 1/R(t) 2 works for more structured solutions as well, see below, when there is a clear separation between toroidal and poloidal fields. There is a number of analytical solution that confirm the B ∝ 1/R(t) 2 scaling: the Prendergast solution (Prendergast 2005, see also §4), and the expanding spheromak solution (Lyutikov & Gourgouliatos 2011).
Thus, multi-dimensional magnetic explosions are qualitatively different from both hydrodynamics explosions and 1D MHD flows: the requirement of conservation of the magnetic flux dominates the overall dynamics of magnetic explosions. This is missed in conventional hydrodynamic models of the ejections as expanding shell with parametrically added magnetic field, as well as onedimensional models of magnetic disturbances. Qualitatively, this is the same issue as the "sigma problem" in Pulsar Wind Nebulae: as Kennel & Coroniti (1984) argued, the presence of the global magnetic field changes the dynamics completely if compared with the fluid case. Blandford (2002) reformulated the problem in terms of the magnetic flux: the σ-problem is the problem of the (properly defined) toroidal magnetic flux conservation. Below we apply those ideas to the launching region of magnetic explosions.
To illustrate the effect of requirement of magnetic flux conservation consider a Solar flare-type event near the surface of the neutron star that produces a magnetically disconnected magnetized cloud with complicated (linked) internal magnetic structure and total internal pressure (magnetic field and pairs), slightly exceeding the local dipole field (Lyutikov 2022). Let the initial energy associated with the blob be E 0 ∼ B 2 N S R 3 0 (where R 0 ≤ R N S is the initial size of the blob; we assume that that the blob originates near the surface). As the magnetic blob expands its pressure quickly becomes larger than of the surrounding dipolar one. Expansion quickly becomes relativistic (within the magnetosphere). Still, the energy within the blob decreases with increasing radius.  (Lyutikov 2022). A magnetic bubble with tangled magnetic field is created near the neutron star; it expands within the magnetosphere and enters the pre-explosion wind when its size reaches the size of the light cylinder. At this point the blob is over-pressurized with respect to the preceding wind. We neglect the asymmetry due to the linear motion of the blob, and the polar-angle dependence of the preceding wind.
By the time the blob expands to the size of the light cylinder it's energy is ∼ E 0 (R 0 /R LC ). This energy should be compared with the energy of the dipolar field, as measured at the light cylinder, by the time the blob expanded to the size of the light cylinder it's energy is still larger than the dipolar magnetic energy at the light cylinder.
What is the ensuing fate of the explosion? As we discuss in the present paper two outcomes are possible. First, as the blob expands into the wind, magnetic field in the blob scales as B b ∝ 1/R 2 , decreasing much faster than the magnetic field in the wind, B wind ∝ 1/R. As a result, the expanding bubble comes into force-balance with the preceding wind near the light cylinder, at (since R 0 ≤ R N S ). Thus, the expelled blob quickly reaches pressure equipartition with the wind flow, at ∼ few light cylinder radii. It is then advected (puffs-up) "quietly" with the wind.
Alternatively, as we demonstrate in this paper, if the explosion is powerful enough the expelled blob detonates, and is totally disrupted. The condition for the explosion is The plan of the paper as follows. In §2 we perform a number of 3D relativistic MHD simulations of highly-over-pressurized spheromak configurations into expanding pre-explosion wind. In §3 we describe corresponding 2D force-free simulations, with external magnetic field decreasing with radius, (either as 1/r or 1/r 2 -this mimics the dependance of the decreasing magnetic field in the accelerating pulsar wind). These two types of simulations probe somewhat different aspects of the problem: MHD is better suited to take account of dissipative plasma effects, while force-free limit studies the extreme magnetic dynamics. In §4 we discuss a compact analytical solution to the Prendergast problem: relativistic self-similar expansion of force-free plasma.
The main part of the simulations' set up is the properties of the external medium. We consistently explore several possibilities. First, for magnetically confined spheromak the static external field is reduced with respect to the internal magnetic field -this mimics over-pressurized magnetic cloud ready to explode, §2.3. Second, in addition to magnetic pressure drop at the edge of a spheromak we introduce a Hubble-like flow outside, §2.4 (this mimics the pre-explosion wind, see Fig.  1). We investigate various combinations of these set-ups, the dependence of the dynamics on the magnetization and the properties of the field jump at the boundary.
The simulations are performed with PLUTO code (Mignone et al. 2007) in 3D relativistic MHD (RMHD) approximation, which used the WENO3 based method (Yamaleev & Carpenter 2009), 3 rd -order TVD Runge Kutta time stepping procedure, and an HLLD Riemann Solver (Mignone et al. 2009). PLUTO is a public modular Godunov-type code entirely written in C intended mainly for astrophysical applications, and high Mach number flows in multiple spatial dimensions.
Energy, and not entropy, is the independent variable, and an Eulerian scheme is adopted when solving the equations of RMHD. Special relativity is sufficient to describe the problem at hand. The duration of the time steps is determined by the adopted CFL number of 0.2, to keep the stability of the simulations. A static uniform smooth grid allows an adequate control of the resolution in the whole computational grid during a whole run. The simulated flow was approximated by ideal gas with polytropic index Γ g = 4/3. In this section we used non-dimensional relativistic units p g,0 = ρ 0 c 2 .

Basic spheromak
The internal magnetic field of a basic spheromak in spherical coordinates r − θ − φ can be expressed in terms of the flux function Ψ in .
where α = 4.49341, a is the radius of a spheromak, Fig. 2. Explicit expressions for magnetic field are given in Appendix A. The spheromak is set as a sphere with radius a = 1. Inside the spheromak we set initial density and pressure ρ = p g = 1.
We start with a model of magnetically confined spheromak with external magnetic field

Magnetically over-pressured basic spheromak in static medium
To initialize an explosion we decrease the normalization of external magnetic field This introduces a jump at the surface so that the spheromak is magnetically over-pressurized.
Outside of the spheromak we set the magnetic field with λ = 10 −3 . The initial velocities were set equal to 0.
We performed small scale simulations of spheromak evolution in rarefied surrounding media ρ = p g = 3 × 10 −3 . We choose magnetic field which corresponds to high magnetization σ = B 2 /4π( g + p g ) = 60, medium magnetization σ = 2 and low magnetization σ = 1/10 3 . The initial total pressure was kept the same in all simulations m + p m + g + p g = const. and for σ = 60 we choose p g = 1. So, reducing the magnetic pressure we increase the plasma pressure respectively keeping total energy density the same. Here we performed so-called small scale simulations with the grid size was taken to be x, y and z ∈ [−2, 2], in initial spheromak radius units "a". The total number of cells in each direction: The examples of velocities profile in X-direction (equatorial plane of spheromak) are shown in Fig. 3. The strongly magnetized solution significantly differ from low magnetized one. Low magnetized solution shows graduate acceleration of flow on its boundary and due to absence of hoop stress the relatively slow acceleration in the central zone. The highly magnetized case shows no expansion inside the spheromak so far rarefaction wave did not reach the given radius. The shock formed on the spheromak's surface accelerated quickly at first, but at some point outflow The middle magnetized case shows very slow growth with time. The strongly magnetized case shows fast initial rise and slow down later on. In late phases the speed of outflow is inversely proportional to the magnetization of spheromak.

Over-pressured spheromak in external Hubble-like flow (the wind)
Next we model interaction of the ejected magnetic blob with the preexisting wind. We assume that a magnetic cloud reached a size of the light cylinder, and starts to interact with the preexisting wind. The wind is accelerating linearly, with v ∝ r. Instead of varying the properties of the wind (e.g., mass loss rate) it is more convenient to vary the velocity of the wind. We use different scalings v = η v cr/a, with η v in the range 0.01 ≤ η v ≤ 0.1. Small values of η v correspond to nearly static external configurations (this allows comparison with results in §2.2) , while large η v → 1 mimic realistic pulsar winds.
In this section we used large scale simulations the grid size was taken to be x, y and z ∈ [−4.5, 4.5], in units of initial spheromak radius. To probe different velocity expansion, we performed simulations with total number of cells in each direction: N x = N y = N z = 384.
In the wind zone spheromak can be described in the comoving rest frame. In the wind zone the toroidal component of magnetic field dominates with radius, so in this section we can choose Z axis direction along large scale toroidal magnetic field (Komissarov et al. 2007(Komissarov et al. , 2009). The magnetic strength of the spheromak was chosen to set σ = 3 at the spheromak center and uniform density and pressure ρ = p g = 1. The current density and pressure distribution are plotted in Fig. 5.
In the absence of the expansion (see §2.4) the highly magnetized over-pressurized spheromak expands up to new equilibrium radius which can be estimated as r eq ≈ a(p tot,in /P tot,out ) 1/4 (Lyutikov & Gourgouliatos 2011) and for chosen parameters r eq ≈ 2a. The radial evolution of spheromak in stationary environment with a drop of the external pressure in 16 times are shown on Fig. 6. The initial expansion in equatorial plane (X-direction) after passing the equilibrium radius slows down, and later switched to shrinking and saturation near 2.5a. In polar direction (Z-axis) expansion was not so strong and saturates near 2a.    evolved ones. In the case of slow expansion (η v = 0.01) magnetic flux of the spheromak preserve its shape for many tens of dynamical times. The explosion like solutions leads to strong redistribution of magnetic field in spheromak. Fig. 9 illustrates the key difference in the expansion process: the Mach number maps for fast magnetosonic speed on the (We calculate the Alfven velocity as v a /c = σ/(σ + 1) and the fast magnetosonic speed can be estimated as v f /c = (σ + c 2 s )/(σ + 1) Lyubarsky 2003). In the explosion-like expansion we clearly see supersonically expanding parts of the flow. If the spheromak initially was in equilibrium even for fast expansion speed (η v = 0.3) the flow remains subsonic, no explosion is generated, in contrast to over-pressured spheromaks (see Fig. 9).
Finally, we investigated evolution of the spheromak in the Hubble expanding wind which initially was in pressure equilibrium. We increase the expansion speed up to η v = 0.3 to increase the effect of expansion. The expansion takes place in the subsonic regime during all the simulation evolution. The Mach number distribution and velocity field near the end of the simulation are presented in Fig. 10. The flow outside of the spheromak is supersonic, but spheromak itself is subsonic. We did not see the change of the magnetic field topology and its significant redistribution in radial direction.

Higher order spheromak
Next, we used the same setup as in the previous Section 2.4, but we change topology of the magnetic field to the 3rd order spheromak with the flux function Ψ 3 ∝ 1 + 5 cos 4 (θ) − 6 cos 2 (θ) × 15a 3 sin αr a α 3 r 3 − 15a 2 cos αr a α 2 r 2 − 6a sin αr a αr + cos αr a where α = 6.98793, Fig. 11. This magnetic field has three changes of magnetic polarity in the θ direction. This configuration mimics expected highly tangled magnetic field of the ejection. The current density of the high order spheromak is presented in Fig. 12 and the 3D rendering of magnetic field lines in Fig. 13 and slices of the magnetic field in Fig. 14. We observe a behavior similar to the basic spheromak: for supersonically expanding case the magnetic field is concentrated near the surface, while subsonically expanding case approximately keeps the initial field structure. The expansion is more isotropic in this case.

Force-free simulations of a magnetic explosion
In order to further verify the results of MHD simulations, we performed 2D force-free simulation of a spheromak using PHAEDRA code (Parfrey et al. 2012). PHAEDRA is a pseudo-spectral code developed specifically to study highly magnetized plasma regime, force-free electrodynamics, the vanishing-inertia limit of magnetohydrodynamics. The grid is defined in axisymmetric spherical coordinates (r, θ), with N r points in r and N θ points in θ direction. For this work, we chose N r = 280 and N θ = 160.
First, we repeat the experiment with magnetically confined over-pressured spheromak, §2.3. Similar to §2, the external field is reduced with respect to the internal magnetic field. We set the external field to be a fraction λ of the equatorial spheromak's field on the surface.

As a new ingredient, we investigated different profiles of the external magnetic field, scaling with radius as
(This is the scaling of the flux function Φ out , to conserve the div-B condition.) Qualitatively, this mimics the decreasing magnetic field in the wind.
The outer boundary of the setup was fixed at r = 20, absorbing layer beginning at r = 15 whereas the radius of the spheromak was chosen to be a = 2 (somewhat larger initial size is chosen to avoid a quickly developing tilting instability Mehta et al. 2020). The distances are in code unit. We tried various grid sizes and found no substantial difference in the results for very high resolutions. We ran our simulations for number of combinations of λ = 1, 1/3 and γ = 0, 1. Results are plotted in Figs. 15-16.
Results of force-free simulations are qualitatively similar to the MHD case: the evolution of spheromak depends on the profile of the external field. We first initialize the over-pressured spheromak by bringing it out of force balance by setting λ = 1/3. The magnetic blob first puffs up, launches Alfvén waves in the external medium, but later finds a new equilibrium (resembling original configuration), Fig. 15. In the case of decreasing external magnetic field ∝ 1/r, Fig. 16, the spheromak detonates, creating an expanding flux tube. (We also did simulations for B ex ∝ 1/r 2 , results are similar.)

Analytical solution to the Prendergast problem
Next we supplement our numerical results with analytical solutions for expanding force-free spheromaks. The mathematical equation describing this problem has been derived by Prendergast (2005) who considered a system of inertial frames of reference without any electric field, but only force-free magnetic field. Then, by transforming back to the origin frame of reference and postulating an appropriate expansion rate, he obtained the partial differential equation describing this system. Solutions that resemble expanding jets and analogues of the magnetic towers in the relativistic regime have been explored in this context (Gourgouliatos & Lynden-Bell 2008;Gourgouliatos & Vlahakis 2010), nevertheless, the spheromak family of these solutions has not been studied thoroughly. Here, we present a more intuitive derivation for the governing equation, where both the magnetic and electric field are accounted for an observer who is at rest with respect to the origin by considering both the magnetic and the electric field (Gourgouliatos 2009).
Consider an axisymmetric time-dependent magnetic field in spherical coordinates (r, θ, φ): By construction the above magnetic field satisfies Gauss' law, and is expressed in terms of two scalar functions Ψ and T for the poloidal and toroidal field. Let the whole system expand with the We note that this velocity is not the drift velocity defined as E×B/B 2 , as the former has components parallel to the magnetic field, whereas the latter is always normal to the magnetic field.
Let us introduce the dimensionless parameter  (8). This mimics magnetic field in the wind. An expanding flux tube is generated: the spheromak first expands, and then detonates generating causally disconnected outflow. The spheromak is completely destroyed. (Motion in the wind ahead of the spheromak-generated Alfvén pulse, bottom left panel, appears due to the fact that the initial configuration is not in force balance.
We demand that the poloidal flux depends on the following combination of variables: and The induced electric field is: We substitute into the induction equation: Ther andθ components are satisfied identically, but theφ component sets the extra constraint: v which is straightforward to integrate. Thus, the expression for T becomes: Next, we take the force equation: For the evaluation of the electric charge and current we use Gauss's law for the electric field and Ampère's law with Maxwell's correction: We set all three components of the force to zero. From the φ component we take the equation: On division by v 2 and we obtain the following equation: The above expression is the Jacobian of Ψ and 1−v 2 v W (v, θ), therefore we can write: with β an arbitrary function of Ψ. We substitute into the r and θ components of the force equation, we substitute µ = cos θ and we obtain the "Prendergast Equation" : We found fully analytic, separable, self-similar solution for a relativistically expanding forcefree flow for a dipolar structure and linear dependence between Ψ and β in the form: Under this condition, the equation for f becomes: The equation admits the analytical solution: Combining equation (17), (23) and (26) we obtain the expression for T : It is then straightforward to evaluate the fields and the current: Even though these equations may seem over-complicated: they closely match the well know spheromak solution (e.g., function f ), yet the independent variable is given in terms of the variable arctanhv. These are fully analytic self-similar solutions for time-dependent axially-symmetric forcefree expanding magnetic bomb. Fig. 17.-The function f of the Prendergast solution for three choices of α = 1 , 2 , 5 shown in red, green and blue lines respectively. The solution for α = 1 has its first root at v 0 = 0.9992, corresponding to maximum Lorentz factor Γ 0 = 25.4, the one for α = 2 has v 0 = 0.9715 and Γ 0 = 4.2 and the one for α = 5 has v 0 = 0.71 and Γ 0 = 1.4. We also plot the solutions for two stationary spheromaks appropriately scaled so that their first roots lie at 0.71 (Sph1) and 1 (Sph2), for comparison.
The location of the first positive root, v 0 > 0, of the function f (v) depends on the value of α. For smaller values of a it is pushed towards v → 1, whereas for a larger α the radius of the spheromak becomes smaller. As α is a factor relating the poloidal current to the poloidal field, a stronger current leads to a field of smaller radius as shown in Fig. 17, where we present the solution of equation (28), for α = 1 , 2 , 5. Indeed, we can evaluate the relation between the outer boundary of the bubble with α for the expanding spheromak, Fig. 18 For example: In the stationary spheromak solution, the change of α is a simple scaling factor as is evident from equation (4). In the relativistic case however, a change in α is not a mere rescaling, but a drastic alteration in the structure of the spheromak. This is evident from Fig. 19, where we show the structure of two relativistic spheromaks for α = 2 and α = 5, a higher value of α leads to a solution that is reminiscent of the stationary one, whereas a smaller α leads to a field where most of the field lines are concentrated near the edge of the expanding bubble and so does the energy and the pressure. We further comment, that the contours of the poloidal field lines given by Ψ do not coincide to the contours of T , which is the case for the stationary spheromak. This difference is related to the fact that in the stationary case it is only the j × B term that contributes to the force, thus the electric current must flow along the field lines. On the relativistically expanding case, the contribution of the electric field and the charge in the dynamics leads to this difference.
Moreover, while, α, the constant of proportionality between Ψ and β is smaller this does not imply that the toroidal field becomes weaker. This is because the maximum of f is pushed towards higher velocities. Given the factor v/(1 − v 2 ) that appears in the form of T , this term dominates, leading to a strong toroidal field. The part of the solution that is beyond the first positive root, shows an oscillatory behavior. This part corresponds to field lines that are not connected to the rest of the spheromak, thus we focus our attention to the inner sphere. Overall, the field is rather strong near the boundary of the bubble, Fig. 20. We note that in the limit of low expansion velocities (v 1), where only the first order term is kept, it reduces to the standard spheromak solution as arctanhv ≈ v.
Even though the field within the spheromak itself is force-free, as the bubble expands the energy within the bubble decreases as ∝ 1/t. Indeed the total energy is The energy decrease can be clearly seen from the requirement of the conservation of the magnetic flux: B ∝ 1/R 2 , so that B 2 R 3 ∝ 1/t. This decrease in energy is the a work done by the expanding blob on the external medium. In the origin frame the surface magnetic field is The corresponding work made by the pressure due to the expansion is Thus, the energy of the explosion evolves according to We note that the magnetic pressure at the surface of the bubble is highly anisotropic as is evident by the form of the meridional component of the field there, Equation (36), which is much stronger at the equator than higher latitudes. Thus, if such an expanding spheromak is embedded within a uniform medium it will expand primarily along the equator. Such a disturbance will affect its spherical boundary and the force-equilibrium at its interior. A possible extension could be the inclusion of a thermal pressure pressure profile, which removed the electric current discontinuity which is responsible for this anisotropy Gourgouliatos & Lyutikov 2012).

Discussion
In this work we investigate the dynamics of relativistic magnetically driven explosions. Our analytical results, MHD and force-free simulations consistently show a complicated, parameterdependent picture, with various important effects. Most importantly, the dynamics differs qualitatively from the fluid case. The most exemplary case, demonstrating the differences between highly a b Fig. 19.-The structure of the spheromak for α = 2 (a) and α = 5 (b). The poloidal field lines are shown in black, and in color we plot log T and T respectively. As T depends on both v and r, equation (29) the snapshot is taken at time t = 1/c. magnetized and fluid cases, is the magnetic explosion in static external medium, §2.3. For example, Fig. 4 demonstrates that low-σ over-pressurized configurations create relativistically expanding flows which accelerated with time, while high-σ configurations shows tendency of de-acceleration and saturation to the constant velocity.
Similarly, just over-pressurized highly magnetized cloud in constant density/constant confining magnetic field would just puff-up -no explosion is generated. Slightly over-pressurized spheromak in a Hubble flow just quickly expand trying to rich pressure balance with preceding wind. Only in the case of highly over-pressurized initial configuration, and the external wind environment the magnetic cloud "explodes" -generates supersonically expanding (causally disconnected) outflow.
In case of over-pressured spheromak expansion into external wind, if the surrounding expansion speed is small (η v = 0.01, or η v = 0.032) the spheromak evolution is similar to the steady case. The initial growth is not so fast and we see bounce from the external matter. The fast expansion of surrounding matter (η v = 0.1) allows fast expansion of spheromak and generation of explosion-like, causally-disconnected expansion. In such the case the spheromak transforms to shell like structure with "hollow" central part.
Thus, the internal structure of the explosion depends crucially on the type of expansion. Significant redistribution of the magnetic field take place in the case of supersonic expansionmagnetic field is concentrated near the surface of the cloud: this is seen in MHD simulations and is analytically predicted, see §4. On the other hand, if expansion is subsonic, the magnetic field keeps approximately the initial configuration.
Puffing-up solutions did not show significant explosion like events. So it is difficult to expect significant observational manifestations. In other hand, explosion like solution can form shocks and particle acceleration. We shown what for the Hubble like expansion and over-pressured spheromaks if expansion speed exceed > 0.1cr/a. We can interpret it as formation of the spheromak at light cylinder with radius about r sh ≈ 0.1r l , so the variability time can be up to t f lare ≈ 2π/0.1 ∼ 100 times smaller compare to a magnetar rotation period.

Acknowledgements
We are appreciated Prof. S. Nagataki and Dr. J. Mahlmann for useful discussion. The 3D RMHD calculations were carried out in the CFCA cluster XC50 of National Astronomical Observatory of Japan. This work had been supported by NASA grants 80NSSC17K0757 and 80NSSC20K0910, NSF grants 1903332 and 1908590. KNG acknowledges funding from grant FK 81641, University of Patras ELKE.

Data availability
The data underlying this article will be shared on reasonable request to the corresponding author.
B. The magnetic field structure for Hubble expansion case.
The over-pressured spheromak evolve differently depends on Hubble expansion speed. To illustrate this point we show vertical component of the magnetic field. If the surrounding expansion speed is small (v r = 0.01 cr/a) Fig. 21 or (v r = 0.032 cr/a) Fig. 22 the spheromak evolution is similar to the steady case. Here the initial growth is not so fast and we see bounce from the external matter (see the velocity field evolution) and settle to quasi-static expansion. The internal structure of the magnetic field kept the same along all simulation. The relation of maximal to minimal values of the vertical component of the magnetic field is preserved (∼ 0.3).
The fast expansion of surrounding matter (v r = 0.1 cr/a) leads to fast expansion of spheromak