Analytical solution of the SIR-model for the temporal evolution of epidemics. Part A: time-independent reproduction factor

We revisit the susceptible-infectious-recovered/removed (SIR) model which is one of the simplest compartmental models. Many epidemological models are derivatives of this basic form. While an analytic solution to the SIR model is known in parametric form for the case of a time-independent infection rate, we derive an analytic solution for the more general case of a time-dependent infection rate, that is not limited to a certain range of parameter values. Our approach allows us to derive several exact analytic results characterizing all quantities, and moreover explicit, non-parametric, and accurate analytic approximants for the solution of the SIR model for time-independent infection rates. We relate all parameters of the SIR model to a measurable, usually reported quantity, namely the cumulated number of infected population and its first and second derivatives at an initial time t = 0, where data is assumed to be available. We address the question of how well the differential rate of infections is captured by the Gauss model (GM). To this end we calculate the peak height, width, and position of the bell-shaped rate analytically. We find that the SIR is captured by the GM within a range of times, which we discuss in detail. We prove that the SIR model exhibits an asymptotic behavior at large times that is different from the logistic model, while the difference between the two models still decreases with increasing reproduction factor. This part A of our work treats the original SIR model to hold at all times, while this assumption will be relaxed in part B. Relaxing this assumption allows us to formulate initial conditions incompatible with the original SIR model.


