# Cosmic ray acceleration by multiple shocks in the jets of Active Galactic Nuclei

Ana Laura Müller<fnref n="1" />,\* and Anabella Araudo<fnref n="1" />,<fnref n="2" />,\*\*

**Abstract.** Active galactic nuclei are one of the most promising sources for accelerating particles up to the highest energies. In this contribution, we present a scenario in which cosmic rays are accelerated in multiple shocks created by the interaction of relativistic AGN jets with the winds of embedded massive stars. We solve the Fokker-Planck equation considering escape and radiative losses as well as the collective effect of the shocks and the reacceleration of the particles. Finally, we calculate the maximum energies that the particles can achieve and discuss the possibility of producing ultra-high energy cosmic rays in this astrophysical situation.

## 1 Introduction

Starburst galaxies and active galactic nuclei (AGNs) have been in the center of the discussion about the origin of ultra-high-energy cosmic rays (UHECRs) since the last years. Arguments in favour and against each of these astrophysical objects have been exposed by several authors based on theoretical and observational studies [1–3].

When proposing cosmic-ray sources, it is not only necessary to explain the high energies, but also the mass composition inferred from observations, which agrees with intermediate mass nuclei dominating at the higher energies [4, 5]. In the case of starburst galaxies, there is a large amount of massive stars and supernova explosions injecting heavy nuclei in the galactic disk, however, their strong photon fields enhance the photo-disintegration of these particles, and therefore, the superwind region in the halo was proposed as the most suitable place in these galaxies to accelerate particles up to the highest energies [1, 3]. Nevertheless, from the theoretical point of view to achieve energies of $10^{19} – 10^{20}$ eV has been demonstrated to be very unlikely there [see e.g., 6–9]. Besides that, starbursts do not satisfied the minimum energetic requirement to be sources of UHECRs [10, 11]. Jets from powerful radio-galaxies, on the other hand, fulfil the energetic requirements, but justifying the presence of atomic nuclei inside them is not trivial [see e.g., 12–18]. Moreover, the mechanism and efficiency of relativistic jets to accelerate particles are continuously a matter of debate.

The situation of a jet colliding with stellar winds was proposed a long time ago as a way to mass-loading jets with baryons [19]. Star-jet interactions were also suggested as a possible scenario to explain the bright X-ray knots observed in Centaurus A jet and as sources of

gamma-rays in AGNs [20–24]. Additionally, the wind of massive stars is very rich in intermediate-mass nuclei. In particular, the presence of oxygen, nitrogen, and carbon makes the atmosphere of Wolf-Rayet (WR) stars to have a composition similar to the one observed in cosmic rays [25, 26]. Motivated by these ideas, we propose in this work a scenario where cosmic rays are accelerated by multiple shocks. Those shocks are produced by the interaction of the relativistic jet with massive stars, as it is sketched in Fig. 1. In Section 2 we describe our scenario. Section 3 contains the description of the procedure that we follow to calculate the particle distribution of the reaccelerated particles. In Section 4 we present our results and finally, Section 5 is dedicated to our conclusions.

## 2 Model

When the jet collides with the wind of a massive star a pair of shock waves is created (see Fig. 1). A non-relativistic shock propagates in the wind with velocity ~ $v_{wind}$ and a relativistic shock propagates backwards into the jet with speed ~ $v_{jet}$. A fraction of the kinetic energy of the jet is dissipated in each interaction and, therefore, the jet slowly decelerates. We assume that its power evolves as $L_{jet}(z) = L_{jet}^0 e^{-\tau}$, where $L_{jet}^0$ is the initial jet kinetic luminosity and

