© The Authors 2025

1. Introduction

Extremely luminous, highly variable, non-thermal, high-energy cosmic sources of radiation are often interpreted as originating from powerful collimated outflows known as relativistic jets (Madejski & Sikora 2016; Blandford et al. 2019). In many supermassive active galactic nuclei (AGNs), narrow jets and their relativistic motions can be observed directly (e.g. Jorstad et al. 2001; Lister et al. 2016). In the subclass of blazars, a small viewing angle allows for a huge (∼10⁵) relativistic boost of non-thermal radiation produced by particles accelerated during dissipative processes in the approaching jet (Urry & Padovani 1995). The most luminous blazars, belonging to the subclass of flat-spectrum radio quasars (FSRQs), produce giga-electronvolt γ-ray flares with apparent luminosities of L_fl, FSRQ ∼ 10⁵⁰ erg s⁻¹ and variability timescales of t_fl, FSRQ ∼ 10 min (e.g. Abdo et al. 2011; Ackermann et al. 2016; Shukla & Mannheim 2020)¹. Another subclass of blazars, high-energy peaked BL Lac objects (HBLs), produce tera-electronvolt γ-ray flares with apparent luminosities of L_fl, HBL ∼ 10⁴⁷ erg s⁻¹ and variability timescales of t_fl, HBL ∼ 3 min (Aharonian et al. 2007). Even after accounting for the relativistic boost, the energetic requirements for the blazar zones can be very demanding in terms of local energy density and radiative efficiency (e.g. Nalewajko et al. 2012; Ackermann et al. 2016; Appendix A).

The toroidal magnetic field is an essential ingredient of relativistic jets launched by magnetised accretion onto spinning black holes ((BHs); e.g. Begelman et al. 1984; Davis & Tchekhovskoy 2020). It is the main carrier of electromagnetic momentum (Poynting) flux (the initial form of jet power) (Blandford 1976). Its pressure acting along the velocity streamlines is the primary agent of acceleration to relativistic speeds (e.g. Camenzind 1989; Komissarov et al. 2007). In jet regions that expand laterally, the toroidal field decreases slower than the poloidal field, likely becoming the dominant component at distances of ∼0.1 − 1 pc typically inferred for emission of blazar flares (Nalewajko et al. 2014).

Toroidal magnetic field is also an essential ingredient of plasmoids (flux ropes) produced during large-scale magnetic reconnection due to the tearing instability of thin current sheets (Furth et al. 1963; Loureiro et al. 2007). Relativistic reconnection is one of the primary candidate mechanisms of energy dissipation considered for relativistic jets to explain the non-thermal particle acceleration and broad-band radiation in blazars (Sironi et al. 2015). In particular, rapid gamma-ray flares of blazars have been modelled in terms of reconnection-driven relativistic Alfvenic outflows – minijets (Giannios et al. 2009; Nalewajko et al. 2011)², slower and denser plasmoids (Giannios 2013; Petropoulou et al. 2016; Christie et al. 2019), head-on plasmoid mergers (Nalewajko et al. 2015; Zhang et al. 2024), or tail-on plasmoid mergers (Ortuño-Macías & Nalewajko 2020). Although produced at vastly smaller scales than relativistic jets, reconnection plasmoids can be described by essentially the same physical principles.

A fundamental property of the magnetic field is that the associated stress tensor has an anisotropic component called magnetic tension (e.g. Parker 2007). For any loop of the toroidal field, tension induces a force acting towards the symmetry axis with the effect of compressing the loop interior. In the absence of other forces, an axisymmetric force-free equilibrium radial profile of the toroidal magnetic field is B_ϕ(r)∝r⁻¹, implying a strongly concentrated energy density profile of u_B(r) = B_ϕ²(r)/8π ∝ r⁻². If such a profile could be maintained between widely separated scales, r₁ ≪ r₂, energy density at the inner radius would be greatly enhanced: u_B, 1/u_B, 2 ∼ (r₂/r₁)².

In relativistic jets, the separation between the macrophysical scale (jet width) and the microphysical scale (particle gyroradius) can in principle reach ∼10¹⁰ (see Appendix B). The potential for energy density enhancement due to magnetic tension is thus enormous. On the other hand, scale separations that can be realised in numerical simulations are ∼10⁴. Moreover, in the presence of other forces acting away from the symmetry axis, such as pressure gradients of plasma or the axial magnetic field (e.g. Ortuño-Macías et al. 2022), the compression of energy density is reduced and divided between the magnetic field and plasma in various proportions.

In Section 2, we consider the case of plasmoids emerging from reconnecting current layers in relativistically magnetised plasma. The internal structure of plasmoids can be investigated by kinetic numerical simulations, in particular using the particle-in-cell (PIC) method (e.g. Sironi et al. 2016; Petropoulou & Sironi 2018; Ortuño-Macías & Nalewajko 2020; Hakobyan et al. 2021; Schoeffler et al. 2023; Hakobyan et al. 2023; Chernoglazov et al. 2023). Previous numerical studies of reconnection plasmoids (summarised in Section 2.1) adopted various configurations of PIC simulations, such as 2D or 3D, in periodic or open boundaries, various background magnetisations, relativistically cold or hot plasma, with or without a guide field, and with or without radiative cooling.

The main result of this work is the analysis of energy density compression in plasmoids produced in a series of PIC simulations of relativistic reconnection, presented in Section 2.3. As is described in Section 2.2, these simulations were initiated from hot Harris-type current layers in periodic boundaries, leading to the formation of relaxed ‘monster’ plasmoids. We checked the effects of the guide magnetic field and the third dimension. While we do not include radiative cooling in those simulations, in Section 2.4 we compare our results with previous simulations presented in Ortuño-Macías & Nalewajko (2020), which included synchrotron radiative cooling and open boundaries.

In Section 3, we consider the lateral structure of relativistic jets and the role of the toroidal magnetic field. The context for this is the possibility of producing luminous and rapid flares observed in blazars powered by plasmoids produced by relativistic reconnection. The key question is whether reconnection can produce plasmoids in the right region across the jet, where the compression of plasmoid cores can be multiplied by the compression of the jet.

Previous studies summarised in Section 3 indicate that tension of the toroidal field competes with other forces (e.g. centrifugal force); while it does not dominate jet collimation, it contributes to a gradual differentiation of the jet structure into a quasi-cylindrical inner core and a paraboloidal outer layer. One of the key effects of toroidal field tension is poloidal field bunching (Tchekhovskoy et al. 2009), which effectively focuses the inner loops of toroidal field (and part of the Poynting flux) around the cylindrical jet core, in an adjacent quasi-cylindrical layer that we call the bunched spine. The focusing of the toroidal field is expected to destabilise the jet core to current-driven m = 1 kink modes, seeding reconnecting current layers and plasmoids around the jet core, where compression of plasmoid cores and high-σ particle acceleration can be combined most effectively with compression of the inner jet and beamed relativistically to produce rapid and bright gamma-ray flares of blazars.

Section 4 presents the discussion and conclusions. We propose a specific scenario for the production of luminous radiation flares from relativistic jets by combining energy density enhancements due to tension of the toroidal magnetic field at the scales of reconnection plasmoids (investigated by means of kinetic PIC simulations) and the global lateral jet structure. High-resolution observations of a central ridge along the nearby misaligned jet of M87 may constrain the proposed jet structure.

2. Reconnection plasmoids

2.1. Previous studies

Physical properties of plasmoids produced during relativistic magnetic reconnection have been studied by means of kinetic (PIC) numerical simulations. Here, we report a selection of previous studies initiated from Harris current layers. These studies adopted various assumptions: 2D or 3D, periodic or open boundaries, cold or hot background plasma, with or without a guide field, and with or without synchrotron cooling.

The guide field in the problem of magnetic reconnection is the field component that is not reversed in the initial configuration. As the anti-parallel field component becomes the toroidal field of individual plasmoids, the guide field becomes the axial field concentrated in the plasmoid cores. In the absence of the initial guide field, the plasmoid cores are supported by plasma pressure, resulting in the Z-pinch configuration.

