<figure><img src="image_1.png" /></figure>

HAL
open science

# Extreme ion acceleration at extragalactic jet termination shocks

Benoît Cerutti, Gwenael Giacinti

► To cite this version:

Benoît Cerutti, Gwenael Giacinti. Extreme ion acceleration at extragalactic jet termination shocks. Astronomy and Astrophysics - A&A, 2023, 676, pp.A23. 10.1051/0004-6361/202346481. hal-04060434

HAL Id: hal-04060434

https://hal.science/hal-04060434

Submitted on 5 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Extreme ion acceleration at extragalactic jet termination shocks

Benoît Cerutti<fnref n="1" /> and Gwenael Giacinti<fnref n="2" />,<fnref n="3" />

Received 22 March 2023 / Accepted 15 May 2023

## ABSTRACT

**Context.** Extragalactic plasma jets are some of the few astrophysical environments able to confine ultra-high-energy cosmic rays, but whether they are capable of accelerating these particles is unknown.

**Aims.** In this work, we revisit particle acceleration at relativistic magnetized shocks beyond the local uniform field approximation by considering the global transverse structure of the jet.

**Methods.** Using large two-dimensional particle-in-cell simulations of a relativistic electron-ion plasma jet, we show that the termination shock forming at the interface with the ambient medium accelerates particles up to the confinement limit.

**Results.** The radial structure of the jet magnetic field leads to a relativistic velocity shear that excites a von Kármán vortex street in the downstream medium trailing behind an over-pressured bubble filled with cosmic rays. Particles are efficiently accelerated at each crossing of the shear flow boundary layers.

**Conclusions.** These findings support the idea that extragalactic plasma jets may be capable of producing ultra-high-energy cosmic rays. This extreme particle acceleration mechanism may also apply to microquasar jets.

**Key words.** acceleration of particles – shock waves – methods: numerical – galaxies: jets

## 1. Introduction

Active supermassive black holes nested in the heart of galaxies produce powerful plasma outflows in the form of relativistic magnetized jets (Blandford et al. 2019; Hardcastle & Croston 2020; Gabuzda 2021). These collimated outflows propagate on extragalactic scales and carry energy and angular momentum away from the central source. They abruptly terminate their journey as bright radio hotspots, where the jet shocks the ambient medium and inflates a giant radio bubble downstream, also known as the lobe. The relativistic speeds of the flow combined with the typical size and magnetic field strength involved make extragalactic jets and their lobes some of the few astrophysical objects capable of confining the most energetic cosmic rays received on Earth (Hillas 1984; Giannios 2010; Matthews et al. 2018), but how these particles are accelerated to such ultra-high energies is still an open question.

The jet termination shock itself is a prime candidate for the location of particle acceleration. Unfortunately, it is well established that relativistic uniform shocks are poor accelerators because the transverse nature of the magnetic field with respect to the shock normal drastically reduces the efficiency of diffusive shock acceleration (Begelman & Kirk 1990; Bell et al. 2018; see however Kirk et al. 2023). Even under the most favorable conditions of an unmagnetized shock, particle acceleration via this mechanism would be too slow to explain the extreme energization rate required to reach the confinement limit of the source (Sironi et al. 2013; Plotnikov et al. 2018). In contrast, relativistic magnetic reconnection is sufficiently fast to alleviate the above difficulties, but particle acceleration up to the highest energies would require a system-size current sheet and an

unrealistically high plasma magnetization (Werner et al. 2016). A large-scale velocity shear provides another way to accelerate particles efficiently (Ostrowski 1998; Rieger 2022). This situation could be found for instance at the interface between the jet and its backflow or the external medium (e.g., Alves et al. 2014; Kimura et al. 2018; Matthews et al. 2019; Sironi et al. 2021; Meli et al. 2023; Wang et al. 2023), but it is unclear how it could operate at the jet hotspots and lobes where in situ particle acceleration must occur.

In this work, we present a novel ab initio plasma model of particle acceleration at the jet termination shock beyond the local uniform field approximation that has been used so far by considering the full radial dependence of the magnetic field across the jet. For this purpose, a new numerical setup inspired by a force-free jet solution is introduced in Sect. 2. Section 3 focuses on the dynamics of the flow, acceleration of the ions, and their escape in the downstream medium. Section 4 summarizes the main findings of this work and their astrophysical implications.

## 2. Numerical methods and setup

We focus our attention on the global structure of the shock front and nonthermal particle acceleration using particle-in-cell simulations performed with the Zeltron code (Cerutti et al. 2013). The plasma is so dilute that it can be considered as collisionless. The electron-ion plasma dynamics was simulated with an artificially low mass ratio $m_i/m_e = 25$ to reduce the computational cost, but it has remained sufficiently large to preserve a safe scale separation (Sironi et al. 2013). The jet was modeled as a cylindrical plasma column of radius $R_j$ moving with a uniform bulk Lorentz factor $\Gamma_j = 1/(1 - V_j^2/c^2)^{1/2} = 10$ along

<footnote n="1">¹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
e-mail: benoit.cerutti@univ-grenoble-alpes.fr</footnote>

<footnote n="2">² Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai 201210, PR China</footnote>

<footnote n="3">³ School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, PR China
e-mail: gwenael.giacinti@sjtu.edu.cn</footnote>
# Extreme ion acceleration at extragalactic jet termination shocks

Benoît Cerutti<fnref n="1" /> and Gwenael Giacinti<fnref n="2" />,<fnref n="3" />

Received 22 March 2023 / Accepted 15 May 2023

## ABSTRACT

**Context.** Extragalactic plasma jets are some of the few astrophysical environments able to confine ultra-high-energy cosmic rays, but whether they are capable of accelerating these particles is unknown.

**Aims.** In this work, we revisit particle acceleration at relativistic magnetized shocks beyond the local uniform field approximation by considering the global transverse structure of the jet.

**Methods.** Using large two-dimensional particle-in-cell simulations of a relativistic electron-ion plasma jet, we show that the termination shock forming at the interface with the ambient medium accelerates particles up to the confinement limit.

**Results.** The radial structure of the jet magnetic field leads to a relativistic velocity shear that excites a von Kármán vortex street in the downstream medium trailing behind an over-pressured bubble filled with cosmic rays. Particles are efficiently accelerated at each crossing of the shear flow boundary layers.

**Conclusions.** These findings support the idea that extragalactic plasma jets may be capable of producing ultra-high-energy cosmic rays. This extreme particle acceleration mechanism may also apply to microquasar jets.

**Key words.** acceleration of particles – shock waves – methods: numerical – galaxies: jets

## 1. Introduction

Active supermassive black holes nested in the heart of galaxies produce powerful plasma outflows in the form of relativistic magnetized jets (Blandford et al. 2019; Hardcastle & Croston 2020; Gabuzda 2021). These collimated outflows propagate on extragalactic scales and carry energy and angular momentum away from the central source. They abruptly terminate their journey as bright radio hotspots, where the jet shocks the ambient medium and inflates a giant radio bubble downstream, also known as the lobe. The relativistic speeds of the flow combined with the typical size and magnetic field strength involved make extragalactic jets and their lobes some of the few astrophysical objects capable of confining the most energetic cosmic rays received on Earth (Hillas 1984; Giannios 2010; Matthews et al. 2018), but how these particles are accelerated to such ultra-high energies is still an open question.