$$ \tau = \int_{1 \text{ pc}}^{z} \left( \frac{\sigma_{\text{sp}}}{\sigma_{\text{j}}} \right) n_{\star}(z') \sigma_{\text{j}} \, dz' \quad (1) $$

where $n_{\star}(z)$ is the density of stars, and $\sigma_{\text{sp}}$ and $\sigma_{\text{j}}$ are the cross sections of the bow shock and the jet, respectively [22]. The lower limit of the integral is given by the minimum distance at which the first star can be located, we adopt $z = 1 \text{ pc}$ for this value. The size of the bow shock around the star is determined by the position of the stagna-

<footnote n="1">¹Extreme Light Infrastructure ERIC,ELI Beamlines Facility, Za Radnicí 835, CZ-25241 Dolní Břežany, Czech Republic</footnote>

<footnote n="2">²Laboratoire Univers et Particules de Montpellier (LUPM) Université Montpellier, CNRS/IN2P3, CC72, place Eugène Bataillon, 34095, Montpellier Cedex 5, France</footnote>

<footnote>\*e-mail: analaura.muller@eli-beams.eu</footnote>

<footnote>\*\*e-mail: anabella.araudo@eli-beams.eu</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Graphic representation of the scenario studied in this work. Several stars are considered inside the AGN jet, represented by the orange spheres. The interaction of the relativistic jet with the stellar winds produces a pair of bow shocks consisting of a pair of shocks, one propagating through the stellar wind and one in the jet. The characteristic size of the bow shocks is given by stagnation point located at a distance $R_{sp}(z)$, which increases with the distance $z$ as indicated in Eq. (3). Particles accelerated in a bow shock can be reaccelerated in another bow shocks at larger $z$.</figcaption></figure>

tion point $R_{sp}$. By equating the jet and wind kinetic pressures we obtain [20]

$$R_{sp}(z) = \sqrt{\frac{\dot{M}_{\star} v_{wind} v_{jet} R_{jet}^{2}(z)}{4 L_{jet}}}, \quad (2)$$

where $\dot{M}_{\star}$ is the mass-loss rate of the star, $R_{jet}(z) = z \tan \theta$ is the jet radius, and $\theta$ is the jet opening angle. Taking into account the effect of the deceleration, $R_{sp}(z)$ can be written in first approximation as

$$R_{sp}(z) = \frac{R_{sp}^{0}(z)}{\sqrt{1 - \left(\frac{R_{sp}^{0}(z)}{R_{jet}(z)}\right)^{2} N_{\star}(z)}}, \quad (3)$$

with $R_{sp}^{0}$ the location of the stagnation point when the jet does not decelerate and $N_{\star}(z) = \int_{1 \text{ pc}}^{z} n_{\star}(z') \pi R_{jet}(z')^{2} dz'$ is the cumulative number of stars into the jet. In Figure 2 we plot $N_{\star}$ using the stellar distribution from [22]

$$\frac{n_{\star}}{\text{pc}^{-3}} \sim 4.6 \times 10^{-4} (497)^{y} \eta_{\text{accr}}^{0.89} \left(\frac{M_{\text{BH}}}{10^{7} M_{\odot}}\right)^{0.89} \left(\frac{z}{\text{pc}}\right)^{-y} \quad (4)$$

with $M_{BH} = 10^{8} M_{\odot}$ the mass of the central black hole, $\eta_{accr} = 0.1$ the AGN accretion rate, and $y = 1$ and 2. The

opening angle of the jet is assumed to be $\theta = 5^{\circ}$. We can see in Fig 2 that for $y = 1$ one star is expected in the first 50 pc of the jet, ten stars at 100 pc, and fifty stars at 300 pc.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Cumulative number of stars inside a jet with opening angle $\theta = 5^{\circ}$ for a black hole of $M_{BH} = 10^{8} M_{\odot}$ and an accretion efficiency $\eta_{accr} = 0.1$. The $0.1 N_{\star}$ are considered to be WR stars.</figcaption></figure>

The stellar distributions in Eq. (4) correspond to the total distribution of massive stars. In the present contribution we are interested only in WR stars given that they have more powerful winds and a chemical composition rich in intermediate-mass nuclei, although the number of WR stars is expected to be lower than the total number of massive stars. We assume that the number of WR stars is $n_{WR} = 0.1 n_{\star}$. We adopt standard parameters for a WR star: $\dot{M}_{\star} = 10^{-4} M_{\odot} \text{yr}^{-1}$, luminosity $L_{\star} = 10^{38} \text{ erg s}^{-1}$, temperature $T_{\star} = 3 \times 10^{4} \text{ K}$, and $v_{wind} = 3000 \text{ km s}^{-1}$. For the jet we adopt $L_{jet}^{0} = 10^{42} \text{ erg s}^{-1}$, Lorentz factor $\Gamma_{0} = 5$, and $\theta = 5^{\circ}$. For calculating the deceleration of the jet we consider all the other stars to be OB stars with typical parameters $v_{wind} = 2000 \text{ km s}^{-1}$ and $\dot{M}_{\star} = 10^{-6} M_{\odot} \text{yr}^{-1}$. In Fig. 3 we plot $R_{sp}(z)$ (and $R_{jet}$ for comparison). Note that $R_{sp}$ increases significantly when the number of WR stars into the jet is much larger than 10 ($z > 180 \text{ pc}$ for $y = 1$ and $z > 63 \text{ pc}$ for $y = 2$). Additionally, the deceleration of the jet makes also the Lorentz factor to be lower than $\Gamma_{0}$. Its change can be found by using the definition $L_{jet}(z) = [\Gamma(z) - 1]\rho_{jet}^{0} c^{2} \pi R_{jet}^{2}(z) v_{jet}(z)$, with $\rho_{jet}^{0}$ the jet density, and solving the transcendental equation

$$(\Gamma_0 - 1) \sqrt{1 - \frac{1}{\Gamma_0^2}} e^{-\tau} = (\Gamma - 1) \sqrt{1 - \frac{1}{\Gamma^2}}. \quad (5)$$

In Fig. 4 we plot $\Gamma$ for the two stellar distributions considered in the present study. We can see that $\Gamma \sim \Gamma_{0}$ until $z \sim 60 \text{ pc}$ for $y = 2$ and $z \sim 200 \text{ pc}$ for $y = 1$, where $N_{WR} = 10$ is large enough to affect the dynamics of the jet and also $R_{sp}(z)$ (See Fig. 3).

The magnetic field in the jet is assumed to be $B_{jet}(z) = B_{0} z^{-1}$, where $B_{0} = \sqrt{4 L_{jet}/(R_{jet}^{2}(z_{0}) c)}$, with $z_{0} = 5 \times 10^{-5} (M_{BH}/10^{7} M_{\odot}) \text{ pc}$ the position where the jet is launched [22, 27, 28]. The magnetic field of the stellar wind is described as in [29]. Using these expressions,
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Comparison of $R_{sp}(z)$, $R_{sp}^0(z)$, and $R_{jet}(z)$. For the distribution of stars we adopt the expression given in Eq. 4 for $M_{BH} = 10^8 M_{\odot}$ and $\eta_{accr} = 0.1$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 5.** Hillas energy at the shocks along the jet considering the stellar population given by Eq. 4 and $Z = 1$ (protons).</figcaption></figure>

along with $R_{sp}(z)$ and $\Gamma(z)$ we estimate the Hillas energy at the shock in the jet as

$$ \frac{E_{\text{Hillas}}(z)}{\text{eV}} = 5.1 \times 10^{17} Z \left(\frac{v_{\text{jet}}(z)}{c}\right) \left(\frac{R_{\text{sp}}(z)}{5 \text{ pc}}\right) \left(\frac{B_{\text{jet}}(z)}{125 \text{ µG}}\right), \quad (6) $$

where Z is the atomic number [30]. Note that $R_{sp} = 5$ pc and $B_{jet} = 125$ µG at $z = 338$ pc for $y = 2$. In Fig. 5 we show that $E_{Hillas}$ increases as a result of the faster increment of $R_{sp}(z)$ with respect to the decrease of $B_{jet}$ and $\Gamma$. At $z \sim 1$ kpc, $E_{Hillas} = 10^{18}$ eV and $2.6 \times 10^{19}$ eV, for protons and iron nuclei, respectively. Nevertheless, $E_{Hillas}$ is an upper limit. The actual maximum energies expected at each shock will be determined by the radiative and escape losses, as we discuss in the next section.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Evolution of the Lorentz factor because of the deceleration of the jet for the stellar distribution given in Eq. 5 with $y = 1$ and $y = 2$.</figcaption></figure>

## 3 Multiple shock acceleration

Particles can be accelerated in the bow shocks created by the interaction of the jet with the winds of WR stars. Acceleration at single bow shocks has been considered before

in order to study the gamma-ray emission from misaligned AGN jets [20–22]. Because a large number of stars can be simultaneously into the jet, particles accelerated in a shock can be reaccelerated at another shock at larger $z$. We follow the next steps to obtain the particle distributions after multiple-shock interactions

1. As we are interested in accelerating nuclei, first we consider the injection of particles by the shock into the wind, because heavy particles are present there. This shock is non-relativistic and strong, with a velocity $v_{sh} \sim v_{wind}$. We assume that particles are accelerated by the Fermi I mechanism whose typical injection is $S(p) = S_0 p^{-\alpha} e^{-p/p_{max}}$, with $\alpha = 4$ for strong shocks, and the maximum momentum $p_{max}$ determined by equalling the acceleration and losses timescales. The only relevant radiative process for the particles are proton-proton inelastic collisions, and we consider also diffusion and convection escape losses. The normalisation is

$$ S_0 = \frac{L_{\text{nt}}}{\int_{p_{\text{min}}}^{p_{\text{max}}} 4 \pi p^2 p^{-\alpha} e^{-p/p_{\text{max}}} p c dp} \quad (7) $$

with $L_{nt}$ the luminosity injected into non-thermal particles. We adopt $L_{nt} = 0.1 \dot{M}_{\star} v_{wind}^2/2$, i.e., 10% of the kinetic energy of the stellar wind. The particle distribution $f(p)$ is obtained by solving the Fokker-Planck equation as in [31].

2. The particles distribution $f(p)$ produced at the first shock propagates to the next one. During propagation, we only consider that particles can lose energy because of the radiative processes or escape the system. To obtain the particle distribution reaching the next shock, we solve the Fokker-Planck equation with temporal and energy dependency and assume that the particles move along the jet convected by the outflow material, thus time and position inside the jet relates as $t = \int_{z_0}^{z} [v_{jet}(z')]^{-1} dz'$.

3. The amount of accelerated particles encountering the next bow shock is $f_{prop} \propto (R_{sp}(z)/R_{jet})^2$. These
cosmic rays will be reaccelerated at the relativistic shocks in the jet. We estimate the distribution after reacceleration in one shock using the following expression [32–34]

$$f_{\text{reacc}}(p) = \alpha_{\text{rel}} \int_{p_0}^{p} \frac{dp'}{p'} f_{\text{prop}}(p') \left(\frac{p}{p'}\right)^{-\alpha_{\text{rel}}} , \quad (8)$$

where $\alpha_{rel} = 3/(1 - \xi^{-1}) = 4.5$ when the compression factor is $\xi = 3$. This expression is formally valid for the case of reacceleration in non-relativistic shocks. We consider it as good approximation in our scenario since $\Gamma_0 = 5$ and the shocks are mildly relativistic [35]. We calculate $p_{max}$ at each shock and introduce a cut-off $e^{-p/p_{max}}$ to the distribution.

4. All these steps are repeated until there are not more bow shocks inside the jet.

# 4 Preliminary results

As a preliminary calculation, we consider the shocks produced by three WR stars arbitrarily located at $z_1 = 100$, $z_2 = 150$, and $z_3 = 200$ pc. At $z_1$ we consider only acceleration in the wind bow shock, because we are interested in accelerating heavy nuclei. At $z_2$ and $z_3$ we consider reacceleration at the shocks in the jet. For the preliminary results presented in this contribution, we do not consider acceleration of particles in the wind shocks at $z_1$ nor at $z_3$.

## 4.1 Maximum energies

Fig. 6 shows the timescales related to the parameters mentioned before for a wind shock at $z_1$ (top panel), and jet shocks at $z_2$ (middle panel) and $z_3$ (bottom panel), respectively. Adopting the Bohm diffusion regime and the escape distance $0.5R_{sp}(z)$, the maximum energy of the particles ($E_{max}$) is determined by the particle escape out of the acceleration region. In the case of the wind shock, it can be seen that $E_{max}$ is much smaller than in the case of the shocks in the jet, mostly because $B_{jet} > B_{wind}$. As mentioned in Section 2, $R_{sp}(z)$ increases faster with $z$ than the decay of the magnetic field of the jet, and therefore $E_{max}$ slowly increases with $z$, obtaining $E_{max}(150 \text{ pc}) = 5 \times 10^{16} \text{ eV}$ and $E_{max}(200 \text{ pc}) = 6 \times 10^{16} \text{ eV}$ for protons.

## 4.2 Particle distributions

We calculate the particle distributions following the procedure described above. In Fig. 7 we exhibit the proton distribution after each step. The distribution injected in one shock arrives to the next one after propagation almost without being modified. This is a consequence of the low jet matter density, which prevents the cooling through proton-proton collisions. After reacceleration, the distribution becomes slightly harder ($\alpha = -3.8$).

Figure 8 displays the injected distribution at $z_1$ and the result after the interaction with two successive shocks at

<figure><img src="image_6.png" /><figcaption>Figure 6. These plots show the acceleration, cooling by proton-proton inelastic collisions, and escape timescales for a shock into the wind of a star located at $z = 100 \text{ pc}$ (top panel), for shocks in the jet of a system at $z = 150 \text{ pc}$ (middle panel), and at $z = 200 \text{ pc}$ (bottom panel). The maximum energies obtained are $10^{14} \text{ eV}$, $5 \times 10^{16} \text{ eV}$, and $6 \times 10^{16} \text{ eV}$, respectively.</figcaption></figure>

$z_2$, and $z_3$. The distribution extends from shock to shock up to higher energies because of the higher $E_{max}$ at each shock. After the interaction with two shocks, the spectral index of the particle distribution becomes $\alpha = -3.6$.

# 5 Conclusions

We investigate the acceleration of particles at multiple shocks inside AGN jets. These shocks are formed by the interaction of the jet with the winds of WR stars. The jet is charged with protons and heavier nuclei which are removed from the stellar winds. In particular, WR stars con-
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Evolution of the particle distribution between two shocks. Protons are injected at the shock into the wind at $z = 100$ pc, propagate up to the next shock located at $z = 150$ pc, where they are reaccelerated.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Particle distribution injected by the wind shock at $z = 100$ pc, and the distribution after reacceleration in the jet shocks at 150 pc and 200 pc.</figcaption></figure>

tribute with oxygen, nitrogen, and carbon, alike the mass composition inferred from the observation of UHECRs.

Particles are accelerated and reaccelerated in the bow shocks by the Fermi I mechanism. We find that the spectrum of the cosmic rays is just slightly modified by propagation between the shocks because of the low jet matter and photon densities. Accounting for the jet deceleration after each interaction, we find that the Hillas energy at the shocks increases with $z$ allowing to reach higher energies. The energy dissipated in the collisions also reduces significantly the Lorentz factor of the jet after the interaction with tens of stars, causing the shocks into the jet to become mildly relativistic. Reacceleration by multiple shocks, on the other hand, causes the spectrum to become flatter than the typical injection slope through first order Fermi acceleration. Overall, we show the potential of this scenario to produce UHECR in AGN jets. Nevertheless, better estimations require to improve the characterisation of the stellar populations and their spatial distribution in the surroundings of AGN jets.

The authors thank the Czech Science Foundation under the grant GAČR 20-19854S. A.T.A thanks the Marie Skłodowska-Curie fellowship.

## References

[1] L.A. Anchordoqui, G.E. Romero, J.A. Combi, Phys. Rev. D **60**, 103001 (1999), astro-ph/9903145

[2] A. Aab, P. Abreu, M. Aglietta, I.F.M. Albuquerque, I. Allekotte, A. Almela, J. Alvarez Castillo, J. Alvarez-Muñiz, G.A. Anastasi, L. Anchordoqui et al., ApJL **853**, L29 (2018), 1801.06160

[3] L.A. Anchordoqui, Phys. Rev. D **97**, 063010 (2018), 1801.07170

[4] The Pierre Auger Collaboration, A. Aab, P. Abreu, M. Aglietta, I.F.M. Albuquerque, J.M. Albury, I. Allekotte, A. Almela, J. Alvarez Castillo, J. Alvarez-Muñiz et al., arXiv e-prints arXiv:1909.09073 (2019), 1909.09073

[5] The Pierre Auger Collaboration, A.A. Halim, P. Abreu, M. Aglietta, I. Allekotte, K. Almeida Cheminant, A. Almela, J. Alvarez-Muñiz, J. Ammerman Yebra, G.A. Anastasi et al., arXiv e-prints arXiv:2211.02857 (2022), 2211.02857

[6] G.E. Romero, A.L. Müller, M. Roth, A&A **616**, A57 (2018), 1801.06483

[7] J.H. Matthews, A.R. Bell, K.M. Blundell, A.T. Araudo, MNRAS **479**, L76 (2018), 1805.01902

[8] A.L. Müller, G.E. Romero, M. Roth, MNRAS **496**, 2474 (2020), 2006.12259

[9] E. Peretti, G. Morlino, P. Blasi, P. Cristofari, MNRAS **511**, 1336 (2022), 2104.10978

[10] M. Kachelriess, *Extragalactic cosmic rays*, in *37th International Cosmic Ray Conference* (2022), p. 18, 2201.04535

[11] A.R. Bell, J.H. Matthews, MNRAS **511**, 448 (2022), 2108.08879

[12] S.S. Kimura, K. Murase, B.T. Zhang, Phys. Rev. D **97**, 023026 (2018), 1705.05027

[13] B. Eichmann, J.P. Rachen, L. Merten, A. van Vliet, J. Becker Tjus, JCAP **2018**, 036 (2018), 1701.06792

[14] J.H. Matthews, A.R. Bell, K.M. Blundell, A.T. Araudo, MNRAS **482**, 4303 (2019), 1810.12350

[15] F.M. Rieger, P. Duffy, ApJL **907**, L2 (2021), 2101.06014

[16] L. Merten, M. Boughelilba, A. Reimer, P. Da Vela, S. Vorobiov, F. Tavecchio, G. Bonnoli, J.P. Lundquist, C. Righi, Astroparticle Physics **128**, 102564 (2021), 2102.01087

[17] R. Mbarek, D. Caprioli, ApJ **921**, 85 (2021), 2105.05262

[18] J.H. Matthews, A.M. Taylor, arXiv e-prints arXiv:2301.02682 (2023), 2301.02682

[19] S.S. Komissarov, MNRAS **269**, 394 (1994)

[20] W. Bednarek, R.J. Protheroe, MNRAS **287**, L9 (1997), astro-ph/9612073
[21] M.V. Barkov, F.A. Aharonian, V. Bosch-Ramon, ApJ **724**, 1517 (2010), 1005.5252

[22] A.T. Araudo, V. Bosch-Ramon, G.E. Romero, MN-RAS **436**, 3626 (2013), 1309.7114

[23] S. Wykes, M.J. Hardcastle, A.I. Karakas, J.S. Vink, MNRAS **447**, 1001 (2015), 1409.5785

[24] B. Snios, S. Wykes, P.E.J. Nulsen, R.P. Kraft, E.T. Meyer, M. Birkinshaw, D.M. Worrall, M.J. Hardcastle, E. Roediger, W.R. Forman et al., ApJ **871**, 248 (2019), 1901.00509

[25] S. Thoudam, J.P. Rachen, A. van Vliet, A. Achterberg, S. Buitink, H. Falcke, J.R. Hörandel, ApJ **595**, A33 (2016), 1605.03111

[26] S. Wykes, B.T. Snios, P.E.J. Nulsen, R.P. Kraft, M. Birkinshaw, M.J. Hardcastle, D.M. Worrall, I. McDonald, M. Rejkuba, T.W. Jones et al., MN-RAS **485**, 872 (2019), 1812.04587

[27] J.H. Krolik, *Active Galactic Nuclei* (Princeton University Press, Princeton, 1999), ISBN 9780691227474, https://doi.org/10.1515/9780691227474

[28] N. Torres-Albà, V. Bosch-Ramon, A&A **623**, A91 (2019), 1902.05008

[29] D. Eichler, V. Usov, ApJ **402**, 271 (1993)

[30] A.M. Hillas, ARA&A **22**, 425 (1984)

[31] B.T. Park, V. Petrosian, ApJS **103**, 255 (1996)

[32] R.D. Blandford, J.P. Ostriker, ApJ **237**, 793 (1980)

[33] P. Schneider, A&A **278**, 315 (1993)

[34] M. Cardillo, E. Amato, P. Blasi, A&A **595**, A58 (2016), 1604.02321

[35] A. Zech, M. Lemoine, A&A **654**, A96 (2021), 2108.12271