Sironi et al. (2016) measured plasmoid profiles (functions of the radius, r) in 2D simulations (2L domain with L ∼ 3600d_e with d_e = c/ω_p the skin depth, ω_p = (4πe²n₀/m)^1/2 the plasma frequency³; using open boundaries) with relativistically cold background plasma (Θ₀ = k_BT₀/mc² = 10⁻⁴) without radiative cooling. They found plasma density to scale like n ∝ r⁻¹, the magnetic energy fraction ϵ_B = u_B/(n₀mc²)∝r^−1.2 (with n₀ the background plasma density), and the kinetic energy fraction ϵ_kin = (⟨γ⟩ − 1)n/n₀ ∝ r^−1.4. The largest plasmoids of half-width w ∼ L/20 ∼ 180d_e reached core densities of n_c/n₀ ∼ 300(σ₀/10)^1/2, magnetic energy fractions of ϵ_B, c ∼ 800(σ₀/10)^3/2, and kinetic energy fractions of ϵ_kin, c ∼ 4000(σ₀/10)^3/2.

Petropoulou & Sironi (2018) analysed comparable simulations (2L domain with L ∼ 1000r_L with a nominal ‘hot’ gyroradius r_L = σ₀^1/2d_e ∼ 3d_e for magnetisation σ₀ = 10; and periodic boundaries) and obtained even larger plasmoids with half-widths of w ∼ 300d_e, and reported for them magnetic energy fractions of ϵ_B ∼ 20. In this and subsequent works of that group, the population of initial drifting particles was excluded from the analysis.

Ortuño-Macías & Nalewajko (2020) investigated evolution of plasmoids in 2D simulations (L ∼ 1500ρ₀ domain with open boundaries) in relativistically hot plasma (Θ₀ = k_BT₀/mc² ∼ 10⁶, hence the background gyroradius ρ₀ = λ_D/σ₀^1/2 is related to the Debye length, λ_D = Θ₀^1/2d_e, with magnetisation of σ₀ = B₀²/(4πΘ₀n₀mc²)∼10), and with synchrotron cooling. In case of slow cooling (Θ₀ = 2 × 10⁵), they obtained plasmoids growing to half-widths of w ∼ 100ρ₀ ∼ 30λ_D with core densities of n_c/n₀ ∼ 60 and magnetic field strengths of B_c/B₀ ∼ 4. In case of fast cooling (Θ₀ = 1.25 × 10⁶), the plasmoids were smaller (w ∼ 60ρ₀ ∼ 20λ_D) with significantly denser cores (n_c/n₀ ∼ 300) and slightly stronger magnetic fields (B_c/B₀ ∼ 6).

Hakobyan et al. (2021) investigated the structure of plasmoids in 2D PIC simulations with background hot magnetisation, σ₀, up to 100, reaching w ∼ 30r_L ∼ 300d_e (with r_L = σ₀^1/2d_e the hot gyroradius, and d_e the skin depth). Considering the force balance between magnetic tension and plasma pressure, they characterised the radial profiles of plasmoids with power laws of ρ ∝ r⁻¹ and B ∝ r^−2/3, noting also strong stratification of the mean particle energy ⟨γ⟩ (resulting in an even steeper profile of plasma energy density). For σ₀ = 100, they achieved an enhancement in the plasma density, ρ, by a factor of ∼30 and in the magnetic field strength, B, by a factor of ∼3.

A comprehensive investigation of compression of plasmoid cores in 3D PIC simulations of relativistic reconnection was reported by Schoeffler et al. (2023). They considered background plasma with hot magnetisation of σ₀ ∼ 25 and a mildly relativistic temperature (Θ₀ = 4), taking into account synchrotron cooling (regulated by the magnetic field strength, B₀) and the guide magnetic field, B_g. Compression was analysed in a parameter space of plasma density n/n₀ and the magnetic field strength B/B₀, by which they showed that the enhancement of n does not always correlate with the enhancement of B, which has complex implications for the luminosity of associated radiation signals. In their reference case (3D, radiative with B₀ = 2 × 10¹¹ G, moderate guide field B_g/B₀ = 0.4), the limit on the compression of both parameters was identified at n/n₀ ∼ 30 (comparable with the initial density of the drifting plasma, n_d/n₀) and B/B₀ ∼ 2. That limit was not significantly different in the non-radiative case (B₀ = 2 × 10⁸ G). It depended somewhat (but not monotonically) on the guide field, B_g, and the aspect ratio, L_z/L_x, of the 3D domain. The density compression increased with increasing magnetisation, σ₀, in proportion to increasing n_d/n₀.

Effects of strong synchrotron cooling on reconnection plasmoids have recently been investigated by Hakobyan et al. (2023) in 2D and Chernoglazov et al. (2023) in 3D. Strong synchrotron cooling reduces plasma pressure, causing further plasmoid contraction and plasma compression. It also severely limits the maximum energy that can be reached by the particles. Nevertheless, plasmoid cores are generally the brightest sources of radiation within sites of relativistic reconnection.

2.2. New simulations

We used a modified version of the public PIC code Zeltron (Cerutti et al. 2013) to perform effectively 2D and 3D simulations of relativistic magnetic reconnection with electron-positron pair plasma without a radiation reaction in Cartesian coordinates (x, y, z) with periodic boundaries. The physical domain was 0 ≤ x ≤ L_x, 0 ≤ y ≤ L_y, 0 ≤ z ≤ L_z, with L_z < L_y = L_x. The numerical resolution was dx = L_x/N_x = L_y/N_y = L_z/N_z = ρ₀/2.56 with the nominal gyroradius ρ₀ = Θ_em_ec²/eB₀, with m_e the electron mass, c the speed of light, and e the electric charge; Θ_e = k_BT_e/m_ec² = 1 the relativistic temperature of the initial Maxwell-Jüttner electron distribution with k_B the Boltzmann constant; with B₀ = 1 G the nominal magnetic field strength. We set two Harris layers: B_x(y < y_1/2) =  − B₀tanh[(y − y_1/4)/δ] and B_x(y ≥ y_1/2) = B₀tanh[(y − y_3/4)/δ] with y_1/2 = L_y/2, y_1/4 = L_y/4, y_3/4 = 3L_y/4. The layer thickness δ = 2ρ₀/u_dr was determined by the drift velocity, β_dr = 0.3 via u_dr ≡ γ_drβ_dr ≡ β_dr/(1 − β_dr²)^1/2 (Kirk & Skjæraasen 2003). The layers were supported by the pressure of a drifting particle population of density n_dr(y < y_1/2) = n_dr, 0cosh⁻²[(y − y_1/4)/δ] (and likewise for y > y_1/2, centred at y_3/4) with n_dr, 0 = γ_drB₀²/(8πΘ_em_ec²). A static background particle population of uniform density n_bg was added to achieve the desired background magnetisation of σ₀ = B₀²/(4πn_bgΘ_em_ec²). In some cases, we added a guide magnetic field component of uniform strength B_z. We did not use any initial perturbation.

The list of performed simulations and their key parameters is presented in Table 1. The cases for N_z = 2 are essentially 2D. Simulations were performed for durations of at least 3L_x/c, over which they developed plasmoid chains that merged into one large relaxed plasmoid per layer. The magnetic energy dissipation efficiency, ℰ_B/ℰ_B, ini, was at least 50% in all cases without a guide field, almost 30% by 5L_x/c for B_z/B₀ = 0.5, and roughly 10% by 4.5L_x/c for B_z/B₀ = 1.

In simulations with a guide field, the total electric energy, ℰ_E, decays at a regular but slow rate of ∼20% per L/c (reflecting elastic oscillations of the plasmoid core). They were thus extended to durations of t ∼ 6L/c, by which point ℰ_E ∼ 10⁻² ℰ_B,0. In contrast, simulations without a guide field (for any considered L/ρ₀ or σ₀) achieved ℰ_E ∼ 10^−2.5 ℰ_B,0 by t ≲ 4L/c.