Introduction
Several recent studies [1][2][3][4] have demonstrated that the normal or Gaussian distribution function for the temporal evolution of the daily number of new cases (deaths, or alternatively infections) at time t due to the COVID-19 pandemic disease provides quantitatively correct descriptions for the monitored rates in many different countries during a single wave. If applied early enough at the beginning of a pandemic wave, within the regime of non-exponential growth, the Gauss model (GM) makes realistic and reliable predictions for the future evolution of the first wave. It has been argued that the assumption of a Gaussian time evolution is well justified by the central limit theorem of statistics [3], an agent-based model [4], a Taylor expansion [4], or as a special case of the SIR (susceptible-infectious-recovered/removed) model [2]. A motivation of the present manuscript was to provide more rigorous arguments in favor of using the GM to estimate the characteristics (peak time, height, and width) of a first epidemic wave well ahead of its climax.
To describe the time evolution of pandemic deceases many models have been developed which are based on the original deterministic compartmental SIR-model pioneered by Mc-Kendrick and Kermack [5] and its later variants such as the SEIR-model (for recent review see Estrada [6]). Important modifications to these mean-field rate equations include stochastic aspects to account properly both, for strong number fluctuations, driving the continuous phase transition of the system near the epidemic threshold, and also for spatially correlated clusters and spreading fronts caused by the disease transmission through nearest-contact infection [7,8]. It is hoped that with these important modifications the so improved SIR models ultimately provide quantitatively correct descriptions for the forecast of pandemics in many countries, even beyond the exponential growth phase of the pandemic. The inclusion of the stochastic effects mentioned above is often done by combining numerical computer simulations with the mean-field rate equations of the SIR-or SEIR-models with constant values of the infection and reaction rates. Therefore these studies will also profit from more general analytical solutions to the SIR-model equations for time-dependent infection and recovery rates which we develop here.
On our way we discovered a new analytical solution of the standard SIR-model [5,[9][10][11][12][13] without vital dynamics describing the temporal evolution of the COVID-19 pandemic disease that applies to the whole range of parameters. This solution allows for an arbitrary time dependence of the infection and recovery rates but assumes that the ratio of the two rates is independent of time. We emphasize that the assumption of a constant value for the ratio k of the recovery to the infection rate is necessary to deduce an analytical solution for the SIR model equations. There is no reason why the ratio of these two rates should be constant, as the recovery rate μ(t) is largely determined by the nature of the disease which is not expected to change strongly over time, whereas on the other hand, the infection rate a(t) is affected by many factors, not the least mitigation measures introduced by the government. Even if we ignore the latter non-pharmaceutical interventions (NPIs) it is not obvious that the infection rate varies proportional to the recovery rate. Nevertheless we are convinced that our calculations based on a constant value of the ratio k are meaningful and important, basically for the following four reasons: first, our case generalizes earlier treatments where time-independent infection and recovery rates were adopted. Secondly, as our solution allows for an arbitrary time dependence of the infection and recovery rates, despite their constant ratio, it is possible for the first time to model analytically the influence of mitigation measures on the time evolution of epidemics as done in reference [14]. Thirdly, our analytical solution can serve as a benchmark for the verification of numerical solutions of the SIR-model and their variants. Fourthly, if the intrinsic time scales of the pandemic event and the taken mitigation measures are rapid enough we may use our analytical results derived for a constant value of k to study slowly time-varying ('adiabatic') time-dependences k(t) by inserting these slow time variations into our final analytical expressions for k. Such a Wentzel-Kramers-Brillouin-approximation has proven useful in many branches of propagation and transport theory of physics including quantum mechanics.
Moreover, in contrast to earlier work we will also calculate analytically the daily differential rate of newly infected persons resulting from the SIR-model which is the key quantity to compare with the monitored data in different countries, as the size of the currently infected, and not yet recovered, compartment is usually not known. Its asymptotic behavior, peak time and peak amplitudes will all be obtained analytically and exactly.
Besides providing analytic expressions for the quantities characterizing the solution of the SIR model, we derive a simple, accurate approximant that can be used in practice, and shares all relevant features with the exact solution, as we will show. This is a significant improvement compared with approaches, where the solution of the SIR model was for example expanded into a divergent but asymptotic series [15,16], or where it had been obtained assuming inequalities that do hold only within a very limited range of SIR parameters, as we will show.
As a side-observation we find that the SIR model exhibits an asymptotic behavior at large times that is qualitatively different from the logistic model, while the difference between the two models still decreases with increasing reproduction factor. Because the SIR model is sometimes used with arbitrary initial conditions, we recall here that initial conditions must be interrelated if the SIR model is assumed to hold at all times. This is the scenario to be investigated here, for reasons to be discussed in detail. The present investigation can be extended to the many variants of the SIR model [17][18][19][20][21][22][23][24][25].

SIR-model
As the dynamics of the COVID-19 pandemic is much faster than the dynamics of births and deaths, the neglect of these demographic factors is well justified. The SIR system is the simplest of the compartmental models used for the mathematical modeling of infectious diseases. The considered population of N 1 persons is assigned to the three compartments S (susceptible), I (infectious), or R (recovered/removed). Persons from the population may progress between these compartments.

Basic equations
In a fixed population of individuals let S(t), I(t) and R(t) denote the susceptible, infected and recovered/removed fractions of persons involved in the infection at time t, so that and because S, I, R are fractions, they must all reside within the interval [0,1]. If a(t) and μ(t) denote the semipositive time-dependent infection and recovery rates, respectively, the SIRmodel is defined with the two dynamical equations [5,9,26] where the dot here and in the following denotes a derivative with respect to t. The equation for the dynamics of R(t) follows from the sum constraint (1), i.e.
where we inserted equations (2) and (3). Equation (3) can be written aṡ so that implying that equation (4) becomeṡ with the potentially time-dependent inverse reproduction factor The solution and its analytic approximant to be derived in this work holds for arbitrary a(t), as long as μ(t) and a(t) remain proportional to each other. In that case k is a constant, usually denoted as inverse basic reproduction number, k = 1/r 0 .

Initial conditions: all-time versus semi-time SIR
There are two qualitative very different cases to be considered that can be regarded as equally valid approaches. We refer to these cases as the all-time case (I) and the semi-time case (II). For the all-time case (I), treated by Kendall [9] and part A of this work, the ratio k = μ(t)/a(t) is regarded as identically constant at all times, from t = −∞ to t = ∞. This implies, that one has to use boundary conditions S(−∞) = 1, I(−∞) = 0, R(−∞) = 0 as the epidemic must not have existed at t = −∞. For the special case of k = 0, μ(t) = 0 and R(t) = 0 at all times according to equation (7), while equations (1)-(3) reduce to I(t) = 1 − S(t) and the simple logistic differential equationṠ(t) = −a(t)S(t) [1 − S(t)] determining all fractions analytically in terms of a(t) (appendix F).
For constant k > 0, equation (7) implies the important relation [9] valid for any choice of t and t . Upon inserting the stated boundary conditions at t = −∞ into equation (9), S(t) and R(t) are simply related via S(t) = e −R(t)/k or R(t) = −k ln S(t). A special case of this relationship is In view of equation (1), there remains the freedom to just specify one single initial condition at a certain t = 0, where data may have become available. One may regard the time t = 0 when the existence of the pandemic wave in the society is realized and monitoring of newly infected persons starts. Besides k there is thus a second parameter of the SIR model, which we denote by the positive ε as it usually represents a small number, but our results are valid for any ε. We define ε via the susceptible fraction of the population at time t = 0, There is no freedom for the remaining R(0) and I(0) if k is known. Inserting the above boundary conditions into equation (10), and making use of equation (1) one has involving a function in f k (ε) that is going to frequently occur in this work Here R(0) and I(0) represent the recovered/removed and infected fractions of the population at time t = 0. For small ε 1 one can write In turn, if two initial values are known at t = 0, the basic reproduction number is given by the initial conditions. We prefer to treat the inverse basic reproduction number k and ε as variables, so that R(0) and I(0) are determined by this set. Any other choice of two variables from the set k, ε, I(0), and R(0) would work equally well, as the two remaining ones are then given by equation (12).
It it important to realize that the all-time case (I) has solutions only for a limited range of k values. While there is no recovery when k = 0 due to equation (7), for k > 1 the number of infected would drop (proven below) rather than grow from the boundary value I(−∞) = 0, which is not possible as I ∈ [0, 1]. For k = 1, I(0) = f 1 (ε) = 1 − ε − e −ε is negative for any ε > 0, and simply remains at its boundary value I(−∞) = I(0) for ε = 0, and thus I(t) = 0 at all times for k = 1. Finally, because I(0) 0 must hold, not all values in k ∈ (0, 1) are compatible with the initial condition (11) at time t = 0. The requirement I(0) 0 is equivalent with k k max with This is why the remaining meaningful range of k values is k ∈ (0, k max ) for case (I). As the SIR model is assumed to hold at all times, its solution allows to calculate all fractions at times before and after t = 0. These two features make case (I) qualitatively very different from case (II). In contrast, for the semi-time case (II) treated by Kermack and McKendrick [5] and in part B of this work, the SIR model is not assumed to hold at all times, but only at times t 0, the time for which feed data may be available. In that case the model cannot be used to calculate fractions at times prior to t = 0, as this would under many circumstances lead to S > 1 or I / ∈ [0, 1] at times prior to the 'observation' time t = 0. Because the boundary conditions of case (I) must not be respected anymore in that case, one can use arbitrary initial conditions incompatible with equation (12). The model in this setup has therefore three independent parameters such as k, S(0) and I(0), while the remaining R(0) is given by equation (1). If one chooses the initial conditions to satisfy equation (12), we are back at case (I). Case (II) therefore chooses initial conditions incompatible with equation (12) such as S(0) = 1 − and I(0) = , implying R(0) = 0. In that case one has moreover the freedom to choose k > 1, as it does not automatically lead to I / ∈ [0, 1]. Case (II) is thus more flexible, but cannot be used to adequately describe the past. While case (I) can be considered as the solution to the SIR model over the whole time domain, case (II) can be used to study the future and the case k > 1 of relevance toward the end of an epidemics. In this work we are going to study case (I) and have to thus assume k ∈ (0, 1) and ε 0.

Three important remarks
First, as has been noted before [9], in the three dynamical equations (2), (3) and (6) there is migration from the S-compartment (susceptible) to the I-compartment (infected) at a rate proportional to SI, and removal from the I-compartment to the R-compartment (recovered, dead or isolated) at a rate proportional to I; there is no exit from the R-compartment and no entry into the S-compartment.
Secondly, it is necessary to start initially at time t = 0 with at least one infected person, i.e., I(0) = f k (ε) > 0 (see the second initial condition (12)), as the dynamical equation (2) implies for the initial changeİ(t = 0) = 0 if I(t = 0) = 0. Nothing can grow out of nothing; a situation very similar to kinetic plasma instabilities where seed electromagnetic fluctuations, often from spontaneous emission [27], are needed for starting the instability. Keeping exact tract of the initial conditions during the analysis will prove to be essential to avoid mathematical singularities in the analysis.
Thirdly, the first term on the right-hand side of the dynamical equation (2) denotes the newly infected population fraction, i.e., the differential rate of the total fraction J(t) of persons that have ever been infected, i.e.
where the second equality results from equation (3). It is this rateJ(t) that can be measured, is usually reported by health agencies along with the cumulative fraction J(t) = t −∞J (ξ)dξ, and that has been modeled by the Gaussian time evolution in our earlier work [3,4]. Because S(0), I(0) are usually not reported as they cannot be measured directly, we will show how to replace SIR parameters k and ε by an initial condition for the typically available J(0) andJ(0).

Results
Throughout this work we use the bar-notation q if we approximate a quantity q.

Reduced time
We are here assuming that the infection and recovery rates have the same time dependence, so that their semipositive ratio is time-independent: corresponding to a SIR-basic reproduction number r 0 > 1. We emphasize that this special case includes the standard case used by most analysis before that the infection and recovery rates are constants with respect to time. Equation (16) still allows us to take into account an arbitrary time-dependence of the infection rate a(t) which, of course, then is identical to the time dependence of the recovery rate μ(t) due to assumption (16). Then it is convenient to introduce for arbitrary time dependence of the infection rate a(t) the new dimensionless time variable τ with τ (0) = 0 via so that Using the reduced time τ , equations (2), (6), and (7) simplify to and respectively. In equation (19) in order to meet the initial condition (12) for I(τ = 0) = f k (ε). Equation (19) includes the well-known threshold [5,9] value k which determines the socalled J-shape or peak-shape [9] of the pandemic wave. Obviously, equation (19) indicates that a peaked pandemic wave cannot emerge if the ratio k > 1: as S 1 the right-hand side of equation (19) in this case is always negative, so that with positive I(τ ) the time derivative dI/dτ < 0, so that the infection rate decreases with time. An epidemic wave emerges if the ratio k < 1; the case k = 1 we already discussed in section 2.2.
The boundary conditions S(t) = 1 and R(t) = 0 for t = −∞ take over to S(τ = −∞) = 1 and R(τ = −∞) = 0 for k > 0, because a positive k implies μ(t) > 0 for all t, and thus also a(t) > 0 for all t. The solution to equation (20) is thus formally equivalent to the relationship we had between R(t) and S(t),

Analytical solution
Since the original pioneering work [5] the procedure to obtain analytical solutions is similar in later work [9,28] and also here: either (1) one expresses two of the interesting variables S(t), I(t) and R(t) in terms of the third one as done in reference [5], or (2) one expresses all three variables in terms of a suitable chosen function as done in reference [28] and here. In both cases one uses equation (1) to calculate the solution. Inserting equations (19) and (21), We emphasize that this equation fulfills the initial conditions (11) for τ = 0. For given rates a(t) and k = μ 0 /a 0 the solution of equation (22) yields the temporal evolution of S(t). Equation (22) can be written as which apart from a different notation corresponds to equations (23) and (27) of Harko et al [28]. The solution to the dimensional version of equation (24) had been expanded into a series S(t) = ∞ n=0 a n t n by Barlow and Weinstein [15], who realized that the convergence radius is rather low and suggested to look for solutions with correct asymptotic behavior, as we will do here.
To this end we follow a slightly different but equivalent approach starting out from equation (22). We introduce the quantity G(τ ) by with the initial condition G(τ = 0) = ε. According to equations (19) and (21) the function G (25) determines The dynamical equation (22) in terms of G(τ ) reads using the function f k defined in (13). Integrating over τ then provides (figure 1) where the integration constant was determined from the initial condition G(τ = 0) = ε. The integrand vanishes at x = 0 and at certain x = G ∞ to be discussed below. The whole τ ∈ [−∞, ∞] range is thus captured by equation (29) upon varying G between zero and G ∞ , and vice versa, The solution (29) generalizes the known analytical solutions in the literature [5,9,28] as it holds for arbitrary time-dependence of the infection rate a(t). The mentioned known solutions can be reproduced with equation (29) by setting τ = a 0 t on its left-hand side resulting from a constant injection rate a 0 . The parametric solution presented in the original SIR work [5] assumed G 1 to simplify the analysis. We will see below that G ∞ 1 holds only within a very narrow range of k values close to k = 1. The contributions within the regime G ∈ [0, ε] are usually not considered in numerical schemes while they can be used to verify the proper boundary conditions stated at the very beginning. Harko et al [28] did show results for t > 0 only, as some of the fractions would get negative or exceed unity at times prior to t = 0, using , according to equation (29). Approximant (86) shown in green. The vertical red lines marks G = G 0 (54), corresponding to peak time τ 0 = τ (G 0 ). The range of G is given by G ∈ [0, G ∞ ] with G ∞ provided by equation (32) or (34). The relative deviation between exact and approximative versions is plotted in figures 16(a) and (b). their analytic solution. We will present an equivalent formulation with τ (J) of the solution (29) in equation (41).
As NPIs during the pandemic wave generate a time-changing infection rate a(t), the generalized solution (29) is highly valuable to assess quantitatively the effect of the NPIs on the time evolution of the disease.

Maximum value G ∞
The solution (29) indicates that the maximum value G ∞ = G(τ = ∞) of the function G is attained when the denominator of the integrand vanishes, i.e. or which is of the form (G1) with c = 1, a 0 = −k and r = 1/k. Consequently, according to equation (G2) we obtain (figure 2) where W 0 is the principal solution of Lambert's equation z = W(z)e W(z) (discussed in appendix G) and α is given in terms of k as The argument α of the Lambert's W 0 function in equation (32) is negative for all k (figure 15(b) in appendix G). Accordingly, W 0 ∈ [−1, 0] is negative as well (figure 15(a)), while G ∞ 0.
In figure 2 we show the maximum G ∞ as a function of k. It can be seen that G ∞ has values significantly greater than unity for k 1, while it vanishes at k = 1. The value G ∞ = 1 is attained at k = 1 − 1/e ≈ 0.63. The dashed black line shows the asymptotic behavior of G ∞ ∼ 1/k at small k 1. Significant deviations from the asymptotic behavior set in at k ≈ 1/3. According to these considerations, G ∞ is well approximated by (figure 2) The relative error of this approximant is below 2.2% for k 0.63; its absolute error is below 0.014 for k 0.63, while it is exact at k = 1 and exhibits the correct asymptotic behavior for k → 0. Similarly, the non-principal solution of Lambert's equation sets the lower bound The lower bound is identical zero because Lambert's equations is trivially solved by W −1 = −1/k for α given by (33). To summarize, . This compatibility is established as a result of the proper boundary conditions.

Differential and cumulative rates of newly infected persons
From the invariantJ(t)dt = j(τ )dτ we obtain with equations (25), (27), and equation (17) in the form dτ/dt = a(t) for differential rate of newly infected persons (figure 3)  (36). The simple approximant for j(τ ), equations (36) and (86), is shown in green. The peak position, located at τ 0 and j max , is highlighted by vertical and horizontal red lines according to equations (59) and (58). The initial value at with the initial value j(τ = 0) = f k (ε)e −ε . The corresponding cumulative distribution picks up all the newly infected individuals, not only the ones that are going to occur at t > 0, i.e., one has (figure 4) In the far past, both quantities j and J emanated from zero, because G −∞ = 0. As indicated we may regard these j and J as functions of G(τ ). Equation (19) in the form dS/dτ = −IS, and dI/dτ = SI − kI imply with j = SI from equation (35) A maximum of ln j and j thus occurs when S − I = k, while (25) and (28) can be used to write S − I in terms of G  (34). Its approximant is given by equation (50).
Hence equation (38) can be written as Having introduced j and J = 1 − e −G , it is worth mentioning that the solution (29) can also be expressed in terms of J as follows.
The range of J values follows from the range of G values, and is thus given by

SIR parameters
The J thus captures all infections up to τ , including those that occurred prior τ = 0, and the usually reported This possibility to determine a parameter does not exist for case (II) where such quantity J is not available from the model equations and can only be manually entered as a new parameter. Only for case (I) it is directly related to S(τ ), as we have just shown. Having determined ε from J(0), and assuming τ = a 0 t with a constant a 0 = a(t = 0), there are two remaining parameters a 0 and k of the SIR model, that we now obtain from the measurable J(0) and its derivativesJ(0) andJ(0) at the initial t = 0. To this end we write down expressions forJ(0) andJ(0) and solve them for k and a 0 . The first relationship is obtained fromJ(t) = a 0 j(τ ) given by equation (36), The second relationship follows fromJ(0) = a 0 (d j/dτ )| τ =0 , using equations (27) and (36). As d j/dτ = (d j/dG) f k (G) we can proceed and expressJ(0) in terms of G(0) = ε, k, and a 0 , where the 2nd line can be used if ε 1. Under such conditions Upon replacing k from equation (43) in equation (45) one can write down an explicit, but lengthy, expression for a 0 in terms of J(0),J(0) andJ(0). We have thus shown how to obtain all three SIR parameters from the cumulative fraction of newly infected persons J and its derivatives at t = 0. In other words, the parameters are obtained from a quadratic fit to the reported J(t), taking into account only data in the vicinity of t = 0. Depending on the available data, other expressions presented in this work may be preferred to extract the coefficients, such as equation (68).

Final values
The knowledge of G ∞ from equation (32) and the solution (29) allow us already to derive a number of important results. We first consider the final values at τ = ∞ of (figure 2) where we have used equations (19), (21) and (25). As a consistency check we note that equation (1) is fulfilled here, i.e. S ∞ + R ∞ + I ∞ = 1. Inserting G ∞ from equation (32) we obtain with the help of equation (32) which both are determined solely by the inverse basic reproduction number k = 1/r 0 . In figure 2 we plot the fractions R ∞ and S ∞ from equation (48) as a function of r 0 > 1 along with their approximants that follow immediately follow from equations (34) and (46).
If the final values were known, one can calculate the corresponding r 0 upon inverting the relationship between S ∞ and r 0 . The result is r 0 = ln(S ∞ )/(S ∞ − 1). The canonical values R ∞ = 2/3 and S ∞ = 1/3 are thus reached for r 0 = 3 ln(3)/2 ≈ 1.65. The final susceptible fraction S ∞ decreases with increasing r 0 , starting with S ∞ = 1 at r 0 = 1, and reaching S ∞ = 0 as r 0 grows.
For these values of the differential and cumulative rates of newly infected persons after infinite time we obtain with equations (35) and (46) figure 2). Its approximant based on equation (34) is With increasing values of k the cumulative fraction of infected persons decreases from 1 at k = 0 to 0 at k = 1.

