Particle acceleration in radio galaxies with flickering jets: GeV electrons to ultrahigh energy cosmic rays

Variability is a general property of accretion discs and their associated jets. We introduce a semi-analytic model for particle acceleration and radio jet/lobe evolution and explore the effect of Myr timescale jet variability on the particles accelerated by an AGN jet. Our work is motivated by the need for local powerful ultrahigh energy cosmic ray (UHECR) sources and evidence for variability in AGN and radio galaxies. Our main results are: i) UHECR and nonthermal radiative luminosities track the jet power but with a response set by the escape and cooling times, respectively; ii) jet variability produces structure in the electron, synchrotron and UHECR spectra that deviates from that produced for a constant jet power - in particular, spectral hardening features may be signatures of variability; iii) the cutoff in the integrated CR spectrum is stretched out due to the variation in jet power (and, consequently, maximum CR energy). The resulting spectrum is the convolution of the jet power distribution and the source term. We derive an approximate form for a log-normal distribution of powers; iv) we introduce the idea of $\sim 10$ GeV 'proxy electrons' that are cooling at the same rate that UHECRs of rigidity 10 EV are escaping from the source, and determine the corresponding photon frequencies that probe escaping UHECRs. Our results demonstrate the link between the history of an astrophysical particle accelerator and its particle contents, nonthermal emission and UHECR spectrum, with consequences for observations of radio galaxies and UHECR source models.