We calculated the volume distributions, F_a(μ) = dF_a/dμ, of the parameter μ_a = log₁₀(u_a/u_B, 0), with u_a the energy density of either the magnetic field, u_B = B²/8π (u_B, 0 = B₀²/8π), or the plasma, u_pl = ⟨γ⟩ nm_ec² (including both electrons and positrons). The distributions are normalised to unity, so that dμ_a∑_iF_a(μ_a, i) = 1.

2.3. New results

Figure 1 compares the functions, f_a(μ_a) = log₁₀(F_a(μ_a)), averaged over the entire durations of each simulation. The peaks of the F_B distributions correspond to the initial background magnetic field strength, μ_B, ini = log₁₀(1 + B_z²/B₀²). The peaks of the F_pl distributions correspond to the initial background plasma population (〈γ〉 ≃ 3.37 for Θ_e = 1): μ_pl, bg ≃ log₁₀(6.7/σ₀). The initial drifting population contributes up to μ_pl, dr ≃ log₁₀(3.4γ_dr²)≃0.57.

Fig. 1. Logarithms, f = log10(F), of volume distributions, F(μ) = dF/dμ, over the argument μ = log10(u/uB, 0), with u the energy density: of magnetic fields, uB = B2/8π (left panel), and of the plasma, upl = ⟨γ⟩nmc2 (right panel). Functions f(μ) were averaged over the duration of each simulation.

The high-energy-density tails of the distributions have a complex structure that can be evaluated at different levels. For the reference case of σ₀ = 10 without a guide field, comparing the results for simulations L1800_σ10 and L3600_σ10, convergence for F_B was achieved to the level of 10⁻⁴, and for F_pl even down to 10⁻⁷. The case of σ₀ = 20 is converged in F_pl down to 10⁻⁵.

2.3.1. Relaxed monster plasmoids

Monster plasmoids are very large (≳0.2L) and slow (v ≪ v_A) plasmoids emerging from hierarchical mergers of plasmoids chains. In simulations performed with periodic boundary conditions, they dominate the final states into which reconnecting plasma relaxes.

The parameter with the most qualitative impact on the structure of monster plasmoids is the strength of the guide field, B_g, relative to B₀. Thus, Figure 2 presents energy density maps for relaxed monster plasmoids in the final states of three simulations for B_g/B₀ = 0, 0.5, 1.

Fig. 2. Maps of log10(u/uB, 0)(x, y) for magnetic energy density, uB = B2/8π (upper panels), and plasma energy density, upl = ⟨γ⟩nmc2 (middle panels), for relaxed monster plasmoids at the end of each simulation. In the lower panels, we compare 1D energy density profiles of (u/uB, 0)(x, y0) measured along the strip indicated in the above maps by the dashed grey lines. A common colour scale for all maps is referenced along the left axes in the lower panels. Here, we present the effect of a guide field for σ0 = 10 and L/ρ0 = 1800. From the left, the columns show simulations: (1) L1800_σ10, (2) L1800_σ10_Bg05, and (3) L1800_σ10_Bg1.

In the case of no guide field, B_g = 0 (L1800_σ10), the plasma energy density has a smooth structure, while the magnetic energy density ripples, especially along the x axis. The plasmoid is dominated by the plasma energy density with u_pl/u_B ∼ 10, which peaks at μ_pl ≃ 0.9 across the central unmagnetised core for 0.45 < x/L < 0.54; in the f_pl distribution (Figure 1, right panel, solid red line), this corresponds to a break at the −1 level. The magnetic energy density peaks at μ_B ≃ −0.05 (just below the initial peak of f_B; Figure 1, left panel, solid red line) along a ring structure crossing the y/L = 0.75 plane at x/L ≃ 0.39 and 0.62. Hence, a relaxed monster plasmoid achieves only a minor compression of u_pl, and no compression of u_B.

In the case of a strong guide field, B_g = B₀ (L1800_σ10_Bg1), the plasmoid is characterised by a roughly uniform magnetic energy density with μ_B ∼ 0.3 (corresponding to the peak of f_B; Figure 1, left panel, solid black line). The plasma energy density shows complex substructures, with high-density arcs tracing the magnetic field lines. Only the plasmoid core (0.4 < x/L < 0.59) is dominated by the plasma with μ_pl ≃ 0.6 (corresponding to a bump in f_pl at the −1.1 level; Figure 1, right panel, solid black line), and hence u_pl/u_B ∼ 2, the plasmoid layer is magnetically dominated with u_B/u_pl ∼ 1.8.

The case of the moderate guide field, B_g = B₀/2 (L1800_σ10_Bg05), shows an intermediate structure, which is qualitatively closer to the case of the strong guide field.

Figure 3 compares the y profiles of magnetic and plasma energy densities measured along the y axis of relaxed monster plasmoids. The first (from the left) panel compares such profiles for three simulation domain sizes, L/ρ₀ = 900, 1800, 3600 (with no guide field, σ₀ = 10). The general characteristic of these profiles (specific to the cases without a guide field) is a very sharp transition between the plasmoid core (0.455 < y/L < 0.54, u_B < u_B, 0/100, uniform u_pl) and the plasmoid layer, supported by a sharp ring structure of electric current density j_z. Across the plasmoid core, the plasma energy density is roughly uniform at the level of u_pl ≃ 8.5u_B, 0. It clearly dominates the magnetic energy density in the plasmoid layer, peaking at the level of u_B, peak ≃ u_B, 0, while remaining in pressure balance – a key factor is that our choice for the relativistic plasma temperature, Θ_e = 1, implies a relation⁴ between the plasma energy density and the pressure, u_pl/P_pl ≃ 3.4. Across the plasmoid layer, the plasma energy density decreases exponentially like u_pl ∝ 10^{−|y − y₀|/(L/7)}, bringing the magnetic field towards the force-free equilibrium u_B ∝ 1/|y − y₀|² for |y − y₀|/L > 0.28. With increasing L/ρ₀, the plasmoid core becomes smaller in units of L, and also larger in units of ρ₀. It scales approximately as ∝(L/ρ₀)^1/2. However, this has little effect on the profiles of u_pl, both in the core and across the layer.

Fig. 3. Profiles of magnetic (green lines) and plasma (red lines) energy densities across (i.e. along the y coordinate) relaxed monster plasmoids. Panels from the left compare: (1) different sizes, L/ρ0, of simulation domain (2D with σ0 = 10, Bz = 0); (2) different magnetisations, σ0 (2D with L/ρ0 = 1800, Bz = 0); (3) 3D and 2D domains (with L/ρ0 = 900, σ0 = 10, Bz = 0); (4) different guide field strengths, Bg/B0 (2D with σ0 = 10, L/ρ0 = 1800).

The second panel of Figure 3 compares the profiles of u_B and u_pl for three values of initial background magnetisation σ₀ = 10, 20, 40 (with no guide field, L/ρ₀ = 1800). The effect of σ₀ on the core and the layer of the plasmoid is rather minor (a clear effect of σ₀ can be seen for y/L < 0.25 and y/L > 0.75, where background plasma still dominates). For σ₀ = 10, the magnetic boundary between the core and the layer is less sharp, which is balanced by a slightly lower u_pl across the core.

The third panel of Figure 3 compares the profiles of u_B and u_pl between an essentially 2D simulation (with N_z = 2) and a 3D simulation with L_z/ρ₀ ≃ 42 (with N_z = 108), other key parameters being exactly the same (L_x = L_y = 900ρ₀, σ₀ = 10, B_g = 0). Profiles in the 3D case were calculated using different statistics along the z coordinate: minimum, mean, maximum. The resulting monster plasmoid shows a broader core with gradual boundaries – a plateau of u_pl ∼ 6u_B, 0 with a total width of Δy ∼ 0.2L, and a base of u_B ∼ 0.07u_B, 0 with a total width of Δy ∼ 0.1L. The plasmoid layer is consistent between the 3D and 2D cases.