The jet termination shock itself is a prime candidate for the location of particle acceleration. Unfortunately, it is well established that relativistic uniform shocks are poor accelerators because the transverse nature of the magnetic field with respect to the shock normal drastically reduces the efficiency of diffusive shock acceleration (Begelman & Kirk 1990; Bell et al. 2018; see however Kirk et al. 2023). Even under the most favorable conditions of an unmagnetized shock, particle acceleration via this mechanism would be too slow to explain the extreme energization rate required to reach the confinement limit of the source (Sironi et al. 2013; Plotnikov et al. 2018). In contrast, relativistic magnetic reconnection is sufficiently fast to alleviate the above difficulties, but particle acceleration up to the highest energies would require a system-size current sheet and an

unrealistically high plasma magnetization (Werner et al. 2016). A large-scale velocity shear provides another way to accelerate particles efficiently (Ostrowski 1998; Rieger 2022). This situation could be found for instance at the interface between the jet and its backflow or the external medium (e.g., Alves et al. 2014; Kimura et al. 2018; Matthews et al. 2019; Sironi et al. 2021; Meli et al. 2023; Wang et al. 2023), but it is unclear how it could operate at the jet hotspots and lobes where in situ particle acceleration must occur.

In this work, we present a novel ab initio plasma model of particle acceleration at the jet termination shock beyond the local uniform field approximation that has been used so far by considering the full radial dependence of the magnetic field across the jet. For this purpose, a new numerical setup inspired by a force-free jet solution is introduced in Sect. 2. Section 3 focuses on the dynamics of the flow, acceleration of the ions, and their escape in the downstream medium. Section 4 summarizes the main findings of this work and their astrophysical implications.

## 2. Numerical methods and setup

We focus our attention on the global structure of the shock front and nonthermal particle acceleration using particle-in-cell simulations performed with the Zeltron code (Cerutti et al. 2013). The plasma is so dilute that it can be considered as collisionless. The electron-ion plasma dynamics was simulated with an artificially low mass ratio $m_i/m_e = 25$ to reduce the computational cost, but it has remained sufficiently large to preserve a safe scale separation (Sironi et al. 2013). The jet was modeled as a cylindrical plasma column of radius $R_j$ moving with a uniform bulk Lorentz factor $\Gamma_j = 1/(1 - V_j^2/c^2)^{1/2} = 10$ along

<footnote n="1">¹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
e-mail: benoit.cerutti@univ-grenoble-alpes.fr</footnote>

<footnote n="2">² Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai 201210, PR China</footnote>

<footnote n="3">³ School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, PR China
e-mail: gwenael.giacinti@sjtu.edu.cn</footnote>
the +z direction. The exact magnetic structure in the jet is not well constrained by observations, but it is most likely dominated by the toroidal component far away from the jet-launching point (Gabuzda 2021). This component stems from the twist of poloidal field lines by the accretion disk orbiting around the central black hole (Blandford & Payne 1982) or by the spin of the black hole itself (Blandford & Znajek 1977).

For a steady, axisymmetric, and Poynting-flux-dominated jet at launching, the toroidal component is obtained from the relativistic Grad-Shafranov equation (Michel 1973; Scharlemann & Wagoner 1973). Assuming an initial uniform poloidal magnetic field, $B_p$, spinning at a constant angular velocity $\Omega$ within the jet radius and gradually cutting off outside, such that $\Omega \cdot B_p > 0$, the solution is

$$B_{\phi} = -\frac{B_0}{2} \left(\frac{R}{R_j}\right) \left[1 - \tanh\left(\frac{R - R_j}{\Delta}\right)\right], \quad (1)$$

where $B_0$ is the fiducial magnetic field strength, $R$ is the cylindrical radius to the jet axis, and $\Delta = 0.5R_j$ is fixed throughout this study. The hyperbolic tangent term has been added here to parametrize a smooth and finite physical edge of the jet within the simulation domain. The jet must transport a current density $J = c\nabla \times B/4\pi$, whose nonvanishing component is along the jet axis,

$$J_z = -\frac{cB_0}{4\pi R_j} \left[ 1 - \tanh\left(\frac{R - R_j}{\Delta}\right) - \frac{R}{2\Delta} \cosh^{-2}\left(\frac{R - R_j}{\Delta}\right) \right]. \quad (2)$$

The first two terms represent a net negative current density uniformly distributed within the core of the jet, while the last term is the positive (return) current that flows along the jet boundary and satisfies the current closure in the system, $\nabla \cdot J = 0$, akin to a coaxial cable (see, Fig. 1). The relativistic bulk motion of the plasma along the jet direction is accompanied by the ideal motion electric field $E = -V_j \times B/c$, whose nonvanishing component is along the radial direction. Thus, the plasma is polarized: the core of the jet is negatively charged, while the sheath is positively charged (and vice versa if $\Omega \cdot B_p < 0$). The motion of this net charge density $\rho$ carries the required current. Although polarized, the plasma remains quasi-neutral: $|\rho|/e \ll n_0$, where $e$ is the electron charge and $n_0$ is the reference plasma density in the jet.

The simulation box is limited to a two-dimensional Cartesian plane (xz coordinates) mimicking the poloidal plane of a cylindrical jet (Rz coordinates). With this choice of coordinates, the expansion of the flow with distance from the jet axis and the drift motion of the plasma associated with the curvature of field lines are missing (Huang et al. 2023). Yet, the Cartesian geometry has the advantage of capturing some of the nonaxisymmetric features that would be missing in a two-dimensional cylindrical simulation. The departure from perfect axisymmetry plays an important role in the escape mechanism of energetic particles as shown below. We defer the study of full three-dimensional simulations of a cylindrical jet to a future work.

The fiducial run is composed of 16 384 × 262 144 cells along the x and z directions, respectively. The box is elongated along the jet direction, it extends between $-4R_j < x < -4R_j$ and $0 < z < 128R_j$. The reference electron plasma skindepth, $d_e = (\Gamma_j m_e c^2 / 4\pi n_0 e^2)^{1/2}$, has been resolved by eight cells and the ion skindepth scale $d_i = d_e (m_i / m_e)^{1/2}$ by 40 cells, meaning that the box size is $410 d_i \times 6554 d_i$. The simulation time step has been fixed at $\omega_{pe} \Delta t = 8.75 \times 10^{-2}$, where $\omega_{pe} = c/d_e$ is the

<figure><img src="image_1.png" /><figcaption>Fig. 1. Transverse magnetic field strength (top panel, $B_y$) and current density (bottom panel, $J_z$) profiles carried by the jet in the upstream medium projected into the xz plane, where $\Delta = 0.5R_j$.</figcaption></figure>

electronic plasma frequency. All boundaries are perfectly reflective for both the fields and the particles. The $z = z_{max}$ boundary mimics the contact discontinuity with the ambient medium that is not simulated here to save on computational resources. Another common optimization procedure is to inject the plasma into the box gradually. Initially, the particle injector is close to the $z_{max}$ wall, and then it recedes away from this boundary at the speed of light. The injected plasma bounces off the wall with no loss of energy such that the beam interacts with itself and the shock forms. The choice of reflecting boundaries on the $\pm x$ directions mimics the presence of a perfectly confining medium surrounding the jet. At injection, the plasma was modeled by eight particles in each cell with a constant density at all radii, including the region outside the jet boundaries. The particles are streaming at the jet bulk Lorentz factor, but a small temperature $T/m_{e,i}c^2 = 0.01$ has been included in addition to spatial filtering of the current densities to limit the growth of spurious Cherenkov radiation (Greenwood et al. 2004).

