© The Authors 2026

1. Introduction

Magnetars are hosts of highly energetic events that release radiation in X-ray and γ-rays. The most energetic of these events are giant flares, which release an isotropic luminosity up to 10⁴⁷ erg/s (Hurley et al. 1999; Palmer et al. 2005; Hurley et al. 2005). These events have been observed from extra-galactic sources (Beniamini et al. 2025; Trigg et al. 2025) and have been associated with strong magnetised explosions (Thompson & Duncan 1995) originating from magnetars where a major reconfiguration of the global magnetic field has occurred (Lyutikov 2006). Magnetar flares also are loaded with mass entrained from the crust and the interior of the neutron star (Frail et al. 1999; Demidov & Lyubarsky 2023; Patel et al. 2025). Thus, a complete view of giant flares requires modelling a mass-loaded magnetised explosion, reaching relativistic velocities. Here, ‘mass-loaded’ refers to finite plasma inertia and pressure within the Magnetohydrodynamical (MHD) framework, rather than explicit source terms in the continuity equation.

The collection of an ample sample of fast radio bursts (FRBs), which are millisecond-duration transients of extragalactic origin, allows us to explore their diverse phenomenology, which indicates coherent emission from compact astrophysical sources (Lorimer et al. 2007; Petroff et al. 2019). Since their discovery, FRBs have challenged theoretical models, demanding mechanisms capable of releasing up to 10⁴⁰ erg in radio waves over millisecond timescales, often with repetition, polarisation, and spectral properties indicative of strong magnetised environments (Marcote et al. 2020). The growing evidence linking at least a subset of FRBs to magnetars, particularly the association of FRB 200428 with the Galactic magnetar SGR 1935+2154 (Bochenek et al. 2020; CHIME/FRB Collaboration 2020), has reinforced the idea that magnetic energy release in highly magnetised neutron stars can power these events.

Despite this association, the detailed physics linking magnetic activity to FRB emission remains uncertain. The proposed mechanisms span a broad range from magnetospheric reconnection and crustal cracking to relativistic shocks in magnetar winds (e.g. Lyubarsky 2014; Beloborodov 2017; Metzger et al. 2019). Most models assume magnetically dominated outflows capable of generating coherent radio emission through synchrotron maser, shocks, or antenna processes, and possibly a combination of these (Lyubarsky 2021). Theoretical and observational developments suggest that magnetar magnetospheres may often be mass-loaded by plasma ejected during starquakes, bursts, or fallback episodes (Thompson & Duncan 1995; Wadiasingh & Timokhin 2019; Yuan et al. 2020). Moreover, the connection between magnetar bursts and flares and FRBs that has recently been confirmed (Bochenek et al. 2020; Tavani et al. 2021; Beniamini et al. 2025) stresses the necessity to consider mass-loaded magnetic detonations that might be related to an event originating from the neutron star interior and propagating to the magnetosphere (Bransgrove et al. 2025). The presence of dense plasma in a strongly magnetised environment can substantially modify the dynamics of magnetic explosions, their ability to accelerate particles, and the conditions for coherent emission.

We explore mass-loaded magnetic explosions as potential origins of magnetar giant flares and FRBs. Using simplified relativistic magnetohydrodynamic (MHD) models, we examine how the interplay between the magnetic field, plasma mass loading, and self-similar expansion shapes the energetics and temporal properties of such events. The proposed structures can provide an appropriate environment that can lead to the generation of short, intense events with energetics and timescales consistent with observed FRBs.

To model the coupled evolution of the magnetic field and plasma during these explosive events, we adopted the formalism developed in Gourgouliatos & Vlahakis (2010), which provides a self-consistent description of self-similar, co-expanding magnetised outflows. In this framework, the plasma and magnetic field evolve together under a homologous (Hubble-type) expansion, preserving the topology of the magnetic field while allowing for a time-dependent redistribution of electromagnetic and kinetic energy.

The underlying concept traces back to the Prendergast (2005) family of self-similar, uniformly expanding, axisymmetric magnetic equilibria. These equilibria, later extended by Gourgouliatos & Lynden-Bell (2008) to include non-linear poloidal currents and toroidal fields, demonstrated that magnetic structures can expand coherently, a key insight for modelling transient magnetic explosions in compact objects. Building on this foundation, Gourgouliatos & Vlahakis (2010) formulated a relativistic configuration that couples the magnetic and velocity fields in a co-expanding MHD flow, enabling the analytical and numerical exploration of outflows that evolve self-similarly while conserving magnetic helicity and flux. In these models, the field originates from an appropriately twisted dipolar located at the origin. Due to the singularity at the origin, the field is studied beyond an inner radius that serves as a lower boundary condition.

More recently, Barkov et al. (2022) applied similar ideas for numerical solutions of magnetic explosions in magnetically dominated and mass-loaded environments, demonstrating that magnetic energy release can drive relativistic outflows and shocks capable of producing bright, short-duration transients. In the same work, analytical solutions of force-free magnetic fields in the absence of thermal pressure and inertia forces were presented. These solutions do not diverge at the origin. It was found that the maximum expansion rate that can be attained by the expanding structure depends on the amount of twist of the magnetic field, with highly twisted magnetic fields not being able to expand at highly relativistic rates. In the absence of any thermal pressure and inertia, the pressure on the surface of the boundary of the spherical explosion is highly anisotropic due to the magnetic field structure, that is, in the dipole case, it is maximum at the equator and becomes zero at the axis of symmetry. When a structure like this is embedded within a uniform pressure environment, it will expand sideways on the equator and will form a highly anisotropic explosion, while it will lose its internal equilibrium and coherence. The inclusion of mass density and thermal pressure leads to configurations in the form of large-scale plasmoids whose surface pressure, including the electromagnetic term, is more isotropic, and under certain conditions, it can even become completely isotropic, with the surface magnetic field completely vanishing, as was shown in the non-relativistic approach (Gourgouliatos et al. 2010). These structures continue their uniform expansion when they are situated in a uniform density medium and maintain their internal equilibrium (Gourgouliatos & Lyutikov 2012).

In the present work, we extend these ideas by embedding mass-loaded magnetic explosions within the relativistically exapanding self-similar MHD formalism of Gourgouliatos & Vlahakis (2010). This allows us to consistently follow the co-evolution of magnetic fields and plasma during an explosive expansion and to explore the interplay between magnetic field strength, mass, and pressure and its effect on the expansion rate. This will allow us to identify the physical conditions under which these events can generate coherent radio emission with energetics and timescales comparable to those of magnetar flares and FRBs.