The last panel of Figure 3 compares the profiles of u_B and u_pl for three values of the guide field strength, B_g/B₀ = 0, 0.5, 1 (with L/ρ₀ = 1800, σ₀ = 10). In the presence of the guide field, both profiles appear less regular, which reflects a much slower relaxation of these plasma configurations. Nevertheless, one can notice that the profiles of u_B are flatter, with plasmoid cores filled with uniform u_B at the levels of ∼u_B, 0 for B_g = B₀/2, and ∼2u_B, 0 for B_g = B₀. Still, the plasmoid cores are dominated by u_pl, in the B_g = B₀ case by a factor of ∼2. The plasmoid layer is roughly in equipartition for B_g = B₀/2 and magnetically dominated for B_g = B₀.

2.3.2. Plasmoid mergers

Figure 4 shows energy density maps for merging plasmoids in the same three simulations as in Figure 2.

Fig. 4. Same as Figure 2, but for plasmoid mergers that maximise the plasma energy density, upl, for the same set of simulations: (1) L1800_σ10, (2) L1800_σ10_Bg05, and (3) L1800_σ10_Bg1.

In the reference case L1800_σ10, the entire system of merging plasmoids including the secondary current layer is dominated by the plasma energy density reaching μ_pl ∼ 1.3 (corresponding to f_pl ≃ −3 in the high-density tail; Figure 1, right panel, solid red line). The magnetic energy density reaches μ_B ∼ 0.6 (corresponding to f_B ≃ −3.7; Figure 1, left panel, solid red line) in the vicinity of the current layer (possibly related to the formation of a minor secondary plasmoid). The electric energy density reaches μ_E ∼ −0.9 in the same region. Compared to the relaxed monster plasmoid, compression of u_pl is stronger by a factor of ∼2.5, but compression of u_B is stronger by a factor of ∼4.5.

In the case of a strong guide field, B_g = B₀ (L1800_σ10_Bg1), we show a merger that forms a secondary current layer where the plasma energy density reaches μ_pl ≃ 0.95 (corresponding to f_pl ≃ −4; Figure 1, right panel, solid black line; higher by a factor of ≃2.2 than in the monster plasmoid) and the magnetic energy density reaches μ_B ≃ 0.55 (corresponding to f_B ≃ −2; Figure 1, left panel, solid black line; higher by factor ≃1.8 than in the monster plasmoid).

The case of the moderate guide field, B_g = B₀/2 (L1800_σ10_Bg05), shows plasma energy densities comparable to the case of no guide field.

2.3.3. Structures in the volume distributions

With this analysis, we can identify the main structures in the volume distributions f_pl(μ_pl) and f_B(μ_B) presented in Figure 1.

The plateau in the f_pl distribution corresponds to relaxed plasmoids. In the absence of a guide field, it extends to μ_pl ≃ 1 or u_pl/u_B, 0 ∼ 10, only weakly depending on the background magnetisation, σ₀. Compared with the initial drifting particle population, this means a modest compression by a factor of ∼2.7. A guide field of B_z = B₀ reduces this compression by a factor of ∼1.8. Such a plateau is absent in the f_B distribution – relaxed plasmoids do not amplify u_B.

The steep high-density tail of the f_pl and f_B distributions, extending to the level of −4, corresponds to plasmoid mergers. Considering the tail of the f_pl distribution at the level of −4, in our reference case (σ₀ = 10 in the absence of a guide field) it reaches μ_pl ≃ 1.5. Compared with the initial drifting particle population, this means compression by a factor of ∼8.5. A higher magnetisation of σ₀ = 40 enhances this compression by a factor of ∼1.6, while a guide field of B_z = B₀ reduces it by a factor of ∼3.5. Considering the tail of f_B at the same level of −4, for σ₀ = 10 in the absence of a guide field it reaches μ_B ≃ 0.75 (compression of u_B, 0 by factor ∼5.5), for σ₀ = 40 compression increases by a factor of ∼1.4, and for B_z = B₀ it decreases by a factor of ∼1.25.

The second high-density tail of the f_pl distribution, extending below the −4 level, can be identified as secondary plasmoids. Considering the level of −6, in the reference case it reaches at least μ_pl ≃ 2.1 (a factor ∼4 higher than for primary merging plasmoids). In the cases of σ₀ = 20, 40, it reaches μ_pl ≃ 2.4 (factor ∼2 higher), and in the guide field case B_z = B₀ it reaches μ_pl ≃ 1.35 (a factor of ∼5.5 lower).

2.4. Comparison with open-boundary radiative simulations

Figure 5 presents volume distributions of magnetic and plasma energy densities for simulations of relativistic reconnection initiated from a single Harris current layer using open boundary conditions and including synchrotron cooling, first presented in Ortuño-Macías & Nalewajko (2020). Open boundary conditions were applied across both ends of the current layer (x = 0 and x = L_x), which allowed for the free escape of plasmoids and prevented the formation of a single slow (due to the cancellation of x momentum) monster plasmoid. The background plasma was extended to a larger volume (L_z = 4L_x), largely preserving it in pristine condition throughout the simulation time. We thus find a stronger dominance of background plasma at u_B, u_pl < u_B, 0 (x_B, x_pl < 0). In the F_pl distribution, a bump for 0 < x_pl < 1.1 corresponding to large plasmoids is still prominent. A major difference from periodic boundary simulations is the presence of extended power law tails in both F_B and F_pl distributions. In the high-magnetisation case (σ₀ = 50), the tail of F_B extends over 0.3 < x_B < 2.8, and the tail of F_pl extends over 1.6 < x_pl < 3.7. The tails are less extended in the low-magnetisation cases (σ₀ = 10), and especially the cut-off of F_B is sensitive to radiative cooling efficiency, shifting to x_B ≃ 1.2 for Θ = 2 × 10⁵ (slow cooling). The power law index is approximately p ≃ −5/3 for both F_B and F_pl, although for σ₀ = 50 it hardens towards the highest x_B, x_pl to p ≃ −1, and for σ₀ = 10 it softens to p ≃ −2.

Fig. 5. Same as Figure 1, but for the open-boundary simulations with synchrotron cooling first presented in Ortuño-Macías & Nalewajko (2020). The black dashed lines indicate power laws of index −5/3.

3. Relativistic jets

Relativistic jets are launched from rotating magnetospheres with large fluxes in the poloidal magnetic field (Blandford 1976). The twist of the poloidal magnetic field generates a toroidal magnetic component and a perpendicular electric field, combining into outflowing electromagnetic momentum – the Poynting flux. In the process of Blandford & Znajek (1977), an electromagnetic version of the Penrose process (Lasota et al. 2014), torque is exerted on the BH, converting its rotational energy to the energy of the toroidal magnetic field and the associated electric field.

A magnetosphere connected to the BH horizon achieves highly relativistic magnetisation (σ = b²/w ≫ 1) and force-free conditions (negligible inertia, ρc², pressure, P, enthalpy, w; sufficient charge density, ρ_e, current density, j) by unloading its plasma onto the BH. This results in a sharply edged cavity of low plasma density, which (in combination with relatively uniform magnetic enthalpy density, b²) corresponds with the σ > 1 magnetosphere (e.g. Nalewajko et al. 2024).

An outflowing BH magnetosphere is gradually collimated by external plasma pressure and accelerated by the pressure of the toroidal magnetic field (Camenzind 1989), converting its electromagnetic energy into kinetic plasma energy. According to the conservation of the Michel parameter (specific energy) μ = Γ(1 + σ_ϕ′) (Michel 1969), magnetisation σ_ϕ′=B_ϕ^′2/(4πw′), based on the co-moving toroidal magnetic field B_ϕ′ and relativistic enthalpy density w′, converts into the bulk Lorentz factor Γ. The efficiency of such acceleration depends critically on the geometry of poloidal magnetic field lines – it is particularly low in the spherical case (Begelman & Li 1994), but can be very high in the parabolic case (Beskin & Nokhrina 2006). In this acceleration-collimation zone (ACZ), extending to distances of r (10³ − 10⁵)M_BH (Komissarov et al. 2007, 2009), the jet is characterised by high magnetisation and largely ordered magnetic fields.

