{'fraction_non_alphanumeric': 0.03957945787935158, 'fraction_numerical': 0.010878952963061387, 'mean_word_length': 4.627009345794392, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical phases in growing populations: understanding recovery from bottlenecks

Since steep declines in a population's size also typically alter its composition, population bottlenecks are considered highly important for evolution. However, despite such significance, the mechanisms governing the impact of a given population bottleneck remain poorly understood. In this context, we show that long-term post-bottleneck outcomes can depend crucially on the rate at which a diminished population grows whilst recovering. That is, two otherwise identical populations, each having undergone the same bottleneck, can fixate on dramatically different demographics due to different rates of post-bottleneck growth. This behaviour is moreover shown to change non-trivially with different levels of mutation. The underlying mechanism can be traced to intrinsic fluctuations whose standard deviation decreases in time with the inverse square root of the growing population size. Crucially, such decreasing fluctuations couple to the underlying dynamics and result in distinct regimes of (non-equilibrium) demographic behaviour, delimited by abrupt transitions at critical population sizes. Seen through this lens, growth and mutation alter evolutionary outcomes by changing the duration and character of dynamical phases; a feature that we speculate may be of generic importance across many systems.


I. INTRODUCTION
Population bottlenecks involve a steep decline in population size due to exogenous events, such as disease, changes to the climate, or population fracture (often referred-to as a founder event) [1,2]. They are widely accepted as an important facet in the modern understanding of evolution, and have been implicated in the reduction of both genetic and phenotypic variation across a variety of organisms, including viruses [3,4], song sparrows [5,6], tropical surgeonfish [7], elephant seals [8] and humans [9].
Despite such conceptual importance, the relationship between the short-term restriction of a population and its subsequent long-term evolutionary outcomes is still poorly understood, with a lack of both representative models and/or classifying phenomenology. This is particularly relevant since post-bottleneck recovery is typically characterised by rapid or sustained population growth, which several pioneering works have now shown can lead to a variety of ostensibly counter-intuitive behaviours. These include driving novel demographic transients [10]; facilitating cooperation [11,12], even when deterministic selection favours defection [13]; and altering the success of invading variants [14].
We therefore set out to disentangle the effects of a bottleneck itself-manifest as the population's initial size and demographic mix following the bottleneck-and those of post-bottleneck characteristics, such as the rates * These authors contributed equally to the manuscript. † To whom correspondence should be addressed: shashi@ncbs.res.in and r.g.morris@unsw.edu.au. of growth and mutation. To do this, we use evolutionary game theory [15][16][17] and, specifically, a growing variant of the otherwise well-studied Iterated Prisoner's Dilemma (IPD) under replicator-mutator dynamics [18,19].
The growing IPD has several appealing characteristics; it is at once minimal, facilitating both analysis and simulation, but it is also sufficiently rich as to act as an effective 'toy' system for understanding broader phenomenology. For example, the IPD includes two stable attractors, one of which is a limit cycle, which speaks to models of stable population diversity, such as that observed in Escherichia coli [20][21][22] and side-blotched lizards [23] that are understood via cyclic dominance [24]. The IPD also crudely captures the effects of mutation, which have so far been overlooked in the literature concerning growth, but are shown here to be crucial.
Ultimately, we demonstrate that the post-bottleneck rates of population growth and/or mutation critically determine how fixation probabilities are conditioned on initial population size and demographic mix. This means that even marginal differences to the rates of growth or mutation can result in dramatically different likelihoods of long-term outcomes, despite starting from the same founding population. For example, a bottleneck that restricts the population's composition to the location of one of the two stable attractors typically results in long-term fixation to the same attractor. However, for certain combinations of the rates growth and mutation, this is turned on its head, and long-term fixation to the other attractor occurs with almost certainty. This counter-intuitive and critical-like behaviour broadly interpolates between known behaviours in the deterministic, large population limit, where fixation is fully determined by initial conditions, and systems at finite population size, that cannot fixate due to persistent fluctuations.
The underlying mechanism, generic across well-mixed populations, is an intrinsic noise whose standard deviation decreases in time with the inverse square root of the growing population size, N . This ensures that, as a population grows, its dynamics will eventually fixate upon the stable attractors that characterise the system in the limit of very large population size. However, it also permits a diverging correlation time, so that long-term outcomes can depend on early events, even as t → ∞, despite the ostensibly randomising effects of intervening fluctuations [25]. Moreover, the rate at which fluctuations are suppressed turns out to be vital to understanding the statistics of long-term fixation.
Crucially, the decreasing scale of the fluctuations couples to the geometry of the deterministic attractors, manifesting in distinct regimes of demographic behaviour. These can be explained in terms of effective nonequilibrium phases, which capture the transient characteristics of fixed-N systems by averaging over large but finite times (hence excluding the disproportionate effect of certain highly rare events). Such phases are moreover shown to be delimited by transitions at critical values of N . This is reminiscent of the classical study of nonequilibrium phase transitions [26,27], which have previously been implicated across a number of areas, including directed percolation [28], self-organisation of particle suspensions [29], surface growth [30], epidemiology [31] and even hard sphere packing [32]. It is within this framework that we understand the importance of the rates of population growth and mutation, which together alter the character and duration of such regimes, and hence determine the long-term impact of bottlenecks.
Our results are organised as follows. In Section II A we introduce the growing IPD as a toy model for understanding recovery from bottlenecks. We demonstrate that there are only three possible outcomes in the longtime limit: fixation on one of the two attractors or extinction. In Section II B we compute the likelihood of such fixation as a function of bottleneck (initial conditions), and the rates of growth and mutation. This illustrates our first key message: growth and mutation can critically dictate the impact of a given bottleneck on evolutionary outcomes (long-term fixation probabilities). To understand these results, Section II C outlines our second key message: dynamical behaviour can always be characterised by one of three effective non-equilibrium phases, dependent on the population size. These are: 1. a stochastically-induced phase, at small population sizes, where state-dependent fluctuations dominate, and crossings between both basins of attraction are permitted; 2. an asymmetric phase at intermediate population sizes, where crossings between the two basins of attraction are overwhelmingly likely to be in one direction only, and; 3. a locked-in phase, at large population sizes, where FIG. 1. The growing frequency-dependent IPD: a toy model for bottleneck recovery via population growth. By de-coupling the rates of birth and death in the frequencydependent IPD [18], we adopt a model that combines the notions of frequency-dependence and mutation with that of population growth (see main text for definitions of rates). We ask: what is the long-term effect of a population bottleneck? Specifically, how does it influence evolutionary outcomes, and how does this change with the rate at which the recovering population grows?
escape from the basins of attraction is extremely unlikely.
Section II D then shows how a growing population exhibits these phases sequentially, in a manner controlled by the rates of growth and mutation. This suggests that our first message can be explained in terms of the second. Finally, Section II E validates this conjecture by satisfactorily reconstructing the outcome statistics of Section II B using a decomposition of conditional probabilities motivated by the observed dynamical phases.

A. Growing frequency-dependent IPD
The frequency-dependent IPD is a well-established evolutionary game involving three strategies. It abstracts Since overall population growth is exponential (with rate b − d) the system converges on canonical deterministic dynamics as t → ∞ (panel b, colours represent the direction of deterministic flows). This results in only three different evolutionary outcomes. Either, the population goes extinct (panel a in black, panel c), or the system fixates on one of the two deterministic attractors: the mixed-strategy limit cycle (panels a in blue, panel d) or the AllD fixed point (panels a in red, panel e). The statistics of these two latter outcomes are dictated not only by the founding population, but by fluctuations. These determine the likelihood of crossing the separatrix that marks the boundary between the two basins of the stable attractors (dashed, panel b). Lower-case font represents population fractions, e.g., n AllD = N AllD /N .
the key facets of biological evolution and population dynamics (e.g., birth, death, and mutation) into a tractable framework whereby the fitness of individuals is determined by their relative successes and/or failures when playing each other in the repeated game, the IPD.
In the classical formulation, each player of the repeated IPD assumes one of three strategies: "always cooperate" (AllC), who cooperate in every round; "always defect" (AllD), who defect in every round, or; "tit-for-tat" (TFT), who default to cooperation for the first round and then, at a small complexity cost, copy their opponents' moves thereafter. In a single repeated game, players accumulate payoff over the rounds according to the standard Prisoner's Dilemma rules: if both players cooperate, they receive a larger pay-off than if they both defect, but if one player cooperates and the other defects, then the defector receives the highest possible pay-off whilst the cooperator gets the lowest payoff. The accumulated payoffs for an m-round repeated game are encoded by the matrix [18]: where T > R > P > S and R > (T + S)/2. We use T = 5, R = 3, P = 1, S = 0.1, and c = 0.2.
This game is then set against a backdrop of birth, death and mutation (Fig. 1). With a rate b, the fraction of players born into a given strategy is proportional to the fraction of the total pay-off accumulated by that strategy when all players play each other. This is given is the average fitness of a given strategy when played against the whole population (including individuals of the same strategy). Here, Latin indices denote the three different strategies, i.e., AllC → 1, AllD → 2 and TFT → 3 such that the a ij are the components of the above payoff matrix and N i is the number of individuals playing each strategy. In our analysis, but not in simulations, we use the simplification that which, although including self-interactions, still retains all the relevant features associated with the IPD [19]. In either case, the mean fitness is ϕ = Rates of growth and mutation critically alter how a given bottleneck translates into fixation probabilities. For a range of initial conditions, growth and mutation rates, we use repeated Gillespie-Itô simulations to compute p AllD , the probability that a growing population will fixate on the AllD fixed point. Choosing the initial demographic mix as the mixed-strategy fixed point, we see that p AllD depends on both the rate of population growth and the rate of mutation (panel a). Conversely, fixing the rates of growth and mutation, p AllD also depends on both the initial demographic mix n0 (panels b, c, e & f ) and the initial population size, N0 (panel d & g). Behaviour can be qualitatively classified into two regions separated by the red dashed line (panel a). In the upper region, behaviour agrees with expectations: higher rates of growth more rapidly suppress O(1/ √ N ) intrinsic fluctuations and therefore increase the likelihood of fixation within the basin of attraction that the system started (panel a). Increasing the initial population size reduces the likelihood of large fluctuations at early times, therefore exacerbating this effect (panels b-d). In the lower region, behaviour is more complex: there are high likelihoods of an AllD outcome at both low (b − d 0.01) and intermediate (0.03 b − d 0.07) rates of growth (panel a). This results in a non-monotonic dependence of p AllD on b − d, which implies that population growth can actually increase the likelihood of crossing the separatrix between the two basins of attraction, despite more rapidly suppressing O(1/ √ N ) fluctuations. There is a similarly non-trivial structure to the effects of initial demographic mix on p AllD where, for certain n0, an increase in N0 actually causes p AllD to increase rather than decrease (representative initial states s and s , respectively, panels e-g). undergo mutation and are assigned a different strategy. Death also occurs at random, with a rate d.
Choosing a population birth rate b that is greater than the death rate d gives rise to unbounded exponential (i.e., Malthusian) population growth (Fig. 2a). This rapidly suppresses O(1/ √ N ) fluctuations and converges to wellstudied deterministic behaviour (represented on the unit simplex in Fig. 2a, inset). In particular, for values of mutation rate µ in the interval 10 −7.5 ≤ µ ≤ 10 −2.5 , there are two stable deterministic attractors [18]: a stable 'AllD' fixed point, where a small fraction of TFT players (who mutually cooperate) are exploited by a large population of AllD defectors, and; a stable limit cycle around an unstable 'mixed-strategy' fixed point ( Fig. 2b & Supplementary Information, Sec. 2). The latter is characterised by a three-phase cycle whose handedness is anti-clockwise in the traditional presentation of the state-space simplex; players of TFT can out-compete those playing AllD due to their capacity for mutual cooperation, however: they are then susceptible to invasion by players of AllC due to the complexity cost, whereby; AllC players can be exploited by those playing AllD, completing the cycle.
As a consequence, only three outcomes are possible as t → ∞.

B. Statistics of evolutionary outcomes
Since the growing IPD is both non-linear and timeinhomogeneous, it resists most standard approaches to probabilistic analysis. Computing the statistics of the aforementioned t → ∞ outcomes-i.e., the fixation probabilities-therefore involves using a high performance computing facility [33] to perform stochastic simulations.
Specifically, we employ a hybrid Gillespie-Itô approach (see Supplementary Information, Sec. 1) to approximate the following fixation probability: }. This is the limiting probability, after long times, that the system converges to the AllD fixed point, given specified post-bottleneck initial conditions, and conditioned on populations that do not become extinct.
The results demonstrate several interesting features. For instance, choosing n 0 to be the unstable mixedstrategy fixed point, we see that p AllD depends on both the rate of population growth and the rate of mutation (Fig. 3a). Similarly, fixing the rates of population growth and mutation reveals sensitivity to initial conditions, where p AllD depends on both initial demographic mix, n 0 , and initial population size, N 0 (Figs. 3b-g).
The parameter space of mutation and growth rates can moreover be divided into two qualitative regions (panel a, red dashed line).
In the upper region, behaviour agrees with the expectation that faster growth rates reduce O(1/ √ N ) intrinsic fluctuations more rapidly and therefore increase the likelihood of fixation within the same basin of attraction that the system started (panel a). Similarly, as the initial population size, N 0 , increases, the system experiences fewer large fluctuations at early times and hence this also increases the likelihood of fixation within the starting basin of attraction (Figs. 3b-d).
In the lower region, however, behaviour is more complex.
There are high likelihoods of an AllD outcome at both low (b − d 0.01) and intermediate (0.03 b − d 0.07) rates of growth (panel a). This results in a non-monotonic dependence of p AllD as a function of population growth. That is, despite more rapidly reducing the fluctuations that are ostensibly required to cross the separatrix between mixed-strategy and AllD basins, growth can actually increase the likelihood of fixating on the AllD fixed point. Similarly, rather than decreasing the likelihood of crossing the separatrix and fixating there, increasing N 0 can actually increase this likelihood for certain initial states, confounding expectations regarding the role of fluctuations at early times when populations remain small (Figs. 3e-g).
In the context of population bottlenecks, this demonstrates that the long-term ramifications of reducing a population to a particular size and demographic mix can depend, critically, on the rates of post-bottleneck growth and mutation.

