Time-dependent branching processes: a model of oscillating neuronal avalanches

Recently, neuronal avalanches have been observed to display oscillations, a phenomenon regarded as the co-existence of a scale-free behaviour (the avalanches close to criticality) and scale-dependent dynamics (the oscillations). Ordinary continuous-time branching processes with constant extinction and branching rates are commonly used as models of neuronal activity, yet they lack any such time-dependence. In the present work, we extend a basic branching process by allowing the extinction rate to oscillate in time as a new model to describe cortical dynamics. By means of a perturbative field theory, we derive relevant observables in closed form. We support our findings by quantitative comparison to numerics and qualitative comparison to available experimental results.

www.nature.com/scientificreports/ of a particle corresponds to the change from detection to no detection 8,14,28 . Branching processes display a phase transition from asymptotic extinction with probability 1 to asymptotic survival with a positive probability, depending on the average number of created particles. At the critical point, the avalanche size of a branching process is power-law distributed 16 .
Recently, there has been evidence that neuronal activity, when compared to branching processes, is not at criticality but in a reverberating regime close to criticality 14,29,30 . The reason for the strong interest in the system's distance to a critical point is the criticality hypothesis, which states that information processing in the brain might be optimized by the cortical network being close to a critical point 10,12,13,31,32 . However, fitting power laws is notoriously difficult 30,33 and other means of verifying the criticality of the neuronal avalanches are essential. For this reason, the avalanche shape, defined as the average temporal profile of the avalanches, has received more attention in recent years. It is debated whether, at criticality, the neuronal avalanche shape takes the universal form of an inverse parabola 16,26,[34][35][36][37][38] , which is the case of a critical branching process. The universality of this shape has been particularly challenged in Ref. 38 by observations of oscillations that modulate the avalanche shape. In Ref. 38 , the oscillations are identified as γ-oscillations, which are a particular frequency band of brain waves between 30 and 100 Hz. Brain waves are electric oscillations spanning the entire brain that can be recorded using electroencephalography (EEG) 39,40 . The oscillations are organized into bands covering frequencies between 0.05 Hz (slow 4 band) and 500 Hz (ultra fast band). Their power spectrum follows approximately a 1/f distribution 41 . Statistically, they are linked to the quiet time distribution between avalanches and correlations between avalanche sizes 42 . Their physiological role and the biological mechanisms that sustain them continue to attract attention in neuroscience [43][44][45][46][47][48] The observation of oscillations in the EEG activity on the one hand, and the observation of scale-free avalanches of LFP activity on the other hand raises the question of how these two seemingly incompatible descriptions of the same phenomena can be reconciled. The experimental data in Ref. 38 indicates that γ-oscillations modulate the avalanche shape. Can branching processes incorporate oscillations as well? Will those oscillations modulate the avalanche shape, widely regarded as a universal feature?
To answer these questions, we extend in this paper the field theory in Ref. 16 by incorporating oscillatory extinction rates, see "Model". We then calculate observables such as "Moments" of the particles number, its "Covariance", "Survival probability", and the "Avalanche shape" (or temporal profile) and compare them qualitatively to experimental results from Ref. 38 . We conclude in "Discussion and conclusion".