Bell-shaped differential rate. GM-like solution of the SIR
Equation (36) can also be used, even without the explicit inversion of equation (29), to derive generally valid expressions for the time τ 0 of maximum, the maximum level j max and the dimensionless width ω of the differential rate j, which correspond to the three important parameters characterizing a Gaussian differential rate, although here the width ω may not be a τ -independent constant. Assuming a constant ω, equation (51) is known as the GM for the time evolution of the daily number of new infections (or also deaths) [4].
directly from equations (38) and (39). Writing this as It is exact at k = 0 and k = 1 and has a maximum absolute error of 0.008. Inserting equation (52) in equation (36) yields an analytic expression for the maximum value where the 1st, 2nd, and 3rd line are obtained with the help of (36), (54), (53), respectively. Since W −1 (α 0 ) −2 for all k ∈ [0, 1], the j max is positive. For k = 1 the W −1 (α 0 ) = −2 and j max vanishes. For k → 0 the W −1 (α 0 ) diverges to −∞ (appendix G), but the k 2 in front causes j max to approach 1/4. Note that j max (58) is not affected by ε but is solely determined by the inverse basic reproduction number k. This is reflected by the results shown in figure 3: regardless of ε the maximum value j max remains unchanged for given k. A simple approximant for j max is provided by equation (77). We recall that the maximum j max in the measurable daily number of new infections is completely different from the maximum population fraction I max of the infected compartment. Because I(τ ) = f k (G) it achieves its maximum at I max = 1 − k + k ln k according to equation (A2).
with the very useful abbreviations and (figure 6(a)) where 1 + W 0 (α) is identical to the crossover x * from equation (A11). We can use this κ with the feature κ ∈ (0, 1) to rewrite the exact expression (32) as The peak time in real units is τ 0 /a 0 , and corresponds to the number of days between the time of maximum daily infections and the time for which initial conditions have been specified. Because G 0 (k = 1) = 0, depending on the initial condition characterized by ε, the τ 0 becomes negative for sufficiently large k (the peak time has already passed). The relative deviation between exact and approximative versions is plotted in figure 16(d).  k), where 1 − k and κ(1 − k) characterize the exponential increase and decay of the differential rate j at early and late times. By determining k from the regime of exponential growth at early times, the exponential decrease at late times is already encoded in κ(1 − k).
(b) Several k-dependent quantities entering the approximant G(τ ), such as b 1,2 and c 1,2 required in equation (87). Also shown is the exact maximum differential rate of newly infected, j max , and maximum fraction of the infected compartment, I max . Figure 3 indicates that the peak time τ 0 varies inversely with ε, confirming equation (59).