The boundary of a relativistic jet can be defined by the trans-relativistic four-velocity u = Γv/c = 1, which corresponds to the Lorentz factor $\mathrm{\Gamma }={(u^2+1)}^{1/2}=\sqrt{2}$. Within the jet is the density cavity, which in the ACZ can be identified by σ = 1 (beyond the ACZ, σ < 1 across the entire jet; whether the cavity remains well defined depends on the efficiency of plasma mixing processes). Recent numerical works (e.g. Davelaar et al. 2018; Salas et al. 2024) describe the inner σ > 1 region as the jet spine, and the outer σ < 1 region as the jet sheath; the original spine-sheath structure was proposed by Ghisellini et al. (2005).

To describe the jet structure, we use spherical coordinates (r, θ, ϕ) and a cylindrical radius of R = r sin θ. Gradual collimation in the ACZ means that the shape of the jet spine, θ_sp(r), can be described as paraboloidal. Beyond the ACZ, it transitions to conical; this has been confirmed by direct measurements of the innermost jets in nearby AGNs with the very-long-baseline interferometry (VLBI) technique (Asada & Nakamura 2012; Kovalev et al. 2020). A flat or slowly declining profile of external (sheath) pressure may even lead to recollimation of the jet spine (Bromberg & Levinson 2007; Porth & Komissarov 2015).

3.1. Lateral structure of jet magnetic fields

Relativistic jets can be robustly described by stationary and axisymmetric models. In such models, one considers two main components of the magnetic and velocity fields that are closely interrelated. The poloidal plasma velocity is parallel to the poloidal magnetic field; hence, poloidal plasma streamlines θ_s(r) coincide with the poloidal field lines. The relation between the poloidal and toroidal magnetic fields is governed by the trans-field momentum (Grad-Shafranov) equation (Beskin 2010). Detailed solutions presented in the literature (e.g. Appl & Camenzind 1993; Lyubarsky 2009, 2010; Beskin et al. 2023) show several characteristic features. One can start by describing the poloidal component, persistently anchored to the central engine, and then discuss the toroidal component as a train of giant electromagnetic waves propagating along the poloidal streamlines.

Close to the BH, the magnetosphere is uncollimated; hence, the streamlines are roughly spherical, θ_s(r)≃θ_s(r_H). In the initial stage of uniform (equilibrium) collimation (Lyubarsky 2010), the entire jet remains in causal contact, and the streamlines are able to adapt to the jet spine boundary, θ_sp(r), to the tension of toroidal magnetic field, and to the centrifugal forces due to both azimuthal rotation and poloidal curvature (Narayan et al. 2009). In the subsequent stage of differential (non-equilibrium) collimation, beyond the fast magnetosonic surface, the inner layer of poloidal field lines bunches towards the jet symmetry axis (Tchekhovskoy et al. 2009; Chatterjee et al. 2019). An important consequence is that the strength of poloidal magnetic field, B_p, evolves very differently along individual field lines.

Close to the jet symmetry axis, magneto-hydrodynamic (MHD) models predict a narrow cylindrical structure called the jet core (e.g. Appl & Camenzind 1993; Beskin & Nokhrina 2009; Lyubarsky 2009; Beskin et al. 2023). Here, the poloidal magnetic field lines are able, more or less, to co-rotate at the rate imposed by the central engine (this corresponds roughly to the light cylinder of pulsars). The characteristic radius of the jet core is R_c ∼ u_c/Ω_B, where u_c = Γ_cβ_c ∼ 1 is the minimal four-velocity of bulk motion along the core⁵. For a ≃ 0.9, this means that R_c ≃ 6M_BH, effectively imposing the BH length scale over great distances. A strictly cylindrical jet core would maintain constant B_p, support only a weak toroidal field (hence little Poynting flux), and inefficient bulk acceleration with Γ_j ≃ B_ϕ/B_p ≃ Ω_BR (Lyubarsky 2009).

Outside the jet core (R > R_c), differential collimation results in poloidal field lines diverging faster than in the spherical case, resulting in B_p declining steeper than 1/r², enhancing the acceleration efficiency due to the magnetic nozzle effect (Begelman & Li 1994). Efficient conversion of magnetic to kinetic energy still requires maintaining equilibrium with the external pressure (Lyubarsky 2010). On the other hand, a steeply decreasing external pressure (e.g. when a gamma-ray burst jet breaks out of a collapsar) enables additional acceleration of the outermost jet layers via the rarefaction mechanism (Tchekhovskoy et al. 2010; Komissarov et al. 2010).

The toroidal magnetic field is induced by the rotation of poloidal field lines at the jet base, where the centrifugal force wins over BH gravity and dominates the magnetic tension. Black hole magnetospheres do not rotate rigidly, and they inject the toroidal field over a broad range of streamlines, well away from the jet core. The relative strengths of the toroidal and poloidal components define the magnetic pitch, 𝒫 = RB^r/B^ϕ (Appl et al. 2000; Bromberg et al. 2019). The initial stage of jet collimation is characterised by roughly uniform 𝒫. At the jet core radius, R_c, one has B^ϕ ≃ B^r; hence, 𝒫 ≃ R_c across the jet. This implies that radial profiles of toroidal field strength B^ϕ(R) peak well outside the jet core, and hence R_c is not their characteristic lateral scale.

(Chatterjee et al. 2019) performed axisymmetric general-relativistic MHD (GRMHD) simulations of strongly collimated relativistic jets extending to distances beyond 10⁵r_g. They demonstrated strong bunching of poloidal field lines (Tchekhovskoy et al. 2009) leading to the formation of a jet core with a radius of R_c ≃ 4R_L (Beskin et al. 2023) (their a = 15/16 implies Ω_B ≃ Ω_H/2 ≃ 0.17; hence, R_L ≃ 5.7R_g). On the one hand, bunching of poloidal field lines is a consequence of the magnetic tension of the toroidal field; on the other hand, bunched poloidal field lines transmit the toroidal field parallel to the jet core. The toroidal magnetic energy gradually converts into the kinetic energy of the plasma, resulting in high bulk Lorentz factors in the vicinity of the jet core, which was also reported by Chatterjee et al. (2019). At the same time, the outer jet regions are loaded with denser proton plasma from the sheath region, apparently due to the interchange (magnetic Rayleigh-Taylor) instability⁶. Mass loading via instabilities is expected to develop a highly inhomogeneous filamentary structure; whether such filaments can reach the jet core is unknown. It has been suggested that filamentary mass loading would result in highly inhomogeneous jet magnetisation with σ_max ≳ 10³ matching the maximum particle energies inferred from non-thermal emission of blazars (Nalewajko 2016). On average, mass loading should reduce not only the bulk Lorentz factor, Γ, but also the Michel parameter, μ = Γ(1 + σ) (maximum achievable Γ, hence the potential for further acceleration). The overall effect of poloidal field bunching is thus the concentration of the jet energy flux along the jet core.

The poloidal pressure of the toroidal magnetic field is the main force accelerating the jet to relativistic speeds. A non-uniform structure of B^ϕ(θ) results in a non-uniform structure of Γ(θ). The acceleration efficiency depends also on the geometry of jet streamlines (the rate of velocity divergence). Thus, acceleration should be less efficient in the jet core, and particularly along the jet axis. The toroidal magnetic field is the main carrier of the poloidal Poynting flux, which is converted into plasma kinetic energy density along each streamline.