Model
A branching process in continuous time t can be regarded as a reaction of a single particle type B that splits into K ∈ N 0 copies of itself, with Poissonian rate s [18][19][20] . The particle number at time t, which is also called population size, is denoted by N(t).
In the event of a reaction, the population size increases by K − 1 particles. The waiting time for an individual particle to undergo any such reaction is exponentially distributed with rate s. The number of offspring K may be a random variable itself, with probability distribution P(K = k) = p k . We call the reaction in (1) a branching event if K > 1 and an extinction event if K = 0 . Throughout the present work we initialize the process with one particle at time t 0 = 0 , so that N(0) = 1 . The totality of a realisation of a branching process, from the initialization until the termination after which the population size remains 0 indefinitely is referred to as an avalanche. Times between avalanches are not modelled 49 and there are no correlations between avalanches 42 .
In Ref. 16 , we introduced a Doi-Peliti field theory for the continuous-time branching process with timeindependent but arbitrary offspring distribution p k . We refer to such a branching process as standard branching process. The action functional of this field theory is with and time-dependent Doi-shifted annihilation and creation fields φ(t) and φ(t) respectively. The parameter r is called the mass in the context of field theories and is, according to Eq. (3), closely linked to the first moment of the offspring distribution E[K] . Traditionally 18 , the branching process is called sub-critical if E[K] < 1 (and thus r > 0 ), it is called critical if E[K] = 1 (and thus r = 0 ), and super-critical if E[K] > 1 (and thus r < 0).
The time-scale of this branching process is set by s, the rate with which any single-particle event takes place. The rate with which any particle spontaneously disappears, the extinction rate, is thus ǫ = sp 0 .
In the following we will focus on binary branching processes, i.e. p k = 0 for all k except k = 0 and k = 2 , so that p 2 and p 0 = 1 − p 2 in Eq. (3) are determined by r/s = 1 − 2p 2 . The rest of parameters follow then immediately, such as E[K] = 2p 2 , q 2 /s = p 2 = 1/2(1 − r/s) and q j = 0 for j ≥ 3 , as well as bounds such as s ≥ r + q 2 ≥ 0 . In particular r = ǫ − q 2 , so that the branching process is critical if q 2 balances the extinction rate. Extensions to branching processes with other time-independent offspring distributions are straight forward 16 .
Our extension to this model (2) consists of a time-dependent, oscillating extinction rate www.nature.com/scientificreports/ where ν is the frequency of the oscillations. The dimensionless amplitude A is the small parameter of the field theoretic perturbation theory about A = 0 . Our model takes the oscillations as given and for simplicity idealizes them as sinus functions 50 . It does not model their emergence from underlying mechanisms. In the following, r ≥ 0 is considered. The magnitude |A| is bounded by p 0 because ǫ(t) , as an extinction rate, has to be non-negative at all times. If the amplitude A is positive, then the extinction rate is initially suppressed, favouring more branching events. In the field theory, this extension amounts to adding the term (details in "Extension of the action"), to the original action A 0 , Eq. (2), where r = sp 0 − sp 2 is still the "bare mass". Our extension can be summarised as a standard branching process where binary branching takes place with rate q 2 = sp 2 and extinction with rate ǫ(t) , Eq. (4). We retain the parameterisation p 0 + p 2 = 1.
It will turn out that the perturbation remains noticeable at all times in the process in all observables considered. In particular, even at criticality, the "Avalanche shape" carries a clear signature of the oscillations despite its expected universality 34 . The model and the analysis in the present work therefore provide an explanation for the shape of the temporal profile of neuronal avalanches recently reported by Miller et al. 38 . Figure 1 shows two example trajectories of the population size N(t), together with the perturbation of the extinction rate, ǫ(t)/s − p 0 = −A sin(νt) . In all figures, the latter is shown in green with the ordinate on the right. All data in this work is presented in dimensionless form, in particular time as st.
In the following, we state the central results, whose field-theoretic derivation is relegated to the Appendix. In particular, we consider the "Avalanche shape" at criticality, which is a common observable in LFP recordings of the brain 5-8 and the "Covariance", which recently gained more attention in neuroscience as a tool to estimate the system's distance to the critical point 14,29,30 . Moments first moment. Since the extinction rate varies periodically, we expect that the first moment varies accordingly. As shown in "Appendix: First moment", the expected particle number is, where �(t) is the Heaviside function reflecting the initialization at t 0 = 0 , henceforth dropped from all expressions. This result is consistent with the known result of an inhomogenous Poisson process with a time-dependent event rate 51 . Figure 2 shows an estimate of the first moment based on Monte-Carlo simulations together with the analytical result Eq. (6) for r = 0 . In addition to this exact result, the dashed line in Fig. 2 shows the first order approximation E[N(t)] = g , Eq. (10), which is discussed below in "n-th factorial moment to first order".
Because the extinction rate is reduced during the first half-period, Eq. (4), branching dominates the process and the population size increases at first. During the second half-period, extinction dominates and the population size decreases, perfectly balancing, on average, the creation of particles in the first half of the period. According to Eq. (6) and Fig. 2, the expected population size E[N(t)] at r = 0 is always equal or larger than unity, which is . Both trajectories were initialized with one particle at t 0 = 0 . The perturbation of the extinction rate, −A sin(νt) with ν/s ∈ {π/4, π/2} , A = 0.5 is shown in green (right ordinate). When the extinction rate is lower, branching is effectively promoted. www.nature.com/scientificreports/ the expected number of particles of a standard branching process at criticality. Thus, while the extinction rate is not shifted on average, the expected particle number is on average larger than in the process without oscillations. Just like in a standard branching process, E[N(t)] converges to 0 for r > 0 and diverges for r < 0 , indicating that the present process is critical at r = 0 . However, unlike a standard branching process, the expected population size never converges for r = 0 as the oscillations never cease. It will turn out that the effective mass acquires no shift due to the perturbation. Second moment. The second moment can be calculated in closed form from a convolution integral involving the first moment only, "Appendix: Second moment". From Eqs. (32b), (33) and (35), At r = 0 the second moment is, to first order in A, which diverges asymptotically linear in t. Because E[N(t)] is bounded, the variance also diverges linearly in t. Figure 2 shows the ratio E[N 2 (t)]/(1 + 2q 2 t) , which illustrates the deviation of the second moment at A > 0 from that at A = 0 . For large t, the second moment shows a linear increase with an amplitude that oscillates mildly with period ν , so that E[N 2 (t)] oscillates around the A = 0 behaviour. For large q 2 t the ratio E[N 2 (t)]/(1 + 2q 2 t) is to leading order 1 − (As/ν)(2 cos(νt) − 1).
Although the extinction is not shifted on average, the second moment is shifted on average due to the oscillations, see n-th moment. In principle, all moments of N(t) can be calculated exactly, following the same lines as in the previous section "Second moment" and in "Appendix: Second moment". However, even the second moment involves an integral that cannot be carried out in closed form and, for higher moments, the procedure quickly becomes unpleasantly complicated. To calculate higher moments of the population size E[N n (t)] , we resort to a perturbative expansion in powers of the amplitude A, which can be systematically expressed in terms of diagrams shown in the appendix.
Following 16 , all moments can be written in terms of factorial moments, which are naturally produced by the diagrammatics of the field theory. We denote the n-th factorial moment by g n (t 0 , t) , so that where n ℓ are the Stirling numbers of the second kind. In the following, the factorial moments are expressed in orders of A,  www.nature.com/scientificreports/ The n-th factorial moment at A = 0 , given by was calculated in closed form from the diagrams as a matter of combinatorics 16 . The function g (0) n is dominated by (q 2 (t − t 0 )) n−1 for small r(t − t 0 ) ≪ 1 , while it is exponentially decaying for large r(t − t 0 ) ≫ 1 . In the present work, the factorial moments acquire a dependence on the initial time t 0 , when N(t 0 ) = 1 . Only to zeroth order in A, at A = 0 , do the factorial moments become time-homogeneous and reduce to those calculated in Ref. 16 . The next order term in the small-A expansion equals whose derivation is explained in "Appendix: n-th factorial moment to first order".
In the subcritical regime r > 0 , in large r(t − t 0 ) ≫ 1 all moments vanish exponentially, because g (0) n vanishes exponentially, see Eq. (11). For small r(t − t 0 ) ≪ 1 and large q 2 (t − t 0 ) ≫ 1 , the moments are dominated by the largest factorial moment, This can be seen by expanding E[N n (t)] in terms of the factorial moments which are asymptotically dominated by Within the small-A expansion of the n-th factorial moment, the terms g The first two moments at r = 0 can be approximated to first order in A by on the basis of Eqs. (12) and (15). The expressions are consistent with the exact expressions Eqs. (6) and (7) [as given in Eq. (8)], respectively. Figure 2 shows a comparison between the exact expressions (solid blue lines) and the approximation Eq. (16) (black dashed lines). Deviations are clearly noticeable only for large amplitudes A.