This paper is organised as follows. In Sect. 2 we outline the physical and mathematical framework for magnetic plasmoid explosions and describe the key ingredients of our model. Section 3 presents the solutions and the exploration of the parameter space. The results are discussed in Sect. 4, focusing on energy release, plasma dynamics, and physical properties of the solution. In Sect. 5 we discuss the implications for astrophysical sources with a main focus on magnetar bursts and FRBs. We conclude in Sect. 6 and outline the limitations of the current model and the prospects for future work.

2. Mathematical formulation

We considered a fluid of rest-mass density ρ and pressure p containing a magnetic field B. The electric and magnetic field obey Maxwell’s equations,

$$\begin{array}{cc}\hfill \mathrm{\nabla }·\mathit{B}& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \mathrm{\nabla }·\mathit{E}& =\frac{4\pi }{c}{j}^0,\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \mathrm{\nabla }\times \mathit{E}& =-\frac{1}{c}\frac{\partial \mathit{B}}{\partial t},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \mathrm{\nabla }\times \mathit{B}& =\frac{1}{c}\frac{\partial \mathit{E}}{\partial t}+\frac{4\pi }{c}\mathit{j},\hfill \end{array}$$(4)

where j⁰ is the electric charge density, j is the electric current density, and c is the speed of light. We further considered ideal MHD, and there is thus no electric field in the co-moving frame of the fluid. For an observer who sees the fluid moving at velocity v, normalised to the speed of light, Ohm’s law becomes

$$\begin{array}{c}\hfill \mathit{E}=-\mathit{v}\times \mathit{B}.\end{array}$$(5)

The continuity equation reads

$$\begin{array}{c}\hfill {\partial }_t\left(\mathrm{\Gamma }\rho \right)+\mathrm{\nabla }·\left(\mathrm{\Gamma }\rho c\mathit{v}\right)=0,\end{array}$$(6)

where Γ is the Lorentz factor,

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{1}{\sqrt{1-v^2}}·\end{array}$$(7)

The momentum equation is

$$\begin{array}{c}\hfill \mathrm{\Gamma }\rho ({\partial }_t+c\mathit{v}·\mathrm{\nabla })\left(\xi \mathrm{\Gamma }c\mathit{v}\right)+\mathrm{\nabla }p-\frac{j^0\mathit{E}+\mathit{j}\times \mathit{B}}{c}=\mathbf{0},\end{array}$$(8)

where p is the pressure, and the relativistic specific enthalpy is given by the following expression:

$$\begin{array}{c}\hfill \xi =1+\frac{\gamma }{1-\gamma }\frac{p}{\rho c^2},\end{array}$$(9)

where γ is the ratio of the specific heats. Furthermore, the entropy equation reads

$$\begin{array}{c}\hfill ({\partial }_t+c\mathit{v}·\mathrm{\nabla })\left(\frac{p}{{\rho }^{\gamma }}\right)=0.\end{array}$$(10)

In what follows, we assume a relativistic degenerate gas where γ = 4/3. We further considered spherical coordinates (r, θ, ϕ) and prescribed a velocity normalised to the speed of light given by the following expression:

$$\begin{array}{c}\hfill \mathit{v}=\frac{r}{\mathit{ct}}\widehat{\mathit{r}}.\end{array}$$(11)

This velocity corresponds to uniform expansion from the centre. Thus, each fluid element moves at constant velocity without any acceleration, and the system remains implicitly at force balance. Here and throughout the paper, v denotes the dimensionless radial expansion velocity v = r/(ct).

It is possible to find self-similar, separable solutions using the polar angle θ and the expansion velocity v as main variables. The full details of this procedure were shown by Prendergast (2005), Gourgouliatos & Lynden-Bell (2008), Gourgouliatos & Vlahakis (2010), and we summarise the basic steps for clarity and consistency below.

We express the magnetic field in terms of two functions: one related to the poloidal field Ψ(θ, v), and the other to the toroidal T(θ, v),

$$\begin{array}{c}\hfill \mathit{B}=\frac{1}{2\pi r^{2}sin\theta }({\partial }_{\theta }\mathrm{\Psi }\widehat{\mathit{r}}-v{\partial }_v\mathrm{\Psi }\widehat{\mathit{\theta }}+T\widehat{\mathit{\varphi }}).\end{array}$$(12)

Combining Ohm’s law, Eq. (5), with the momentum equation, Eq. (8), and taking the toroidal component, we deduce that the poloidal and toroidal functions are related to each other through the following expression:

$$\begin{array}{c}\hfill T={\mathrm{\Gamma }}^{2}vH(\mathrm{\Psi }),\end{array}$$(13)

where H(Ψ) is an arbitrary function of Ψ.

The pressure and the density are related through the entropy equation, Eq. (10), which is satisfied for

$$\begin{array}{c}\hfill \frac{p}{{\rho }^{4/3}}=Q(v,\theta ),\end{array}$$(14)

where Q(v, θ) is an arbitrary function of v and θ. Upon integration of the continuity equation, Eq. (6), we obtain that

$$\begin{array}{c}\hfill \rho =\stackrel{\sim }{\rho }(v,\theta )\frac{{\mathrm{\Gamma }}^3{v}^3}{r^3},\end{array}$$(15)

and combining this with Eq. (14), we find that the pressure has the following form:

$$\begin{array}{c}\hfill p=\stackrel{\sim }{p}(v,\theta )\frac{{\mathrm{\Gamma }}^4{v}^4}{r^4},\end{array}$$(16)

where the quantities Q$\stackrel{\sim }{p}$, $\stackrel{\sim }{\rho }$, all depending on the separation variables v and θ, are related to each other through the following expression:

$$\begin{array}{c}\hfill Q=\frac{\stackrel{\sim }{p}}{{\stackrel{\sim }{\rho }}^{4/3}}·\end{array}$$(17)

Next, we evaluate the current and charge densities through Maxwell’s equations, Eqs. (2) and (4), and substitute them into the momentum equation, Eq. (8). We find that $\stackrel{\sim }{p}=\stackrel{\sim }{p}(\mathrm{\Psi })$ and that the system should satisfy the following equation:

$$\begin{array}{cc}& v^2\frac{{\partial }^2\mathrm{\Psi }}{\partial v^2}-\frac{2v^3}{1-v^2}\frac{\partial \mathrm{\Psi }}{\partial v}+\frac{sin\theta }{1-v^2}\frac{\partial }{\partial \theta }\left(\frac{1}{sin\theta }\frac{\partial \mathrm{\Psi }}{\partial \theta }\right)\hfill \\ \hfill & +\frac{v^2}{{(1-v^2)}^2}H\frac{\mathit{dH}}{d\mathrm{\Psi }}+16{\pi }^3{sin}^2\theta \frac{v^4}{{(1-v^2)}^3}\frac{d\stackrel{\sim }{p}}{d\mathrm{\Psi }}=0.\hfill \end{array}$$(18)

Separable solutions in v and θ are feasible under the following assumptions:

$$\begin{array}{c}\hfill \mathrm{\Psi }=g(v){sin}^2\theta ,H\frac{\mathit{dH}}{d\mathrm{\Psi }}=c_0\mathrm{\Psi },\frac{d\stackrel{\sim }{p}}{d\mathrm{\Psi }}=\frac{c_1}{16{\pi }^3}·\end{array}$$(19)

The function g(v) determines the radial distribution of the poloidal magnetic flux under self-similar expansion. Its first non-zero root, v₀, defines the outer boundary of the magnetised plasmoid, beyond which all field components vanish or become disconnected from the inner region. The overall shape of g(v) therefore controls both the spatial extent of the configuration and the presence or absence of surface currents and pressure anisotropies.

Upon integration, and considering the fact the B_ϕ component of the field must remain finite,

$$\begin{array}{c}\hfill H(\mathrm{\Psi })=\pm \sqrt{c_0}g{sin}^2\theta ,\end{array}$$(20)

with c₀ a non-negative constant, which governs the intensity of the toroidal field, or equivalently, the poloidal current. The pressure term becomes

$$\begin{array}{c}\hfill \stackrel{\sim }{p}={\stackrel{\sim }{p}}_0+\frac{c_1}{16{\pi }^3}g{sin}^2\theta .\end{array}$$(21)

In the above expression, ${\stackrel{\sim }{p}}_0$ is the pressure on the axis and the spherical surface of the plasmoid; and it has to be a non-zero constant. There is no constraint on the sign of c₁, but the combination of ${\stackrel{\sim }{p}}_0$ and c₁ must ensure that $\stackrel{\sim }{p}$ remains non-negative throughout the solution.

Finally, substituting the expressions from Eq. (19) into Eq. (18), we obtain an ordinary differential equation for g(v),

$$\begin{array}{c}\hfill v^{2}g″-\frac{2v^3{g}^{\prime }}{1-v^2}-\frac{2g}{1-v^2}+\frac{c_0{v}^{2}g}{{(1-v^2)}^2}+\frac{c_1{v}^4}{{(1-v^2)}^3}=0.\end{array}$$(22)

The first three terms describe the poloidal magnetic structure under relativistic expansion, the fourth term represents the contribution of the toroidal magnetic field (poloidal current), and the final term accounts for pressure and inertia effects. The solution of the above equation provides the radial structure of the configuration.

3. Solutions

The behaviour of the solution near the origin can be inspected for v ≪ 1. The leading terms of the equation are the following:

$$\begin{array}{c}\hfill v^{2}g″-2g=0.\end{array}$$(23)

This admits solutions in following form:

$$\begin{array}{c}\hfill g(v)=A_1{v}^2+A_2{v}^{-1},\end{array}$$(24)

where A₁ and A₂ are constants. Thus, there are two main classes of solutions of Eq. (22). A class where g ∝ v⁻¹, leading to solutions that reduce to a magnetic dipole near the origin, and a class where g ∝ v², where the field near the centre is approximately uniform. While the former case was explored in detail in Gourgouliatos & Vlahakis (2010), the latter has not been studied. In this work, we focus on the latter class of solutions, which are not singular at the origin, by exploring appropriate combinations of the main parameters c₀ and c₁.

3.1. Non-relativistic limit

In the non-relativistic limit, while keeping the terms proportional to c₀ and c₁, Eq. (22) reduces to

$$\begin{array}{c}\hfill v^{2}g″-2g+c_0{v}^{2}g+c_1{v}^4=0.\end{array}$$(25)

The solution of the above equation is

$$\begin{array}{cc}\hfill g(v)& =B_1(\frac{sin\left(\sqrt{c_0}v\right)}{\sqrt{c_0}v}-cos\left(\sqrt{c_0}v\right))\hfill \\ \hfill & +B_2(\frac{cos\left(\sqrt{c_0}v\right)}{\sqrt{c_0}v}+sin\left(\sqrt{c_0}v\right))-\frac{c_1}{c_0}{v}^2,\hfill \end{array}$$(26)

where B₁ and B₂ are constants. The term multiplied with B₁ is zero at the origin, whereas the term multiplied with B₂ is singular. Therefore, we kept the B₁ terms and set B₂ = 0. This solution is mathematically identical to the solutions found in Gourgouliatos et al. (2010, Eq. (6)), where the role of v is taken by r. In the present case, the sphere instead of being static, it expands uniformly.

3.2. Relativistic force-free limit

In the absence of pressure and inertia, we have c₁ = 0. The results of the integration for c₁ = 0 and various values of c₀, ranging from 1 to 100, are shown in Fig. 1. They are fully consistent with the analytical solutions reported by Barkov et al. (2022), with the note that the value of α there corresponds to $\sqrt{c_0}$ in this work. Clearly, when the value of c₀ increases, the maximum velocity attained by the expanding field becomes lower. This velocity corresponds to the first root, beyond the origin v = 0, as g(0) = 0, of the function g(v₀) = 0. Thus, by increasing the ratio of the poloidal current over the poloidal field, or equivalently, the ratio of the toroidal field over the poloidal field, the maximum expansion velocity is lower. We note that beyond v₀, the solutions become oscillatory, with the flux taking negative and positive values. However, this part of the solution is completely disconnected from the inner part in terms of the poloidal and the toroidal fields, and thus, no field lines from the central part are connected to the outer ones. Thus, we halted the integration at v₀. These solutions are denoted as Z-type solutions, as the pressure and density are zero.

Fig. 1. Z-family solutions for g(v) for various values of c0 when c1 = 0.

3.3. General solution

When all terms are included so that c₀ ≠ 0 and c₁ ≠ 0 and the system is relativistic, the equation does not permit general analytical solutions, and we proceed with a numerical solution. Given the singularities at v = 0 and v = 1, we integrated within an interval v ∈ [v_min, v_max], where v_min = 10⁻³ and v_max < 1. We set g(v_min)′ = 0 and g(v_min) = g₀, and we chose g₀ so that the maximum value of g was set to unity. This is straightforward for the homogeneous version of the equation with c₁ = 0, as a simple division with the maximum value of a trial solution with arbitrary g₀ leads to a normalised solution. However, in the non-homogeneous case (c₁ ≠ 0), it requires an iterative processes as simple division of the function does not lead to a solution. Thus, we repeated the process by dividing g₀ by the maximum value of the solution and integrated until convergence. We used a tolerance limit of 10⁻⁴ for the normalisation, and we typically found that the solution converged to a normalised form after fewer than 100 iterations. For the integration of the ordinary differential equation, we used the Runge-Kutta-Fehlberg method, applying the built-in routines of numpy and scipy in Python (Harris et al. 2020; Virtanen et al. 2020).