The lateral structure of a relativistic jet spine in the differentiated ACZ stage is summarised schematically in Figure 6. Across such a jet spine, we distinguish three zones: the innermost cylindrical jet core, the intermediate quasi-cylindrical bunched spine, and the outermost paraboloidal main spine. The bunched spine is distinguished as the location of the peak bulk Lorentz factor, Γ, and the toroidal magnetic field, B_ϕ, which means the peak energy (Poynting) flux density. The main spine zone is the most extended and carries the bulk of the energy flux. It can be polluted by protons from the jet sheath via boundary instabilities (Chatterjee et al. 2019), but this reduces the bulk Lorentz factor and dilutes the energy flux density. The jet core is prone to current-driven instability (CDI) (Kruskal & Schwarzschild 1954; Freidberg 1982; Begelman 1998), which can accelerate particles via reconnection in high-σ conditions (Alves et al. 2018; Davelaar et al. 2020; Ortuño-Macías et al. 2022). However, the reduced bulk Lorentz factor of the jet core would reduce the relativistic boost of associated radiation.

Fig. 6. Proposed lateral structure (not to scale) of a relativistic jet differentiated due to poloidal field bunching. The ‘bunched spine’ zone is introduced as the region maximizing the jet energy flux density. The introduction of plasmoids from reconnection layers created due to CDI in the jet core allow the energy density enhancement factors of the plasmoids and of the bunched jet spine to multiply.

Efficient poloidal field bunching would result in a narrow bunched spine with a radius of peak toroidal field strength R_peak ≪ R_j. It also implies the enhancement of the peak Poynting flux density, S, compared with the mean value across the jet spine. In Appendix C, we consider a generic model for radial profiles S(R) and calculate the dependence between the enhancement factor f_S = S_peak/⟨ S⟩ and R_peak/R_j. We find that f_S ∝ (R_peak/R_j)^−3/2, for example enhancement by a factor of f_S ∼ 10, requires R_peak/R_j ∼ 0.15.

In numerical simulations of the CDI, perturbations tend to spread from the most unstable toroidal field core to the outside (Bromberg & Tchekhovskoy 2016; Bromberg et al. 2019; Ortuño-Macías et al. 2022). Hence, we propose that plasmoids from reconnection layers induced by CDI in the jet core may spread to the bunched spine zone, as is indicated in Figure 6. During that process, they would be need to be carried by the poloidal bulk flow of the jet in order to relativistically boost the emitted radiation. Carrying large plasmoids by background high-σ flow is not expected to be a highly efficient process, because of their significant inertia, which in a slightly different context results in an anti-correlation between plasmoid growth and bulk acceleration (Sironi et al. 2016; Ortuño-Macías & Nalewajko 2020). Moreover, the poloidal velocity shear is rarely accounted for in studies of CDI (Nalewajko & Begelman 2012), and its effect is still unclear. Alternatively, strongly sheared magnetic perturbations may induce reconnection layers directly in the bunched spine zone.

4. Discussion and conclusions

We considered an analogy between the lateral structures of relativistic jets and plasmoids produced by magnetic reconnection. In both cases, a key role is played by the toroidal component of magnetic field, the tension of which introduces an anisotropic stress that in principle can compress the energy density of both plasma and the magnetic field.

Compression of the energy density in reconnection plasmoids was investigated by kinetic PIC numerical simulations in 2D and 3D periodic domains, starting from hot Harris layers and evolving to relaxed monster plasmoids. We considered the effects of domain size, background magnetisation, and the guide magnetic field. In our simulations without a guide field, the plasma energy density achieves values of u_pl ∼ 10²u_B, 0 (relative to the initial background magnetic energy density), higher than the magnetic energy density, u_B ∼ 10u_B, 0. The highest energy densities correspond to the secondary plasmoids emerging between merging primary plasmoids. However, the initial plasma energy density extends to ∼4u_B, 0 for the hot drifting particles. In the presence of a guide field or third dimension, the plasma energy densities are reduced by a few factors.

We have also re-analysed the results of previous PIC simulations, first presented in Ortuño-Macías & Nalewajko (2020), which used open boundaries and synchrotron cooling. Those results show power law tails in the distributions of both u_pl and u_B with indices close to −5/3. The cooling efficiency affects mainly the extent of the u_B distribution, which may explain the apparent lack of such a power law tail in the results of our new simulations.

In relativistic jets, the effect of the toroidal magnetic field has been controversial. Eventually, it was understood that its tension is insufficient to collimate such jets. However, toroidal field tension has the subtle effect of gradually differentiating the lateral jet structure. The key effect of poloidal field bunching has been identified by Tchekhovskoy et al. (2009). We propose that this effect forms an intermediate layer across the jet spine, referred to as the bunched spine, adjacent to the central jet core, corresponding roughly to the peaks of both the toroidal field strength, B_ϕ(R), and the bulk Lorentz factor, Γ(R). The bunched spine is characterised by an enhanced Poynting flux density, S. For generic radial profiles of S(R), the enhancement factor, f_S = S_peak/⟨ S⟩, is sensitive to the relative peak radius, R_peak/R_j. Poloidal field bunching was demonstrated in 2D GRMHD global simulations of relativistic jets (Chatterjee et al. 2019); however, it appears to be less efficient in 3D simulations.

An additional enhancement of the luminosity of short flares of energetic radiation can be provided by the effect of kinetic beaming (Cerutti et al. 2012), also known as pitch-angle anisotropy (Comisso & Sironi 2021; Sobacchi et al. 2023). This effect is not expected to be important in relaxed reconnection plasmoids, but rather during the formation of plasmoid chains, and likely also during plasmoid mergers. Our PIC simulations confirm that plasmoid mergers produce the highest energy densities in reconnection layers, and the resulting radiation intensity could be further enhanced by kinetic beaming.

The possibility of locally concentrating the momentum flux density of relativistic jets is important also for alternative models of rapid blazar flares; for example, the obliteration of stars wandering into the jet (e.g. Bednarek & Protheroe 1997). The efficiency of bulk acceleration can be increased in a sharply inhomogeneous jet due to impulsive acceleration (Granot et al. 2011), which can also be described as a relativistic whip.

4.1. The jet of M87

The proposed scenario can be confronted with observations of relativistic AGN jets resolved by VLBI. In the best-studied case of M87 (for review see Hada et al. 2024), the approaching jet (seen at the viewing angle of ≃18°) has long been known to be limb-brightened at scales⁷ of ∼50 mas (Reid et al. 1989), ∼1 mas (Junor et al. 1999), and ∼0.2 mas (Kim et al. 2018). The < 0.6 mas limbs were modelled by Punsly (2022, 2023, 2024) in terms of the thick tubular jet sheath. Recently, a third central ridge has been reported on scales of ∼0.5 mas (Lu et al. 2023) and ∼8 mas (Tazaki et al. 2023), and possibly also at ∼300 mas⁸. The 0.5 mas ridge was interpreted by Lyutikov & Ibrahim (2024) in terms of strong beaming due to the instant bulk acceleration of a pair-dominated hollow-cone jet. The central ridge appears more likely to be a signature of the jet core weakly beamed due to a reduced bulk Lorentz factor. The highly relativistic jet spine would be de-beamed at such a large viewing angle.

According to the AGN unification paradigm (Urry & Padovani 1995), a Fanaroff-Riley type I radio galaxy such as M87 should be a misaligned counterpart of a BL Lac type blazar. Beamed non-thermal blazar emission should be produced in a highly dissipative blazar zone. Models of BL Lac type blazars do not provide strong constraints on the location. A de-projected distance of 1 pc corresponds to ≃4 mas, indicating that the ∼8 mas ridge is potentially adjacent to a blazar zone in the form of bunched spine.

A well known dissipative zone in the M87 jet is the HST-1 knot observed at ≃860 mas (Giroletti et al. 2012). It shows evidence of apparently superluminal motions (Biretta et al. 1999; Cheung et al. 2007), and produced a spectacular X-ray/optical/radio outburst over a timescale of a few years (e.g. Harris et al. 2006, 2009). It has been suggested to be the blazar zone of the M87 jet (Barnacka et al. 2014), which would pose severe problems of sufficient energy density and a short variability timescale at a de-projected distance of ∼220 pc.

4.2. Conclusions