further observables
In the following, we analyse observables that are somewhat more involved to derive in the present framework.
In particular, the avalanche shape and covariance are more of immediate interest to experimentalists because these are more accessible from LFP recordings of the brain 5-8,14,29,30 . covariance. The autocorrelation function is the covariance of N(t 1 ) and N(t 2 ) and it is a common way to quantify how strongly correlated data are at different points in time. If the population size at different times was independent, the autocorrelation function at times t 1 and t 2 , t 1 = t 2 , would be zero (the converse does not hold in general 52 ). As shown in "Appendix: Covariance", the covariance can be calculated in closed form up to an integral, where t min = min{t 1 , t 2 } and t max = max{t 1 , t 2 }.  Fig. 3a,b. For A = 0 , the maximum of Cov(N(t 1 ), N(t 2 )) occurs at t 1 = t 2 = t and equals (2q 2 /r + 1)(exp −rt − exp −2rt) 16 . Furthermore, its maximum rises with increasing t, provided that rt < log(2) . When A = 0 , these results are only approximations and the exact position of the maximum depends on the phase and amplitude of the oscillations. In Fig. 3a, all st 2 are chosen to fulfil that condition, while in Fig. 3b, the larger st 2 don't meet this condition and show a decreasing maximum.
As an autocorrelation function, the covariance quantifies how the activity in the system at one instance influences the system at later instances. This has recently proven to be a valuable tool for determining the distance of neural networks from the critical point 14,29,30 . What the imposed oscillations imply for this tool is discussed in "Discussion and conclusion". Survival probability. The "First moment" shows that the expected number of particles is on average larger if A > 0 , i.e. if the extinction rate drops before growing in every period, compared to the case without oscillations or reversed order of rise and fall of extinction, A < 0 . As noted in "First moment", the sign of the mass r still determines whether E[N(t)] ultimately vanishes or diverges, even when E[N(t)] oscillates indefinitely for r = 0 . This observation raises the question whether the survival probability P s (t 0 , t) , that is the probability of N(t) > 0 at a given time t after initialization at t 0 , displays a corresponding behaviour.
Based on the derivation in "Appendix: Probability of survival to first order", to leading order in A we find The first term, which is independent of A, is the probability of survival of the critical branching process with constant extinction rate 16 (Fig. 4). The second term is the first-order correction and indicates a shift of the probability of survival. For positive A (leading to an initial decrease of the extinction rate) the survival probability increases compared to the system without oscillations. For negative A (corresponding to an initial increase of the extinction rate) it decreases. For A > 0 , the initial push into the supercritical phase seems to dominate the entire survival probability, even many oscillations later, see Fig. 4. Despite the shifted survival probabilities, Fig. 4, and shifted average number of particles, Fig. 2, the avalanches do not survive indefinitely. They die eventually with probability 1 at r = 0 . As the survival probability at criticality is shifted, Eq. (19), it can be expected that the avalanche duration distribution P T is affected by the oscillations. Considering that the avalanche duration is equal to the time of death T of an avalanche that started at t = 0 , it can be derived from the survival probability P T (t) = − dP s (t) dt 16 : It shows that the distribution of durations still follows the ∼ t −2 power law of conventional critical branching with constant extinction rate 16 , Fig. 5 (left panel). However, this power law has oscillations superimposed, which means that, strictly, scale invariance at criticality is broken. Yet, in a data-analysis based on binning the avalanche durations, these oscillations will be averaged out and thus not be visible. Similarly, simulations of the avalanche size distribution P S show that for small oscillation amplitudes, the power law distribution with exponent 3/2 appears to be maintained and the effect of oscillations is averaged out, Fig. 5 (right panel).
To appreciate better the effect of the extinction oscillations on the ultimate survival lim t→∞ P s (t 0 , t) , we also consider the asymptotics of large t for r > 0 or equivalently p 0 > p 2 , lim r→0 P T (t) = q 2 (1 + q 2 t) 2  www.nature.com/scientificreports/ to leading order. For r < 0 , or equivalently p 0 < p 2 , the limit is positive, where we have to rely on Eq. (56) being the analytic continuation for the result obtained at positive mass r. Eq. (21) implies that P s (t 0 , t) , with or without oscillations, vanishes in the limit of large times as extinction prevails, since r > 0 . Eq. (22) indicates that the linear increase in −r of the ultimate survival probability is present with or without oscillations, however, that the amplitude of that increase depends on A. As far as the frequency ν is concerned, the effect of the oscillations on the ultimate survival is most pronounced for ν = ±r, with the minimum attained if ν = r and the maximum for ν = −r . It is noteworthy that Eq. (23) no longer vanishes as r → 0 . The constraints mentioned above, such as As ≤ r + q 2 , do not affect this result, as As = q 2 /2 still produces lim r→0 − lim t→∞ P s (t 0 , t)| r<0 = 1/4 . Together with Eq. (19) this seems to suggest the possibility of a sudden onset survival, whereby lim t→∞ P s (t 0 , t) jumps from 0 at r → 0 to a finite value. However, it is crucial in which order limits are taken. Eq. (19) remains valid if r → 0 (tying ν = ±r ) before t → ∞ . When taken in this order, the ultimate survival probability is zero at the critical point. We therefore conclude that the ultimate survival increases linearly and continuously from 0 for r ≥ 0 , Eqs. (19) and (21), to a finite value at r < 0 , Eq. (22). If the critical point of the present process is defined as the  After suitable normalisation 16 the time-homogeneous version of the process displays a universal, parabolic shape. Not least because of its universality, the shape has gained some popularity to serve as a fingerprint of a process 34,54,55 . However, as shown in Figs. 6 and 7, the periodic extinction imposes characteristic humps on the shape. They are of course rooted in the periodic extinction, as an analytical calculation shows, see "Appendix: Avalanche shape to first order" for details. These oscillations remain visible in the shape V(t, T) at all T and all t, but become particularly vivid whenever the termination time T is commensurate with the period of the oscillations, 2π/ν . As we have used our field-theoretic scheme only to first order, the slight mismatch with simulation results at larger A, such as those shown in Fig. 7, is not surprising. However, at such large amplitudes, the resulting shape resembles that of recent experimental results, where γ-oscillation modulated the average shape of neuronal avalanches 38 . Furthermore, for larger amplitudes A, the avalanche shapes are asymmetric, which can be seen particularly well in Fig. 7 for sT = 16.
Comparison to experiments. The mean avalanche shapes, such as in Figs. 6 and 7 can be qualitatively compared to avalanche profiles recorded in the brain. We reproduce in Fig. 8 plots from Miller et al. 38 (Figure 4(a) in Ref. 38   www.nature.com/scientificreports/ tion 4.0 International License. The data was collected through multielectrode arrays implanted in three adult nonhuman primates (see Ref. 38 for details). The plots show qualitatively a good agreement between the data and our model. The general shape and periodicity is captured extremely well. Small qualitative disagreement occurs in the third and fourth plot in the relative heights of consecutive maxima and minima. This disagreement may be down to higher order correction terms in the amplitude A.
The figure illustrates that branching processes, which are commonly used to explain the statistical properties of avalanches recorded in the brain, can be extended to also incorporate neuronal oscillations. Although the observation of non-universal avalanche shapes questioned the criticality hypothesis of the brain 10,32,38 , our analytical results clearly show that criticality is compatible with avalanche profiles that are modulated by oscillations.
Universal parabolic shape. The universal parabolic shape may be recovered by suitably averaging across different termination times T. Devising such a scheme in a given numerical or experimental setting may not always be feasible 56 . Figure 9 shows the T-averaged and T-rescaled expected avalanche shape, Varying small amplitudes do not seem to alter this shape and appear to converge to the universal parabola shape when criticality is approached r → 0.

Discussion and conclusion
Our discussion is two-fold: first, we focus on aspects of interest to research in stochastic processes. We then follow with a discussion on the implications for research in the area of neuronal avalanches.   Fig. 7 above. On y-axis, s is the mean profile in Ref. 38 , corresponding to our V(t, T); on x-axis, duration T from Ref. 38 is denoted by time t in our article. Data was collected through multielectrode arrays implanted in adult nonhuman primates (see Ref. 38 for details). Grey lines, single electrode data, black line mean of array, red mean size-per-timestep. www.nature.com/scientificreports/ In this paper, we extend the standard branching process by including time-dependent, deterministic variations of a reaction rate. While we focus on varying the extinction rate only, our approach is applicable to any reaction rate in a Doi-Peliti field theory.
Neuronal avalanches can show oscillating behaviour, which we capture with a branching process model with oscillating extinction rate. These oscillations can be observed in any moment, of which we present the zeroth (the survival probability), first and second, and in any correlation function, of which we present the two-time covariance. All of these observables are calculated exactly. Furthermore, we introduce an approximation scheme for factorial moments to first order in the amplitude of the extinction rate oscillations. This allows approximating more complicated observables such as the survival probability and the avalanche shape in the limit of small oscillation amplitudes.
All analytically calculated observables are compared with simulation results from Monte Carlo simulations. While the exact analytic results match perfectly, we also evaluate the first-order approximation scheme and find good agreement in the limit of small oscillation amplitudes.
Although the extinction rate is unchanged on average, its oscillation leads on average to a shift in the moments. For example, if the oscillations start with a decrease of the extinction rate, the expected particle number is always greater compared to the process without oscillations. Conversely, if the oscillations start by increasing the extinction rate, the expected number of particles is always lower than without oscillations. Despite the average shift of particle numbers, the onset of indefinite survival is not shifted by the oscillations. We therefore conclude that the critical point remains unchanged.
Both subcritical branching processes and neuronal activity show exponentially decaying auto-correlation functions whose decay rate is independent of the spatial sub-sampling of the neural network 14 . In particular, if the neural system is close to the critical point, in the reverberating regime, the exponential decay is slow and can be observed over more than 100ms even in a single neuron's activity 14 . As the oscillations in the neuronal activity can have periods well below 100 ms (e.g. in the spindle band), they should be visible in the data.
Similarly, both the oscillating branching process presented here and recordings of neuronal avalanches show oscillating avalanche shapes 38 . In our model and in the data, the oscillations modulate the shape and are most pronounced when the avalanche duration is an integer multiple of the period of the oscillations. Surprisingly, these results defy the assumption that the avalanche shape is universal at criticality. However, when the shapes are also rescaled and averaged over the avalanche duration, our simulations indicate that the universal parabola shape is recovered.
In future research, the oscillating branching process should be compared quantitatively to local field potential recordings of oscillating neuronal avalanches. Thus, the presented model can contribute to the understanding of neuronal avalanches by rejecting or supporting the branching process picture of signal propagation in the brain. Furthermore, it would be interesting to find out whether other types of variation of the extinctions rate, such as random fluctuations, result in similar behaviour.

Appendix: extension of the action
The extension of the model consists of making the extinction rate time-dependent. Hence, in the following derivation, we focus on introducing the new extinction rate and ignore the branching process.
Let P(n, t) be the probability that there are n particles in the system at time t, i.e. N(t) = n . Then, the master equation for the extinction process only is The first step in the derivation of the field theory is the introduction of bra-ket vectors �n| and |n� , respectively, and ladder operators a and a † . The ket vector |n� represents that there are n particles in the system, while the bra vector �n| projects out the state of n particles, �n|m� = δ n,m , the Kronecker function. The ladder operator a † creates an additional particle, i.e.a † |n� = |n + 1� , and the operator a annihilates a particle (with the additional prefactor n) a|n� = n|n − 1� . Thus, the entire, probabilistic state of the system at time t can be represented by a ket vector |M(t)�, which is nothing else than a representation of the probability generating function.
In the next step, following 57 , the master Eq. (25) is used to derive a corresponding equation for the state vector |M(t)�, In the final step, based on work by Peliti 58 , equations such as (27) can be solved by a path integral with an action A ǫ , where in every normal-ordered compound operator, a is replaced by a field φ(t) and a † is replaced by a field φ † (t) , and where an additional time-derivative term is added: P(n, t) = ǫ(t) (n + 1)P(n + 1, t) − nP(n, t) .
Scientific RepoRtS | (2020) 10:13678 | https://doi.org/10.1038/s41598-020-69705-5 www.nature.com/scientificreports/ In general, it is useful to shift the creation field φ † = φ + 1 . An explanation of implications of the shift, as well as a more detailed derivation of the path integral can be found in Ref. 59 . Thus, a simpler action is obtained The total action enters the path integral in the form exp A so that field theoretic expectations become and are calculated by performing the path integral over the bilinear part of the action A 0 , Eq. (2), and expanding perturbatively for the term with coupling sp 0 − ǫ(t) , which is unaccounted for in A 0 . In connection with the observation of oscillating neuronal avalanches, we choose ǫ(t) = sp 0 − As sin(νt) and the perturbative part of the action becomes Eq. (5), giving rise to a perturbation theory in small A.

Appendix: Derivation of moments
In order to keep a readable, consistent notation, we use the mathematical notation, e.g. E[N(t)] , when we characterize the statistical properties of N(t). When referring to field-theoretic calculations, we use the notation common in physics, i.e. O for the expectation of an observable O . Of course, the two are closely related, e.g. Appendix: Second moment. The second moment E[N 2 (t)] is determined field-theoretically on the basis of the square of the particle number operator, (a † a) 2 = a † 2 a 2 + a † a , which is easily expressed in terms of fields φ † and φ once in normal-ordered form, Diagrammatically, the first term is which is a convolution in direct time, Appendix: n-th factorial moment to first order. We define the n-th factorial moment of N(t) 16 as the function where the additional extinction −A sin(νt) , Eq. (4), is dealt with perturbatively about A = 0 . The zeroth-order correction g n (t 0 , t) corresponds to the sum of binary tree diagrams with one in-coming leg and n out-going legs. Since the zeroth-order correction does not include any perturbation in the extinction rate, the process is homogeneous in time and, therefore, it does not depend on the initial time t 0 , but rather on the amount of time t − t 0 .
In Ref. 16 , the function g n (t 0 , t) was derived from the following forest (meaning it contains a sum of tree-like diagrams) This zeroth-order correction satisfies the identity which can be obtained by choosing an arbitrary time t ′ ∈ [t 0 , t] when a diagram has ℓ legs, cf. Eq. (40a), and reconstructing the entire forest. As the left-hand side of Eq. (38) does not depend on t ′ it might be tempting to differentiate with respect to t ′ , but for our purposes the sum will prove particularly useful.
We proceed with the construction of the first order correction Ag (1) n (t) of g n (t) . This correction contains the external perturbation dt ′ As sin(νt ′ ) once in each diagram. It needs to be attached to every tree diagram at time t ′ , as indicated by the vertical dotted line in the diagram below, where the last identity is obtained by the substitution t ′′ = t − t ′ and using the time-homogeneity g  ℓ (0, t − t 0 ) . Defining the sum the first order perturbation can be written as n (t 0 , t) , Eq. (38), was it not for the factor ℓ in front of g (0) ℓ (0, t − t 0 − t ′′ ) , corresponding to the ℓ legs at time t ′ to choose from to insert the external perturbation. Nevertheless, the summation Eq. (41) can be carried out by using g ℓ (t 0 , t) and multiples of g ℓ (t 0 , t) according to Eq. (44) allows the summation to be carried out, using Eqs. (38) and (39), which eventually gives and therefore Eq. (12). For what follows, this is best re-written as with Appendix: further observables Appendix: covariance. The covariance is defined as where all expectations are conditioned to N(0) = 1 . In the field theory, the term E[N(t 1 )N(t 2 )] is calculated as where t max = max{t 1 , t 2 } and t min = min{t 1 , t 2 } . The first term of Eq. (49) is a convolution similar to Eqs. (34) and (35), while the second term is a product of two first moments, At t 1 = t 2 , when φ(t max ) φ(t min ) = 1 , Eq. (49) recovers Eq. (33).
exp −r(t − t 0 )(cos(νt) − cos(νt 0 )) + 1 r 2 + ν 2 r sin(νt) − ν cos(νt) + exp −r(t − t 0 )(ν cos(νt 0 ) − r sin(νt 0 )) . n (t 0 , t) , Eq. (37), may be written as g n (t 0 , t) = n!αβ n−1 with α = exp −r(t − t 0 ) and β = q 2 (1 − exp −r(t − t 0 ))/r, as shown in Eq. (42) in Ref. 16 . Using the first-order corrected g n (t 0 , t) , Eq. (46), the summation does not change significantly as far as the term Au(t 0 , t) is concerned, but A(n − 1)v(t 0 , t) requires some additional work. Recognising that Eq. (53) is effectively a generating function, the required sum is easily obtained from Eq. (53), To reduce clutter, we expand u(t 0 , t) and v(t 0 , t) only for the specific case of t 0 = 0 . Using for the zeroth order approximation, the survival probability to first order reads Appendix: Avalanche shape to first order. In Refs. 16,59 , it is shown that the avalanche shape (temporal profile) of a branching process can be related to the incompletely normalised average which is the average population size N at time t taken over the joint probability P(N(t) = N, N(T) = 0|N(t 0 ) = 1) of having size N at time t and size 0 at time T, conditioned to having size 1 at time t 0 . This quantity is easily captured in the field theory, Using identities (B2) and (B3) of Ref. 16 in the form and and writing g n (t 0 , t) = n!αβ n−1 and g n (t, T) = n!ab n−1 , Eq. (58) can be written at A = 0 after some tedious algebra To zeroth order in A, we recover the result in Ref. 16 with α = exp −r(t − t 0 ) , β = q 2 r 1 − exp −r(t − t 0 ) , as used in "Probability of survival to first order", as well as a = exp −r(T − t) and b = q 2 r 1 − exp −r(T − t) . To incorporate the first-order correction, we adjust α , β , a and b so that g n (t 0 , t) and �N(t), N(t 0 ) = 1, N(T) = 0� can be calculated order by order via generating functions. Defining (52) P s (t 0 , t) = 1 − exp −φ(t)φ † (t 0 ) = − ∞ n=1 (−) n n! φ n (t) φ(t 0 ) = ∞ n=1 (−) n−1 n! g n (t 0 , t) .