Next, we explored solutions with non-zero values of c₁. Because of the very many possible combinations, we decided to keep the toroidal field parameter constant (c₀ = 1) and then explore positive values of c₁ ranging from 0.1 to 100. The results are shown in Fig. 2. The inclusion of pressure and density leads to a structure that is less extended in the velocity space and reaches a lower maximum velocity than in the c₁ = 0 case. We denoted these solutions P-type solutions, as the pressure and density increase with increasing magnetic flux function.

Fig. 2. Characteristic P-family solutions of g(v) for various values of c1 while keeping c0 = 1.

Finally, it is possible to explore negative values of c₁. These cases correspond to plasmoids where the pressure decreases as the magnetic flux function increases. Despite this negative relation, the freedom of an additive constant ${\stackrel{\sim }{p}}_0$ in Eq. (21) prevents negative values of the pressure and density. These solutions have the remarkable property that for an appropriate combination of c₀ and c₁, not only the flux function reaches 0 at v₀, but its derivative also becomes zero there, and thus, g′(v₀) = 0, as shown in Fig. 3. This has the important consequence that all three components of the magnetic field are zero on the surface of the sphere. We note that when the values of c₁ are smaller algebraically than those corresponding to the zero derivative on the boundary, the solutions become unphysical, with the field lines not closing, but increasing infinitely. The resulting structures correspond to underdensities, where the drop of thermal plasma pressure and mass density within the sphere is compensated for by the magnetic and electric field pressure and tension. We denoted these as N-type solutions as in these solutions the density and pressure decrease with increasing flux function.

Fig. 3. N-family of solutions showing g(v) for various values of c0 and negative c1, so that the derivative is zero at the maximum velocity.

The parameters provided in Table 1 were employed for the three families of solutions: c₁ = 0 Z-models, without mass or thermal pressure, c₁ > 0 P-models, with varying mass and pressure with the same sign as the flux function, and c₁ < 0 N-models, where, in addition to the negative sign of c₁, the derivative of g at v₀ vanishes.

4. Discussion

4.1. General properties of the solutions

We describe the main characteristics of the solutions depending on the values of the parameters below. For all solutions, lower values of c₀ lead to a v₀ that approaches unity, and thus, the system becomes more relativistic. This implies that the more twisted the solution, or equivalently, the higher the poloidal current, the smaller the sphere containing the plasmoid. This property is related to the fact that by assuming a constant c₀, a particular length-scale is introduced. In this case, the units of $\sqrt{c_0}$ are inversely proportional to the length. Consequently, by increasing c₀, this length scale becomes smaller. This is consistent with the well-known non-relativistic linear force-free axisymmetric solutions in spherical geometry (Chandrasekhar & Kendall 1957), whose radius shifts to lower values for a higher toroidal-to-poloidal magnetic field ratio. This property is well preserved and reflected on its relativistic analogue.

In the Z and P solutions, a non-vanishing meridional magnetic field exists on the surface v = v₀. Because of this magnetic field discontinuity, an azimuthal electric current sheet flows on the surface of the plasmoid. Moreover, there also exists a volume electric charge because the meridional component of the electric field diverges. Following the convention of a positive value for the flux function, the charge is positive in the northern and negative in the southern hemisphere. While this structure is in equilibrium at its interior, the magnetic pressure is anisotropic on its surface, stronger at the equator, and weaker at the poles. Therefore, when it is embedded within a medium of uniform pressure, it is expected to expand in the equatorial direction, while the surface current might trigger resistive instabilities and lead to reconnection in a layer that propagates relativistically.

Considering the solutions with c₁ > 0, we find that higher values of c₁ lead to plasmoids with smaller v₀. This is due to the inclusion of additional pressure and density that prevent the system from extending to a v₀ as high as in the purely electromagnetic case (c₁ = 0). The characteristics of the overall structure are similar to those of the previous solution regarding the surface azimuthal current sheet and pressure profile, and the structure is prone to similar instabilities.

Solutions of the N family with c₁ < 0 lead to different qualitative properties. In general, choosing a negative c₁ leads to v₀ closer to unity than in solutions with zero or positive c₁ and the same c₀. Below a critical value of c₁, the solution has no roots and thus becomes unphysical. In the N-class of solutions, we determined these critical values for various choices of c₀. The solutions corresponding to these critical values have the remarkable property that at v₀, not only is the function g zero, but its derivative also vanishes. This implies that all three components of the magnetic field and the two components of the electric field are zero on the spherical surface defined by v = v₀. Because of this, there is no anisotropy or a surface current due to the meridional component. Thus, when spherical plasmoids obeying these solutions are embedded within a medium of uniform pressure, the system can maintain its equilibrium, as opposed to structures corresponding to the other families of solutions. Moreover, the lack of surface current sheets makes the appearance of current-driven instabilities less likely.

4.2. Physical quantities

Based on these solutions, it is possible to calculate the corresponding physical quantities arising from these calculations. In what follows, we focus on the electromagnetic field energy, the mass of the sphere, the work of the thermal pressure, and the temperature of the plasmoid. A set of characteristic configurations and the corresponding physical properties are shown in Figs. 4 and 5.

Fig. 4. Contour plots of the solutions with c0 = 1 and c1 = 1 (top row) and c0 = 1 and c1 = 10 (bottom row). The black contours correspond to the poloidal magnetic field lines, and we plot in colour the Bϕ component of the magnetic field (panels a and e), the rest-mass density ρ (panels b and f), the temperature T (panels c and g), and the electric charge density j0 (panels d and h). The solution remains invariant in the v-space. We present a snapshot at some time t, so that r = 1 corresponds to v = 1.

Fig. 5. Contour plots of the solutions with c0 = 10 and c1 = −1.448 (top row) and c0 = 100 and c1 = −211.3 (bottom row). The black contours correspond to the poloidal magnetic field lines, and we plot in colour the Bϕ component of the magnetic field (panels a and e), the rest-mass density ρ (panels b and f), the temperature T (panels c and g), and the electric charge density j0 (panels d and h). As in the plots of Fig. 5, the snapshot is taken at some time t, so that r = 1 corresponds to v = 1.

4.2.1. Electromagnetic field

The electric and magnetic fields are given from Eqs. (5) and (12) with appropriate substitution of the solutions and the velocity profile. The electromagnetic field is contained within a sphere in the v-space. The poloidal field is approximately uniform close to the axis, and it closes within an expanding sphere of radius R = cv₀t. The B_ϕ field peaks at the equatorial plane. However, this peak does not coincide with the maximum of the flux function Ψ, which corresponds to the null point of the closed poloidal loops (Figs. 4 and 5 panels a and e) because B_ϕ is proportional to Ψ, but is also multiplied by a factor $\frac{v{\mathrm{\Gamma }}^2}{r^2}$, which for the highly relativistic case dominates and pushes the maximum to more distant radii. Therefore, in the highly relativistic case, a thin layer containing a toroidal field expanding relativistically is expected, while in the non-relativistic case v₀ ≪ 1, the toroidal field is contained closer to the centre.