Exact behavior in the vicinity of the maximum.
With the help of equations (28) and (40) one has so that a 2nd order Taylor expansion of ln j(τ ) around its maximum can be evaluated, where ln j max is given by equation (58). Making use of e −G 0 as given by (53), equation (63) evaluated at τ = τ 0 , or equally, G = G 0 , simplifies to with ( figure 7(b)) We thus have for the Taylor expansion of ln j(τ ) around its maximum value at τ 0 to second order, ln j(τ ) = ln j T (τ ) with ω 0 given by equation (66). 3.7.4. Exact asymptotic behavior of the SIR model. Next we show that the rate j(τ ) exhibits exponential behavior at early and late times with respect to the peak time and evaluate the growth coefficients. Making use of the limiting values G(τ = −∞) = 0 and G(τ = ∞) = G ∞ in equation (40) we find the growth coefficients characterizing exponential growth and decrease at early and late times, with κ from equation (61) and τ 0 from equation (59). Note that the exponent in (68) is positive, and the exponent in (69) is negative since κ > 0 for all k ∈ (0, 1).
To calculate the important prefactors in these expressions we start from equation (B2). At small times, where g = G/G ∞ approaches zero, we use the Taylor expansion ln(1 − g) = −g to obtain for g 1 as the solution of equation (B2). Since the argument of Lambert's principal solution W 0 gets small in this limit, cf equation (G10), because j = (1 − k)G + O(G 2 ), and with the help of τ 0 (59) and A(x) (60) we obtain for the regime of early times τ τ 0 In the opposite limit, where g = G/G ∞ approaches unity, we can use the Taylor expansion for 1 − g 1 as the solution of equation (B2). Since the argument of Lambert's function W gets again small in this limit, because j = κ (1 − k) , we obtain for the regime at late times τ τ 0 We have thus shown that the initial condition S(0) enters only the peak time τ 0 and that the prefactors j early and j late characterizing the amplitudes of the asymptotic behaviors depend on k alone. Their dependence on k is shown in figure 8. While j early reaches its maximum at k → 0, where the SIR reproduces the SI model (section 3.6), j late becomes very small in this limit. As a matter of fact, the asymptotic behavior of the SI model, corresponding to a value unity in figure 8, is all that can be observed up to huge times τ ∼ k −1 (denoted as τ c in appendix F), as not only the prefactor j late becomes extremely small in the limit k → 0, but also the contribution j late /κ(1 − k) from the asymptotic exponential to the cumulative J. These expressions (68), (69), (71) and (73) exactly capture the short and long time behavior of j(τ ) shown in figure 3. The coefficients ε and k can thus also be read off from the regime of exponential growth, if a 0 is already known. It is worthwhile noticing that the initial condition S(0) does only affect the peak time, but not the prefactors j early and j late in the asymptotic behavior of j (figure 9).