An initially Poynting-flux-dominated jet gradually converts its magnetic energy into plasma bulk motion, such that if internal dissipation along its way to the termination shock is weak as in Fanaroff-Riley class II jets (Fanaroff & Riley 1974), the magnetic field and the particles could be close to equipartition (Li et al. 1992; Komissarov et al. 2007). In this work, we assume perfect equipartition which translates into the magnetic enthalpy density to the particle enthalpy density ratio, $\sigma \equiv B_0^2 / (4\pi\Gamma_j n_0 m_i c^2) = 1$ (the electronic contribution is neglected in this expression due to the mass ratio difference). We also ran the same numerical setup with $\sigma = 0.1$ and found similar results. We focus this study on the case where the poloidal magnetic component is negligible, $B_z = 0$. In reality, it is likely that a small but finite poloidal component still exists and stabilizes the jet against current driven instabilities such as the kink mode which could lead to its disruption akin to a Z-pinch configuration
the +z direction. The exact magnetic structure in the jet is not well constrained by observations, but it is most likely dominated by the toroidal component far away from the jet-launching point (Gabuzda 2021). This component stems from the twist of poloidal field lines by the accretion disk orbiting around the central black hole (Blandford & Payne 1982) or by the spin of the black hole itself (Blandford & Znajek 1977).

For a steady, axisymmetric, and Poynting-flux-dominated jet at launching, the toroidal component is obtained from the relativistic Grad-Shafranov equation (Michel 1973; Scharlemann & Wagoner 1973). Assuming an initial uniform poloidal magnetic field, $B_p$, spinning at a constant angular velocity $\Omega$ within the jet radius and gradually cutting off outside, such that $\Omega \cdot B_p > 0$, the solution is

$$B_{\phi} = -\frac{B_0}{2} \left(\frac{R}{R_j}\right) \left[1 - \tanh\left(\frac{R - R_j}{\Delta}\right)\right], \quad (1)$$

where $B_0$ is the fiducial magnetic field strength, $R$ is the cylindrical radius to the jet axis, and $\Delta = 0.5R_j$ is fixed throughout this study. The hyperbolic tangent term has been added here to parametrize a smooth and finite physical edge of the jet within the simulation domain. The jet must transport a current density $\mathbf{J} = c\nabla \times \mathbf{B}/4\pi$, whose nonvanishing component is along the jet axis,

$$J_z = -\frac{cB_0}{4\pi R_j} \left[ 1 - \tanh\left(\frac{R - R_j}{\Delta}\right) - \frac{R}{2\Delta} \cosh^{-2}\left(\frac{R - R_j}{\Delta}\right) \right]. \quad (2)$$

The first two terms represent a net negative current density uniformly distributed within the core of the jet, while the last term is the positive (return) current that flows along the jet boundary and satisfies the current closure in the system, $\nabla \cdot \mathbf{J} = 0$, akin to a coaxial cable (see, Fig. 1). The relativistic bulk motion of the plasma along the jet direction is accompanied by the ideal motion electric field $\mathbf{E} = -\mathbf{V}_j \times \mathbf{B}/c$, whose nonvanishing component is along the radial direction. Thus, the plasma is polarized: the core of the jet is negatively charged, while the sheath is positively charged (and vice versa if $\Omega \cdot B_p < 0$). The motion of this net charge density $\rho$ carries the required current. Although polarized, the plasma remains quasi-neutral: $|\rho|/e \ll n_0$, where $e$ is the electron charge and $n_0$ is the reference plasma density in the jet.

The simulation box is limited to a two-dimensional Cartesian plane (xz coordinates) mimicking the poloidal plane of a cylindrical jet (Rz coordinates). With this choice of coordinates, the expansion of the flow with distance from the jet axis and the drift motion of the plasma associated with the curvature of field lines are missing (Huang et al. 2023). Yet, the Cartesian geometry has the advantage of capturing some of the nonaxisymmetric features that would be missing in a two-dimensional cylindrical simulation. The departure from perfect axisymmetry plays an important role in the escape mechanism of energetic particles as shown below. We defer the study of full three-dimensional simulations of a cylindrical jet to a future work.

The fiducial run is composed of 16 384 × 262 144 cells along the x and z directions, respectively. The box is elongated along the jet direction, it extends between $-4R_j < x < -4R_j$ and $0 < z < 128R_j$. The reference electron plasma skindepth, $d_e = (\Gamma_j m_e c^2 / 4\pi n_0 e^2)^{1/2}$, has been resolved by eight cells and the ion skindepth scale $d_i = d_e (m_i / m_e)^{1/2}$ by 40 cells, meaning that the box size is $410 d_i \times 6554 d_i$. The simulation time step has been fixed at $\omega_{pe} \Delta t = 8.75 \times 10^{-2}$, where $\omega_{pe} = c/d_e$ is the

<figure><img src="image_3.png" /><figcaption>Fig. 1. Transverse magnetic field strength (top panel, $B_y$) and current density (bottom panel, $J_z$) profiles carried by the jet in the upstream medium projected into the xz plane, where $\Delta = 0.5R_j$.</figcaption></figure>

electronic plasma frequency. All boundaries are perfectly reflective for both the fields and the particles. The $z = z_{max}$ boundary mimics the contact discontinuity with the ambient medium that is not simulated here to save on computational resources. Another common optimization procedure is to inject the plasma into the box gradually. Initially, the particle injector is close to the $z_{max}$ wall, and then it recedes away from this boundary at the speed of light. The injected plasma bounces off the wall with no loss of energy such that the beam interacts with itself and the shock forms. The choice of reflecting boundaries on the $\pm x$ directions mimics the presence of a perfectly confining medium surrounding the jet. At injection, the plasma was modeled by eight particles in each cell with a constant density at all radii, including the region outside the jet boundaries. The particles are streaming at the jet bulk Lorentz factor, but a small temperature $T/m_{e,i}c^2 = 0.01$ has been included in addition to spatial filtering of the current densities to limit the growth of spurious Cherenkov radiation (Greenwood et al. 2004).

An initially Poynting-flux-dominated jet gradually converts its magnetic energy into plasma bulk motion, such that if internal dissipation along its way to the termination shock is weak as in Fanaroff-Riley class II jets (Fanaroff & Riley 1974), the magnetic field and the particles could be close to equipartition (Li et al. 1992; Komissarov et al. 2007). In this work, we assume perfect equipartition which translates into the magnetic enthalpy density to the particle enthalpy density ratio, $\sigma \equiv B_0^2 / (4\pi\Gamma_j n_0 m_i c^2) = 1$ (the electronic contribution is neglected in this expression due to the mass ratio difference). We also ran the same numerical setup with $\sigma = 0.1$ and found similar results. We focus this study on the case where the poloidal magnetic component is negligible, $B_z = 0$. In reality, it is likely that a small but finite poloidal component still exists and stabilizes the jet against current driven instabilities such as the kink mode which could lead to its disruption akin to a Z-pinch configuration
<table><caption>**Table 1.** All the simulations reported in this work.</caption><thead><tr><th>Run</th><th>Cells</th><th>Size (in d<sub>i</sub>)</th><th>σ</th><th>B<sub>z</sub>/B<sub>0</sub></th></tr></thead><tbody><tr><td>Fiducial</td><td>16 384 × 262 144</td><td>410 × 6554</td><td>1</td><td>0</td></tr><tr><td>S01</td><td>8192 × 65 536</td><td>205 × 1638</td><td>0.1</td><td>0</td></tr><tr><td>Bz0</td><td>8192 × 65 536</td><td>205 × 1638</td><td>1</td><td>0</td></tr><tr><td>Bz001</td><td>8192 × 65 536</td><td>205 × 1638</td><td>1</td><td>0.01</td></tr><tr><td>Bz01</td><td>8192 × 65 536</td><td>205 × 1638</td><td>1</td><td>0.1</td></tr></tbody></table>