We then calculated the energy of the electric and the magnetic field,

$$\begin{array}{cc}\hfill E_{\mathrm{mag}}& =\frac{2\pi }{8\pi }{\displaystyle {\int }_0^{\pi }{\int }_0^R}\frac{1}{{(2\pi r^{2}sin\theta )}^2}[{\left({\partial }_{\theta }\mathrm{\Psi }\right)}^2\hfill \\ \hfill & +{\left(v{\partial }_v\mathrm{\Psi }\right)}^2+{\left({\mathrm{\Gamma }}^{2}v\sqrt{c_0}\mathrm{\Psi }\right)}^2]r^{2}sin\theta drd\theta \hfill \\ \hfill & =\frac{1}{16\pi ct}{\displaystyle {\int }_0^{\pi }{\int }_0^{v_0}}\frac{1}{v^2{sin}^2\theta }[{\left({\partial }_{\theta }\mathrm{\Psi }\right)}^2\hfill \\ \hfill & +{\left(v{\partial }_v\mathrm{\Psi }\right)}^2+{\left({\mathrm{\Gamma }}^{2}v\sqrt{c_0}\mathrm{\Psi }\right)}^2]sin\theta dvd\theta .\hfill \end{array}$$(27)

It follows from Eq. (27) that the magnetic field energy scales inversely with time. This behaviour can be understood from simple scaling arguments: under self-similar expansion, the magnetic flux is conserved, implying that the characteristic magnetic field strength scales as B ∝ R⁻². Since the magnetic energy density scales as B² ∝ R⁻⁴, integration over the expanding spherical volume (∝ R³) yields a total magnetic energy E_mag ∝ R⁻¹. Because the outer radius of the configuration evolves as R = cv₀ t, this directly leads to E_mag ∝ t⁻¹.

Similarly, the energy stored in the electric field is given from the following expression:

$$\begin{array}{cc}\hfill E_{\mathrm{el}}& =\frac{2\pi }{8\pi }{\displaystyle {\int }_0^{\pi }{\int }_0^R}\frac{v^2}{{(2\pi r^{2}sin\theta )}^2}[{\left(v{\partial }_v\mathrm{\Psi }\right)}^2+{\left({\mathrm{\Gamma }}^{2}v\sqrt{c_0}\mathrm{\Psi }\right)}^2]\hfill \\ \hfill & \times r^{2}sin\theta drd\theta \hfill \\ \hfill & =\frac{1}{16\pi ct}{\displaystyle {\int }_0^{\pi }{\int }_0^{v_0}}\frac{1}{{sin}^2\theta }[{\left(v{\partial }_v\mathrm{\Psi }\right)}^2\hfill \\ \hfill & +{\left({\mathrm{\Gamma }}^{2}v\sqrt{c_0}\mathrm{\Psi }\right)}^2]sin\theta dvd\theta .\hfill \end{array}$$(28)

As the electric field is the θ and ϕ components of the magnetic field multiplied by the velocity, it follows a similar behaviour, with the electric field energy being inversely proportional to time. Thus, as the system expands, it decreases. We show the electric and magnetic field energy contained within a sphere as a function of c₀ in Fig. 6, where the solution for the Z (c₁ = 0) and N-models (c₁ < 0) is shown. For lower and higher values of c₀, the magnetic field energy is higher and has a minimum for intermediate c₀. We remark that we normalised the maximum value of g and consequently Ψ to unity in the solutions. Because of this, while the magnetic flux is the same in all solutions, its energy content is related to its distribution in the sphere. In the case of a smaller sphere and large c₀, all components of the magnetic field are large since they scale inversely with r. For spheres where v₀ approaches unity, while the radius is higher, g develops a sharp derivative at v₀, leading to a strong B_θ component and energy content. In contrast, the electric field energy has a monotonic dependence on c₀ as its intensity is proportional to v, and therefore, its overall contribution is smaller for smaller spheres in velocity space, even though the magnetic field energy is higher.

Fig. 6. Electric and magnetic energy for the Z and N solutions as functions of c0.

4.2.2. Mass

The density profile of the plasmoid is related to the flux function Ψ, Eqs. (17) and (21), multiplied by the prefactor Γ³v³/r³, as demonstrated in Eq. (15). The two families of solutions we studied here for positive and negative c₁ have qualitatively different density profiles. For c₁ > 0, the minimum is located at the centre of the sphere (r = 0) and peaks close to the location of the maximum of Ψ with some modulation due to the Γ³ prefactor that pushes the maximum closer to the surface. In contrast, for c₁ < 0, the minimum is located close to the maximum of Ψ.

We first considered the rest mass contained within the expanding plasmoid. We integrated the rest-mass density expression given by Eq. (15) within a sphere of given radius R, so that R = cv₀t at some time t, which corresponds to the physical location of the non-zero root of g. By accounting for the axial symmetry of the system, integration over ϕ provides a factor of 2π,

$$\begin{array}{cc}\hfill M& =2\pi {\displaystyle {\int }_0^{\pi }}{\displaystyle {\int }_0^R}{({\stackrel{\sim }{p}}_0+c_1\frac{g(v){sin}^2\theta }{16{\pi }^3})}^{3/4}Q{(v,\theta )}^{-3/4}\hfill \\ \hfill & \times {\left(\frac{\mathrm{\Gamma }v}{r}\right)}^3{r}^{2}sin\theta drd\theta \hfill \\ \hfill & =2\pi {\displaystyle {\int }_0^{\pi }}{\displaystyle {\int }_0^{R/(ct)}}{({\stackrel{\sim }{p}}_0+c_1\frac{g(v){sin}^2\theta }{16{\pi }^3})}^{3/4}Q{(v,\theta )}^{-3/4}\hfill \\ \hfill & \times {\mathrm{\Gamma }}^3{\left(\frac{r}{\mathit{ct}}\right)}^{2}sin\theta d\left(\frac{r}{\mathit{ct}}\right)d\theta \hfill \\ \hfill & =2\pi {\displaystyle {\int }_0^{\pi }}{\displaystyle {\int }_0^{v_0}}{({\stackrel{\sim }{p}}_0+c_1\frac{g(v){sin}^2\theta }{16{\pi }^3})}^{3/4}Q{(v,\theta )}^{-3/4}\hfill \\ \hfill & \times {\mathrm{\Gamma }}^3{v}^{2}sin\theta dvd\theta .\hfill \end{array}$$(29)