Gaussian width.
While j is not a perfect Gaussian, as it is not perfectly symmetric (to be discussed below), we can calculate a width via three routes. Route (A) fits j(τ ) by the Gaussian function (51) with constant w. Route (B) determines w by the known GM relationship w = J ∞ / √ π j max using our above results for J ∞ and j max . The third route is to estimate a width from the behavior in the vicinity of the peak times. The latter leads to width ω 0 (66). The former two routes give rise to very compatible widths ω as a function of k (all widths shown in figure 7(b)).  Making use of our above results (49), (46) and (32), we can thus express a τ -independent characteristic width of j(τ ) analytically as where α 0 = 2α/e and α expressed in terms of k in (33), and where we have used G 0 from equation (54). Note, that the width is unaffected by ε as well, it is solely determined by the inverse basic reproduction number k. Taking into account all limiting behaviors our approximant for ω is given by implying, with the help of our approximant G ∞ (34) Both are compared with the exact results in figure 7. The two frequencies ω 0 and ω are actually very comparable which lead us to the GM-like approximant to be discussed next. We recall that the dimensionalJ max = a 0 j max for a constant a(t), and that the dimensional width is then w = ω/a 0 .

GM-like approximant for j(τ ).
We have by now calculated the behavior of ln j(τ ) around its maximum including peak height (denoted as j T (τ ) in equation (67)), and the asymptotic behavior of j(τ ) at small and large τ . Using a continuity requirement we can construct an approximant j(τ ) for j(τ ) that captures the behavior of j(τ ) qualitatively at all times. The continuity requirement is, in light of equations (68) and (69), where j early and j late are given by equations (71), (73), and where the two times τ early and τ late have to be determined by the continuity conditions. The condition (78), with j T (τ ) from equation (67), yields Similarly, condition (79) yields The full GM-like approximant for j(τ ) therefore reads (figure 10) with j early and j late given by equations (71)   Not under all conditions does the central Gaussian part meet the asymptotic exponential tail. Both D early and D late can become complex-valued as soon as the argument of the square root becomes negative. This is the case for D late if τ 0 is sufficiently large and negative. Under such conditions the relevant crossover times are the real parts of τ early and τ late , and the asymptotic branches are only shown up to the level that is reached by j T (τ ).
Our GM-like approximant (84) suggests that the GM is a suitable approximation for times between τ early and τ late . The ω from equation (74) and ω 0 from equation (66) are found to be very similar, which implies that J ∞ calculated for the GM can serve as an estimate for the exact J ∞ , as done in reference [4]. The difference between the GM and the SIR can also be visualized by eliminating time, and plotting j/ j max versus J/J ∞ (figure 11).
The value of τ early can be most easily determined from the differential and cumulative doubling times [29] of the GML-approximant (84), as both are constants below t early . t early coincides with the time when the monitored doubling times start to increase with time.
Similarly, τ late can be determined when the differential decay time (although this is seldomly reported) or the cumulative doubling time become constants at late times of a pandemic wave.
The cumulative as well as the final J ∞ corresponding to the GM-like j(τ ) then follow by integration. Analytic expressions are provided in appendix H.