in full 3D (e.g., Kruskal & Schwarzschild 1954; Begelman 1998; Lyubarskii 1999; Mignone et al. 2010). We performed two additional simulations with $B_z/B_0 = 0.01$ and $B_z/B_0 = 0.1$, but we did not find significant differences with respect to the fiducial run presented below. Table 1 summarizes all the runs performed in this work.

# 3. Results

## 3.1. Shock structure and evolution

Figure 2 shows the final state of the fiducial simulation. After the initial reflection of the beam, the flow quickly reconfigured and produced a shock and the downstream medium. For $|x| > R_j$ where $\sigma \approx 0$, a Weibel-mediated shock formed, which is in agreement with previous studies of unmagnetized relativistic shocks (Nishikawa et al. 2003; Sironi et al. 2013). This solution is characterized by the growth of multiple thin plasma filaments near the shock, in both the upstream and the downstream media, generated by particles reflected at the shock front. Within the jet radius, the flow behaves very differently. The plasma collapses on the axis where the magnetic field vanishes due to the partial loss of the support from the jet electric force. Near the edges, the shock becomes gradually more magnetized ($\sigma \sim 1$) and therefore weaker (Kennel & Coroniti 1984). The most striking feature is the formation of an over-pressured cavity at the shock front along the jet axis that is drilling through the upstream medium. A similar pattern has recently been reported in the context of pulsar-wind termination shocks in pair plasmas (Cerutti & Giacinti 2020), but this phenomenon differs in this work by its extreme magnitude. The lateral size of this cavity does indeed increase with time up to the point where it fills the full radial extent of the jet. The cavity is under-dense (~0.1–0.2$n_0$), weakly magnetized, and filled with energetic particles (see Fig. 2). At this stage, the Weibel-mediated shock that prevailed in the early phases of the simulation has nearly completely disappeared because of the overwhelming presence of the cavity. In fact, it becomes even difficult to define a clear shock front in the simulation. It can be approximately located at the base of the cavity where the flow significantly decelerates and the plasma density increases by a factor $\approx 3$ at $z \approx 105R_j$ in Fig. 2.

Another important property is the presence of a strong backflow in the cavity moving relativistically toward the upstream medium against the jet at a velocity of $V_z \approx -0.75c$. A sharp tangential velocity discontinuity between the edges of the cavity and the incoming flow drives the formation of strong Kelvin-Helmholtz vortices with a size on the order of the width of the cavity itself at any given time. The downstream medium close to the shock front is dominated by the dynamics of an oscillating wake trailing behind the cavity. The wake periodically carves part of the cavity out, which is then advected further downstream in the form of a series of hot, low-density bubbles (Fig. 2), resembling a von Kármán vortex street. Their position is nearly

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Structure of the jet termination shock in the final state of the fiducial simulation at $\omega_{pi}t = 2500$. From left to right, this figure shows the following: the plasma density $n/n_0$, the out-of-plane magnetic field component ($B_y/B_0$), the axial plasma bulk velocity ($V_z/c$) and streamlines (solid contour lines with arrows), as well as the mean individual ion Lorentz factor $\langle \gamma_i \rangle$. The lower panels show the radial profile of each quantity across the large plasma cavity at $z \approx 95R_j$ (blue solid line), which can be compared with the upstream solution (red dashed line).</figcaption></figure>

at rest with respect to the simulation frame (see top panels in Fig. 3), but the flow inside the bubbles has a strong velocity vorticity which in turn leads to efficient mixing and its progressive obliteration in the downstream medium. We identify this phenomenon as the main escape mechanism which allows energetic particles momentarily stored in the cavity to move away from the acceleration zone. The shock front cavity and vortices also form in the presence of a weak vertical magnetic field component, $B_z/B_0 = 0.01, 0.1$; we observe a stronger plasma compression in the downstream medium with increasing $B_z$ (top panels in Fig. 4).

## 3.2. Particle acceleration

In contrast to plane-parallel uniform magnetized shocks, the strong variations of the transverse magnetic field strength and, in particular, the presence of a magnetic null line along the jet axis
<table><caption>**Table 1.** All the simulations reported in this work.</caption><thead><tr><th>Run</th><th>Cells</th><th>Size (in d<sub>i</sub>)</th><th>σ</th><th>B<sub>z</sub>/B<sub>0</sub></th></tr></thead><tbody><tr><td>Fiducial</td><td>16 384 × 262 144</td><td>410 × 6554</td><td>1</td><td>0</td></tr><tr><td>S01</td><td>8192 × 65 536</td><td>205 × 1638</td><td>0.1</td><td>0</td></tr><tr><td>Bz0</td><td>8192 × 65 536</td><td>205 × 1638</td><td>1</td><td>0</td></tr><tr><td>Bz001</td><td>8192 × 65 536</td><td>205 × 1638</td><td>1</td><td>0.01</td></tr><tr><td>Bz01</td><td>8192 × 65 536</td><td>205 × 1638</td><td>1</td><td>0.1</td></tr></tbody></table>

in full 3D (e.g., Kruskal & Schwarzschild 1954; Begelman 1998; Lyubarskii 1999; Mignone et al. 2010). We performed two additional simulations with $B_z/B_0 = 0.01$ and $B_z/B_0 = 0.1$, but we did not find significant differences with respect to the fiducial run presented below. Table 1 summarizes all the runs performed in this work.

# 3. Results

## 3.1. Shock structure and evolution

Figure 2 shows the final state of the fiducial simulation. After the initial reflection of the beam, the flow quickly reconfigured and produced a shock and the downstream medium. For $|x| > R_j$ where $\sigma \approx 0$, a Weibel-mediated shock formed, which is in agreement with previous studies of unmagnetized relativistic shocks (Nishikawa et al. 2003; Sironi et al. 2013). This solution is characterized by the growth of multiple thin plasma filaments near the shock, in both the upstream and the downstream media, generated by particles reflected at the shock front. Within the jet radius, the flow behaves very differently. The plasma collapses on the axis where the magnetic field vanishes due to the partial loss of the support from the jet electric force. Near the edges, the shock becomes gradually more magnetized ($\sigma \sim 1$) and therefore weaker (Kennel & Coroniti 1984). The most striking feature is the formation of an over-pressured cavity at the shock front along the jet axis that is drilling through the upstream medium. A similar pattern has recently been reported in the context of pulsar-wind termination shocks in pair plasmas (Cerutti & Giacinti 2020), but this phenomenon differs in this work by its extreme magnitude. The lateral size of this cavity does indeed increase with time up to the point where it fills the full radial extent of the jet. The cavity is under-dense (~0.1–0.2$n_0$), weakly magnetized, and filled with energetic particles (see Fig. 2). At this stage, the Weibel-mediated shock that prevailed in the early phases of the simulation has nearly completely disappeared because of the overwhelming presence of the cavity. In fact, it becomes even difficult to define a clear shock front in the simulation. It can be approximately located at the base of the cavity where the flow significantly decelerates and the plasma density increases by a factor $\approx 3$ at $z \approx 105R_j$ in Fig. 2.