In the multi-scale system of a relativistic jet hosting plasmoids produced by magnetic reconnection, we consider it more likely for a significant (factor ≳10) energy density enhancement to be achieved due to tension of the toroidal magnetic field on the local scale of plasmoid cores, rather than on the global scale of a bunched jet spine (the layer adjacent to the central jet core). Even a moderate (∼3) global-scale enhancement combined with a stronger (∼10) local-scale enhancement would boost the luminosity of rapid blazar flares by a factor of ∼30 at standard jet bulk Lorentz factors of Γ_j ∼ 10 − 20, strongly reducing the extreme energetic requirements of such flares.

Acknowledgments

We gratefully acknowledge Poland’s high-performance Infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH, PCSS, CI TASK, WCSS) for providing computer facilities and support within computational grants no. PLG/2023/016444 and PLG/2024/017356; and the Nicolaus Copernicus Astronomical Center for providing the Chuck cluster. This work was supported by the Polish National Science Centre grant 2021/41/B/ST9/04306.

References

Appendix A: Energetics of blazar flares

Rapid and luminous blazar flares require efficient conversion of very high energy densities. For a typical estimate, consider a compact flaring region of radius R_fl′ propagating with speed β_fl = v_fl/c and Lorentz factor Γ_fl = (1 − β_fl²)^−1/2 (defining reference frame 𝒪′), and observed at viewing angle θ_obs ≃ 1/Γ_fl corresponding to a Doppler factor 𝒟_fl = [Γ_fl(1 − β_fl cos θ_obs)]⁻¹ ∼ Γ_fl resulting in a luminosity boost⁹L_fl, obs = 𝒟_fl⁴L_fl′ and shortening of variability timescale t_fl = t_fl′/𝒟_fl. In blazars, relativistic motions with bulk Lorentz factors Γ_fl ∼ 20 can be deduced from apparently superluminal proper motions with speeds v_app ∼ 20c (with c the speed of light) along the jet (e.g. Jorstad et al. 2001; Lister et al. 2016). The intrinsic parameters can be approximated as t_fl′∼R_fl′/c and L_fl′∼4πR_fl′²u_fl′c with u_fl′ the intrinsic radiation energy density. The total power carried by the flare photons can be estimated as P_fl ∼ πR_fl²Γ_fl²u_fl′c ∼ [Γ_fl²/(4𝒟_fl⁴)]L_fl, obs ∼ L_fl, obs/(4Γ_fl²) with R_fl = R_fl′∼𝒟_flct_fl. That power can be compared with the Eddington luminosity of the supermassive BH L_Edd ≃ 1.5 × 10⁴⁷(M_BH/10⁹ M_⊙) erg/s with M_BH ∼ 10⁹ M_⊙ the BH mass in units of the solar mass M_⊙. For flaring FSRQs this gives P_fl, FSRQ ∼ 0.4(Γ_fl/20)⁻²L_Edd, which is dangerously close to L_Edd, and this is one reason to consider significantly higher values of Γ_fl (Ackermann et al. 2016). For flaring HBLs, P_fl, HBL ∼ 4 × 10⁻⁴(Γ_fl/20)⁻²L_Edd, more comfortable at least in terms of L_Edd.

Consider that this compact blazar zone is only a small part of a relativistic jet with a total power P_j and Lorentz factor Γ_j. At distance scale r the jet has lateral radius R_j = θ_jr with θ_j the half-opening angle. Interferometric radio observations of pc-scale blazar jets suggest that Γ_jθ_j ∼ 0.2 (Clausen-Brown et al. 2013). The jet can be subdivided laterally into streamtubes, of which we are interested only in the streamtube of radius R_st = R_fl that completely and minimally contains the flaring region. Neglecting the poloidal line bunching and assuming that power per unit solid-angle is roughly uniform across the jet, the fraction of jet power P_st in that streamtube would scale as P_st/P_j ≃ (R_st/R_j)². We can ask, at what distances is the jet power passing through a flaring region sufficient to power the observed flare, P_j → fl > P_fl. We obtain the following estimate:

$$\begin{array}{c}\hfill r<c{t}_{\mathrm{fl}}\frac{2{\mathcal{D}}_{\mathrm{fl}}^3}{{\mathrm{\Gamma }}_{\mathrm{fl}}}\frac{{\mathrm{\Gamma }}_{\mathrm{j}}}{{\mathrm{\Gamma }}_{\mathrm{j}}{\theta }_{\mathrm{j}}}\sqrt{\frac{P_{\mathrm{j}}}{L_{\mathrm{fl},\mathrm{obs}}}}.\end{array}$$(A.1)

In the case of FSRQs we set Γ_j = 20 and obtain r < 0.02(Γ_fl/20)²(P_j/L_Edd)^1/2 pc. This is shorter from constraints using other methods, which typically find r_FSRQ ∼ 0.1 pc (Nalewajko et al. 2014). This is another reason to increase the Lorentz factor of the flaring region. Setting Γ_fl = 50 (Ackermann et al. 2016) would result in P_fl, FSRQ ∼ 0.07L_Edd and r < 0.11 pc. This is particularly challenging for FSRQs producing ∼100 GeV flares like PKS 1222+216, for which pair-production absorption by optical/UV broad emission lines would be important at distances r ≲ 0.5 pc (Nalewajko et al. 2012). In the case of HBLs, r < 0.17(Γ_fl/20)²(P_j/L_Edd)^1/2 pc.

Appendix B: Scale separation in relativistic jets

Relativistic jets are characterised by very large separation between their macroscopic global scales (e.g. the lateral jet radius R_j) and the microscopic plasma scales (e.g. the gyroradius R_L′, formally in the jet co-moving frame). A typical relativistic AGN jet with Lorentz factor Γ_j ∼ 20 and half-opening angle θ_j ∼ 0.2/Γ_j ∼ 0.01 (Clausen-Brown et al. 2013) at the distance of r ∼ 0.1 pc (typical for blazar zones; Nalewajko et al. 2014) would have a radius of R_j ∼ 10⁻³ pc ≃ 3 × 10¹⁵ cm and may contain magnetic fields of co-moving strength B′∼2.7 G¹⁰, and electrons of mean random Lorentz factor γ_e ∼ 100 having a gyroradius of R_L′=γ_em_ec²/eB′∼6 × 10⁴ cm, resulting in a huge scale separation of R_j/R_L′∼5 × 10¹⁰. This would mean an enormous potential for energy density enhancement, if it could be fully utilised.

The presence of compact jet core would introduce an intermediate scale – the core radius R_c ∼ 6M_BH (for high BH spin a ≃ 0.9). For a BH mass of M_BH ∼ 10⁹ M_⊙ ∼ 1.5 × 10¹⁴ cm, one could expect R_c ∼ 10¹⁵ cm (and R_j ∼ 20M_BH at r ∼ 2000M_BH), implying R_j/R_c ∼ 3. This ratio could be increased by considering a lower M_BH or a broader, less relativistic jet, but it can hardly be much larger than ∼10, and thus large energy density enhancement should not be expected to result from a compact jet core. Global GRMHD numerical simulations of jets reaching distances of ∼10³M_BH with well resolved jet spine of radius R_sp, both in 2D (Chatterjee et al. 2019) and in 3D (Ripperda et al. 2022; Salas et al. 2024), show bulk Lorentz factors Γ_j monotonically increasing across the jet from its axis until at least ∼0.4R_sp, which is consistent with toroidal magnetic field peaking at ≳R_sp/2, which implies no significant energy density enhancement due to tension of global toroidal field.