C. Effective non-equilibrium phases
To understand the non-trivial behaviour in Fig. 3, we repeatedly simulate the fixed population size IPD (Supplementary Information, Secs. 1 & 7), computing a socalled empirical distribution [34]. Specifically, we calculate the mean fraction of time spent in the AllD basin, where 1 is the indicator function, N denotes the size of the ensemble over which the average is taken, and the integration is understood in the Itô sense. The time T > (log N max ) /(b − d) exceeds the entire duration of our growing simulations (which we stop at N max = 10 10 ) and represents an effective cutoff, so that the statistics of fixed size simulations are not skewed by events that are highly unlikely to occur in the growing simulations (i.e., with characteristic rates 1/T ). Such a large but finite T therefore aims to capture the average transient behaviour of a growing population at a particular N .
Our results (Fig. 4a) are suggestive of a large deviation principle, such that p(F | N ) exp [−N I µ (F N )], where I µ is a convex rate function. Although determining the precise functional form of I µ is considered out-of-scope for this article, our data suggests that it has only three zeros, despite varying N over ten orders of magnitude: This means that, depending on the population size, there are three statistically distinct types of characteristic demographic behaviour. Due to the finite size of T , we call these effective non-equilibrium phases. The three effective phases are characterised as follows: 1. Stochastically-induced phase-Demographic trajectories at small population sizes are characterised by large intrinsic fluctuations and an intermediate value of F ( Fig. 4a in magenta, Figs. 4b & c). Fluctuations are both correlated and statedependent; features that are captured by the symmetric 3 × 3 correlation matrix, B † ij , that can be obtained by performing a Van Kampen system-size expansion [35] and projecting the results onto the unit simplex using Itô's lemma ( Fig. 4b & c, orange crosses, and Supplementary Information, Sec. 5). The frequency-dependent nature of births means that fluctuations at the centre of the simplex are large and isotropic, whilst the components normal to the boundaries decrease rapidly as the edges and corners are approached. This gives rise to stochastically-induced effects [36,37], where fluctuation gradients bias stochastic trajectories, driving them towards the simplex edges and corners, on average ( Fig. 4b & c and Supplementary Information Secs. [3][4][5]. Despite such overall behaviour, B † ij is not symmetric under the interchange of n AllC , n AllD and n TFT , and stochastic trajectories retain characteristics encoded by the payoff matrix; including a bias for anti-clockwise motion, and a comparatively low likelihood of reaching the AllC corner (when compared to AllD and TFT corners). These behaviours are crucial to understanding the precise µ-dependent mechanisms that underpin the value F ≈ 0.5. . Beyond this regime, behaviour becomes increasingly deterministic (green half-arrows, panels e-d). For intermediate population sizes, behaviour is asymmetric: demographic trajectories can cross the separatrix from the mixed-strategy limit cycle, but not from the AllD fixed point (panel d). The upper critical population size of this regime scales as ∼ 1/µ (white, panel a). For large enough N , behaviour becomes increasingly deterministic, and trajectories remain locked-in to the mixedstrategy limit cycle (beige, panel e). At low µ, the separatrix is typically crossed where it intersects either the TFT or AllD edges, which is where the magnitude of stochastic effects, B † F , are largest relative to deterministic flows, |Ai|, with · F and | · | denoting Frobenius-and 2 -norms, respectively (panels c, d & f). Representative trajectories and population fractions (bars) only show a fraction of the total time simulated (see bar legend) in order to aid visualisation. Crossing statistics (circles) are taken from single simulations lasting 10 6 s. try of the underlying attractors becomes increasingly important. In particular, the system enters an asymmetric regime at populations above N ≈ 10 2.7 , for which F ≈ 1 (Fig. 4a, yellow). Here, state-dependent fluctuations permit the system to cross from the mixed-strategy limit cycle to the AllD basin, but not from the AllD fixed point to the basin of the limit cycle (Fig. 4d). In other words, once the separatrix has been crossed, trajectories are extremely unlikely to come back within the time T .
3. Locked-in phase-Once N is sufficiently large, fluctuations are small and demographic trajectories are effectively locked into the mixed-strategy basin for times < T , implying F ≈ 0 (Fig. 4e, dark gray).
The aforementioned effective phases moreover couple to mutation, which alters both the character of the stochastically-induced phase (Figs. 4b & c), and the population size at which the system transitions from asymmetric to locked-in phases (Fig. 4a, white dashed line).
The origin of the former is that the standard devia- The critical scaling, at which the footprint occupies the entire mixed-strategy basin, is N ∼ 1/µ, which agrees with ensemble statistics (panel e). Since this corresponds to onset of the locked-in phase, the corollary is that the asymmetric phase has a duration that is µ-dependent (cf. Fig. 4). . This not only exacerbates stochastically-induced effects, but also increases the mean residence times associated with the corners. In particular, at comparatively high levels of mutation, residence times are less than T , which results in a stochastic cycling between TFT and AllD corners (recall that there is a lower likelihood of finding the AllC corner) (Fig. 4b). Since these two corners have comparable mean rates of escape to the opposite basin of attraction (Supplementary Information, Sec. 8), F takes a value of approximately 0.5. By contrast, lower rates of mutation imply dwell times greater than T (Fig. 4c). On average, therefore, the system will find either the AllD or TFT corner and then remain there. Here, the value of ≈ 0.5 results from the position of the mixed strategy fixed point, and the correspondingly equal probability that trajectories are expelled to either the AllD or TFT corners (typically, via the AllC-TFT edge, see Fig. 4f ).
For the latter, the ∼ 1/µ dependence of the asymmetric to locked-in transition can be understood in terms of the stochastic 'footprint' of evolutionary trajectories in the locked-in regime-i.e., those that repeatedly (and stochastically) navigate the limit cycle (Fig. 5). In particular, due to the finite nature of T , the convex hull of this footprint is well-defined, reflecting the shape of the limit cycle at a different value of µ ( Fig. 5a-d). The result is a family of scaling relations, where the footprint of small populations with high levels of mutation (i.e., high noise, small limit cycle, Fig. 5c) is approximately equivalent to that of large populations with low levels of mutation (i.e., low noise, large limit cycle, Fig. 5d). The critical scaling that defines the onset of the regime occurs when the stochastic footprint fills the mixed-strategy basin ( This behaviour can be recast as a size-dependent antagonistic relationship between mutation and intrinsic noise, where the former favours population heterogeneity (expelling towards the simplex boundaries) and the second homogeneity (attracting towards the centre of the simplex). This is particularly important in growing populations, since the balance between the two factors Growth and mutation control the relative duration of dynamical regimes. By virtue of the sizedependent effective non-equilibrium phases (Fig. 4), a growing population exhibits three distinct dynamical regimes in sequence. Changing the overall rate of growth therefore rescales the duration of the stochastically-induced and asymmetric phases. Changing the rate of mutation, by contrast, alters only the asymmetric phase. The duration of the asymmetric phases is important to overall fixation probabilities, since a longer duration implies a higher likelihood of crossing the separatrix, from which there is no return. changes over time.

D. Dynamical regimes in growing populations
As a population grows, it exhibits the three effective phases in sequence, with transitions occurring abruptly at critical population sizes (Fig. 6). Thus, despite the population size steadily increasing, dynamical behaviour manifests in three distinct regimes.
Seen through this lens, it is clear that changing the rate of growth changes the duration of the stochasticallyinduced and asymmetric phases by an overall factor (the locked-in phase is, in principle, open-ended). By contrast, the rate of mutation changes only the duration of the asymmetric phase (by virtue of the µ-dependent transition to the locked-in phase).
The asymmetric phase is especially important for determining fixation probabilities. If the separatrix is crossed from the mixed-strategy basin during this phase, then the system remains in the basin of the AllD fixed point until the onset of the locked-in phase, where it remains for all until the population reaches N max = 10 10 . Therefore, the longer the duration of the asymmetric phase, the greater the likelihood of crossing the separatrix.
However, despite the appeal of this heuristic, we must also account for the fact that this likelihood is also con-ditioned on the state at which the system enters the asymmetric regime. This is determined by the duration of preceding stochastically-induced regime, which is set by the rate of growth. It is also determined by the character of that regime, which is set by mutation. Despite the presence of two µ-dependent mechanisms in the fixed-N ensemble (see previous Section) it should be stressed that, in a growing system the duration of the stochastically-induced regime is sufficiently short in comparison to corner dwell times that, for all but pathologically slow growth rates or very large µ, the separatrix is only ever likely to be crossed once. In turn, this puts greater emphasis on the role of the system's initial conditions.
To test our understanding of this complex interplay, we construct an approximation to the full outcome statistics. This is based on a decomposition in terms of conditional probabilities associated with each phase, and four 'equivalence classes' of states at the start of the asymmetric phase.

E. Equivalence classes
Consider the conditional probability p asy AllD (s) = Pr n t lock ∈ AllD-b | n tasy = s , where avoidance of extinction is now assumed implicitly. That is, the likelihood of being in the AllD basin at the onset of the locked-in phase, t lock , given that the system was in a state s at the onset of the asymmetric phase, t asy . Computing this probability via stochastic simulation demonstrates the existence of three approximate equivalence classes, S i , such that p asy AllD (s i ) is agnostic to the demographic mix s i ∈ S i , at the onset of the asymmetric phase ( Fig. 7a and Supplementary Information, Fig. 11). These are: those states in the AllD basin, S 1 ; those along the AllD-TFT edge that stretch from the separatrix to the TFT corner, S 2 , and; those along the TFT-AllC edge (excluding the TFT corner) on the mixed-strategy side of the separatrix, S 3 . The remaining states of the mixed-strategy basin are labelled S 4 . Whilst these do not form an equivalence class, we assume (and later show) that they only minimally contribute to overall fixation probabilities.
The existence of equivalence classes prompts the following simplification (Supplementary Information, Sec. 9): where p stoc Si = Pr n tasy ∈ S i | N 0 (again, with avoidance of extinction assumed implicitly), and states s i can be chosen arbitrarily from S i . This approximation allows us to verify our heuristic understanding of how the rates of growth and mutation impact fixation probabilities by controlling the duration and stochastic character of dynamical phases and therefore the likelihood of (stochastic) behaviours, such as crossing the separatrix, that are crucial in dictating long-term outcomes. It also dramatically reduces the computational time needed to calculate p AllD for a range of different initial conditions, since we only need to recalculate the four likelihoods that the stochastically-induced regime finishes in each of the equivalence classes, respectively. The conditional probability p asy AllD , by contrast, does not depend on the specific initial conditions (but rather in which equivalence class the system is at time t asy ).

Growth and mutation
Fixing the mixed-strategy fixed point as the initial demographic mix, we can use Eq. (6) to deconstruct the dependence of p AllD on the rates of growth and mutation (Fig. 7).
In the stochastically-induced regime, fluctuation gradients 'drive' trajectories from the centre of the mixed-strategy basin towards the simplex edges and then the corners. As a result, the growth-dependent (average) duration of the regime, t asy = log 10 2.7 /N 0 /(b − d), dictates in which class the trajectories are likely to start the asymmetric phase: rapid growth rates are required to confine trajectories that end in the mixed-strategy basin to the S 4 region, whilst intermediate and slow growth rates suffice for the S 3 and S 2 regions, respectively (Figs. 7c-e). Although the average time spent in the two basins by such trajectories during the stochasticallyinduced regime is independent of µ (Fig. 4a), the likelihood that the system is in S 1 at time t asy actually increases with decreasing µ (Fig. 7b). The reason is that decreasing µ changes the shape of the separatrix, therefore reducing the the size of the S 3 region ( Supplementary  Information, Sec. 10).
The asymmetric regime, by contrast, has a duration that is both growth-and µ-dependent: τ µ = t lock −t asy = − log µ  Fig. 3a-fixation probabilities can rely critically on initial conditions. That is, initial demographic mixes that are 'close' (and not necessarily near the separatrix) can give rise to dramatically different long-term fixation probabilities. This is particularly striking at low levels of mutation (log µ = −6.5) and comparatively high initial population size (N0 = 250) (cf. Fig. 3f). The reason for this is twofold. First, larger initial populations reduce the effect of the stochastically-induced phase, which is to expel populations to the simplex boundaries and the AllD and TFT corners (panels c & d). Second, the duration of the asymmetric phase is significant and cannot be neglected (middle row). This implies that trajectories finishing the stochastically-induced phase in classes S3 or S4 can still cross the separatrix and fixate on the AllD fixed point (panels g & h). Whilst this provides a heuristic understanding, the contribution to p AllD from the S4 region violates one of the assumptions of (6), and therefore the error increases with N0 for the stated values of µ and b − d (panel j and Supplementary Information, Figs. 19 & 20).
the separatrix hinges, principally, on the likelihood of avoiding the TFT corner and its associated large confinement times (Figs. 7g-i). For example, τ asy µ must be extremely long in order to permit crossings from the S 2 region, since populations starting the asymmetric regime from this region encounter the TFT corner with almost certainty. Crossings from the S 4 region, however, occur at more modest τ asy µ -achieved by either low growth and high µ, or modest growth and low µ-reflecting the possibility that trajectories might avoid the TFT corner. Those from the S 3 region can happen at the smallest τ asy µ , since there is a high likelihood that trajectories will avoid the TFT corner (recall the anti-clockwise dynamics) and the ratio of the magnitude of stochastic effects to the magnitude of the deterministic flow is large in S 3 (see Fig. 4f ).
Combining these conditional probabilities using Eq. (6) satisfactorily reproduces the overall statistics of demographic outcomes (cf. Fig. 7n and Fig. 3a): the difference between our approximation and the full simulations have a mean value of 0.027, when averaged over growth and mutation rates, and a maximum value of 0.099 (Fig. 7o). This also confirms our assertion that the trajectories that start the asymmetric phase from the S 4 region do not impact long-term outcomes. The reason is that this only happens with significant likelihood when growth rates are high, which simultaneously ensures that such trajectories never cross the separatrix (Figs. 7e, i & m). The same decomposition is shown for initial demographic mixes other than the mixed-strategy fixed point in Supplementary, Sec. 11.

Initial conditions
The approximation in Eq. (6) further provides insight into the founder-like dependence of fixation probabilities on initial conditions ( Fig. 8 and Supplementary Information, Figs. [16][17][18]. At high values of µ, the asymmetric phase has negligible duration and behaviour is trivial (Figs. 3b & c, Supplementary Information, Figs. 16 & 17). At low values of µ, however, the asymmetric regime cannot be ignored and has a significant bearing on fixation probabilities.
In this case, if the initial population size is small (Fig. 3e, Supplementary Information, Fig. 18), then the stochastically-induced phase is sufficiently long as to expel trajectories to the AllD or TFT corners from the AllD or mixed-strategy basins, respectively. Only the former trajectories impact fixation probabilities, however, since for all but the slowest growth rates, trajectories stuck in the TFT corner have a residence time longer than the duration of the asymmetric regime.
By contrast, if the initial population size is large (Fig. 3f and Fig. 8), then the duration of the stochastically-induced phase is not long enough to expel trajectories to the boundaries and/or corners, resulting in a non-zero likelihood of starting the asymmetric regime from either the S 3 or S 4 regions. For the former, there is a high probability of crossing to the AllD basin during the asymmetric regime, since the deterministic flows direct demographic trajectories towards the separatrix. For the latter, this probability is much lower, since trajectories are more likely to be entrained to the limit cycle. Nevertheless, the small contribution that results from the S 4 region breaks one of the assumptions on which Eq. (6) is based, which also explains why the error in our decomposition increases with N 0 for certain values of µ and b − d (the mean value of the error in Fig. 8j is 0.079, while the maximum value is 0.285). Of note, the values shown in Fig. 8 are the worst case of those we have simulated ( Supplementary Information, Sec. 12).

III. DISCUSSION
Using a growing variant of an iconic model of evolutionary game theory, we have demonstrated the existence of non-trivial fluctuation-mediated effects in growing populations, whereby the rates of growth and mutation critically determine how fixation probabilities are conditioned on initial population size and demographic mix. This has ramifications for the understanding of population bottlenecks and their long-term impact. The implication is that the population growth commonly associated with post-bottleneck recovery can, in fact, be as important as the effects of the bottleneck on population size and demographic mix.
Our findings apply to populations that are well-mixed, and as such pertain to systems that combine short-range interactions with a mechanism for mixing that is fast on the timescales of the population dynamics, or those that otherwise have effectively long-range interactions, either explicitly or via the mutual interaction with public resource (although some of these assumptions have recently been brought into question [38]).
In this context, we follow several pioneering works that have characterised various effects of growth in well-mixed populations [10][11][12][13][14]. Although they do not explicitly examine the role of initial conditions, the behaviours appearing in two of these works [13,14] are related to the those reported here. As are several studies concerning growing systems of binary 'spins' [25,39,40].
The shared mechanism in all these cases is the pres-ence of O(1 √ N ) intrinsic fluctuations that decrease as the system grows. This has two main ramifications: first, ergodicity is broken, implying that a population fixates even in the presence of non-zero mutation, and; second, there is a decreasing scale by which fluctuations can couple dynamically to the geometry of underlying deterministic attractors of the system. It is the combination of these effects that gives rise to our headline behaviour: fixation probabilities that depend non-trivially on initial conditions (or otherwise stochastic events at early times), even as t → ∞.
Such fluctuation-mediated effects are further characterised by a dependence on the rates of growth and mutation, which we explain by putting forward the notion of effective non-equilibrium phase transitions, and showing that these delimit distinct demographic regimes in our model. The rate of growth controls the overall rate at which intrinsic fluctuations decrease, and therefore also the relative duration of such demographic regimes. Mutation, by contrast, has two related effects. First, it changes the nature of the underlying deterministic attractors (the size of the mixed-strategy limit cycle, in our case) and hence the structure to which decreasing fluctuations couple. Second, it also changes the statedependence of the fluctuations. Whilst the former impacts the character of the initial stochastically-induced phase, it is the combination of both of these effects that sets the ∼ 1/µ dependence of the critical transition between the intermediate asymmetric and the final lockedin phases.
The latter behaviour is an example of a potentially interesting and unexpected antagonistic relationship between two sources of stochasticity: mutation and the intrinsic effects of finite-sized populations. In our model and other studies, mutation promotes heterogeneity, whilst intrinsic fluctuations typically drive the system towards homogeneity [39]. This is particularly important in the context of population growth, since the effects of mutation do not depend on the population size, whereas the stochastically-induced forcing due to intrinsic noise decreases with increasing population size. As a result, we speculate the other growing systems may also exhibit mutation-dependent critical transitions, where the effects of intrinsic fluctuations are balanced by those of mutation.
More generally, for growing well-mixed populations with all but the simplest of deterministic attractorsi.e., fixed points, limit cycles and stable manifolds etc.the implication is that growth may be synonymous with effective non-equilibrium phase transitions. Understanding whether and how these fit into the existing literature is an open question. Of potential interest is the ongoing challenge to classify non-equilibrium phase transitions by their universality classes [27]. Work has been undertaken to describe single absorbing state transitions (of the directed-percolation type) [26,41,42] and also symmetric absorbing state transitions [43], but there appears to be very little literature on asymmetric absorbing state transitions.
Theoretical considerations aside, we put forward that our ideas may be examined within the context of directed evolution [44,45]. Here, another, albeit direct, interplay between growth and mutation has already been demonstrated: mutations occurring at the genetic loci associated with growth-control promote so-called genetic instabilities [46]. Our results also appear relevant to stateof-the-art in silico representations of directed evolution, where the role of intrinsic fluctuations during growth has so far been overlooked [47]. A further setting that may prove relevant is that of quasispecies such as viruses.
Viruses have both high rates of mutation and growth (cf. Fig. 3a, upper region), and appear to be largely unaffected by the many population bottlenecks involved in host-to-host transmission, as well as its intra-host (e.g., plaque-to-plaque) spreading [3,4,48].
Nevertheless, bridging the gap between the abstract setting of the present work and the aforementioned applications will undoubtedly involve significant work. The extent to which this will be possible remains an open question, and may hinge on features that are not included in our model, such as spatial structure and/or other physical constraints [49][50][51][52]. So-called "patch" models of interacting locally well-mixed subpopulations is one avenue that may prove promising. Exploring how the ideas set out here translate across a wider class of systems is therefore an important avenue of future research, and we welcome further work in the area.

Gillespie-Itô simulations
The Gillespie algorithm [53,54] allows the exact simulation of the stochastic dynamics of the growing IPD when N is small, i.e., when such dynamics cannot be approximated by simply integrating the SDEs (see Secs. 3 & 4). The downside of this algorithm is that its computational time scales linearly with N , becoming impractical as the population grows. We therefore adopt a hybrid approach: when N is smaller than a chosen threshold the system's dynamics are simulated with the Gillespie algorithm, and when the threshold is exceeded they are simulated by numerically integrating the SDEs (Euler-Maruyama). Since the latter method results in values of N i along the real line, rounding is required if stochastic fluctuations trigger a switch back to the Gillespie algorithm.
Our choice for the algorithm switch threshold is N = µ −1 . This is motivated by the existence of the locked-in phase, which begins when N ≈ µ −1 (see Fig. 4). During this phase, the risk that the SDE approximation would lead to an 'accidental' crossing of the separatrix is extremely small. A minimum threshold of N = 10,000 is used for simulations with large µ.
All results were obtained using the High Performance Computing facility Katana [33]. The results in Fig. 3a  The fixed size dynamics were also simulated using the hybrid Gillespie-Itô approach. The Gillespie algorithm was modified by setting d = 0 and 'killing' a randomly selected player at every birth. The Itô part of the algorithm involved replacing the noise correlation matrix (see details in Sec. 7). The results in Fig. 4 were obtained by running the fixed size Gillespie-Itô algorithm 1000 times for 20,000 seconds for each combination of µ and N .

The deterministic limit
In the N → ∞ limit, the behaviour of the growing Iterated Prisoner Dilemma described in the main text is equivalent to that of the following continuous, deterministic system: a b c d log 10 μ=-6.9 log 10 μ=-5.1 log 10 μ=-3.9 log 10 μ=-2. where n i = N i /N . These equations can be easily derived from Fig. 1. As shown in Fig. 9 (see also [18]), the mutation rate has a clear effect on the dynamics of the continuous system. For 10 −7.5 ≤ µ ≤ 10 −2.5 , the system exhibits two stable attractors: a fixed point characterised by a population of almost entirely AllD players and a stable mixed-strategy limit cycle. The most notable effect of varying µ within this range is the resizing of the limit cycle: the smaller the mutation rate, the larger the limit cycle (panels a-c and e-g). For large values of µ (approximately 10 −3.1 ≤ µ ≤ 10 −2.5 ), the limit cycle collapses into a point, although the separatrix still exists (panels d and h). Varying µ also slightly changes the shape of the separatrix, especially in the proximity of the AllC-TFT edge of the simplex.

System-Size Expansion
Following [35,55], we derive a set of coupled SDEs that approximate the dynamics of the underlying protocol when N 1. In terms of step-operators, the Master-equation has the form: Expanding the step-operators in the usual fashion, and retaining only the leading and next-to-leading order terms, gives where r = ±1. Substituting Eq. (9) into Eq. (8) results in an equation of the Fokker-Planck type [55]: After some manipulation, it can be shown that and Equation (10) implies an SDE for the variables N i : The deterministic part of Eq. (13) is equivalent to Eq. (7). The noise sources that appear in Eq. (13) have zero mean (i.e., ξ i = 0) and are correlated according to: Summing over index i in Eq. (13), we have where ξ is a single source of zero-mean Gaussian white noise.

Multiplicative delta-correlated noise
Eqn. (13) can be re-written in terms of delta-correlated noise and multiplicative pre-factors. To do this we must choose a matrix b that satisfies B = b T · b. For b to be square, this requires a Cholesky decompostion, and ensures no more independent noise sources than there are variables in the system. However, a more natural approach is to decompose B according to the rules set out in Gillespie's CLE approach [54]. This results is SDEs of the following form: where the index α runs from 1 . . . 12 and correspond to the reactions in Fig. 1. The twelve independent noise sources each have mean zero, and are delta-correlated-i.e., η α = 0 and η α (t)η β (t ) = δ αβ δ(t − t ). The matrix b T is given by

Projected dynamics
We wish to project the dynamics onto the unit simplex-i.e., in terms of variables n i = N i / i N i . For this, we require the multivariate form of Itô's lemma [55]. For finite N , we have where the shorthand ∂ i = ∂/∂N i has been used. First, we deal with the deterministic parts. Using alongside (11), it can be shown that Similarly, using Notice that, when summed over i, both (20) and (22) are zero by virtue of the fact that i f i n i /ϕ = 1. Also, trivially, (22) goes to zero in the deterministic N → ∞ limit. For the stochastic part of (18), define a new matrix b † αi = j b αj ∂ j n i such that (18) can be re-cast in terms of correlated noise sources-i.e., To understand the impact of these correlations, (24) can be computationally decomposed into an eigenbasis for different values of n i . This always results in one zero-eigenvalue eigenvector pointing perpendicular to the simplex. The remaining in-simplex eigenvalues reveal that populations towards the centre of the simplex experience large uncorrelated fluctuations whilst, closer to the simplex boundary, correlations suppress fluctuations in the direction normal to the boundary (Fig. 4b-c, orange crosses). Moreover, the magnitude of along-boundary fluctuations decrease as a corner is approached.

Boundary effects
Equation (23) can be evaluated at the simplex edges, where we are particularly interested in the both deterministic drift and fluctuations in the direction of the bulk, which is captured by the dynamics of the strategy who's concentration is zero along a given edge.

Langevin equation for the fixed size model
The dynamics of the system were were simulated using the hybrid Gillespie-Itô approach (see Sec. 1). A fixed size version of the system was also used in our analysis (see Figs. 4 & 5): The Gillespie algorithm was modified by setting d = 0 and 'killing' a randomly selected player at every birth, while the Itô part of the algorithm involved replacing the matrix b T in Eq. (17)

Escape time from corners in the stochastically-induced regime
Here we show the statistics of the time t † that is needed by the fixed-size system to cross the separatrix for the first time (Fig. 10). We compare the cases of the systems starting from the AllD corner (panels a and c) and the TFT corner (panels b and d). In both cases, t † is proportional to µ: in the order of thousands of seconds when log 10 µ = −3.5 (panels a-b) and in the order of millions of seconds when log 10 µ = −6.5 (panels c-d). However, the influence of the population size N 0 on t † is different in the two cases. On the one hand, the crossing time starting from the AllD corner increases very rapidly when N 0 approaches the onset of the asymmetric regime at around N 0 = 10 2.7 (panels a and c), indeed showing that crossing the separatrix from the AllD basin to the mixed-strategy basin becomes extremely unlikely in the asymmetric regime. On the other hand, the crossing time starting from the TFT corner decreases with N 0 .

Equivalence classes decomposition
Consider a trajectory of the growing IPD that starts from an arbitrary state in the stochastically-induced phase (i.e., N 0 < 10 2.7 ) at time t = 0, enters the asymmetric phase at time t = t asy and the locked-in phase at time t = t lock . Due to the nature of the locked-in phase, the evolutionary outcome at t → ∞ is already known at time t lock , and thus Eq. (3) can be simplified: where avoidance of extinction is implicitly assumed and AllD-b indicates the AllD basin of the simplex. Moreover, since the growing IPD is a Markovian process, p AllD can be decomposed as follows: where U is the set of all possible states of the system at time t asy . Eq. (30) is impractical, since it requires the numerical estimation of the probabilities Pr{n t lock ∈ AllD-b | n tasy = s} and Pr{n tasy = s | N 0 } for all s in the very large set, U . However, the approximation given in Eq. (6) can be made by inspecting Pr{n t lock ∈ AllD-b | n tasy = s} in Fig. 11a A first area of the simplex, S 1 (outlined in red), corresponding to the AllD basin, can be immediately identified as it is characterised by a homogenous probability of AllD outcomes of approximately 1. If the system is anywhere within S 1 at time t asy , then it is expected never to move to the mixed-strategy basin, since the crossing of the separatrix in this direction is extremely unlikely. Any point s 1 would be good candidate for representing the entire area S 1 , however, we chose the AllD corner (red dot in Fig. 11b) as the system is in the proximity of such point for the vast majority of the time spent in S 1 .
A second area can similarly be identified: S 2 (orange), corresponding to the mixed-strategy part of the AllD-TFT edge (including the TFT corner) and characterised by a homogenous probability of AllD outcomes of approximately 0. If the system is in this area at time t asy , then it is very unlikely to cross the separatrix because of the long time (inversely proportional to µ) spent in the TFT corner while growing (i.e., as the fluctuations become smaller and smaller). Again, any point s 2 can represent the area S 2 , but we chose the TFT corner (orange dot in Fig. 11b).
A third area, S 3 (yellow), corresponding to the mixed-strategy part of the AllC-TFT edge with exclusion of the TFT corner, also stands out: here we see the probability of AllD outcomes quickly increases with the fraction of AllC players. Due to the anti-clockwise dynamics, in S 3 the system directed towards the point where the separatrix a b meets the AllC-TFT edge, which is where a move from the mixed-strategy basin to the AllD basin is most likely (see Fig. 4f ). Since the probability of AllD outcomes is less homogeneous, an approximation must be made for s 3 . The dynamics along the AllC-TFT edges are slow, and they get even slower close to the corners. We observe that a good proxy for the average position of the system over time during a climb of the AllC-TFT edge from the TFT corner to the separatrix is around n AllC = 0.25 and n TFT = 0.75 (yellow dot in Fig. 11b).
Within the remainder of the simplex, S 4 (blue), Pr{n t lock ∈ AllD-b | n tasy = s} can be very heterogeneous, however, as discussed in the main text, this contributes very little to the outcome statistics for small values of N 0 . We chose s 4 = MS-fp to represented this area (blue dot in Fig. 11b).  (Fig. 3). The conditional probability decomposition is always very accurate except for the case of low mutation rate (e.g., log 10 µ < −6), growth rates around b − d ≈ 0.05 and large initial populations (e.g., N 0 > 200), for which it begins to produce less accurate results (Fig. 8 panel o).

Statistics of the stochastically-induced regime in detail
In this section we describe the effect of the mutation rate on the statistics of the stochastically-induced regime outcomes more in detail. Fig. 12a shows the same probabilities p stoc Si shown in Figs. 7b-e, but for a single value of b − d = 0.026 (we remind that n 0 = MS-fp and N 0 = 128). We can see that all probabilities depend on µ for high mutation rates (i.e., log 10 µ > −4), but only p stoc S1 and p stoc S3 depend on µ for lower mutation rates. For a more detailed analysis, we decompose p stoc S1 into p stoc , the probability of n tasy being anywhere else within the AllD basin. We can now see that p stoc S AllC-TFT 1 decreases with µ, in contrast with p stoc S3 , which instead increases with µ. This is explained by the shape of the separatrix changing with µ: the point where the separatrix meets the AllC-TFT edge moves towards the TFT corner as µ decreases (see panels c & d).
At the same time, we also see that the probability of the system being in one of the two corners at time t asy , p S AllD-c 1 and p S AllC-c 1 , are independent on the mutation rate for log 10 µ < −5. However, for higher mutation rates it becomes more and more likely for the system to be in the AllD corner rather than in the AllC corner (see also panels c & d).
Finally, for mutation rates log 10 µ > −5 the probability of the system being away from the boundaries of the simplex p stoc S 1 into cases of the AllD corner, the AllC corner, the AllD basin side of the AllC-TFT edge (excluding the AllC corner) and the rest of the AllD basin. Panels c and d show the location of nt asy within the simplex for the same 1000 stochastic trajectories used for estimating the probabilities in a and b, for the cases of log 10 µ = −6.5 and log 10 µ = −3.5 respectively. The size of the circles is proportional to the occurrences of nt asy in the location and the colours reflect the different equivalence classes Si (see panel a).

Outcome statistics over growth/mutation rates: more initial mixes
Here we show that our decomposition of the outcome statistics based on equivalence classes holds not only for n 0 = MS-fp but in general. We consider other three cases that are very different: n 0 = centre, i.e., n AllC = n AllD = n TFT = 1/3 (Fig. 13); n 0 = TFT-c, i,e., n AllC = n AllD = 0, n TFT = 1 (Fig. 14); and n 0 = AllD-c, i,e., n AllC = n TFT = 0, n AllD = 1 (Fig. 15). In all figures, comparing the results obtained via decomposition (panels n) against those obtained via full simulation (panels p) yields a low error (panels o), demonstrates the validity of our method. The average error over all the considered values of growth and mutation is approximately 0.029 for the n 0 = centre, 0.003 for n 0 = TFT-c, and 0.001 for n 0 = AllD-c. The statistics of evolutionary outcomes can be decomposed in terms of conditional probabilities that are based on four equivalence classes of states (panels a-i). The landscape of p AllD (AllD-c) for different rates of growth (b − d) and mutation (µ) can be reconstructed with good agreement (panels j-n and panel o). The initial population size is N0 = 128. tasy and t lock are derived from the critical population sizes identified in Fig. 4.

Outcome statistics over initial population's size and mix
In this section we illustrate how the outcome statistics in Fig. 3 panels b, c, d, e and g can be reconstructed using our probability decomposition method.
At high mutation rate (Figs. 3b-c), the asymmetric phase is very short. This is reflected in the probability decomposition in Figs. 16 and 17, where the probabilities in panels f-h are all zeros, i.e., the statistics of the outcome are given exclusively by the fixation of the system into one of the two basins during the stochastically induced phase (panels a). The accuracy of the probability decomposition is very high for both small (Fig. 16) and large (Fig. 17) initial populations, as illustrated in panels j.
At low mutation rate the asymmetric regime is very long, however, such regime has little effect on the outcome statistics for small initial populations. This is explained by the probability decomposition in Fig. 18. For some starting points N 0 the system can be in the equivalence class S 2 at the end of the stochastically-induced regime (panel b), however, from this area of the simplex it is extremely unlikely for the system to cross the separatrix towards the AllD basin, even if the asymmetric regime is long (panel f ). Moreover, it very unlikely for the system to be in S 4 at the of the stochastically induced regime, since with a small initial population such regime is longer and fluctuations are likely to drive the systems towards the edges and corners (panels d). Thus, even if the system can cross the separatrix from S 4 during the asymmetric regime, the combined probability in panel h is approximately zero for every initial condition. Finally, we can see that for some initial state s the system has a small chance of being in S 3 at the end of the stochastically-induced regime (panel c). From this area the system is very likely to cross the separatrix during the asymmetric regime, leading to the small outcome probability contribution in panel g.
Since the probabilities in panel h are very low, the accuracy of the probability decomposition is very high (panel j) also in this case. The accuracy gets worse only for middle values of growth rate (around b − d = 0.05), very small mutation rates (e.g., log 10 µ < −6) and large initial populations (e.g., N 0 > 200). This case is reported in Fig. 8, which shows that the probability of AllD outcomes is underestimated for some initial starts around the mixed-strategy unstable fixed point.  Fig. 3. Panels a-d show the probability of the system being in the four equivalence classes Si at the end of the stochastically-induced regime. The probabilities of AllD outcomes given that the system is in the representative points si at time beginning of the asymmetric regime are reported in the second row. Panels e-i show how the probabilities are combined to approximate the outcome probability and panel j shows accuracy of such approximation (cf. Fig. 3g).