Another important property is the presence of a strong backflow in the cavity moving relativistically toward the upstream medium against the jet at a velocity of $V_z \approx -0.75c$. A sharp tangential velocity discontinuity between the edges of the cavity and the incoming flow drives the formation of strong Kelvin-Helmholtz vortices with a size on the order of the width of the cavity itself at any given time. The downstream medium close to the shock front is dominated by the dynamics of an oscillating wake trailing behind the cavity. The wake periodically carves part of the cavity out, which is then advected further downstream in the form of a series of hot, low-density bubbles (Fig. 2), resembling a von Kármán vortex street. Their position is nearly

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Structure of the jet termination shock in the final state of the fiducial simulation at $\omega_{pi}t = 2500$. From left to right, this figure shows the following: the plasma density $n/n_0$, the out-of-plane magnetic field component ($B_y/B_0$), the axial plasma bulk velocity ($V_z/c$) and streamlines (solid contour lines with arrows), as well as the mean individual ion Lorentz factor $\langle \gamma_i \rangle$. The lower panels show the radial profile of each quantity across the large plasma cavity at $z \approx 95R_j$ (blue solid line), which can be compared with the upstream solution (red dashed line).</figcaption></figure>

at rest with respect to the simulation frame (see top panels in Fig. 3), but the flow inside the bubbles has a strong velocity vorticity which in turn leads to efficient mixing and its progressive obliteration in the downstream medium. We identify this phenomenon as the main escape mechanism which allows energetic particles momentarily stored in the cavity to move away from the acceleration zone. The shock front cavity and vortices also form in the presence of a weak vertical magnetic field component, $B_z/B_0 = 0.01, 0.1$; we observe a stronger plasma compression in the downstream medium with increasing $B_z$ (top panels in Fig. 4).

## 3.2. Particle acceleration

In contrast to plane-parallel uniform magnetized shocks, the strong variations of the transverse magnetic field strength and, in particular, the presence of a magnetic null line along the jet axis
<figure><img src="image_3.png" /><figcaption>Fig. 3. Zoomed-in view and time evolution of the von Kármán vortex street. The top panels shows the formation of a vortex that was carved out of the shock front cavity, and its later evolution in the downstream medium. The bottom panel shows the evolution of the ion spectrum within the vortex encapsulated in the box drawn in the upper panels. The total final ion spectrum is shown for comparison with the black dashed line.</figcaption></figure>

leads to an extremely efficient acceleration of particles. Figure 5 shows that the total particle spectrum is composed of a low-energy thermal bump centered around the upstream jet Lorentz factor $\Gamma_j$, and a high-energy power-law tail with an index $\approx -2.2$ extending up to the accelerator confinement limit defined by the Hillas criterion (Hillas 1984),

$$ \gamma_H \equiv \frac{2R_j e B_0}{m_i c^2} \approx 1000, \qquad (3) $$

in the fiducial simulation. A clear sign that the limit of the accelerator has been reached is the saturation of the maximum ion spectrum cut-off Lorentz factor, $\gamma_{i,max}$, as it approaches $\gamma_H$. This saturation is visible as a small bump in the spectrum just below the high-energy cutoff. In the early phases, $\gamma_{i,max}$ grows approximately linearly with time, which is compatible with an ideal Bohm-like acceleration regime. The full evolution is well

<figure><img src="image_4.png" /><figcaption>Fig. 4. Termination shock structure (top panels: plasma density) and total ion spectrum (bottom panel) as a function of the strength of the vertical magnetic field component $B_z/B_0 = 0, 0.01$ and $0.1$ at $\omega_{pi}t = 1575$.</figcaption></figure>

approximated by the following empirical formula

$$ \gamma_{i,max}(t) \approx 1 + \gamma_H \tanh\left(\frac{2\omega_{pi}t}{3\gamma_H}\right). \qquad (4) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Zoomed-in view and time evolution of the von Kármán vortex street. The top panels shows the formation of a vortex that was carved out of the shock front cavity, and its later evolution in the downstream medium. The bottom panel shows the evolution of the ion spectrum within the vortex encapsulated in the box drawn in the upper panels. The total final ion spectrum is shown for comparison with the black dashed line.</figcaption></figure>

leads to an extremely efficient acceleration of particles. Figure 5 shows that the total particle spectrum is composed of a low-energy thermal bump centered around the upstream jet Lorentz factor $\Gamma_j$, and a high-energy power-law tail with an index $\approx -2.2$ extending up to the accelerator confinement limit defined by the Hillas criterion (Hillas 1984),

$$ \gamma_H \equiv \frac{2R_j e B_0}{m_i c^2} \approx 1000, \qquad (3) $$

in the fiducial simulation. A clear sign that the limit of the accelerator has been reached is the saturation of the maximum ion spectrum cut-off Lorentz factor, $\gamma_{i,max}$, as it approaches $\gamma_H$. This saturation is visible as a small bump in the spectrum just below the high-energy cutoff. In the early phases, $\gamma_{i,max}$ grows approximately linearly with time, which is compatible with an ideal Bohm-like acceleration regime. The full evolution is well

<figure><img src="image_6.png" /><figcaption>Fig. 4. Termination shock structure (top panels: plasma density) and total ion spectrum (bottom panel) as a function of the strength of the vertical magnetic field component $B_z/B_0 = 0, 0.01$ and $0.1$ at $\omega_{pi}t = 1575$.</figcaption></figure>

approximated by the following empirical formula

$$ \gamma_{i,max}(t) \approx 1 + \gamma_H \tanh\left(\frac{2\omega_{pi}t}{3\gamma_H}\right). \qquad (4) $$
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Particle acceleration up to the confinement limit of the system. Top panel: time evolution of the total ion spectrum $\gamma_i dN/d\gamma_i$ (color-coded solid lines). The final total electron spectrum is shown for comparison as a function of $\gamma_e \beta_e \times (m_e/m_i)$ to compare the energy scales rather than the Lorentz factors between the two species (black dashed line). A -2.2 power-law slope is shown for comparison (dotted line). Bottom panel: time evolution of the maximum ion Lorentz factor $\gamma_{i,max}$ (blue solid line) and the latitudinal width of the shock front cavity normalized to the full jet width ($W_c$, red circles). The jet confinement limit, $\gamma_H$, is marked by the horizontal dotted line. The blue dashed line shows a simple empirical formula that captures the time evolution of $\gamma_{i,max}$.</figcaption></figure>