INTRODUCTION
Multi-wavelength variability in the observed fluxes of accreting sources is observed on a range of timescales. This flux variability has a number of near-universal characteristics. Often, the variability is consistent with "flicker" noise, also known as 1/ noise or pink noise because the power spectrum follows a 1/ slope that is steeper than white noise but shallower than red. A pedagogical discussion of flicker noise in astronomy is given by Press (1978). A linear 'rmsflux relation' has been observed in X-ray binaries (Uttley & McHardy 2001;Uttley et al. 2005;Heil et al. 2012), accreting white dwarfs (Scaringi et al. 2012; Van de Sande et al. 2015), blazars (Biteau & Giebels 2012) and other X-ray AGN (Uttley et al. 2005). A related phenomenon is a log-normal distribution of observed fluxes, seen in both disc-dominated (Gaskell 2004;Alston 2019) and jet-dominated sources (H. E. S. S. Collaboration et al. 2010Collaboration et al. , 2017Chevalier et al. 2019; see also Morris, Chakraborty, & Cotter 2019). In addition to flickering-type variability, sources also undergo episodic bursts of accretion activity, with prominent examples being the outburst cycles of X-ray binaries (Fender et al. 2004), restarting radio galaxies (Kaiser et al. 2000;Brienza et al. 2018;Konar et al. 2019) and, possibly, the optical 'changing-look' phenomenon in quasars (LaMassa et al. 2015;MacLeod et al. 2016;Runnoe et al. 2016). ★ matthews@ast.cam.ac.uk The accretion disc is intimately connected to the outflows produced by the system, which can take the form of winds or relativistic jets. If the jets are launched by the Blandford & Znajek (1977) mechanism, then for a black hole (BH) with gravitational radius , the power of the jet depends on the BH spin parameter ( * ) and magnetic flux (Φ ) threading the event horizon as BZ ∝ ( * Φ / ) 2 . The precise relationship between the accretion rate, , and jet power is complicated; the production of powerful jets might require special conditions such as a magnetically arrested disk (MAD) state (Narayan et al. 2003;Tchekhovskoy et al. 2011;Liska et al. 2018) and/or the presence of a disc wind to collimate the flow (e.g. Globus & Levinson 2016;Blandford et al. 2019). If we assume that Φ 2 is proportional to ( 2 ) and that * changes on relatively long timescales, then BZ ∝ 2 . This proportionality is also expected on general energetic grounds (e.g. Chatterjee et al. 2019;Davis & Tchekhovskoy 2020). It is therefore reasonable that the jet variability is in some sense a filtered version of the accretion-induced variability, a concept which has been successfully applied to jet modelling in X-ray binaries (Malzac 2014;Malzac et al. 2018). What effect does this jet variability have in AGN and radio galaxies? How do multiple episodes of accretion, which are themselves variable, affect the jet propagation, feedback and any observable radiative or UHECR signatures?
In AGN, there are also longer-term aspects of the accretion process that can imprint variability in the observed accretion flux and the power in the jet. For example, on long timescales ( Myr), the accretion rate is likely to be determined by the supply of cold gas to the central region of the galaxy. In the chaotic cold accretion model proposed by Gaspari (2016), the accretion rate on to the BH is predicted to follow a log-normal distribution with a pink/flicker noise power spectrum, just as is observed on shorter timescales in accretion discs. Simulations by Yang & Reynolds (2016) also show flickering type variability in the jet power caused by self-regulated accretion on to the central AGN, with jet power consequently varying over a large range, (> 2 dex; see also Beckmann et al. 2019). Long-term variability in jet power therefore seems inevitable from a fuelling perspective.
Astrophysical jets accelerate nonthermal particles (see Matthews et al. 2020 for a review), which radiate as they interact with magnetic fields or radiation. Our primary way of learning about radio galaxies, is, as the name suggests, through radio emission from synchrotronemitting electrons. As well as these nonthermal electrons, AGN jets can accelerate high-energy protons and ions, which we refer to as cosmic rays (CRs). The origin of the highest energy CRs, known as ultrahigh energy cosmic rays (UHECRs), is not known. The maximum energy attainable in a particle accelerator is given by the Hillas energy = , where = / is the characteristic velocity, is the magnetic field strength and is the size of the acceleration region. The Hillas energy is a general constraint that can be understood in terms of moving a particle of charge a distance through an optimally arranged − × electric field. The Hillas criterion states that any accelerator must be a factor −1 larger than the Larmor radius of the highest energy particles it accelerates, i.e.
> −1 ( ), where denotes the Larmor radius. Calculating with some characteristic numbers for radio galaxies reveals they are one of the few sources capable of reaching >EeV energies and accelerating UHECRs. For this reason, they have long been discussed as potential UHECR sources (e.g. Hillas 1984;Norman et al. 1995;Hardcastle et al. 2009;Eichmann et al. 2018;Matthews et al. 2018Matthews et al. , 2019a, along with other classes of AGN jets and alternative sources such as gamma-ray bursts (GRBs), starburst winds and cluster-scale shocks.
The Hillas energy can be used to obtain a power requirement for UHECR production, first derived using a dynamo model by Lovelace (1976), and discussed further by various authors (e.g. Waxman 1995;Blandford 2000;Nizamov & Pshirkov 2018;Eichmann et al. 2018). This power requirement is sometimes referred to as the Lovelace-Hillas or magnetic luminosity condition. If a particle reaches the Hillas limit, then the magnetic power must satisfy −1 10 43 erg s −1 for acceleration to a rigidity of / = 10EV. For reasonable departures from equipartition and accounting for a likely maximum energy some factor below , UHECR acceleration probably requires kinetic jet powers in the region j 10 44 erg s −1 . There are relatively few nearby radio galaxies (within the various UHECR horizons) capable of reaching these powers, if the current jet power estimates are taken at face value (Massaglia 2007b,a;Matthews et al. 2018Matthews et al. , 2019a. Any variability in the jet power will inevitably affect the morphology and dynamics of radio galaxies, as well as any feedback processes at work. There are now many examples of restarting or 'double-double' radio galaxies (e.g. Kaiser et al. 2000;Brienza et al. 2018;Konar et al. 2019), in which there appear to be distinct, discrete episodes of jet activity, but there is also subtler evidence for variability in the jet power. For example, Fornax A displays evidence for a varied and complex history (Iyomoto et al. 1998;Lanz et al. 2010;Maccagni et al. 2020), while Centaurus A shows a distinction between different lobe structures on scales of ∼ 2 kpc and ∼ 300 kpc (Morganti et al. 1999;Croston et al. 2009), which appear to be connected with different episodes of activity. The merger and gas fuelling history of the sources is complex in both Fornax A (e.g. Iyomoto et al. 1998;Mackie & Fabbiano 1998;Iodice et al. 2017) and Centaurus A (e.g. Stickel et al. 2004;Neff et al. 2015). In a more general sense, it is clear that there are many aspects of radio galaxies that are time integrated -for example, the total energy input into the surroundings, or the total energy stored in synchrotron-emitting electrons -and variability creates a disconnect between these integrated quantities and the instantaneous jet properties.
In this paper, our aims are threefold: (i) to introduce a numerical method capable of modelling, in a simple parameterised fashion, the morphology, radiation and UHECR signatures from radio galaxies with variable jet power; (ii) to study the effect of jet variability on observational signatures such as the synchrotron luminosity and broadband spectral energy distribution; (iii) to study the acceleration and escape of UHECRs in radio galaxies with variable jet powers. In our modelling, we focus on jets with a flicker noise power spectrum and a log-normal power distribution. We make a number of further simplifying assumptions and the model is unlikely to provide quantitative matches with real astrophysical sources. Indeed, our approach is mostly heuristic -we aim to demonstrate some key principles regarding particle acceleration in variable jets that can be used to study UHECR and electron acceleration. We begin by describing our method (section 2) and present the results from a single simulation in section 3. In section 4, we introduce the concept of 'proxy electrons', before discussing extragalactic CR propagation, observational applications and limitations of the model. We conclude in section 5.

A SIMPLE VARIABLE JET AND PARTICLE INJECTION MODEL
We begin by introducing our model for the evolution of a flickering jet and the nonthermal particles it accelerates. The jet has a kinetic power that behaves stochastically according to an input power spectrum and probability density function; we begin by describing the process of generating this jet power time series, before describing the dynamics of the jet and the treatement of particle acceleration.

Synthetic Jet Power Histories
To create a synthetic jet power time-series, we use an input power spectral density (PSD) and jet power probability density function (PDF) to generate a time series for the jet power for a given set of PSD and PDF parameters. We use the method described by Emmanoulopoulos et al. (2013) and implemented in python by Connolly (2015). The algorithm is similar to the widely used Timmer & Koenig (1995) method, except that it allows for more flexibility in specifying the underlying The PSD is specified as a power-law model of the form where is the temporal frequency. We set = 1 for a "pink" or flicker noise spectrum and use bins of 0.1 Myr when generating the synthetic time series. We use a log-normal distribution of jet power, j , as our input PDF, given by where ln 0 is the natural logarithm of the median jet power and is the standard deviation of ln 0 . The mean and mode of the distribution are given by exp(ln 0 + 2 /2) and exp(ln 0 − 2 ) respectively. These quantities are shown for comparison in Fig. 1.
As mentioned above, log-normal flux distributions are common in both disc-and jet-dominated accreting systems. In using a lognormal jet power distributions we make an implicit assumption that some form of multiplicative process imprints variability on the jet on long timescales. The mean of a log-normal distributions is larger than both the median and the mode. This asymmetry is important for UHECR production, since AGN jets must reach high powers ( 10 44 erg s −1 ) in order to accelerate CRs to energies beyond 10EeV. Furthermore, the jets that have higher values of (more variable) also have a greater difference between the mean and median values of .
We model the propagation of a jet with a constant mass density, j . The mass density is parameterised via a density contrast with the surrounding medium, , which is defined as the density relative to density of the ambient medium at a radius of = 0.1 kpc. The jet nozzle has a constant radius, j . We assume that the power of the jet is variable and driven by accretion rate fluctuations, and that the jet is kinetically dominated (negligible jet pressure). The kinetic power of a one-sided relativistic jet can be derived from the equations of relativistic hydrodynamics (e.g. Taub 1948;Landau & Lifshitz 1975;Wykes et al. 2019) and is given by where j is the jet mass density, Γ j is the bulk Lorentz factor, j = 2 j is the cross-sectional area of the jet nozzle and j is the jet velocity. We assume that the jet variability is entirely accounted for by variation of the jet velocity/Lorentz factor; thus for a given value of the two-sided jet power j = 2 j,1 we invert the above equation to solve for j (or equivalently, Γ j ). The jet is transrelativistic, so the velocity is allowed to transition between non-relativistic and relativistic regimes. The jet power ultimately determines the jet advance speed, the energy input into the lobes, the maximum CR energy and the normalisation of the particle source term (see subsequent sections). The mass input rate from the jet is j = j j Γ j j , which, together with the rate of change of volume, determines the density in the lobes.

Jet and lobe dynamics
A number of analytic and semi-analytic models for jet propagation have been developed (e.g. Begelman & Cioffi 1989;Falle 1991;Kaiser & Alexander 1997;Turner et al. 2018;Hardcastle 2018). Generally, these models concern themselves with relatively powerful jets, that are well confined and FRII-like in their morphology; our approach is similar, with the main differences being that we allow for a flickering jet power and model the expansion of the lobe by discretising into a series of cylindrical cells along the direction of propagation (the -axis). We consider the propagation of a light jet, with typical density contrasts ∼ 10 −4 . The speed of the jet varies depending on the power, with typical Lorentz factors ranging from non-relativistic to Γ j ≈ 10. We assume that the advance of the jet is determined by equating, in the rest frame of the working surface, the momentum flux of the jet with that of the ambient medium (Marti et al. 1997), with an additional geometric factor (see Appendix A). We refer to the bubble inflated by the jet as a 'lobe' (rather than, say, a cocoon) throughout this paper for simplicity. The sideways expansion of the lobe is set by the bow-shock jump conditions, based on the difference in pressure between the lobe and surrounding envi- Figure 1. The relationship between various properties of a log-normal jet power distribution, plotted on a logarithmic (base 10) scale. The blue solid and orange dashed lines show the mean and mode of the distribution, given by ln 0 + 2 /2 and ln 0 − 2 respectively. The mean increases as , the natural logarithm of the standard deviation, increases. The black dotted line is shown on a different (linear) -axis, and shows the fraction of time the jet has a power exceeding 10 0 . This is equivalent to the survival function of the distribution (or one minus the CDF) evaluated at 10 0 erg s −1 . The plot illustrates how a wider log-normal distribution (a more variable jet power) produces more favourable conditions for UHECR acceleration; as increases, the jet spends more time in a 'high' state and has a higher integrated energy output.

Figure 2.
Schematic illustrating the main aspects of our method. A light, collimated jet propagates into an isothermal ambient medium with decreasing density and pressure, and inflates a cocoon or lobe. In the process it accelerates radiating electrons and CR ions, which gradually cool, and, in the case of CR ions, escape from the lobe. The dynamic model is described in Section 2.2 and Appendix A. The evolution of nonthermal particles is described in Section 2.4, and the escape of CR ions is described in Section 2.5. . ronment. The lobe pressure is calculated in a self-consistent manner from the energy injection and volume of the lobe.
The jet propagates into an ambient medium with a density and pressure calculated from the 'universal pressure profile' described by Arnaud et al. (2010). This pressure profile is characterised by a single variable, 500 , which is the mass contained within the radius where the density is 500 times the critical density of the Universe. The ambient medium is assumed to be isothermal, with the temperature, , set by the − 500 relation from Arnaud et al. (2005), which gives ≈ 2.3 keV for 500 = 10 14 . The density then follows from = ( / ), where = 0.62 is the mean particle mass (we assume a fully ionized solar abundance plasma). This choice of parameters leads to a particle number density of ≈ 2 × 10 −3 cm −3 at 100 kpc. The approximate functional dependence is a smooth broken power-law of the form ( ) ∝ ( / ) − 1 [1 + ( / )] ( 1 − 2 ) , with 1 = 0.305, 2 = 5.70 and = 715 kpc for our adopted 500 = 10 14 . The domain is discretised in the direction of propagation, , with cells, and the evolution of the width of the lobe is calculated for each of the cells containing jet material. We compute a dynamic time step and a particle evolution time step, with the latter evolved during subcycles within the main dynamics loop. Our model for the jet dynamics is described further, with the relevant governing equations and more detail regarding the numerical scheme, in Appendix A. A schematic diagram describing the main features of the method is shown in Fig. 2.

Energy partitioning
To conduct our simulations, we have to decide what fraction of the jet's kinetic power is transferred to the various forms of energy. Following Hardcastle (2018), we assume that a fraction = 0.5 of the jet's power goes into doing work on the surroundings (see also Blundell et al. 1999;Bourne & Sĳacki 2020); the remaining half is stored as internal energy in the lobes and used to calculate the lobe pressure and sideways lobe expansion (see Appendix A). The use of this work factor is an approximation and while it is fairly accurate in estimating the overall work done it does not account for the fact that the work done varies over time and can in principle be higher than the energy input in periods of low jet power (see Appendix A for further details). The energetic particle populations comprise nonthermal electrons ( ) and nonthermal protons and ions -which we refer to as CRs ( ). We assume that a fraction of the jet power is transferred into nonthermal electrons, such that / = j . For simplicity, we assume the same amount of energy is transferred into CRs, giving / = j where = . The equations and source terms for the particle populations are given in sections 2.4.2 and 2.4.1 and account for adiabatic losses. For the magnetic field, we set / = j . Our energy partioning factors are defined at injection, so the fact that the different components cool at different rates means that these quantities differ to the instantaneous equipartition factors often used in the literature. In addition, our quantities are defined as a fraction of the total jet power rather than relative to other components. For the equation of state of lobe plasma, we use an adiabatic index of 4/3, although a four-fluid 'effective adiabatic index' could be adopted, in a manner similar to Pfrommer et al. (2017). The partitioning factors adopted in our simulations are given in Table 1. We experimented and chose values that gave reasonable radio luminosities and typical magnetic field strengths in the lobe of 10s of G, in line with the results of Croston et al. (2004).

Nonthermal particles
AGN jets transfer energy to nonthermal particles. We model this process by evolving populations of both electrons and various CR ion species, with a source term that depends on the jet power. We do not model the particle acceleration physics, and we remain agnostic about the details of the process. We assume that the particles are accelerated in e.g, shocks close to the jet head, over a short acceleration time, acc ∼ ( / ), where is the so-called 'gyrofactor (e.g. Aharonian 2000) such that ≈ 1 in the Bohm regime. We then allow the particles to cool in (and escape from) the lobe of the radio galaxy. Hereafter we use the subscript to denote electrons and to denote CR ion species.

Nonthermal electrons
We evolve the electrons in the lobes according to the continuity equation where ( ) = / is the differential spectrum of electrons at energy . The first term accounts for synchrotron, inverse Compton and adiabatic losses through the total cooling timescale ( ), and ( , ) is the electron source term, given by where is a normalisation constant and = ln( max, / 0, ) for = 2. The injection energy is set to 0, = 10 2 . The maximum electron energy max, is kept constant at 100 TeV. We have verified that this maximum energy is achievable for the range of input jet powers, under the assumption that the maximum particle energy is limited by synchrotron cooling in the jet hotspot. However, the exact value is not well constrained. We discuss the reasons for this and possible improvements further in Section 4.5. We solve equation 4 following the method described by Chang & Cooper (1970) and Chiaberge & Ghisellini (1999), in which the equation is discretised and written in tridiagonal matrix form and then solved using a tridiagonal matrix algorithm (TDMA; also known as a Thomas algorithm). The total radiative cooling rate due to inverse Compton and synchrotron processes for electrons with energy is Here = 2 /8 is the magnetic field energy density, rad is the radiation field energy density and crit is the energy density of the Schwinger field, crit = 4.41 × 10 13 G. We set rad = 4.17 × 10 −13 erg cm −3 , corresponding to the energy density of the cosmic microwave background (CMB) as reported by Fixsen (2009). The energy density of the CMB is equivalent to that of a 3.24 G magnetic field. We assume that the CMB is the dominant radiation field in our simulations, which is not true close to the host galaxy, but is likely to be a good approximation once the lobe has reached ∼ 10kpc in length. Electrons and CR ions both also lose energy adiabatically as the lobe expands at the rate / | ad = 1/3( / ) / . The adiabatic loss timescale is then ad = /( / | ad ), and the total cooling timescale is given by the reciprocal sum such that

Nonthermal ions (cosmic rays)
In the same manner as the electrons, we evolve a CR distribution for each CR ionic species, , with a continuity equation that differs slightly to the electrons, given by where ( ) = / is the differential spectrum of CR species at energy , cr esc is the CR escape time (see Section 2.5) and loss is the energy loss timescale due to photopion, photodisintegration and pair production losses. The total differential CR spectrum is then . We do not model the coupling between the different ion species (see below). We again use a power-law with an exponential cutoff for the source term, ( / , , ), given by (for where again is a normalisation constant and = ln( max, / 0 ) for = 2. The injection energy is assumed constant over time for each species and set to 0 = 10 2 (a CR Lorentz factor of 10). The maximum energy max at each time step is set according to the power requirement for UHECR production based on the Hillas energy, which is discussed in the introduction, and is given by where H = max,i / denotes the fraction of the Hillas energy attainable. The maximum CR energy therefore accounts for the partitioning of magnetic energy via the term . In addition, the use of this maximum energy condition ensures that the UHECR acceleration time is always short compared to j / . It does also implicitly assume that particle acceleration is taking place close to the Bohm regime, but this must be the case for almost any UHECR accelerator (Hillas 1984). Our method could feasibly be made more complicated by incorporating more detailed aspects of shock acceleration physics. However, our approach here captures the general behaviour of the maximum CR energy depending on time through its dependence on the square root of the jet power, which has an important impact on the results.
The quantity controls the relative fraction of ion in the CR population. An enhancement of heavy ions relative to intrinsic abundance is expected on theoretical grounds (e.g. Ellison et al. 1981;Caprioli et al. 2011;Marcowith et al. 2016;Matthews et al. 2020), and there is empirical evidence for a heavy composition in both Galactic CRs (e.g. Aglietta et al. 2004;Blasi & Amato 2012a) and UHE-CRs beyond the ankle (Pierre Auger Collaboration 2017; de Souza 2017). We take a phenomonological approach for charge and mass dependent injection that approximately fits the TeV-range Galactic CR spectrum, proposed by Wykes et al. (2017). In this approach, the spectrum in energy per nucleon is scaled by 2 / , where is the atomic mass number of the ion and is the solar abundance. This is equivalent to setting = 2 −2 . We note that alternative scalings for chemical enhancements have been suggested, for example by Caprioli et al. (2017). However, the enhancement factor they propose is rather mild, scaling only as ( / ) 2 . We include all ions that are at least as abundant by number as Fe in the Sun (H, He, C, N, O, Ne, Mg, Si, Fe).
For CR losses, we consider the normal loss mechanisms at ultrahigh energy ( 1 EeV), where protons undergo losses due to the GZK effect (Greisen 1966;Zatsepin & Kuz'min 1966), ions with > 1 undergo photodistintegration, and all species undergo pair production losses. For each of these processes, we tabulate rates from (Alves Batista et al. 2016) at = 0 for the CMB radiation field and the Gilmore et al. (2012) model for the extragalactic background light (EBL). We only include the CMB and EBL radiation fields and we do not allow for the cascade in decreasing A and Z for heavy nuclei as they photodistintegrate and lose nucleons; instead, we treat the process as an energy loss in equation 9 via the / loss term. We do not include proton synchrotron losses in our model, because, although they can be important for 10 19 eV protons in 1mG fields (e.g. Aharonian 2002), the field strengths we consider instead range from 10 − 250 G, such that the proton synchrotron timescale is longer than the simulation time even for the highest energy protons.

CR Escape Time
The escape time is difficult to estimate. Ultra-high energy protons and ions will gradually escape out of the lobe via a combination of drifts, streaming, advection and diffusion. We expect the particles to be accelerated in shocks relatively close to the jet head, although we remain agnostic about the details of the acceleration mechanism. Particles might be advected fairly quickly into the jet lobe by hydrodynamic backflow, but this backflow slows and dissipates its energy through shocks or is broken up by vorticity and shear; as such, the flow gradually becomes slower and subsonic or transonic (e.g. Falle 1991; Reynolds et al. 2002;Matthews et al. 2019b). The lobe is therefore broadly characterised by transonic turbulence with bulk flows of ∼ where is the sound speed. Flows can reach a significant fraction of (Reynolds et al. 2002;Matthews et al. 2019b), implying a minimum advective timescale adv ∼ 100 kpc/ ∼ 1 Myr. However, in practice, bulk flows are slower, and less uniform, at larger distances from the jet head, and even if the CRs are advected deep into the lobe it is not clear that this facilitates quicker CR escape. This argument depends on the details of the CR escape and the exact magnetic field structure in the lobe, particularly the degree of connectivity between the lobe field and surroundings, motivating further study. Generally, we expect the UHECRs to diffuse out of the lobes more quickly than they are advected away, and even in the case of being advected long distances they must still escape the turbulent cluster magnetic field. For our purposes, we assume that the CRs escape diffusively from the lobe.
We make a rough estimate of the escape time by assuming Bohm diffusion and using the magnetic field estimate described above. If Bohm diffusion applies, the escape time for ultra-relativistic particles of energy and charge from a sphere of radius is given by where = /3 is the Bohm diffusion coefficient and is the Larmor radius. Our lobes are not spherical, so to quickly estimate the characteristic escape distance we approximate them as ellipsoids. The centre of mass of the lobes is at the origin, since we are assuming symmetrical jet propagation. The average distance to the lobe edge, which we take as a characteristic distance of escape, can then be calculated as esc = 2 K ( )/ . Here is the lobe length and K is the complete elliptic integral of the second kind (which ranges from 1 to /2 in this case).
= 1 − 2 / 2 is the square of the eccentricity and is the lobe width, defined at the widest point, which is generally close to = 0. Our final UHECR escape time estimate is We do not account for the magnetic field of the ambient medium or time delays introduced by diffusive propagation to Earth, but we discuss this and other limitations of our approach in Section 4.2.

UHECR and Radiated Luminosities
One of our aims is to examine the effect of jet variability on the particle populations in the lobe, as well any observational consequencies for radio galaxies and UHECR origins. It is therefore useful to consider a few different (multimessenger) luminosites that are relevant for observers. We record the escaping UHECR luminosity and the radio luminosity due to synchrotron-emitting electrons, but also calculate the full broadband SED at some time steps The UHECR luminosity at a given time can be calculated by integrating over the particle distribution for each species, or equivalently, by counting the number of CRs with > that leave the system in the solution of equation 9. We record the CR luminosity above 8 EeV, in order to compare to various studies of CR anisotropies, spectrum and composition by the Pierre Auger Collaboration. We also record the full CR spectrum over time.
We calculate the synchrotron power per unit solid angle, ( ), where is the electron Lorentz factor, by assuming an isotropic pitch-angle distribution (e.g Crusius & Schlickeiser 1986; Ghisellini et al. 1988) and a single magnetic field strength throughout the lobe volume at a given time. The total synchrotron luminosity is then calculated by integrating over the electron population to find the synchrotron emissivity per unit solid angle, , and multiplying by the volume of the lobe times 4 , such that the specific luminosity is given by where 0 and max are the minimum and maximum electron Lorentz factors. We record at 144 MHz and 1.4 GHz for comparison with typical observation frequencies. We also present full broadband SEDs in Section 3.4, calculated using (Hahn 2015). We tested our synchrotron calculation against and the code from Hardcastle et al. (1998), finding excellent agreement.
Protons can produce gamma-rays via and collisions. As above, we assume that the CMB and EBL radiation fields are dominant in our simulations. This assumption means that gamma-ray radiation from interactions is negligible. The energy loss timescale to collisions is pp ≈ 600 (10 −4 cm −3 / ) Gyr (Sikora et al. 1987), where is the target proton density. The typical densities is our modelling are 10 −5 cm −3 , so we expect hadronic gammarays to be negligible and do not consider their contribution. We also checked, using , that the luminosity from collisions was much lower than the IC contribution. However, we note that hadronic gamma-rays may be important either near the base of the jet or when the lobe densities and/or lobe energy contents are higher. In addition, it may be that densities are higher outside the lobe, or inside the lobe where there is significant mixing with the external medium due to, e.g., instabilities at the contact discontinuity. Hadronic gamma-rays could therefore feasibly be detectable from CRs in these regions.

RESULTS FROM A SINGLE, ILLUSTRATIVE JET MODEL
We produce results from a single jet history. We choose to model a light ( = 10 −4 ) jet with a median (two-sided) jet power of 0 = 10 45 erg s −1 and a variability parameter of = 1.5. This median jet power is fairly typical of estimates for FRII radio galaxies, if  Fig. 4. The insets show the PSD used to generate the jet power time series and the resulting histogram of jet powers (with logarithmically spaced bins). On the right-hand y-axis, the corresponding value of 2 j max / is shown, for = 0.1 and H = 0.3 (the parameters used in the simulation), with dotted horizontal lines marking 0.5 dex intervals. The right-hand axis is shown to give the reader a feel for the range of maximum energies protons can attain from a jet with these parameters, in the rough range from 10 18 eV to 3×10 19 eV.
slightly towards the lower end (Godfrey & Shabala 2013;Ineson et al. 2017). We chose these values of 0 and to give reasonable overall energetics and dynamics for our source, but we also discuss the sensitivity to these parameters (as well as the PSD slope, ) in Section 4.5. The full set of parameters for this simulation, which we refer to as our reference model, is given in Table 1. We generate the synthetic jet power time series according to the procedure in Section 2.1 and evolve the jet until it is 300 kpc long. As the jet and lobe evolve, we record the various luminosities, lobe dimensions, lobe physical quantities and timescales, as described in the previous section. The jet power time series is shown in Fig. 3, with an inset showing the resulting histogram of jet powers.

Time evolution of the system
In Fig. 4, we show the evolution of the important timescales, dimensions and physical quantities ( , , , j Γ j ) from our reference simulated radio source. The UHECR escape time and adiabatic cooling time are defined in the previous section, the synchrotron cooling time is denoted sync , and is the age. We calculate sync for 1 GeV electrons and for electrons emitting with a characteristic frequency of 5 GHz, and give the CR escape time for 10 EeV protons. We also define the sound crossing time, cs = / where is the lobe length and ≈ 700 km s −1 is the sound speed in the (isothermal) ambient medium. Our reference model has a supersonic jet and its advance is also supersonic, so the longest timescale in the system is generally cs . We do not show the IC cooling time, which, for 1 GeV electrons, is approximately 1.2 Gyr. This timescale exceeds the longest synchrotron cooling time, as expected since > CMB is always satisfied in this case. The particle acceleration time acc is not shown since this is always less than a few hundred years.
At early times ( 5 Myr), the magnetic field is high (∼ 100 G) before gradually dropping to around 10 G at later times ( 60 Myr). The strength of the magnetic field is important because it sets the hierarchy of the UHECR escape time, cr esc , and electron synchrotron cooling time, sync . Stronger (weaker) magnetic fields  are better (worse) at confining UHECRs and cause quicker (slower) synchrotron cooling (see the next subsection). The UHECR escape time also depends on the size of the radio lobe. The strength of the magnetic field is proportional to the square root of the internal energy of the lobe; thus it is determined by a competition between the rate of change of volume and rate of change of energy. The former depends on the dynamics of the lobe, in particular the momentum flux, and the latter is equal to the jet power. Early on in the lobe's evolution, the environment is denser, so the relative change of volume for a given jet power is smaller. As a result, the magnetic field tends to decrease over time, because the lobe inflated by the jet finds it easier to expand. This behaviour is broadly consistent with the results of Croston et al.
(2018), who find higher mid-lobe pressures (for which we expect correspondingly higher B fields) in shorter lobes. We find that the length of the lobe increases faster than the width, meaning that the lobe aspect ratio -or axial ratio, both defined as length divided by width -increases over time, consistent with Blundell et al. (1999).
The jet lobe responds in a number of ways to the variability of the jet. Each time the jet power increases, the advance speed also increases, energy is injected more quickly into the lobe and the particle populations are re-energised. The exact impact on a given quantity is complicated, because it depends on the relative importance of a given period of high activity compared to the overall history of the jet up to that point. Generally, each spike in jet activity causes a cor- Figure 5. The evolution of the various luminosities over time in the reference simulation. The jet power is shown in black, matching Fig. 3. The UHECR luminosity is calculated above 8EeV for comparison with studies from the Pierre Auger Observatory, and the radio luminosity is given in units for observation frequencies of 144MHz and 1400MHz. responding increase in synchrotron and UHECR luminosity, which then decays away with a characteristic response time set by the synchrotron cooling and UHECR escape time, respectively. This can be seen in Fig. 5. The effect of the cooling and escape times on the light curve is therefore to act as a 'low-pass filter'; high-frequency variability is smoothed out and the power spectrum steepens above a temporal frequency = 2 / where is the relevant timescale. This phenomenon has been discussed in the context of blazars by, e.g., Finke & Becker (2014) and Chen et al. (2016). Fig. 6 shows a number of different representations of CR spectra, in units of 2 ( ) normalised to the maximum within the energy range 10 16 −10 21 eV. The CRs that have not yet escaped from the lobe follow the ( ) ∝ −2 injection spectrum up to a characteristic break energy at which the escape time equals an effective source age. In contrast, the escaping CRs follow ( ) ∝ −1 , which peaks around the same characteristic break energy. The low-energy slope for the escaping CRs, which is dominated by the protons, is a result of the assumed rigidity-dependence of the diffusion coefficient, which we have set to the Bohm diffusion coefficient such that esc ∝ / . Both the internal and escaping spectra have a cutoff at higher energies. The rightmost panel of the Fig. 6 shows CR spectra from 30 random times throughout the simulation. There is significant diversity in the CR shape, particularly at high energies where the CRs can escape easily and are particularly sensitive to the recent activity of the jet. Below the spectral break, both the escaping and internal spectra behave in the same way as the time-averaged spectrum.

CR spectra and composition
The CR are split into different species, and a clear trend of heavier composition with increasing energy can be seen in both the internal and escaping CR spectrum. To illustrate this further, we show the value of ln , the mean of the natural logarithm of the atomic mass number contributing to a given energy bin, in Fig. 7. Both CR acceleration and CR escape are rigidity-dependent rather than energy-dependent; as a result the value of ln naturally increases with energy. In the case of the CRs internal to the lobe, higher energies will exceed the CR energy at which light species have already escaped the lobe. In the case of the escaping CR, higher energies will exceed the maximum CR energy for light species. There is therefore a subtle difference between the increasing value of ln for the escaping and internal CRs. The fact that the escaping CRs have a lighter composition can be understood in two ways. Firstly, that at a given , CRs with lower have a higher rigidity and so find it easier to escape. Alternatively, one can think of the energy difference between the curves in Fig. 7 as encoding the difference between the characteristic maximum energy, and the characteristic energy of escaping CRs. The difference between these is maximised when max / → ∞ and minimised when the source lifetime is comparable to the escape time for CRs of rigidity max / .
The abundances of species outside the source environment itself is dictated by the transport properties in this external region. In the case where energy dependent transport dominates (as may be expected if particles propagate diffusively), the resulting composition will change. Specifically, there will be an ehancement of heavy species outside the source; this effect is discussed further in Section 4.2.

The CR cutoff
In our modelling we apply an exponential cutoff with a maximum energy that is proportional to √︃ j / j (equation 11), as derived from the Hillas condition. As the jet power varies, so does this maximum energy (see Fig. 3); in more powerful episodes, both the maximum energy and normalisation associated with the source term increases, as can be seen in the range of cutoff energies in the spectra at random times shown in Fig. 6. To gain more insight into the problem, we consider the total integrated spectrum (both internal and escaped CRs) and neglect the / loss term. In this limit, the total spectrum is the integral of the source term, i.e. ( ) = ∫ Δ 0 ( / , , ) , where Δ is the outburst time. In the limit of large Δ , and setting j = 1 for simplicity, it can be shown that (see Appendix B) where = H √ is the constant of proportionality in the maximum energy equation. We have written this using a general PDF ( j ) so that it can be applied to other situations where the PDF of input powers takes a different form. The specific form of equation 16 for our adopted log-normal ( j ) is This integral is not analytically straightforward, but a very good approximation can be obtained by considering only the peak of the integrand (see Appendix B). The general behaviour can already be appreciated from equation 16: a spread in jet powers acts to stretch out the exponential cutoff in the spectrum, such that the cutoff function is the convolution of the jet power PDF and the source term cutoff. The effect is similar in nature to the convolution of a CR source term with a luminosity function, and a comparable discussion in the context of supernova remnants is given by Shibata et al. (2010). The fact that the source term normalisation is proportional to j biases the spectrum towards higher energies, even for a symmetric PDF. This effect is even more pronounced if ( j ) is positively skewed as is the case for a log-normal distribution. The result is that the most important episode of activity for determining the CR maximum energy cutoff is likely to be the most energetic episode that has occured within a UHECR loss-time. Figure 6. Cosmic-ray spectra from the reference simulation, in 2 ( ) units where ( ) = / is the differential spectrum; an −2 CR spectrum appears as a horizontal line. Left The CR spectra inside the lobe averaged over the jet history. The thick black line shows the total CR spectrum, and the coloured lines show the contributions of individual species. The spectral shape matches the injected spectrum up to a characteristic break energy, beyond which the spectrum steepens due to the escape of high-energy CRs. Centre As in the left panel, but for the CRs escaping the lobe. The CR spectra behaves in an opposing manner to the left panel, with an inverted spectrum below a characteristic energy, and a flat ∼ −2 spectrum above this energy. The spectrum steepens above a break energy determined by the range of maximum energies during the simulation due to the variable jet power; the shape of this cutoff is discussed further, for an idealised case, in Section 3.3. Right CR spectra from 50 random times during the simulation. Broadly speaking, the spectra are similar to the averaged spectrum, except that jet variability causes the location of the maximum energy cutoff for the escaping CRs to jump around over-time. In addition, the value of cr esc increases over time as the lobe becomes larger, which changes the break energy for the CRs inside the lobe and the peak in the escaping spectrum. Figure 7. The value of ln as a function of energy in the simulation. Solid and dashed lines show the mean value across the jet history, and the shaded region shows the standard deviation at each value of . The CR escape time is a rigidity ( / ) dependent quantity, so at a given energy lighter species find it easier to escape due to their lower . This means that the escaping CRs (orange) have a lighter composition than the CRs remaining inside the lobes (blue).

Spectral Energy Distribution
In Fig. 8, we show broadband spectral energy distributions (SEDs) at 10 Myr intervals throughout the jet history. The SEDs are calculated using as described in section 2.6. Each curve is colourcoded according to the elapsed time, and the inset axis shows the corresponding electron spectra. The top panel shows the jet power over time, as in Fig. 3, but with the time intervals marked with vertical dashed lines to aid interpretation of the plot. We also show characteristic observing frequencies for radio and X-ray, as well as for the Fermi Large Area Telescope (Fermi LAT) with coloured horizontal bands.
The SED is characterised by a classic double-humped shape, with the low energy bump caused by synchrotron emission and the highenergy bump caused by inverse Compton scattering off the CMB. We also computed the spectrum from collisions but found it was insignificant due to the low density of target protons in the lobes. Although a double-humped SED shape is often associated with blazars (e.g Fossati et al. 1998;Ghisellini et al. 1998), it is a generic feature of a population of electrons interacting with magnetic fields and radiation fields when the energy densities of the two fields are comparable. For IC scattering in the Thomson regime, the relative contribution of the synchrotron and IC processes is given by the ratio of the relevant energy densities ( and CMB in this case). In the reference model, > CMB at all times, and the decrease of over time causes the synchrotron bump to be more dominant at early times compared to late times.
There are clear kinks and inflection points in both the SED and electron spectrum (inset). A correspondence can be seen between features in the synchrotron hump and both the IC hump and the electron spectrum, as expected, although some of the smaller features in the electron spectrum are smoothed out in the SED since each electron energy range produces radiation over a broader range of frequencies.
The kink and inflection features are caused by the variability of the jet; effectively, the result is a 'many-populations' model where the electron spectrum is a series of injected populations that have each been subjected to a different cooling history. The behaviour of the electron spectrum can be understood further as follows. If a given short period of high activity dominates over the integrated historical energy input of the jet, the system behaves like a system in "outburst", and the spectrum resembles that of a single power-law spectrum with a cooling-break. However, in some circumstances, the period of high activity can only really be seen at high frequencies or particle energies; this is because the low energy particles haven't had time to cool and so still reflect the integrated activity, whereas the high energy particles (above the cooling break) from previous activity have all cooled, and the new powerful episode injects fresh particles in this regime that can radiate. As a result, a characteristic spectral shape is seen, with inflection points and a secondary peak and break associated with the recent "outburst", leading to a spectral hardening in the spectrum. A similar result is found by Turner (2018), who refers to the spectral hardening as a 'steep-shallow' spectrum produced by a source with short distinct outbursts. This spectral hardening ef-fect could be an important signature of variability, so we discuss the observational perspective in Section 4.3.

Which photon frequencies track the UHECR escape time?
From Fig. 5, we can see that the radio and UHECR luminosities respond to impulses from the jet with a decay time set by the synchrotron cooling and CR escape times, respectively. We can therefore calculate at which electron energy the electrons are cooling at the same rate that the UHECRs are escaping; we refer to the electrons at this energy as 'proxy electrons' for the UHECRs. These proxy electrons emit their synchrotron and inverse Compton radiation at characeristic frequencies in radio and higher-energy bands. These radiation frequencies can be determined by comparing the UHECR escape time to the relevant cooling times. The characteristic frequency of synchrotron emission from an electron with Lorentz factor can be conveniently written in terms of the Schwinger field such that Similarly, we can write the synchrotron cooling time as We can combine these two expressions and eliminate the dependence so that Finally, equating this with equation 13 and rearranging for frequency gives proxy ≈ 20 GHz This is the characteristic frequency of synchrotron radiation emitted by the proxy electrons, which cool at the same rate that UHECRs of rigidity / escape from the lobe. Unfortunately, this equation is extremely sensitive to many of the parameters, particularly the magnetic field. Decreasing the magnetic field by a factor of 10 changes proxy by a factor of 10 5 , which shifts the observing window from the radio to the far infrared (∼ 15 m). The corresponding IC emission from such "proxy electrons" is, We have marked both proxy and IC /ℎ with vertical dot-dashed lines on Fig. 8. For these specific parameters, the IC emission from the proxy electrons falls in the ∼MeV range, at lower energies than probed by current -ray telescopes such as Fermi LAT. The similarity of the features in the SED at proxy and IC /ℎ is apparent -whenever a kink or feature occurs in the synchrotron curve at proxy a matching one is seen in the IC bump, confirming that the two frequencies correspond to electrons of the same energy.
The results of this section are summarised in Fig. 9, where we show Figure 9. A comparison of the synchrotron cooling time, sync , and CR escape time, cr esc , as a function of . We show sync at two different frequencies, compared to estimates of cr esc for two different escaping diffusion coefficients. cr esc is calculated for esc = 100 kpc and / = 10EV. The Bohm diffusion estimate intersects with the 20GHz synchrotron curve at ≈ 10 . We also show, in a grey band, the interquartile range (IQR) of observational estimates of from Croston et al. (2005, C05, their  a comparison of synchrotron cooling times and UHECR escape times as a function of for some appropriate parameter values. To give a feel for typical magnetic fields in radio galaxy lobes, we also show the interquartile range (IQR) of observational estimates of from table 9 of Croston et al. (2005). Our reference model has higher magnetic fields than these observational estimates at early times but at late times the ∼ 10 G fields fall within this range. Fig. 9 can be thought of as a graphical representation of the simple derivation above; longer CR escape times (from larger lobes, higher or slower transport) are tracked by proxy electrons with longer synchrotron cooling times, so the value of proxy decreases as cr esc increases. The above derivation is carried out for a given CR energy, but an alternative approach is to assume that the UHECRs have reached the Hillas energy, = . However, it is important to draw a distinction here between the conditions in the accelerator and the conditions when the CRs are escaping. We might expect the acceleration site to be close to the hotspot, for example in secondary shocks near the jet head (e.g. Matthews et al. 2019b), where the characteristic velocities and magnetic field strengths are higher than the average values in the lobe. Bell et al. (2019) showed that acceleration by multiple shocks in radio galaxy backflows provides a mechanism to reach close to the Hillas energy ( max ≈ 0.6 ). One advantage of these backflows is that they are not highly relativistic, thus the difficulties associated with UHECR acceleration at relativistic shocks (Lemoine & Pelletier 2010;Reville & Bell 2014;Bell et al. 2018) can be avoided. Let us denote the magnetic field strength, velocity and size scale of the acceleration site as acc , acc and acc , respectively. We can then use the Hillas energy to write the escape time in terms of the ratios of these quantities to the lobe values, that is where we have used some characteristic values assuming the UHE-CRs are accelerated in non-relativistic ( acc ∼ 0.1) backflow shocks that are smaller than esc by a factor of 50 but also have stronger magnetic fields than in the lobes by a factor of 10, based on results from Matthews et al. (2019b).

Extragalactic CR transport
We have studied how CR diffusion out of the lobes affects the escaping composition, spectrum and luminosity over time, but we have not accounted for the propagation time of the CRs through extragalactic and Galactic magnetic fields, or the resulting effect on CR composition. Deflection in the extragalactic magnetic field (EGMF) introduces a time delay relative to rectilinear propagation or light travel. This delay is in addition to the escape time from the source, such that the total delay relative to the radiative signatures is tot = cr esc + prop − ( / ) for a source at distance , where prop is the propagation time (outside the lobe) and we assume esc . The propagation time depends on the diffusion coefficient as prop ∼ 2 / . For the rest of this discussion, we assume that the scattering length is small compared to the source distance, > / . The rigidity dependence of the diffusion coefficient depends on the ratio of the Larmor radius to the coherence length of the magnetic field, (see Guedes Lang et al. 2020, their eq. 1). For low rigidity particles, < , the propagation is diffusive with ∝ where we have chosen values such that the < regime applies. This propagation time here can be extremely long, but since in a nGstrength EGMF the Larmor radius of a 1 EV rigidity CR is 1.08Mpc, the above regime is only realised for very large coherence lengths or particle rigidities below 1 EV. For rigidities above 1 EV we generally expect > in the EGMF, and we instead obtain prop ∼ 0.28 Gyr 1 Mpc where the field is assumed to be random and turbulent. For smaller , such that < / , an expression of similar form to equation 26 can be obtained in the small-angle scattering limit. The above estimates vary depending on the exact numerical prefactors adopted. Time delay effects have been dicussed in the context of variable UHECR sources with reference to transient sources such as gammaray bursts; Miralda-Escude & Waxman (1996) show that the observed UHECR spectrum of a transient or 'bursting' source tends to peak at the UHECR energy where prop is equal to the time since photon arrival. The basic principle is similar to the idea of 'proxy electrons' we introduced in Section 4.1, except that in the case of a flickering source the particles are accelerated over many episodes rather than in a transient burst or impulse. A further time delay can be introduced by the escape from the environment surrounding the source -for example, the cosmic rays might have to escape from a cluster. Cosmic-ray escape from cool-core clusters has been studied by Kotera et al. (2009), with typical magnetic field strengths of 10 G in the cluster leading to escape times comparable to those considered here. The estimates here are uncertain, but they show that the propagation times are of larger or comparable magnitude to the escape time from the lobes for CRs with >EV rigidities, propagating in ∼ 1 nG fields for sources within ∼ 10s of Mpc. The propagation time is significantly shorter if the source distance becomes comparable to the scattering length. Any additional delay modifies the equations in Section 4.1 and implies that lower energy proxy electrons may be more appropriate as tracers of the propagating UHECRs. At lower CR rigidities and for larger source distances, the propagation times quickly become restrictive with respect to typical source lifetimes.
Diffusive propagation in the extragalactic region between the source and the Milky Way will lead to the CR density to accumulate due to the increased CR residence time in this region. Assuming the propagation is purely diffusive, once the steady-state level is reached, the CR density from a single source of distance will sit at a factor of ( / ) larger than the CR flux level from the source. Assuming that the rigidity dependence of the diffusion takes the form ∝ ( / ) , with > 0, rigidity dependent build up will occur. Subsequently, the steady-state density for CRs of the same energy will have heavier species abundance enhanced relative to light species, when compared to the escaping species abundance ratio from the source.

Observational Applications
Observations of radio galaxies commonly reveal complex radio morphologies and evidence for restarting or episodic activity. In addition to these morphological signatures, our work demonstrates the potential of spectral signatures as probes of variability. In particular, inflection points where the spectrum is a sum of two or more components with different cooling breaks can be a signature of variability and a recent 'outburst' (section 3.4). Related behaviour is seen in Fornax A; Maccagni et al. (2020) find that the central unresolved component has a spectral break at a higher frequency than for the outer lobes, suggesting a recent injection of fresh nonthermal electrons consistent with a jet power varying on Myr timescales. Spectral hardening in gamma-rays could also be a potential signature of variability, as has been observed at GeV photon energies in Centaurus A (Abdalla et al. 2020). Variability in radio galaxies could also be important when determining the ages of the sources. There is often a discrepancy between the measured spectral age of a source (as measured from the synchrotron break frequency) and the dynamical age, known as the 'spectral age problem' (Eilek 1996;Harwood et al. 2013Harwood et al. , 2015. This problem can be alleviated if the magnetic field is below equiparition (Mahatma et al. 2020). Jet power variability can also help solve the spectral age problem in a few different ways. Firstly, a recent injection of fresh nonthermal electrons could lead to the spectral age being measured as much smaller than the true age of the source. Secondly, a series of powerful episodes can cause the source to grow faster (advance more quickly) than if the source had a steady jet with a power equal to the median or mode of the jet power history, resulting in a measured dynamical age much longer than both the true age and typical synchrotron age of the nonthermal electrons.
Based on our modelling, we can also identify a number of interesting areas of the synchrotron and inverse Compton spectra that may provide further clues about particle acceleration and jet variability. In the gamma-rays, we suggest that characterising, in detail, the steepening of the spectrum associated with the cooling of electrons could be used to search for spectral hardening indicating recent jet activity. Observations close to the frequencies associated with proxy electrons could also be useful for predicting UHECR luminosities from individual objects although, as we have shown, the appropriate frequency for the proxy electrons is strongly dependent on the CR transport.
Additionally, we note that our work has some implications for AGN more generally, in the context of accretion modes and radio properties of quasars. Long-term accretion variability is likely to be important for explaining why quasars with nearly identical emission line properties can have very different radio properties (Richards et al. 2011, Rankine et al., submitted). Such a characteristic can be partly explained by a disconnect between the properties of the quasar accretion disc and the nonthermal electrons, as expected if the synchrotron cooling time is longer than the timescale over which the quasar changes its ultraviolet/optical properties. A similar principle might also be applied to the high-/low-excitation radio galaxy (HERG/LERG) dichotomy, which classifies radio galaxies based on the presence or otherwise of strong optical emission lines (Laing et al. 1994;Tadhunter et al. 1998;Best & Heckman 2012). The HERG/LERG distinction is thought to be fundamentally driven by Eddington fraction and accretion mode (Best & Heckman 2012), but the two classes are found across a wide range of luminosities. The optical emission lines would be expected to respond on shorter timescales than the radio, so it is plausible that some HERGs and LERGs are similar objects, but with the disc caught in different states; for example, HERGs could correspond to objects that have recently (within a synchrotron cooling time) transitioned into a radiatively efficient state. A refinement to our model would be needed to investigate this in more detail; we create a synthetic jet power time series, but one could determine the jet power from an underlying accretion model with different modes of fuelling. Future work in this area might focus on incorporating particle acceleration physics into numerical simulations that simultaneously model the fuelling of the AGN as well as the resulting jet (e.g. Yang & Reynolds 2016;Beckmann et al. 2019).

Parallels with Galactic accelerators
There are important parallels to draw between radio galaxy models for UHECRs and the phenomenology of Galactic CRs. Galactic CRs up to the knee are thought to be accelerated in the shocks of supernova remnants (SNRs), although reaching to PeV energies is a challenge and requires (at least) significant CR-driven magnetic field amplification (Blasi 2013;Bell et al. 2013;Bell 2014). Jets in X-ray binaries might also contribute, particularly at the high-energy end (Heinz & Sunyaev 2002;Fender et al. 2005;Cooper et al. 2020), and the discovery of TeV gamma-rays from SS433 indeed suggests that X-ray binary jets might be good high-energy CR accelerators (Abeysekara et al. 2018). Within the SNR paradigm for Galactic CRs, the most energetic galactic CRs are thought to come from young SNRs with fast-moving shocks, with the bulk of the lower energy CRs originating from older, more numerous SNRs and forming a near isotropic 'soup'. As a result, it is natural that stochasticity becomes important at CR energies close to the knee (e.g. Blasi & Amato 2012a,b), because the number of sources capable of reaching these energies is relatively small. Can similar behaviour be invoked at ultrahigh energy? Nearby radio galaxies are compelling sources for explaining the UHECR spectrum and anisotropies at extreme ( 10 EeV) energies Eichmann et al. 2018;Eichmann 2019b;Guedes Lang et al. 2020). At lower energies (0.1 to a few EeV) the energy loss length-scales from the GZK effect and photodistintegration are larger. Thus, the CRs in this energy range could instead come from an ensemble of background radio galaxies as discussed by Eichmann (2019a) and Matthews et al. (2019b), forming a roughly isotropic component with a spectrum largely determined by the integrated radio galaxy luminosity function. This isotropic component would still have to be produced in relatively local sources, since the propagation times impose a magnetic horizon at low rigidities, as can be seen from inspection of equation 25. The overall importance of stochasticity in determining the UHECR spectrum is emphasized by the fact that the very same nearby radio galaxies that are often invoked as UHECR sources also show evidence for variable jet powers, a complex merger history and unusual radio lobe morphologies.

Model limitations, parameter sensitivity, and future work
Our modelling approach is relatively simple and heuristic, and there are a number of physical characteristics of lobes and jets that are not well modelled by our scheme. We have already discussed cosmic ray transport. Another limitation is the treatment of the lobe as a single, uniform bubble, when in reality there is a complex interplay between the jet head, the jet itself and the lobe (e.g. Falle 1991), and the hotspot is significantly higher pressure than the rest of the lobe. Our assumption affects the observed integrated spectrum, but also the cooling history of the electrons. A future model might include two or more 'zones', with particles gradually moved from the "hotspot" population to the lobe at an appropriate advective rate. A technique like this has been applied in the RAiSE semi-analytic models . Alternatively, magnetohydrodynamic simulations using tracer particles to model ensembles of electron could be used to account for the mixing of populations and non-uniformity of the magnetic field; Vaidya et al. (2018) describe a recent implementation in the hydrodynamics code (Mignone et al. 2007). We have chosen to present results from one jet model, and the specific jet power time series chosen will clearly affect the results. Even the random number seed changes the exact jet power history, but it is more relevant to focus on the fundamental parameters of the time series: the variability parameter, , the PSD slope, , and the median jet power, 0 . Increasing the median jet power increases the typical maximum CR energy attainable in the accelerator (from equation 11). As a result, the spectra shift to the right in Fig. 6. The effect on the magnetic field is more nuanced. The magnetic energy being injected into the lobes increases with 0 , but the energy density of the lobes also depends on the rate of change of lobe volume, and is expected to increase as 0 raised to an exponent < 1. We therefore expect the magnetic field to be stronger in the lobes of more powerful sources, leading to slower CR escape and faster synchrotron cooling. Adopting a higher makes the jet power PDF more skewed, such that the mean jet power is higher compared to 0 . However, increasing also increases the amplitude of variability, which has the effect of exacerbating variability signatures in the time evolution (Figs 4 and 5) and broadband SED (Fig. 8). Changing the shape of the PDF would also have a profound effect, with more symmetric distributions (e.g. a Gaussian) qualitatively mimicking the behaviour of a low , and more asymmetric or skewed distributions increasing the difference between the mean and median j . Finally, we consider the PSD slope, . The pink/flicker noise PSD we used ( = 1) puts equal power in each logarithmic frequency range. A steeper PSD slope, > 1, leads to more power in low-frequency modes. For = 2, the time series behaves as red noise and follows a random walk. As a result, once the jet enters a high-state, it tends to stay there for longer. This behaviour could in principle lead to more dramatic signatures in the UHECR luminosity or calculated spectra. The spectrum would be more sensitive to the exact location in the time series, with longertimescale trends apparent in the luminosities and composition, and the opposite true for a flatter PSD slope.
There are numerous other parameters that are important, such as the energy partitioning factors ( , , ) and the particle spectral index, . The particle spectral index is particularly important for the UHECR luminosity. Our chosen value of = 2, the canonical value for shock acceleration (Bell 1978), spreads out the total energy equally in each decade of energy, but there are many reasons to ex-pect steepening of the CR spectrum (e.g. Bell et al. 2019). Even a mild steepening to = 2.1 causes (for a single specie) the fraction of total CR energy contained in UHECRs above 8EeV to drop to ≈ 2.5% (compared to ≈ 10% for = 2). Finally, we note that our assumptions about the maximum particle energy could be relaxed in future work. In particular, we assumed a constant maximum electron energy, max, , but it would be interesting to investigate the impact of a power-dependent max, , obtained by balancing the electron acceleration time with the synchrotron time in the amplified magnetic field in the jet hotspot. The value of max, primarily affects the shape of the SED in Fig. 8. It may also be necessary to account for more detail plasma physics for both electrons and CR ions, since, particularly in relativistic shocks, the growth rate of CR-driven turbulence at the shock can be severely limiting in terms of the maximum particle energy (e.g. Reville & Bell 2014;Araudo et al. 2016Araudo et al. , 2018Bell et al. 2018).

CONCLUSIONS
We have presented a study of particle acceleration in radio galaxies with jet powers that vary over time according to a flicker noise power spectrum, and examined the effect of this flickering on the UHECR and electron populations accelerated by the source. Our main conclusions are as follows: • A log-normal distribution of jet powers results in a mean jet power that is significantly higher than the median, and especially, the mode (Fig. 1). If such a distribution is realised on long timescales ( Myr) in nature, it provides a situation where UHECRs can be accelerated in the most powerful episodes, even when the most commonly observed jet power is significantly lower. This is true more generally for positively skewed probability distribution functions for the jet power and could also apply to systems that undergo distinct outbursts.
• We present a semi-analytic model for studying the evolution of radio jets and lobes and the particle populations they accelerate and contain (section 2). Our model is influenced by other semi-analytic approaches (e.g. Turner et al. 2018;Hardcastle 2018), but explicitly accounts for a variable power, transrelativistic jet and also models the non-radiating CR ion populations and their diffusion out of the lobe. The lobe is vertically discretised into a series of discs and the advance of the jet head, expansion of each disc and resulting lobe properties are calculated self-consistently in an iterative manner, also feeding into the particle solver. This method involves a series of assumptions, and various limitations of the method together with possible improvements are discussed in Section 4. We present one single jet model with a flickering jet power and use it to illustrate some general principles relating to the problem of particle acceleration in a variable AGN jet source.
• We find that the synchrotron and UHECR luminosities track the jet power, but with a characteristic response that is determined by the radiative cooling time and UHECR escape time, respectively (Fig. 5). The cooling and escape times act like a low-pass filter on the jet power time series, smoothing out fast variability. The response time changes throughout the jet history, as, typically, the lobe size increases and the magnetic field decreases over time.
• The variable jet power creates a series of electron populations with different normalisations, which are in turn each subjected to a different cooling history. Clear kinks and inflection points are observed in the electron spectrum and resultant SED (Fig. 8). Jet variability therefore produces complexity and curvature in the electron and radiated spectrum that clearly differs from the expected spectrum for constant jet power, or for a single electron population.
• The CR spectrum inside the lobe follows the spectrum at acceleration, but with a break energy determined by the CR escape time (Fig. 6). This break is smoothed out due to the varying magnetic field and lobe size. For our assumed energy-dependence of the diffusion coefficient, the escaping spectrum follows an − +1 spectrum up to a cutoff at a maximum CR energy dictated by the jet power.
• The total integrated CR spectrum (the sum of the escaping and internal spectrum) displays a cutoff that is similar to a stretched exponential. In the absence of CR losses, the shape of this exponential is equal to the expectation value of the cutoff function, such that the differential spectrum obeys (see Section 3.3 and Appendix B) In situations where the maximum particle energy is determined by the source power, it may be possible to search for evidence of jet variability by studying the cutoff shape. In addition, the cutoff is likely in general to be smoother and more gradual than a pure exponential, and the energy of the cutoff is more reflective of the most powerful outbursts than the mean or median maximum energy.
• We introduce the idea of 'proxy electrons' -nonthermal electrons that are cooling at the same rate that UHECRs are escaping from the source (section 4.1). We use this concept to derive an optimum observing frequency, the frequency of radiation that best tracks the UHECR luminosity. For our assumptions about CR escape, this frequency is typically 20 GHz for a 10 G field, with corresponding inverse Compton emission at ∼MeV energies. These values suggest radio luminosity may be reasonable UHECR luminosity proxy, but this is strongly dependent on physical parameters in the individual sources as well as the uncertain CR transport physics.
• Some local radio galaxies (such as Cen A and Fornax A) display characteristics of variable jet powers in a morphological sense. These same radio galaxies make for compelling UHECR sources. We discussed parallels with supernova remnants and propose a scenario somewhat analogous to Galactic CR production, in which the highest energy CRs are accelerated in recent outbursts in local sources, whereas UHECRs around the ankle would come from a near-isotropic background of radio galaxies as discussed by, e.g., Eichmann (2019a).
Our work highlights the influence of variability in the jet power on particle acceleration in radio galaxies and their resulting spectra and morphologies. It also emphasizes the importance of stochasticity for UHECR source models, particularly for local sources and at the highest energies.

APPENDIX A: DYNAMIC MODEL FOR THE JET AND LOBE
We consider a light jet propagating in the direction into an isothermal ambient medium with a decreasing density and pressure. Our jet model involves the discretisation of the -domain into a series of disc-shaped cells with height Δ = 0.01 kpc. The length of the jet after time steps is then = Δ where the dynamic time-step Δ = Δ / ℎ (note that there is a separate time step for the particle solver, see Section A1 of this appendix). The advance speed ℎ of the jet head is chiefly governed by ram pressure balance, such that ℎ depends on the jet speed and density constrast with the surrounding medium and is given by (e.g. Marti et al. 1997) where * follows the notation of Marti et al. (1997) and represents the relativistic generalisation of the density contrast. For nonrelativistic internal energies in the jet beam and ambient medium, as assumed here, the specific enthalpy can be neglected and * = Γ 2 j j / ( ), where ( ) is the ambient medium density at distance . Note that * tends to the density ratio j / as Γ j → 1. = 1/4 is a geometric factor, which we chose to give reasonable advance speeds of ≈ 0.01 . The physical motivation for the factor is that the jet termination shock has a smaller area than the bow shock, so the pressure acts over a wider area and slows the advance. Each disc-shaped cell has zero width until the jet material enters it, at which point the width of cell is evolved according to +1, = , + ⊥ Δ . This sideways expansion is governed by the expansion speed of the bow-shock. If the shocked gas between the bow-shock and the contact discontinuity is in pressure equilibrium with the lobe, then, from momentum conservation, we can write where is the angle between the local bow shock normal and the direction of jet propagation and is the position of the edge of the cell. This equation states that the transverse expansion of the lobe is driven by the pressure difference between the lobe and the ambient medium. ⊥ is a function of position since and both decrease away from the origin. The ambient pressure is modelled according to the pressure profile described in Section 2.2. We obtain by calculating the internal energy in the lobes and converting to pressure where = 4/3 is the adiabatic index and the volume is obtained by integrating over each cylindrical cell of width , that is This set of equations is evolved numerically as an initial value problem and thus requires an initial condition for the starting lobe geometry -for this purpose, we assume that the base of the lobe initially has a width ( 1 )/2. The outline of the lobe at 5 Myr intervals is shown in Fig. A1.

A1 Particle Solver and Time-stepping
The particle solver operates inside the main dynamic loop and undergoes a series of subcycles per dynamic time step. The time step is calculated as where = 0.4 is a Courant-type number and is an index denoting the energy bin, and we only include bins with > 10 −15 ,max so as to avoid small time-steps from fast-cooling bins with negligible numbers of particles. The TDMA solver and synchrotron emissivity calculation form part of a Python package called . We have tested the TDMA code using single injections of powerlaw and delta function electron distributions and comparing to analytic formulae. We have also tested the synchrotron code against the code from Hardcastle et al. (1998) and the code (Hahn 2015). The tests result in extremely close agreement and are documented within the code repository 1 .

APPENDIX B: THE SHAPE OF THE CUTOFF IN THE COSMIC RAY SPECTRUM
The variation in jet power causes a corresponding variation in the maximum CR energy. The maximum CR energy is given by equation 11, and can be written more simply as max,i = √︃ j / j where = H √ . The variation in maximum energy means that the source term has a cutoff function that varies over time, with a corresponding effect on the overall spectrum. To illustrate this, we consider the continuity equation 9 governing the evolution of CR ions. We focus on one CR specie (protons) and ignore both the cooling and escape loss terms, such that ( )/ = ( / , , ). The source term ( / , , ) is given by equation 10, and, if we neglect the (weak) variation in the ln( max / 0 ) term, we can write it as ( / , , ) ∝ j ( ) − − / max,i .
We can integrate ( )/ so that, after a time interval Δ the spectrum is (B2) 1 The code is publicly available at https://github.com/ jhmatthews/msynchro