Further reduction of solution (29)
In order to derive explicit analytical expressions G(τ ) it is necessary to calculate, at least approximately, the remaining integral on the right-hand side of equation (29). Whereas in the references [9,28] the analytical evaluation of this remaining integral is not detailed, we recall that in the pioneering work [5] this was done for small values of G less than unity by expanding the exponential function in the denominator to second order in x about x = 0, i.e. e −x 1 − x + x 2 /2. However, this approximation cannot be applied here for all values of G as figure 2 indicates that G ∞ attains values significantly greater than unity. Therefore more generally, without referring to the regime x 1, an analytic approximant τ for τ can be derived, as shown in appendix A that is valid for all k ∈ (0, 1) and all ε 1. The κ is defined in terms of k by equation (61). We note that equation (86) can alternatively be written as by equation (60). Obviously, τ (G = ε) = 0 is fulfilled for this approximant. Upon replacing G by G 0 in (86), we arrive at the expression (59) for the peak time. In figure 1 we compared the numerically evaluated integral (29) with the approximation (86) for various k. We notice the excellent agreement.
Even more important in practice, we were able to invert the relationship (86) to obtain an approximant G(τ ). This is done in appendix (B). There we use this approximant to confirm that it exhibits the correct asymptotic behaviors, that it reproduces the exact j max , and shares the peak time τ 0 with τ (G) exactly. To be specific we obtain where the k-dependent parameters b 1,2 and c 1,2 ( figure 6(b)) are specified in appendix B, τ 0 given by equation (59), and W 0 is the principal solution of Lambert's equation (appendix G).
With G(τ ) at hand, one can express all results in terms of τ , and thus also in terms of time t-for arbitrary a(t)-without preparing a parametric plot. The approximant (87) is implemented in appendix J.

Time-dependent infection rate a(t)
As we have demonstrated, J ∞ and all other final values are insensitive to the precise timedependency of a(t). It is possible to flatten the curve of daily infectionsJ max = a(t max ) j max by interventions that reduce a(t max ) at peak time t max given by τ 0 = τ (t max ). HereJ max is the maximum population fraction that is newly infected within a single day, if we choose day as the time unit in which we also express the infection rate a(t).
While J(t) = J(τ ) is invariant, the differential rates are related asJ(t) =τ j(τ ) = a(t) j(τ (t)). Inserting τ (t) for an arbitrarily chosen a(t) into the expression for j = f k (G)e −G with G = G(τ ) from equation (87), the daily number of new infections is available analytically as a function of time.
Starting from equations (36) and (37) one has G = − ln(1 − J) and S = 1 − J, so that we can write down a relationship between the differential j(τ ) and the cumulative J(τ ) number of infected persons, The same relationship follows from the derivative of equation (41) with respect to τ . Since j(τ ) = a(t)J(t), the infection rate a(t) can be extracted from the measuredJ(t) and J(t) via the right-hand side of equation (88), divided byJ(t), as long as k is time-independent.

Summary and conclusions
We demonstrated that the exact parametric solution S(t), I(t) and R(t) to the SIR model with time-dependent a(t), time-independent k and arbitrary initial conditions is determined by the function G(τ ) uniquely specified equation (29). The t-dependency of G is given by τ (t) from equation (17). With G t = G(τ (t)) as a function of t therefore at hand, and solve the SIR equations at all times, and satisfy the boundary and initial conditions exactly. The measurable differential rate of newly infectedJ and its cumulative counterpart J are then available fromJ(G t ) = a(t)S(t)I(t) and J(G t ) = 1 − S(t), according to equations (36) and (37). Alternatively, one can go over the J(τ ) route instead of G(τ ) as explained in section 3.4. While the SIR model is usually solved numerically, and the relationship between G and τ given by equation (29) can also be calculated numerically, we provide analytic approximants G(τ ) (87) as well as for the inverse τ (G) (86) that both have the correct asymptotic behaviors and may be use instead of the exact result, as they are sufficiently accurate for any practical purpose. Using this approximant (87), that works for any k and any initial condition, the usually reported daily number of new infections we have expressed explicitly as a function of time t for an arbitrarily time-dependent infection rate a(t) in section 5. The presented solution of the alltime SIR model can thus be used to replace the numerical solution of the SIR model. Having an analytical solution should be also advantageous to understand the effect of interventions reflected by a(t) much easier, and in a more transparent fashion.
The whole shape of the differential rate j depends on k only. The shape is characterized by the maximum j max (58), width ω (74), asymptotic exponents 1 − k and −κ(1 − k) as well as their prefactors j early (71) and j late (73), while the initial condition affects the peak time τ 0 (59). All these quantities have been determined analytically in this work for the all-time SIR model. In addition, we have provided approximants for all quantities characterizing the exact solution of the SIR model, to appreciate their qualitative behavior, or in case the Lambert functions are unavailable. For the special case of time-independent a(t) = a 0 , and thus τ = a 0 t, the SIR parameters k, ε, and a 0 can all be determined from initial conditions J(0),J(0), andJ(0) by equations (42), (43) and (45). This is important because the I and R fractions are usually not reported at a certain time to provide an initial condition.
Our derivation not only generalizes existing solutions, but may appear very compact compared with existing ones, while it includes them as a special case. We have relaxed the assumptions of t 0 and G 1 from previous approaches, and present the unique solution that captures consistently both the future and the past, if the SIR model is assumed to hold in the past as well. The classical assumption G 1 we have shown to fail except at inverse basic reproduction numbers close to unity.
The evolution of the usually reported differential rates assumes a bell-shaped curve characterized by a peak time τ 0 , height j max , and width ω. We have provided analytic expressions for all these parameters by equations (59), (58), and (74), respectively. They involve the two solutions W 0 and W −1 of Lambert's equation, usually available in scientific software. To appreciate the validity of our arguments and derivations, we have collected the required exact features of W 0 and W −1 in appendix G.
We have inspected the limiting case of k → 0 of the SIR model to find that the SI model, corresponding to k = 0, does not trivially emerge from the SIR model in this limit. While the asymptotic behavior of the SIR model is qualitatively different from the SI model, the quantitative difference between them still decreases with decreasing k. Insofar is the SI model compatible with the limiting case of an SIR model, but at the same time asymptotically different.
The semi-time SIR model, which can be considered as the SIR model subject to initial conditions incompatible with the all-time SIR model, and thus useful to predict only the future time-evolution, will be elaborated in part B of this work. and (b) k = 0.6 (G ∞ ≈ 1.13) (black). The asymptotic behavior is captured at small x by 1/h 1 (x) (dashed gray) and at large x by 1/h 2 (x) (dashed gray). Our approximant (A5) (yellow) appears to basically coincide with the integrand, and has the correct asymptotic behavior by construction.
with κ from equation (61). Using these asymptotes we obtain the following approximant Examples for the cases of k = 0.3 and k = 0.6 are shown in figure 12. Both partial integrals can be performed analytically as and τ (G) = τ 1 (G) + τ 2 (G). While τ 1 dominates the integral at small G, at G = ε and beyond, τ 2 is responsible for the sharp increase at large G, when G approaches G ∞ (see figure 1). With A(x) defined by equation (60), τ (G) can alternatively be written as The exact asymptotic behaviors at small and large G are where the intermediate expressions are also required, e.g., to calculate j late in equation (73). As already noted, h 1 (x) and h 2 (x) dominate the behavior of the integrand in distinct x regimes. To specify the regimes, we solve h 1 (x * ) = h 2 (x * ) to obtain with the help of (32), and where we recall that W 0 (α) ∈ [−1, 0) is negative for all k ∈ [0, 1]. For x < x * the dominating contribution is h 1 (x). While x * = 1 for k = 0 and x * ≈ 1 for k < 0.2, it depends almost linearly on k for larger k and reaches x * = 0 for k = 1. From equation (A6) we infer and G 0 given by equation (54) we can derive from equation (A11) the inequality for all k, while G 0 = x * at k = 1 follows from the mentioned properties of Lambert's functions.