We confirm, as expected, that the rest mass contained within a spherical configuration does not evolve with time or the physical radius of the fluid, as the integral is on the v-space, where the solution does not evolve. However, the amount of mass contained inside the sphere depends on the solution. In what follows, we assume that the function Q is constant and sets the scale for the mass. Furthermore, the total mass is sensitive to the choice of the additive constant, ${\stackrel{\sim }{p}}_0$. We made the following choices: for c₁ > 0, we set ${\stackrel{\sim }{p}}_0=0$, whereas for c₁ < 0, we chose the minimum acceptable value of ${\stackrel{\sim }{p}}_0$ so that the minimum of the density is exactly zero, and thus, ${\stackrel{\sim }{p}}_0=-\frac{c_1}{16{\pi }^3}$. The resulting masses from the integration are shown in Fig. 7. We note that for higher absolute values of c₁ the mass is higher. We note that the assumption of an additive constant affects the direct comparison of the mass: the P-type solutions have zero density on their surface, where the N-type solutions have a non-zero density on their surfaces. To produce a reasonable comparison, the P-type solutions have a higher mass load than the N-type solutions for the same density on the surface v = v₀.

Fig. 7. Mass M and work W for the P and N solutions as functions of |c1|. While M and W are plotted on the same axis, the units differ. Moreover, W evolves with time, which is inversely proportional to t.

4.2.3. Thermal pressure work

Following a similar procedure as for the density, we evaluated the work W of the thermal pressure. It can be calculated through the integral of pdV. Upon integration for a sphere of radius R at time t, we obtain the following expression:

$$\begin{array}{cc}\hfill W& =2\pi {\displaystyle {\int }_0^{\pi }}{\displaystyle {\int }_0^R}({\stackrel{\sim }{p}}_0+c_1\frac{g(v){sin}^2\theta }{16{\pi }^3})\hfill \\ \hfill & \times {\left(\frac{\mathrm{\Gamma }v}{r}\right)}^4{r}^{2}sin\theta drd\theta \hfill \\ \hfill & =\frac{2\pi }{\mathit{ct}}{\displaystyle {\int }_0^{\pi }{\int }_0^{v_0}}({\stackrel{\sim }{p}}_0+c_1\frac{g(v){sin}^2\theta }{16{\pi }^3})\hfill \\ \hfill & \times {\mathrm{\Gamma }}^4{v}^{2}sin\theta dvd\theta .\hfill \end{array}$$(30)

We note that W is inversely proportional to t, signifying the expansion of the system. This work is associated with the spatial redistribution of the plasma mass and pressure as the system expands. We further note that the thermal pressure on the surface of the sphere decreases with time, and thus, as the system expands, less work in done against the external medium. The results of this integration for various choices is shown in Fig. 7 with dotted lines.

4.2.4. Temperature

Given the relation between the density and the pressure, we evaluated the temperature of the plasmoid through the relation $T=\frac{mp}{k_B\rho }$, where m is the particle mass, and k_B is the Boltzmann constant. Therefore, the temperature is given by the form:

$$\begin{array}{c}\hfill T=\frac{m}{k_B}\frac{\stackrel{\sim }{p}}{\stackrel{\sim }{\rho }}\frac{\mathrm{\Gamma }v}{r}=\frac{m}{k_B}\frac{\stackrel{\sim }{p}}{\stackrel{\sim }{\rho }}\frac{\mathrm{\Gamma }}{\mathit{ct}}·\end{array}$$(31)

Therefore, the local temperature is inversely proportional to t, with the expansion of the system leading to cooling.

As $\stackrel{\sim }{p}$ and $\stackrel{\sim }{\rho }$ are related to each other through Q, we followed the assumption of a constant value of Q(v, θ) in the evaluation of the temperature, as we already chose in the calculation of mass through the density. Therefore, we expect a temperature profile for the spherical structure that decreases inversely proportionally with time.

For positive values of c₁, the maximum value of the temperature approximately coincides with the maximum of the flux function Ψ, while it is cooler at the boundaries of the sphere and the axis (Fig. 4 panels c and g). For negative values of c₁, the minimum of the temperature occurs close to the maximum of Ψ, while the surroundings are hotter, as illustrated in Fig. 5 panels c and g. For the normalisation values assumed here, we refer to the caption of Table 1, and when we assume a proton-electron plasma, the units of the temperature correspond to 10¹¹ K.

5. Astrophysical implications

The solutions derived here provide a framework for the study of astrophysical systems of relativistically expanding configurations of magnetised plasma. It has been proposed that these structures might arise from magnetars expelling magnetic plasmoids, which are associated with magnetar flares (Thompson & Duncan 2001; Lyutikov 2006), and possibly to FRBs (Bransgrove et al. 2025).

5.1. Magnetar flares

The magnetic fields of magnetars are highly active and evolve due to MHD effects in the crust, such as the Hall drift (Goldreich & Reisenegger 1992; Pons et al. 2009; Gourgouliatos et al. 2016; Kojima et al. 2021; Igoshev et al. 2021), and plastic flows (Levin & Lyutikov 2012; Lyutikov 2015; Lander 2016; Gourgouliatos & Lander 2021), and to effects the core due to ambipolar diffusion and to the evolution of the superconducting magnetic field (Graber et al. 2015; Passamonti et al. 2017a,b; Igoshev & Hollerbach 2023; Skiathas & Gourgouliatos 2024; Igoshev et al. 2025). These effects lead to the formation of multipolar magnetic structures in the lower magnetosphere and may release them to the magnetosphere. When an arcade forms in the magnetosphere, it might become unstable through high twist (Parfrey et al. 2013; Ntotsikas & Gourgouliatos 2025) and can disconnect from the stellar magnetic field (Lander et al. 2024; Bransgrove et al. 2025). Such configurations evolve dynamically and give rise to a magnetar flare Fig. 8. Magnetar flares are rare explosive events that occur in magnetars. A flare can release immense amounts of energy that can reach 10⁴⁶ erg in the form of X-rays, shining within a fraction of a second. Moreover, flare ejecta expel masses reaching 10²⁴ − 10²⁶ g (Gaensler et al. 2005; Granot et al. 2006), which subsequently power scaled-down versions of kilonovae through r-process radioactive decay (Cehula et al. 2024).

Fig. 8. Schematic illustration of plasmoid formation and emission regimes associated with the different families of self-similar solutions. (a) Magnetic reconnection in the magnetar magnetosphere leads to the formation of a magnetically confined plasmoid embedded within the closed-field region. (b) Emission properties of plasmoids corresponding to different plasmoid regimes: P-type solutions (overdense) predominantly emit high-energy radiation (X-rays and γ-rays), whereas Z- and N-type solutions (underdense magnetically dominated configurations) correspond to low-density plasmoids capable of producing coherent radio emission.