Figure 5 also demonstrates the perfect coevolution of the maximum particle energy with the width of the cavity, establishing an unambiguous connection between particle acceleration and the growth of the cavity whose dynamics is governed by the cosmic-ray pressure. The final electron spectrum has a very similar shape as the ion spectrum, but it cuts off at slightly lower energies, $\gamma_{e,max} \times (m_e/m_i) \approx 300$. The asymmetry between high-energy electrons and ions is explained by the slightly different acceleration patterns followed by each species and may depend on the polarization of the jet (i.e., the sign of $B_p \cdot \Omega$, Giacinti & Kirk 2018; Cerutti & Giacinti 2020). This effect may also be weaker at higher scale separation. We observe a decrease in the particle acceleration efficiency with increasing strength of the vertical field component (Fig. 4, bottom panel).

To elucidate the origin of this extreme particle acceleration, we tracked a sample of 20 000 simulation particles for each species. They have been randomly chosen throughout the simulation box to avoid any selection biases. The high-energy ions that reach the confinement limit all follow a very similar acceleration history as shown in Fig. 6. The time evolution of the ion energy is characterized by discrete and intense acceleration

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Acceleration history of four representative high-energy ion trajectories tracked from their injection at $\omega_{pi}t = 175$ up to $\omega_{pi}t \gtrsim 2000$ where they have all reached the jet confinement limit. Top: full ion trajectory superimposed onto the final plasma density map (gray scale). The color indicates the particle Lorentz factor $\gamma_i$, while the marker size codes for the acceleration rate $\dot{\gamma}_i/\omega_0$, where $\omega_0 = eB_0/m_i c$, a large marker means a high positive rate while a small marker is for negative rates. Bottom: time evolution of the ion Lorentz factor (blue solid line) and acceleration rate (red dashed line).</figcaption></figure>

episodes during which the particle Lorentz factor grows almost linearly. Far from the saturation, the energy gain during each step is close to $\Delta\gamma \sim \gamma$. The accelerate rate $\dot{\gamma}$ is also independent of the particle energy. We find that each of these acceleration events are unambiguously related to the crossing of velocity shear layers. The first energy boost occurs when the particle crosses the main shear layer at the interface between the upstream medium and the shock front cavity. Once they are captured by the cavity, they are further accelerated in the wake trailing behind.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Particle acceleration up to the confinement limit of the system. Top panel: time evolution of the total ion spectrum $\gamma_i dN/d\gamma_i$ (color-coded solid lines). The final total electron spectrum is shown for comparison as a function of $\gamma_e \beta_e \times (m_e/m_i)$ to compare the energy scales rather than the Lorentz factors between the two species (black dashed line). A -2.2 power-law slope is shown for comparison (dotted line). Bottom panel: time evolution of the maximum ion Lorentz factor $\gamma_{i,max}$ (blue solid line) and the latitudinal width of the shock front cavity normalized to the full jet width ($W_c$, red circles). The jet confinement limit, $\gamma_H$, is marked by the horizontal dotted line. The blue dashed line shows a simple empirical formula that captures the time evolution of $\gamma_{i,max}$.</figcaption></figure>

Figure 5 also demonstrates the perfect coevolution of the maximum particle energy with the width of the cavity, establishing an unambiguous connection between particle acceleration and the growth of the cavity whose dynamics is governed by the cosmic-ray pressure. The final electron spectrum has a very similar shape as the ion spectrum, but it cuts off at slightly lower energies, $\gamma_{e,max} \times (m_e/m_i) \approx 300$. The asymmetry between high-energy electrons and ions is explained by the slightly different acceleration patterns followed by each species and may depend on the polarization of the jet (i.e., the sign of $B_p \cdot \Omega$, Giacinti & Kirk 2018; Cerutti & Giacinti 2020). This effect may also be weaker at higher scale separation. We observe a decrease in the particle acceleration efficiency with increasing strength of the vertical field component (Fig. 4, bottom panel).

To elucidate the origin of this extreme particle acceleration, we tracked a sample of 20 000 simulation particles for each species. They have been randomly chosen throughout the simulation box to avoid any selection biases. The high-energy ions that reach the confinement limit all follow a very similar acceleration history as shown in Fig. 6. The time evolution of the ion energy is characterized by discrete and intense acceleration

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Acceleration history of four representative high-energy ion trajectories tracked from their injection at $\omega_{pi}t = 175$ up to $\omega_{pi}t \gtrsim 2000$ where they have all reached the jet confinement limit. Top: full ion trajectory superimposed onto the final plasma density map (gray scale). The color indicates the particle Lorentz factor $\gamma_i$, while the marker size codes for the acceleration rate $\dot{\gamma}_i/\omega_0$, where $\omega_0 = eB_0/m_i c$, a large marker means a high positive rate while a small marker is for negative rates. Bottom: time evolution of the ion Lorentz factor (blue solid line) and acceleration rate (red dashed line).</figcaption></figure>

episodes during which the particle Lorentz factor grows almost linearly. Far from the saturation, the energy gain during each step is close to $\Delta\gamma \sim \gamma$. The accelerate rate $\dot{\gamma}$ is also independent of the particle energy. We find that each of these acceleration events are unambiguously related to the crossing of velocity shear layers. The first energy boost occurs when the particle crosses the main shear layer at the interface between the upstream medium and the shock front cavity. Once they are captured by the cavity, they are further accelerated in the wake trailing behind.
The coherent oscillating motion of the wake allows for multiple crossings of the velocity shear layer. An acceleration episode is often preceded by a brief deceleration event associated with the sudden large-amplitude change in the direction of the momentum that the particle experiences as it encounters an oppositely directed flow.

The particles shown here reach an energy saturation and remain in the cavity until the end of the simulation without significant further energization. Looking at other high-energy particles shows that, eventually, they are captured in a bubble and advected in the downstream medium without a significant loss of energy, which allows them to escape the acceleration site. In fact, the spectrum measured within a von Kármán vortex presents two components soon after its formation: a low-energy bump and a narrow high-energy component cutting off near the confinement limit. At later times when the vortex is further downstream, the spectrum evolves to a single power-law distribution with an index ~-2.2 without significant changes in the high-energy cut-off (see bottom panels in Fig. 3).

Particles are accelerated by experiencing the large-scale ideal electric field associated with the varying bulk motion of the flow across the velocity shear layers, such that the acceleration rate is $\dot{\gamma} \sim e\mathbf{E} \cdot \mathbf{V}/m_i c^2 \lesssim 0.5\omega_0$, where $\omega_0 = eB_0/m_i c$. This explains the high acceleration efficiency, and is in stark contrast with the microscopic turbulent nature of the scattering centers that is expected in relativistic Weibel-dominated shocks. Particle acceleration at a shear layer can be understood as a Lorentz boost of the particle energy as it is scattered by plasma turbulence frozen into the flow approximately moving at the $\mathbf{E} \times \mathbf{B}$-drift velocity. For an isotropic particle distribution near the shear layer, the mean energy gain per crossing is $\langle \Delta\gamma/\gamma \rangle = \Gamma_s - 1$, where $\Gamma_s = 1/(1 - V_s^2/c^2)^{1/2}$ and $V_s$ is the relative velocity across the layer (Ostrowski 1990; Rieger & Duffy 2004). With typical flow velocities $V_\pm \sim \pm 0.5c$ across the boundary gives $V_s = (V_+ - V_-)/(1 - V_+V_-/c^2) \sim 0.8c$, so a net energy gain per crossing $\langle \Delta\gamma/\gamma \rangle \sim 1$ in agreement with our findings. The $\approx -2$ slope of the nonthermal power-law tail in the particle spectrum can also be well accounted for by a simple Fokker-Planck model assuming a particle scattering time on the order of the Larmor gyration time (Ostrowski 1998; Rieger & Duffy 2006).