Appendix B. Approximant G(τ ) for G(τ )
Starting from the high quality approximant τ (G) derived in the previous section, we are looking for an explicit expression for the inverse function, G(τ ). Unfortunately, τ (G) cannot be inverted analytically over the whole G domain, but we are going to find approximants for G separately over the two disjoint intervals G ∈ [0, G 0 ] and G ∈ [G 0 , G ∞ ] so that with g 1 , g 2 given in terms of k and ε by equations (B4), (B6) to be derived in this section, and compared with the exact G(τ ) in figure 13. We recall that G ∞ is known in terms of k by equation (32). The two partial approximants meet exactly at τ = τ 0 specified by equation (59).
To derive all quantities appearing in equation (B1) we start from τ = τ 1 + τ 2 in terms of G, provided by equations (A6) and (A7). Upon introducing g = G/G ∞ and the known g 0 = G 0 /G ∞ , using the abbreviation κ already introduced by equation (61), the governing equation for g is written involving the g-independent and semipositive constant A(ε) defined by equation (60). Equation (B2) leads to the exact transcendental equation that cannot be solved analytically. To proceed we approximate, for the two cases of g > g 0 and g < g 0 separately, corresponding to τ τ 0 and τ τ 0 , respectively, the smaller of the two logarithmic terms. In each case we make sure that the approximation is correct in all limiting cases, so that not only the asymptotic behavior of G and j will be correctly captured, but also g 1 and g 2 will coincide at τ = τ 0 .
To obtain g 1 we approximate ln(1 − g) −C 1 g in equation (B2) with a constant C 1 in such a way that the approximation is exact at both g = 0 and g = g 0 . This yields C 1 = −g −1 0 ln(1 − g 0 ). If we furthermore introduce x via g = e −x , the corresponding equation for x is of the form (G1) with c = 1, r = A(ε)/κ − (1 − k)τ , and a 0 = κ/C 1 . Its solution is thus x = r + W 0 (e −r /a 0 ), involving W 0 and not W −1 because the argument e −r /a 0 0 (figure 15). With x at hand, we have found g 1 = ln x which can be written with the help of Lambert's equation in two equivalent ways both for a different purpose. In equation (B4) we have introduced the abbreviations b 1 = κ/C 1 = −κg 0 / ln(1 − g 0 ), with g 0 = G 0 /G ∞ , so that x 1 carries the dependency on time, and is unity at τ = τ 0 with τ 0 according to (59), G 0 , κ and A explicitly given by equations (54), (61), and (60). This expression for τ 0 is just another way of writing the expression (59) for τ 0 we obtained from τ (G) in the previous section. The equivalence follows already from equation (B2). The right-hand side of equation (B4) proves that g 1 (τ 0 ) = g 0 .
For the opposite case of large g close to unity. We use the same strategy as before and now approximate ln(g) C 2 (g − 1) with a constant C 2 in such a way that the approximation is exact at both g = 1 and g = g 0 . This yields C 2 = −(1 − g 0 ) −1 ln g 0 . If we proceed in a way analogous to the case of g 1 , we obtain Here we have introduced the abbreviations b 2 = (κC 2 ) −1 = (g 0 − 1)(κ ln g 0 ) −1 , so that x 2 carries the dependency on time, is unity at τ = τ 0 , where τ 0 was already given by equation (59). The two half approximants (B4) and (B6) ( figure 13) perfectly match at G 0 by construction, i.e., by the proper choice of the approximants for ln(1 − y) and ln y. The calculations leading to the present form of the approximants had to make use of properties of Lambert's function, in particular the identity W(z ln z) = ln z from equation (G6) for any positive z. Note that the dependency of G(τ ) on ε is adsorbed by x 1 and x 2 , while all other coefficients are positive and depend on k only. The coefficients are not independent of each other, as holds. To summarize, in this section we have provided an explicit expression for the approximant G(τ ) in terms of k and ε, and also an approximant for the peak time τ 0 in terms of the same parameters. Because the assumptions were asymptotically exact, the approximant has the features as expected, and G(τ ) is compared with the exact G(τ ) in figure 13. To arrive at equation (B9), we used property (G10) of the Lambert's function. Approximants τ (G) and G(τ ) are compared with each other in figure 14.