In the model developed here, the relative content of the spherical configuration in mass, thermal pressure, magnetic and electric energy, temperature profile, and the maximum Lorentz factor attained are all related through the choices of the parameters c₀ and c₁. When we assume that a blob has a radius of 1 km (10⁵ cm) and the maximum value of the magnetic flux corresponds to 10²⁵ G cm², so that the magnetic field scales with 10¹⁵ G, the magnetic and electric field energies E_mag and E_el range from 10⁴⁴ erg to 10⁴⁶ erg. The pressure work W ranges from 10⁴³ erg to 10⁴⁵ erg, and the temperature is at a level of 10¹¹. We note that the temperature and energetics evolve with time as the plasmoid expands, with the decrease in electromagnetic and thermal energy of the sphere related to the work done against the external medium into which the sphere is expanding. The mass load M ranges from 10²⁴ g to 10²⁵ g and does not evolve with time.

5.1.1. Giant flare energetics

We interpret magnetar giant flares as large-scale relativistic magnetic reconfigurations described by self-similar axisymmetric MHD solutions with finite mass loading. For all explored solution families (Z, P, and N), the energy budget at early times is dominated by the electromagnetic component, while the kinetic and thermal contributions depend on the mass loading and the pressure term, which are quantified through parameter c₁. This is a necessary condition for producing a bright, rapidly varying high-energy spike because an outflow dominated by kinetic or thermal energy would radiate inefficiently at X-ray and γ-ray energies and primarily power delayed shock-driven emission. The electromagnetic dominance of the early-time evolution is therefore consistent with the observation that the prompt hard X-ray or γ-ray spike contains the majority of the radiated energy in magnetar giant flares (e.g. Hurley et al. 2005; Palmer et al. 2005).

An immediate observational implication is that giant flares can exhibit significant diversity in their late-time signatures despite comparable prompt luminosities. Solutions with higher mass loading (P-type) predict enhanced interaction between the ejecta and the surrounding medium, potentially giving rise to detectable radio afterglows or nebular disturbances, as observed following the 2004 flare of SGR 1806–20 (Gaensler et al. 2005). In contrast, low-mass solutions (Z- and N-type) correspond to clean flares, in which little kinetic energy is available to power long-lived afterglow emission, consistent with the weaker radio signatures of SGR 1900+14 (Frail et al. 1999).

5.1.2. Timescales and light-curve morphology

The self-similar expansion fixes the characteristic radial scale of the system to

$$\begin{array}{c}\hfill r(t)=c{v}_{0}t,\end{array}$$(32)

with v₀ of order unity. Dissipation is expected to occur preferentially near the outer boundary of the expanding structure, where strong current sheets and pressure gradients develop. This naturally leads to sub-millisecond rise times for the initial hard X-ray or γ-ray spike, consistent with the observed rise times of Galactic giant flares (Palmer et al. 2005; Hurley et al. 2005).

The model predicts temporally asymmetric light curves, characterised by a very rapid rise phase followed by a longer decay. The rise time reflects localised dissipation at the expanding boundary, whereas the decay traces the subsequent global magnetic reconfiguration of the system. Since these timescales primarily depend on the expansion speed and not on the total energy release, similar rise times are expected across giant flares of different luminosities. While the presence of strong currents and discontinuities lead to dissipation and current-driven instabilities, a detailed study is required to confirm and follow the development of these effects in detail.

5.1.3. Angular structure and anisotropy

The meridional structure of the solution leads to highly anisotropic profiles with magnetic field, pressure, and mass loading different on the equator and at the poles. Consequently, the emitted radiation is expected to be anisotropic, with the apparent energy of a flare depending on the observer’s viewing angle relative to the magnetic field symmetry axis.

This anisotropy provides a natural explanation for the wide dispersion in inferred isotropic-equivalent energies among observed giant flares (Hurley et al. 1999), without requiring intrinsic energy releases spanning several orders of magnitude.

5.1.4. Mass loading and flare diversity

The explored parameter space spans ejecta masses ranging from effectively negligible values (Z-type solutions) to ∼10²⁵–10²⁶ g. This range is consistent with estimates of the ejecta mass inferred from modelling the radio afterglow of SGR 1806–20 (Gaensler et al. 2005); Granot:2006. Low-mass events, associated with Z- and N-type solutions, correspond to predominantly magnetospheric reconfigurations with minimum baryon contamination. These events are expected to produce sharp, luminous high-energy spikes and weak or undetectable afterglows.

In contrast, P-type solutions involve substantial mass loading, potentially including baryonic material originating from the stellar crust or inner magnetosphere (Thompson & Duncan 2001; Bransgrove et al. 2025). Such events are expected to exhibit stronger interaction with the ambient medium, leading to detectable radio afterglows or long-lived nebular emission.

5.1.5. Spectral and polarisation signatures

As the system expands self-similarly, the magnetic field strength decreases rapidly with radius, implying a fast spectral evolution of the emitted radiation. The model predicts an initially hard spectrum during the prompt spike, followed by rapid hard-to-soft evolution on millisecond timescales, in agreement with observed giant flare spectra (Feroci et al. 2001; Olive et al. 2004).

Because the magnetic field remains globally ordered during the early stages of the expansion, the prompt emission is expected to be highly linearly polarised, with a relatively stable polarisation angle during the rise phase. This is consistent with expectations for magnetically dominated emission regions (Lyutikov 2006). Significant depolarisation might only occur at later times, after magnetic reconnection and dissipation disrupted the large-scale field structure.

5.2. Fast radio bursts

Within the present framework, FRBs are interpreted as the coherent radio counterparts of magnetically dominated explosions with low baryon loading Lyubarsky (2021). While the total energy release and large-scale dynamics can resemble those of magnetar giant flares, the conditions required for efficient coherent emission impose strong constraints on the ejecta mass and magnetisation. In particular, FRB production is favoured in solutions where the electromagnetic energy dominates the kinetic and thermal components throughout the early expansion (Lyutikov 2017; Beloborodov 2017).

In heavily mass-loaded explosions, such as P-type solutions, the inertia of the ejecta and the rapid formation of a strong compression shock suppress coherent emission by reducing the magnetisation and increasing the plasma frequency. These events are therefore expected to produce luminous high-energy flares, but weak or absent radio bursts. In contrast, Z- and N-type solutions correspond to minimally mass-loaded outflows in which the plasma remains highly magnetised and optically thin, allowing coherent radio emission to escape. This qualitative distinction is consistent with the observed diversity of magnetar activity, including radio-quiet giant flares and radio-loud but energetically modest bursts (Kaspi & Beloborodov 2017; Mereghetti et al. 2020).

The self-similar expansion naturally leads to the formation of strong current sheets and charge-starved regions near the expanding boundary of the magnetised structure. In these regions, large-amplitude electromagnetic disturbances can be generated on millisecond timescales, providing favourable conditions for coherent curvature radiation or synchrotron maser emission at relativistic shocks (Lyubarsky 2014; Plotnikov & Sironi 2019). The characteristic duration of the resulting radio burst is set by the causal timescale t ∼ r/c, which is consistent with the millisecond durations of observed FRBs.