## 4. Conclusion

The toroidal nature of the magnetic field within the jet leads to the formation of an over-pressured cavity filled with cosmic rays and a relativistic shear flow localized near the shock front. These features play a key role in the extreme particle acceleration mechanism reported in this work. Particles are primarily accelerated at the boundaries of the cavity, where the tangential velocity discontinuity is strongest, up to the confinement limit of the system. The cavity acts as an obstacle to the incoming flow and generates a strong von Kármán vortex street-like structure that allows energetic cosmic rays to escape in the downstream medium. The magnetic field inside the cavity is weak such that the high-energy particles it contains should not radiate much via synchrotron emission, in contrast to the downstream medium. Therefore, the cavity may be observable as a dark region in the vicinity of the jet hotspots. Interestingly, a hole in the radio and X-ray images of Cygnus A jets has been recently reported (Snios et al. 2020) and, possibly, also in Pictor A jets (Thimmappa et al. 2022), it would be worthwhile to check whether this feature is ubiquitous among other similar jets as predicted by the model.

Scaling our results up to the typical physical conditions of extragalactic jet hotspots and lobes, assuming a $B_0 \sim 10$ µG, a $L \sim 10$ kpc size yields a maximum proton energy at saturation on the order of $\mathcal{E}_{p,max} = \gamma_{p,max} m_p c^2 \sim eB_0 L \sim 10^{20}$ eV. Thus, the findings of this work support the idea that not only are extragalactic jets able to confine ultra-high-energy cosmic rays, but they may also be capable of generating them. The configuration explored in this study is generic, and therefore the conclusions are applicable to other astrophysical systems involving an organized large-scale toroidal magnetic field. In another work, we argue that a similar mechanism could operate in relativistic pair plasmas at the shock front of pulsar wind nebulae, and explain extreme particle acceleration up to PeV energies (Cerutti & Giacinti 2020). The termination shock of stellar-mass black hole jets such as in the microquasar SS433 (Fabrika 2004) is another relevant astrophysical environment where this mechanism could take place. The recent discovery of $>10$ TeV gamma rays from the lobes spatially coincident with compact X-ray knots suggests that extreme in situ particle acceleration close to PeV energies must take place (Abeysekara et al. 2018), which is an outstanding challenge that may be overcome by the mechanism highlighted in this work.

**Acknowledgements.** We thank B. Crinquand and M. Ostrowski for their comments on an early draft, and the referee Luke Drury for his comments and supportive report. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 863412). Computing resources were provided by TGCC and CINES under the allocation A0130407669 made by GENCI.

## References

Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2018, Nature, 562, 82

Alves, E. P., Grismayer, T., Fonseca, R. A., & Silva, L. O. 2014, New J. Phys., 16, 035007

Begelman, M. C. 1998, ApJ, 493, 291

Begelman, M. C., & Kirk, J. G. 1990, ApJ, 353, 66

Bell, A. R., Araudo, A. T., Matthews, J. H., & Blundell, K. M. 2018, MNRAS, 473, 2364

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433

Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467

Cerutti, B., & Giacinti, G. 2020, A&A, 642, A123

Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147

Fabrika, S. 2004, Astrophys. Space Phys. Res., 12, 1

Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P

Gabuzda, D. C. 2021, Galaxies, 9, 58

Giacinti, G., & Kirk, J. G. 2018, ApJ, 863, 18

Giannios, D. 2010, MNRAS, 408, L46

Greenwood, A. D., Cartwright, K. L., Luginsland, J. W., & Baca, E. A. 2004, J. Comput. Phys., 201, 665

Hardcastle, M. J., & Croston, J. H. 2020, New A Rev., 88, 101539

Hillas, A. M. 1984, ARA&A, 22, 425

Huang, Z.-Q., Reville, B., Kirk, J. G., & Giacinti, G. 2023, MNRAS, 522, 4955

Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 694

Kimura, S. S., Murase, K., & Zhang, B. T. 2018, Phys. Rev. D, 97, 023026

Kirk, J. G., Reville, B., & Huang, Z.-Q. 2023, MNRAS, 519, 1022

Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, MNRAS, 380, 51

Kruskal, M., & Schwarzschild, M. 1954, Proc. R. Soc. London Ser. A, 223, 348

Li, Z.-Y., Chiueh, T., & Begelman, M. C. 1992, ApJ, 394, 459

Lyubarskii, Y. E. 1999, MNRAS, 308, 1006

Matthews, J. H., Bell, A. R., Blundell, K. M., & Araudo, A. T. 2018, MNRAS, 479, L76

Matthews, J. H., Bell, A. R., Blundell, K. M., & Araudo, A. T. 2019, MNRAS, 482, 4303

Meli, A., Nishikawa, K., Köhn, C., et al. 2023, MNRAS, 519, 5410

Michel, F. C. 1973, ApJ, 180, L133

Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010, MNRAS, 402, 7

Nishikawa, K. I., Hardee, P., Richardson, G., et al. 2003, ApJ, 595, 555
The coherent oscillating motion of the wake allows for multiple crossings of the velocity shear layer. An acceleration episode is often preceded by a brief deceleration event associated with the sudden large-amplitude change in the direction of the momentum that the particle experiences as it encounters an oppositely directed flow.

The particles shown here reach an energy saturation and remain in the cavity until the end of the simulation without significant further energization. Looking at other high-energy particles shows that, eventually, they are captured in a bubble and advected in the downstream medium without a significant loss of energy, which allows them to escape the acceleration site. In fact, the spectrum measured within a von Kármán vortex presents two components soon after its formation: a low-energy bump and a narrow high-energy component cutting off near the confinement limit. At later times when the vortex is further downstream, the spectrum evolves to a single power-law distribution with an index ~-2.2 without significant changes in the high-energy cut-off (see bottom panels in Fig. 3).

Particles are accelerated by experiencing the large-scale ideal electric field associated with the varying bulk motion of the flow across the velocity shear layers, such that the acceleration rate is $\dot{\gamma} \sim e\mathbf{E} \cdot \mathbf{V}/m_i c^2 \lesssim 0.5\omega_0$, where $\omega_0 = eB_0/m_i c$. This explains the high acceleration efficiency, and is in stark contrast with the microscopic turbulent nature of the scattering centers that is expected in relativistic Weibel-dominated shocks. Particle acceleration at a shear layer can be understood as a Lorentz boost of the particle energy as it is scattered by plasma turbulence frozen into the flow approximately moving at the $E \times B$-drift velocity. For an isotropic particle distribution near the shear layer, the mean energy gain per crossing is $\langle \Delta\gamma/\gamma \rangle = \Gamma_s - 1$, where $\Gamma_s = 1/(1 - V_s^2/c^2)^{1/2}$ and $V_s$ is the relative velocity across the layer (Ostrowski 1990; Rieger & Duffy 2004). With typical flow velocities $V_\pm \sim \pm 0.5c$ across the boundary gives $V_s = (V_+ - V_-)/(1 - V_+V_-/c^2) \sim 0.8c$, so a net energy gain per crossing $\langle \Delta\gamma/\gamma \rangle \sim 1$ in agreement with our findings. The $\approx -2$ slope of the nonthermal power-law tail in the particle spectrum can also be well accounted for by a simple Fokker-Planck model assuming a particle scattering time on the order of the Larmor gyration time (Ostrowski 1998; Rieger & Duffy 2006).