Appendix C. Approximant S(τ ) for S(τ )
S(τ ) had been defined by equation (25). Using the identity W(z) = ln(z) − ln W(z) = ln[z/W(z)] implies aW(z) = ln[z/W(z)] a and e −aW(z) = [W(z)/z] a . We thus have, starting from the result of the previous section, where S 1 (τ ) captures the regime before, and S 2 (τ ) the regime after the peak time τ 0 . As these approximants have correct asymptotic behavior, the exact asymptotic behavior of S(τ ) follows from these expressions with equation (G10).
with G and S given by equations (B1) and (C2).
With G(τ ) from equation (B1) at hand we can calculate the approximant for j(τ ) via equation (36).
with τ 0 according to equation (59). with the approximants for I(τ ) and S(τ ) in terms of G already written down in equations (D3), (C1), and (C2). Obviously, j 1 meets j 2 at τ = τ 0 , because G 1 (τ 0 ) = G 2 (τ 0 ). Note that this is not the case if equation (35) is used, as the slopes of S 1 and S 2 are not exactly identical at τ = τ 0 . The prefactors of this approximant in the asymptotic limits of early and late times differ from the exact equations (71) and (73) only by a factor C 1 and C 2 , respectively, introduced in section B. Both factors approach unity for small ε, and start to deviate from unity only at large ε 0.2. The maximum of the approximant coincides with the exact j max , because we have determined j using equation (36), and because the definition of G 0 implies j max = (e −G 0 − k)e −G 0 , which remains unchanged, as the maximum of the approximant occurs at τ = τ 0 , where G = G 0 .

Appendix F. The case of vanishing k
For the extremal case of k = 0 (so-called SI model), which is shown here to emerge nontrivially from the SIR model in the limit k → 0, the differential equation (27) for G(τ ) with G(0) = ε is solved by since R SI = 0 at all times for k = 0. The differential rate is given by G SI (τ ) via equation (36) or via equation (35), or from dJ SI /dτ , as The same result is obtained with j SI (τ ) = I SI (τ )S SI (τ ) using the solution of the logistic equation, mentioned in section 2.2. Its maximum with amplitude j max SI = 1/4 thus occurs at τ 0 SI , the asymptotic behavior is given by lim τ→−∞ j SI (τ ) = e τ −τ 0 SI and lim τ→∞ j SI (τ ) = e −(τ −τ 0 SI ) , and the GM-like behavior in the vicinity of τ = τ 0 SI follows from the Taylor-expansion as j SI (τ ) [1 − (τ − τ 0 SI ) 2 ]/4, corresponding to ω 0 = 2. Within the SI model, every person gets infected in the course of time, J ∞ SI = 1, but nobody will ever get uninfected again (through recovery/removal).
Importantly, while some quantities like G ∞ SI = ∞, G 0 SI = ln(2), j max SI = 1/4 obviously agree with the values one obtains for the SIR in the limit k → 0, the exponential decay of j SI (F3) of the logistic SI model (k = 0) is qualitatively different from the SIR model at times after peak time, as the SIR exponent −κ(1 − k) tends to reach zero for k → 0. In fact, the deviations between SI and the k → 0 limit of the SIR get smaller with decreasing k. This can be seen by inspecting τ in the regime where it is large, and G close to G ∞ . For small k one has κ k and κ(1 − k) ≈ k, as well as G ∞ k −1 and thus, according to equation (A7), as k approaches zero. This relationship can be inverted to write G k −1 (1 − e −kτ ). Since e −kτ = 1 − kτ + O(k 2 ) for any finite τ , the G approaches τ , and j thus approaches, albeit very slowly with decreasing k, the limiting e −(τ −τ 0 ) , in agreement with the SI model (F3). To be more precise, the SIR approaches the SI in the sense while for any small but finite k, the SIR j still enters the regime j ∝ e −κ(1−k)τ e −kτ at some finite crossover τ = τ c < ∞ that can be estimated by equating this asymptotics with equation (F3). The asymptotics of the SIR at large times is thus different from the SI, but the numerical difference between the SI and SIR results still decreases with decreasing k. To estimate τ c in the limit of small k, one can safely use the mentioned approximations G ∞ k −1 and κ(1 − k) k. By equating the two differential rates (69) with (73) and (F3) the governing equation for the crossover τ c is thus e −1/k (kε) −k e −kτ c = e −τ c /(e ε − 1), or equivalently, where we have used another relationship, k(1 − k) k, also valid for small k. Further assuming k ε 1, this latter expression simplifies to The crossover τ c thus increases not only strongly with decreasing k but moreover slowly with decreasing ε, and this makes equation (F5) to hold. Another signature of qualitatively different setting between the SI and SIR at k → 0 is the relationship J ∞ = R ∞ for the SIR, which breaks down at k = 0, as R does not depend on time for k = 0, and thus R ∞ SI = 0 and I ∞ SI = 1, while the limiting behavior of the SIR model is lim k→0 R = lim k→0 kG ∞ = 1.
Irrespective these observations, for an arbitrary infection rate a(t), the daily number of new infections as a function of time t is given for the SI model by equation (F3) with τ = t 0 a(t )dt .

Appendix G. Lambert's equation. Solutions and their properties
Lambert's W function [30,31] with several applications in physics [32][33][34][35][36][37] solves the equation z = W(z)e W(z) . It can be used to calculate the solution of transcendental equations of the form e −cx = a 0 (x − r) ( G 1 ) in explicit form as The two real-valued solutions are denoted by W 0 (principal) and W −1 (non-principal) as shown by figure 15(a). These two functions are available in scientific software such as lambertw (Matlab R ) or ProductLog (Mathematica R ).
For the purpose of this manuscript we are interested in W(α) and W(α 0 ) as a function of k ∈ [0, 1], where α = −e −1/k /k (32) and α 0 = 2α/e (54). The two α's and their ranges are depicted as a function of k in figure 15 that implies G 0 G ∞ .
Another important inequality is as it implies G 0 G −∞ = 0.