An important prediction of the model is that FRBs need not be associated with the most energetic magnetar flares. Instead, they preferentially accompany low-mass explosions in which a large fraction of the released energy remains in electromagnetic form and is not degraded into kinetic energy at shocks (Lyutikov 2022). This provides a natural explanation for the observed association of FRBs with magnetars (Popov & Postnov 2013), including the detection of a bright radio emission coincident with an X-ray flare from SGR 1935+2154 (Cameron et al. 2005; Bochenek et al. 2020; CHIME/FRB Collaboration 2020), and the absence of FRBs during the most extreme giant flares.

6. Conclusions

We have presented a class of self-similar axisymmetric MHD solutions describing relativistic magnetically dominated explosions with mass loading, motivated by energetic transient phenomena in magnetar magnetospheres. The solutions encompass a broad range of behaviours, characterised by different degrees of baryon loading and pressure support, while remaining analytically tractable.

A central result of this work is that the early-time evolution of these explosions is generically dominated by electromagnetic energy, with kinetic and thermal components becoming important only as a consequence of mass loading and interaction with the external medium. This regime naturally produces rapid high-energy emission on millisecond timescales, followed by a longer decay associated with magnetic energy depletion and eventual deceleration by a compression shock. The resulting light curves are strongly asymmetric, with rise times controlled primarily by the expansion speed and geometry and not by the total energy release.

Applied to magnetar giant flares, the model reproduces key observational properties, including sub-millisecond rise times, hard prompt spectra, low radiative efficiencies, and the presence or absence of afterglow emission, depending on the ejecta mass. Heavily mass-loaded solutions are associated with efficient energy transfer to shocks and long-lived radio nebulae, while lightly loaded solutions correspond to clean magnetospheric reconfigurations with weak or undetectable afterglows.

In this context, FRBs arise naturally from low-mass highly magnetised explosions in which coherent radio emission can escape before significant baryon loading or shock formation suppresses it. The model therefore predicts that FRBs need not accompany the most energetic magnetar flares, but instead preferentially occur in events where the electromagnetic energy remains dominant throughout the early expansion. This provides a unified explanation for the observed diversity of magnetar activity, including radio-quiet giant flares, radio-loud but energetically modest bursts, and rare events exhibiting both high-energy and radio emission.

The present model is based on a semi-analytical approach, considering ideal MHD, axial symmetry, and self-similarity. These simplifications allow the consistent formation of models and the thorough exploration of the parameter space that illustrates the main qualitative characteristics. The qualitative behaviour of the solutions we presented is consistent with the numerical simulations of magnetically driven explosions in mass-loaded environments by Barkov et al. (2022). In particular, both approaches show that increasing the toroidal-to-poloidal magnetic field ratio leads to more compact configurations and a reduction in the maximum expansion velocity. Our semi-analytical solutions complement these numerical results by providing self-consistent equilibrium plasmoid structures that explicitly relate magnetic twist, plasma pressure, and inertia to the global dynamics of relativistic magnetic explosions. However, as these solutions are force-free while in our solutions, pressure and mass are central, a direct comparison is not yet possible.

The key effects that are associated with observable properties of magnetar flares and FRBs require the relaxation of the assumed idealisations we used here. Namely, the expansion of a spherical magnetised plasmoid, as it propagates away from the neutron star, is not guaranteed to remain self-similar. Moreover, these spherical large-scale plasmoids can deviate from the spherical shape and axial symmetry. Finally, resistive instabilities and the possible formation of shocks due to the interaction of the mass with the external medium shape the radiative properties of the plasmoid. These deviations can be traced with a numerical simulation that will allow us to determine the limitations of the approach made here and their consequences for the observables.

Overall, the results presented here demonstrate that a single physically motivated class of magnetically dominated mass-loaded explosions can account for a wide range of observed transient phenomena associated with magnetars. By explicitly linking the dynamics to a small number of physical parameters, the model provides a coherent framework for interpreting current observations and for guiding future multi-wavelength and polarisation measurements of magnetar-driven transients.

Data availability

The datasets generated during the current study are available from the author upon reasonable request.

Acknowledgments

The author acknowledges funding from grant FK 81641 “Theoretical and Computational Astrophysics”, ELKE. This work was supported by computational time granted by the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility – ARIS – under project ID pr017008/simnstar2.

References

All Tables

All Figures

	Fig. 1. Z-family solutions for g(v) for various values of c0 when c1 = 0.
In the text

	Fig. 2. Characteristic P-family solutions of g(v) for various values of c1 while keeping c0 = 1.
In the text

	Fig. 3. N-family of solutions showing g(v) for various values of c0 and negative c1, so that the derivative is zero at the maximum velocity.
In the text

Fig. 4. Contour plots of the solutions with c0 = 1 and c1 = 1 (top row) and c0 = 1 and c1 = 10 (bottom row). The black contours correspond to the poloidal magnetic field lines, and we plot in colour the Bϕ component of the magnetic field (panels a and e), the rest-mass density ρ (panels b and f), the temperature T (panels c and g), and the electric charge density j0 (panels d and h). The solution remains invariant in the v-space. We present a snapshot at some time t, so that r = 1 corresponds to v = 1.

In the text

Fig. 5. Contour plots of the solutions with c0 = 10 and c1 = −1.448 (top row) and c0 = 100 and c1 = −211.3 (bottom row). The black contours correspond to the poloidal magnetic field lines, and we plot in colour the Bϕ component of the magnetic field (panels a and e), the rest-mass density ρ (panels b and f), the temperature T (panels c and g), and the electric charge density j0 (panels d and h). As in the plots of Fig. 5, the snapshot is taken at some time t, so that r = 1 corresponds to v = 1.

In the text

	Fig. 6. Electric and magnetic energy for the Z and N solutions as functions of c0.
In the text

	Fig. 7. Mass M and work W for the P and N solutions as functions of |c1|. While M and W are plotted on the same axis, the units differ. Moreover, W evolves with time, which is inversely proportional to t.
In the text

Fig. 8. Schematic illustration of plasmoid formation and emission regimes associated with the different families of self-similar solutions. (a) Magnetic reconnection in the magnetar magnetosphere leads to the formation of a magnetically confined plasmoid embedded within the closed-field region. (b) Emission properties of plasmoids corresponding to different plasmoid regimes: P-type solutions (overdense) predominantly emit high-energy radiation (X-rays and γ-rays), whereas Z- and N-type solutions (underdense magnetically dominated configurations) correspond to low-density plasmoids capable of producing coherent radio emission.

In the text