## 4. Conclusion

The toroidal nature of the magnetic field within the jet leads to the formation of an over-pressured cavity filled with cosmic rays and a relativistic shear flow localized near the shock front. These features play a key role in the extreme particle acceleration mechanism reported in this work. Particles are primarily accelerated at the boundaries of the cavity, where the tangential velocity discontinuity is strongest, up to the confinement limit of the system. The cavity acts as an obstacle to the incoming flow and generates a strong von Kármán vortex street-like structure that allows energetic cosmic rays to escape in the downstream medium. The magnetic field inside the cavity is weak such that the high-energy particles it contains should not radiate much via synchrotron emission, in contrast to the downstream medium. Therefore, the cavity may be observable as a dark region in the vicinity of the jet hotspots. Interestingly, a hole in the radio and X-ray images of Cygnus A jets has been recently reported (Snios et al. 2020) and, possibly, also in Pictor A jets (Thimmappa et al. 2022), it would be worthwhile to check whether this feature is ubiquitous among other similar jets as predicted by the model.

Scaling our results up to the typical physical conditions of extragalactic jet hotspots and lobes, assuming a $B_0 \sim 10$ µG, a $L \sim 10$ kpc size yields a maximum proton energy at saturation on the order of $E_{p,max} = \gamma_{p,max}m_p c^2 \sim eB_0 L \sim 10^{20}$ eV. Thus, the findings of this work support the idea that not only are extragalactic jets able to confine ultra-high-energy cosmic rays, but they may also be capable of generating them. The configuration explored in this study is generic, and therefore the conclusions are applicable to other astrophysical systems involving an organized large-scale toroidal magnetic field. In another work, we argue that a similar mechanism could operate in relativistic pair plasmas at the shock front of pulsar wind nebulae, and explain extreme particle acceleration up to PeV energies (Cerutti & Giacinti 2020). The termination shock of stellar-mass black hole jets such as in the microquasar SS433 (Fabrika 2004) is another relevant astrophysical environment where this mechanism could take place. The recent discovery of $>10$ TeV gamma rays from the lobes spatially coincident with compact X-ray knots suggests that extreme in situ particle acceleration close to PeV energies must take place (Abeysekara et al. 2018), which is an outstanding challenge that may be overcome by the mechanism highlighted in this work.

**Acknowledgements.** We thank B. Crinquand and M. Ostrowski for their comments on an early draft, and the referee Luke Drury for his comments and supportive report. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 863412). Computing resources were provided by TGCC and CINES under the allocation A0130407669 made by GENCI.

## References

Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2018, Nature, 562, 82

Alves, E. P., Grismayer, T., Fonseca, R. A., & Silva, L. O. 2014, New J. Phys., 16, 035007

Begelman, M. C. 1998, ApJ, 493, 291

Begelman, M. C., & Kirk, J. G. 1990, ApJ, 353, 66

Bell, A. R., Araudo, A. T., Matthews, J. H., & Blundell, K. M. 2018, MNRAS, 473, 2364

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433

Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467

Cerutti, B., & Giacinti, G. 2020, A&A, 642, A123

Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147

Fabrika, S. 2004, Astrophys. Space Phys. Res., 12, 1

Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P

Gabuzda, D. C. 2021, Galaxies, 9, 58

Giacinti, G., & Kirk, J. G. 2018, ApJ, 863, 18

Giannios, D. 2010, MNRAS, 408, L46

Greenwood, A. D., Cartwright, K. L., Luginsland, J. W., & Baca, E. A. 2004, J. Comput. Phys., 201, 665

Hardcastle, M. J., & Croston, J. H. 2020, New A Rev., 88, 101539

Hillas, A. M. 1984, ARA&A, 22, 425

Huang, Z.-Q., Reville, B., Kirk, J. G., & Giacinti, G. 2023, MNRAS, 522, 4955

Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 694

Kimura, S. S., Murase, K., & Zhang, B. T. 2018, Phys. Rev. D, 97, 023026

Kirk, J. G., Reville, B., & Huang, Z.-Q. 2023, MNRAS, 519, 1022

Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, MNRAS, 380, 51

Kruskal, M., & Schwarzschild, M. 1954, Proc. R. Soc. London Ser. A, 223, 348

Li, Z.-Y., Chiueh, T., & Begelman, M. C. 1992, ApJ, 394, 459

Lyubarskii, Y. E. 1999, MNRAS, 308, 1006

Matthews, J. H., Bell, A. R., Blundell, K. M., & Araudo, A. T. 2018, MNRAS, 479, L76

Matthews, J. H., Bell, A. R., Blundell, K. M., & Araudo, A. T. 2019, MNRAS, 482, 4303

Meli, A., Nishikawa, K., Köhn, C., et al. 2023, MNRAS, 519, 5410

Michel, F. C. 1973, ApJ, 180, L133

Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010, MNRAS, 402, 7

Nishikawa, K. I., Hardee, P., Richardson, G., et al. 2003, ApJ, 595, 555
Ostrowski, M. 1990, A&A, 238, 435

Ostrowski, M. 1998, A&A, 335, 134

Plotnikov, I., Grassi, A., & Grech, M. 2018, MNRAS, 477, 5238

Rieger, F. M. 2022, Universe, 8, 607

Rieger, F. M., & Duffy, P. 2004, ApJ, 617, 155

Rieger, F. M., & Duffy, P. 2006, ApJ, 652, 1044

Scharlemann, E. T., & Wagoner, R. V. 1973, ApJ, 182, 951

Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54

Sironi, L., Rowan, M. E., & Narayan, R. 2021, ApJ, 907, L44

Snios, B., Johnson, A. C., Nulsen, P. E. J., et al. 2020, ApJ, 891, 173

Thimmappa, R., Stawarz, Ł., Neilsen, J., Ostrowski, M., & Reville, B. 2022, ApJ, 941, 204

Wang, J.-S., Reville, B., Mizuno, Y., Rieger, F. M., & Aharonian, F. A. 2023, MNRAS, 519, 1872

Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8
Ostrowski, M. 1990, A&A, 238, 435

Ostrowski, M. 1998, A&A, 335, 134

Plotnikov, I., Grassi, A., & Grech, M. 2018, MNRAS, 477, 5238

Rieger, F. M. 2022, Universe, 8, 607

Rieger, F. M., & Duffy, P. 2004, ApJ, 617, 155

Rieger, F. M., & Duffy, P. 2006, ApJ, 652, 1044

Scharlemann, E. T., & Wagoner, R. V. 1973, ApJ, 182, 951

Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54

Sironi, L., Rowan, M. E., & Narayan, R. 2021, ApJ, 907, L44

Snios, B., Johnson, A. C., Nulsen, P. E. J., et al. 2020, ApJ, 891, 173

Thimmappa, R., Stawarz, Ł., Neilsen, J., Ostrowski, M., & Reville, B. 2022, ApJ, 941, 204

Wang, J.-S., Reville, B., Mizuno, Y., Rieger, F. M., & Aharonian, F. A. 2023, MNRAS, 519, 1872

Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8