In magnetic reconnection involving macroscopically large (∼R_j) and microscopically thin (∼R_L′) current layers, plasmoids will be produced and grown to a broad range of radii R_pl, potentially spanning the entire available range of scales R_L′< R_pl < R_j. For plasmoids produced during relativistic magnetic reconnection, our PIC simulations with initial scale separation up to L/ρ₀ = 3600 produced relaxed monster plasmoids with scale separation of R_out/R_c ≃ 13 and with contrast of total energy density of u_tot(R_c)/u_tot(R_out)≃16, much less than the idealised expectation of (R_out/R_c)² ≃ 170. However, in a hypothetical case of large-scale relativistic reconnection producing a plasmoid reaching a size approaching the jet radius, e.g. R_out ≃ L/2 ∼ R_j/10 and initial scale separation of L/R_L′∼10¹⁰, assuming that the plasmoid core radius scales like R_c/L ∼ 5.6 × 10⁻⁴(L/R_L′)^1/2, this would imply a scale separation of R_out/R_c ∼ 8 × 10⁷, and idealised energy density enhancement up to ∼6 × 10¹⁵. Such large plasmoids were considered by Morris et al. (2019) to explain rapid gamma-ray flares of blazars in terms of synchrotron self-Compton emission, assuming uniform energy density, which in the context of this work appears rather unfortunate.

Appendix C: Model for bunching momentum flux across jet

As discussed in Section 3, bunching (internal collimation) of poloidal magnetic field towards the jet axis (Tchekhovskoy et al. 2009) is expected to bunch the Poynting (electromagnetic momentum) flux into the bunched spine layer. Shifting inwards the radius R_peak of peak momentum flux density relative to the jet (spine) radius R_j affects the mean momentum flux density across the jet. By considering a radial profile S(R) of momentum flux density, one can calculate the ratio of peak value S_peak = S(R_peak) to the mean value ⟨ S⟩ = (2/R_j²)∫₀^R_jRS(R) dR as the jet momentum enhancement factor f_S = S_peak/⟨ S⟩. In the highly magnetised jet core (R < R_peak), the momentum flux is dominated by the Poynting flux of density S(R)≃(c/4π)v_z(R)B_ϕ²(R). The toroidal field should satisfy B_ϕ(R ≪ R_j)∝R in order to assure a finite poloidal current density j_z(R) = (c/4πR) d(RB_ϕ)/dR and a uniform magnetic pitch 𝒫(R) = RB_z(R)/B_ϕ(R) for B_z(R ≪ R_j)≃const. Even in the case of uniformly accelerated core with v_z(R ≪ R_j)≃const, the innermost momentum flux density should scale like S(R ≪ R_j)∝R².

Consider a mathematical model of momentum flux density profile:

$$\begin{array}{c}\hfill S(R)=S_0{(R/R_{\mathrm{j}})}^p{(1-R/R_{\mathrm{j}})}^q\end{array}$$(C.1)

The peak radius is given by R_peak/R_j = p/(p + q), corresponding to S_peak = S₀p^pq^q/(p + q)^p + q. We choose p = 2 to satisfy the inner scaling, and q ≥ 2 in order to have R_peak ≤ R_j/2. Substituting r = R/R_j, one can integrate the mean value:

$$\begin{array}{c}\hfill \frac{\langle S\rangle }{2S_0}={\displaystyle {\int }_0^1}{r}^{p+1}{(1-r)}^q\mathrm{d}r=\frac{\mathrm{\Gamma }(p+2)\mathrm{\Gamma }(q+1)}{\mathrm{\Gamma }(p+q+3)}\end{array}$$(C.2)

The peak enhancement factor is thus:

$$\begin{array}{c}\hfill f_{\mathrm{S}}=\frac{S_{\mathrm{peak}}}{\langle S\rangle }=\frac{p^p{q}^q\mathrm{\Gamma }(p+q+3)}{2{(p+q)}^{p+q}\mathrm{\Gamma }(p+2)\mathrm{\Gamma }(q+1)}.\end{array}$$(C.3)

Examples of profiles S(R)/⟨ S⟩ and the function f_S(R_peak) for p = 2 are presented in Figure C.1. The case without poloidal field bunching (q = 2) corresponds to R_peak = R_j/2 and is characterised by small enhancement of f_S ≃ 1.9. However, as R_peak/R_j decreases due to stronger bunching, f_S increases, e.g. to ≃6.6 for R_peak = R_j/5, to ≃22 for R_peak = R_j/10. The enhancement function for p = 2 and 0.1 < R_peak/R_j < 0.5 can be well approximated as f_S ≃ 0.66(R_peak/R_j)^−3/2.

As will be shown in a different article, in the case without bunching (q = 2) such profile makes a good approximation to the radial distribution of poloidal Poynting flux in relativistic jets obtained in "extreme-resolution" GRMHD simulation of (Ripperda et al. 2022) reaching distances of ∼10³R_g in 3D. The cases with bunching need to be compared with GRMHD simulations reaching larger distances, e.g. ∼10⁵R_g in 2D (Chatterjee et al. 2019).

Fig. C.1. Solid coloured lines show examples of model radial profiles across a cylindrical jet of radius Rj of momentum flux density S(R) (defined by Eq. C.1) normalised by the mean value ⟨ S⟩ (Eq. C.2) for p = 2 and different values of q (corresponding to Rpeak/Rj = 1/10, 1/5, 1/3, 1/2). The thick dashed line shows the enhancement factor fS = Speak/⟨ S⟩ (Eq. C.3) as function of Rpeak/Rj.

All Tables

All Figures

	Fig. 1. Logarithms, f = log10(F), of volume distributions, F(μ) = dF/dμ, over the argument μ = log10(u/uB, 0), with u the energy density: of magnetic fields, uB = B2/8π (left panel), and of the plasma, upl = ⟨γ⟩nmc2 (right panel). Functions f(μ) were averaged over the duration of each simulation.
In the text

Fig. 2. Maps of log10(u/uB, 0)(x, y) for magnetic energy density, uB = B2/8π (upper panels), and plasma energy density, upl = ⟨γ⟩nmc2 (middle panels), for relaxed monster plasmoids at the end of each simulation. In the lower panels, we compare 1D energy density profiles of (u/uB, 0)(x, y0) measured along the strip indicated in the above maps by the dashed grey lines. A common colour scale for all maps is referenced along the left axes in the lower panels. Here, we present the effect of a guide field for σ0 = 10 and L/ρ0 = 1800. From the left, the columns show simulations: (1) L1800_σ10, (2) L1800_σ10_Bg05, and (3) L1800_σ10_Bg1.

In the text

Fig. 3. Profiles of magnetic (green lines) and plasma (red lines) energy densities across (i.e. along the y coordinate) relaxed monster plasmoids. Panels from the left compare: (1) different sizes, L/ρ0, of simulation domain (2D with σ0 = 10, Bz = 0); (2) different magnetisations, σ0 (2D with L/ρ0 = 1800, Bz = 0); (3) 3D and 2D domains (with L/ρ0 = 900, σ0 = 10, Bz = 0); (4) different guide field strengths, Bg/B0 (2D with σ0 = 10, L/ρ0 = 1800).

In the text

	Fig. 4. Same as Figure 2, but for plasmoid mergers that maximise the plasma energy density, upl, for the same set of simulations: (1) L1800_σ10, (2) L1800_σ10_Bg05, and (3) L1800_σ10_Bg1.
In the text

	Fig. 5. Same as Figure 1, but for the open-boundary simulations with synchrotron cooling first presented in Ortuño-Macías & Nalewajko (2020). The black dashed lines indicate power laws of index −5/3.
In the text

Fig. 6. Proposed lateral structure (not to scale) of a relativistic jet differentiated due to poloidal field bunching. The ‘bunched spine’ zone is introduced as the region maximizing the jet energy flux density. The introduction of plasmoids from reconnection layers created due to CDI in the jet core allow the energy density enhancement factors of the plasmoids and of the bunched jet spine to multiply.

In the text

	Fig. C.1. Solid coloured lines show examples of model radial profiles across a cylindrical jet of radius Rj of momentum flux density S(R) (defined by Eq. C.1) normalised by the mean value ⟨ S⟩ (Eq. C.2) for p = 2 and different values of q (corresponding to Rpeak/Rj = 1/10, 1/5, 1/3, 1/2). The thick dashed line shows the enhancement factor fS = Speak/⟨ S⟩ (Eq. C.3) as function of Rpeak/Rj.
In the text