The exact solution of the differential Eq. 4 is given in inverse form by (Appendix A) τ = ∫ 1 − e − ε J   d y ( 1 − y ) [ y + k ln ( 1 − y ) ] , (8) which can be integrated numerically (subject to numerical precision issues), replaced by the approximant presented in paper A (involving Lambert’s function), or semi-quantitatively captured by the simple approximant to be presented next. The solution J ( Δ ) as a function of the relative reduced time Δ = τ − τ 0 , with the reduced peak time approximated by τ 0 = 1 − k f m ( k ) [ ln J 0 1 − J 0 − ln ( e ε − 1 ) ] , (9) corresponding to J = J 0 in Eq. 8, and where f m ( k ) = 1 − k + ln k , is reasonably well captured by (Appendix C) J ( Δ ) = 1 2 [ 1 + tanh Y 1 ( Δ ) ] Θ [ Δ s ( k ) − Δ ] + { 1 − 1 − J ∞ 2 [ 1 + coth Y 2 ( Δ ) ] } Θ [ Δ − Δ s ( k ) ] (10) with the Heaviside step function Θ ( x ) = 1 ( 0 ) for x ≥ ( < ) 0 . In Eq. 10 Y 1 = 1 2 [ f m ( k ) ( Δ − Δ s ) 1 − k + ln 1 − k k ] ,     Y 2 = 1 2 [ E 0 ( k ) ( Δ − Δ s ) + ln k ( 1 − k ) κ ] , (11) with Δ s = 1 − k f m ( k ) ln ( 1 − k ) ( 1 − J 0 ) k J 0 ,     E 0 ( k ) = [ k ( 1 − k ) κ − 1 ] f m ( k ) (12) also tabulated in Table 1. We note that Δ s ( k ) is always positive. Figure 2 shows the approximation (Eq. 10) for the cumulative number as a function of the relative reduced time Δ for different values of k. For a comparison with the exact variation obtained by the numerical integration of Eq. 8 see Appendix C. The agreement is remarkably well with maximum deviations less than 30 percent. The known limiting case of k = 0 is captured exactly by the approximant (Appendix D).

For the corresponding reduced differential rate j ( Δ ) in reduced time we use the right hand side of Eq. 4 with J = J ( Δ ) from Eq. 10, cf. Figure 3. Note, that this j is not identical with the one obtained via j = d J / d τ , because J does not solve the SIR equations exactly. The peak value j max in the reduced time rate occurs when J = J 0 and is thus determined by j max = ( 1 − J 0 ) ( 1 − J 0 − k ) , also tabulated in Table 1.

As can be seen in Figure 3 the rate of new infections (Eq. 12) is strictly monoexponentially increasing j ( Δ ≪ 0 ) ≃ e Γ 1 Δ with Γ 1 ( k ) = f m ( k ) / ( 1 − k ) well before the peak time, and strictly monoexponentially decreasing well above the peak time j ( Δ ≫ 0 ) ∝ e − Γ 2 Δ with the Γ 2 = ( 1 − J ∞ ) Γ 1 / κ . These exponential rates exhibit a noteworthy property and correlation in reduced time: Γ 2 Γ 1 = 1 − J ∞ κ (13) The SIR parameter k affects several key properties of the differential and cumulative fractions of infected persons. If the maximum J ˙ ( t m ) of the measured daily number of newly infected persons has passed already, we find it most convenient to estimate k from the cumulative value J ( t m ) at this time t m . While the maximum of J ˙ ( t ) must not occur exactly at τ ( t m ) = τ 0 (Appendix F), we can still use J 0 as an approximant for the value of J ( t m ) and the relationship between J 0 and k can be inverted to read (Appendix E) k = 2 ( 1 − J 0 ) − 1 1 + ln ( 1 − J 0 ) = 1 − J 0 − J 0 2 2 + O ( J 0 3 ) (14) The dependency of k on J 0 is shown in Figure 1C. With the so-obtained value for k at hand, the infection rate a ( t m ) at peak time can be inferred from a ( t m ) = J ˙ ( t m ) / j max ( k ) . It provides a lower bound for a 0 .

A major advantage of the new analytical solutions in paper A and here is their generality in allowing for arbitrary time-dependencies of the infection rate a ( t ) . Such time-dependencies result at times greater than the observing time t = 0 from non-pharmaceutical interventions (NPIs) taken after the pandemic outbreak [20] such as case isolation in home, voluntary home quarantine, social distancing, closure of schools and universities and travel restrictions including closure of country borders, applied in different combinations and rigor [21] in many countries. These NPIs lead to a significant reduction of the initial constant infection rate a 0 at later times. It is also important to estimate the influence of a later lifting of the NPIs on the resulting increase in the case numbers in order to discriminate this increase from the onset of a second wave. Especially in the papers by Dehning et al. [17], Flaxman et al. [22] and those reviewed by Estrada [4] the influence of NPIs on the time evolution of the Covid-19 pandemics has been studied using numerical solutions of the SIR-model equations. Our analytical study presented here is superior to these results from numerical simulations as its predictions are particularly robust for the late forecast of the pandemic wave.

The corresponding daily rate J ˙ ( t ) and cumulative number J ( t ) of new infections in real time t for given time-dependent infection rates a ( t ) are J ( t ) = J ( τ ( t ) ) and J ˙ ( t ) given by Eq. 2. From a medical point of view the daily rate J ˙ ( t ) is most important as it determines also i) the fatality rate [23] d ( t ) ≃ f J ˙ ( t − t d ) with the fatality percentage f ≃ 0.005 in countries with optimal medical services and hospital capacities and the delay time of t d ≃ 7  days, ii) the daily number of new seriously sick persons [24] NSSPs = 2 f J ˙ ( t − t d ) needing access to breathing apparati, and iii) the day of maximum rush to hospitals t r = t m + t d . In countries with poor medical and hospital capacities and/or limited access to them the fatality percentage is significantly higher by a factor h which can be as large as 10.

To calculate the rate and cumulative number in real time according to Eq. 2 we adopt as time-dependent infection rate the integrable function known from shock wave physics a LD ( t ) = a 0 2 [ 1 + q − ( 1 − q ) tanh t − t a t b ] ≃ { a 0 for  t ≪ t a q a 0 for  t ≫ t a , (15) which implies τ LD ( t ) = a 0 2 [ ( 1 + q ) t − ( 1 − q ) t b ln ( cosh t − t a t b cosh t a t b ) ] ≃ { a 0 t for  t ≪ t a q a 0 t for  t ≫ t a (16) The time-dependent lockdown infection rate (Eq. 15) is characterized by four parameters: i) the initial constant infection rate a 0 at early times t ≪ t a , ii) the final constant infection rate a 1 = q a 0 at late times t ≫ t a described by the quarantine factor q = a 1 / a 0 ≤ 1 , first introduced in Refs. 21 and 24, iii) the time t a of maximum change, and iv) the time t b regularizing the sharpness of the transition. The latter is known to be about t b ≃ 7 –14 days reflecting the typical 1–2 weeks incubation delay. Consequently, the parameter q mainly affects the amplitude J ˙ shown in the left columns of Figures 5 and 6 (note that we also plotted the case of no NPIs taken (i.e., q = 1 ) for comparison). Alternatively, the transition time t b controls the rapidness of the transition in the fraction of infected persons per day and therefore the widespread.

Moreover, the initial constant infection rate a 0 characterizes the Covid-19 virus: if we adopt the German values a 0 ≃ 58  days⁻¹ and t b ≃ 11 determined below, with the remaining two parameters q and t a we can represent with the chosen functional form Eq. 15 four basic types of reductions: 1) drastic (small q ≪ 1 ) and rapid ( t a small), 2) drastic (small q ≪ 1 ) and late ( t a large), 3) mild (greater q) and rapid ( t a small), and 4) mild (greater q) and late ( t a large). The four types are exemplified in Figure 4.

In countries where the peak of the first Covid-19 wave has already passed such as e.g. Germany, Switzerland, Austria, Spain, France and Italy, we may use the monitored fatality rates and peak times to check on the validity of the SIR model with the determined free parameters. However, later monitored data are influenced by a time varying infection rate a ( t ) resulting from non-pharmaceutical interventions (NPIs) taken during the pandemic evolution. Only at the beginning of the pandemic wave it is justified to adopt a time-independent injection rate a ( t ) ≃ a 0 implying τ = a 0 t . Alternatively, also useful for other countries which still face the climax of the pandemic wave, it is possible to determine the free parameters from the monitored cases in the early phase of the pandemic wave. We illustrate our parameter estimation using the monitored data from Germany with a total population of 83 million persons ( P = 8.3 × 10 7 ).

In Germany the first two deaths were reported on March 9 so that ε = 4.8 × 10 − 6 corresponding to about 400 infected people 7 days earlier, on March 2 ( t = 0 ). The maximum rate of newly infected fraction, J ˙ max ≃ 380 / f P , occurred t m = 37  days later, consistent with a peak time of fatalities on 16 April 2020. At peak time the cumulative death number was D m = 3820 / P corresponding to J m = D m / f = 200 D m = 0.009 . This implies k ≈ 1 − J 0 ≈ 0.991 according to Eq. 14 (and not far from the value k = 0.989 to be determined from the fit shown in Figure 5). From the initial exponential increase of daily fatalities in Germany we extract Γ 1 ( k ) a 0 ≃ 0.28 , corresponding to a doubling time of ln ( 2 ) / Γ 1 a 0 ≃ 2.3  days, as we know Γ 1 ( k ) = f m ( k ) / ( 1 − k ) ≃ 0.0046 already from the above k. The quantity a m we can estimate from the measured J ˙ max , as J ˙ max = a m j max and j max ( k ) ≈ 4.2 × 10 − 5 . Using the mentioned value for J ˙ max , we obtain a m ≈ 22 / days as a lower bound for a 0 .

With these parameter values the entire following temporal evolution of the pandemic wave in Germany can be predicted as function of t b and q. To obtain all parameters consistently, we fitted the available data to our model without constraining any of the parameters (Figure 5). This yields for Germany k ≃ 0.989 , t a ≃ 21  days, q ≃ 0.15 , a 0 ≃ 58  days, and t b = 11  days. The obtained parameters allow us to calculate the dimensionless peak time τ m ≃ 1353 , the dimensionless time τ 0 ≃ 1390 , as well as J m ≃ 0.009 , J 0 ≃ 0.011 , Γ 1 ≃ 0.0056 and Γ 2 ≃ 0.0055 .

We note that the value of k = 0.989 implies for Germany that J ∞ ( 0.989 ) = 0.022 according to Figure 1, so that at the end of the first Covid-19 wave in Germany 2.2% of the population, i.e., 1.83 million persons will be infected. This number corresponds to a final fatality number of D ∞ = 9146 persons in Germany. Of course, these numbers are only valid estimates if no efficient vaccination against Covid-19 will be available.

An important consequence of the small quarantine factor q = 0.15 is the implied flat exponential decay after the peak. Because Γ 1 ≃ Γ 2 for k = 0.989 , the exponential decay is by a factor q smaller than the exponential rise prior the climax, i.e., J ˙ ( t ≫ t m ) ∝ e − Γ 2 a 0 q t = e − Γ 1 ( 1 − J ∞ ) a 0 q t κ ≃ e − t / 21.8   days (17) Equation 17 yields a decay half-live of ln ( 2 ) × 21.8  days ≃ 15  days to be compared with the initial doubling time of ≃ 2.3  days. For Germany we thus know that their lockdown was drastic and rapid: the time t a ≃ March 23 is early compared to the peak time t m ≃ April 8 resulting in a significant decrease of the infection rate with the quarantine factor q ≃ a m / a 0 = 0.15 . In Figure 5 we calculate the resulting daily new infection rate as a function of the time t for the parameters for Germany, and compare with the measured data. In the meantime, the strict lockdown interventions have been lifted in Germany: this does not effect the total numbers J ∞ and D ∞ but it should reduce the half-live decay time further.

We also performed this parameter estimation for other countries with sufficient data. For some of them data is visualized in Figure 6, parameters for the remaining countries are tabulated in Tables 2 and 3. Most importantly, with the exception of the six countries ARM, DOM, IRN, PAN, PER, SMR we found values of k > 0.9 for all other countries investigated corresponding to basic reproduction numbers R 0 = 1 / k < 1.11 . These values are significantly smaller than the estimates of R 0 ∈ [ 2.4 , 5.6 ] in the mainstream literature on Covid-19 [4, 22]. Part of these significant differences may be explained by the different definitions of R 0 .

While the inverse basic reproduction number k = 1 / R 0 = μ ( t ) / a ( t ) in the SIR-model is clearly defined as the ratio of the recovery to infection rate, there are alternative definitions of the basic reproduction number R 0 using the effective reproduction factor R ( t ) . As discussed in detail in Sect. 4 of Ref. 25 R ( t ) has to be calculated from the convolution R ( t ) = c ( t ) ∫ ​ 0 ∞ d s   W ( s ) c ( t − s ) (18) of the number of daily cases c ( t ) with the serial interval distribution W ( t ) describing the probability for the time lag between a person’s infection and the subsequent transmission of the virus to a second person. As different choices of the serial interval distribution are used in the literature this leads to differences in the calculated associated effective reproduction factors R ( t ) . As R 0 is identical to the value R ( t 0 ) at the starting time of the outbreak it is not clear in the moment that this R 0 will be identical to the 1 / k of the SIR model [26].

In this work we derived for the first time an analytical approximation for the solution for the SIR-model equations with an accuracy better than 30 percent. The explicit approximation refers to the fraction of newly infected persons per day J ˙ as a function of the relative reduced time with respect to the reduced peak time. This closed form of the analytical solution only depends on a single parameter k, the ratio of infection to recovery rates. We assume that this ratio is independent of time. As J ˙ can be directly compared with the monitored death and infection rates in different countries, we see no advantage in using the more complicated SEIR-model which currently does not allow for a closed analytical solution. An analytic solution of the SIR model with an accuracy better than 5% is available as well from our yet unpublished work where we did not consider time-dependent a ( t ) , but it has the disadvantage that it involves Lambert’s function.

For the first time in the history of SIR-research (these equations have been discovered 93 years ago!) we thus have derived an analytical solution which can be applied successfully to all accumulated data of virus diseases in the world. Being of analytic form it is superior to all existing numerical simulations and results in the literature. We also discovered for the first time how to extract the value of k from the monitored data which is highly nontrivial. We applied this new method to the data taken for the Covid-19 pandemic waves in many countries. Our work includes an estimate on the effects of non-pharmaceutical interventions in these countries. This is possible as our analytical solution holds for arbitrary but given time dependencies a ( t ) of the infection rates.

An example, on how lockdown lifting can be modeled is described in Appendix H. The situation is depicted in Figure 7. The lifting will increase a ( t ) from its present value up to a value that might be close to the initial a 0 . While the dynamics is altered, the final values remain unaffected by the dynamics, except, if the first pandemic wave is followed by a 2nd one. The values for J ∞ provided in Tables 2 and 3 provide a hint on how likely is a 2nd wave. These values correspond to the population fraction that had been infected already. While this fraction is extremely large in Peru (53%), it is still below 1 % in several of the larger countries. The tables also report the unreported number of infections per reported number (column “dark”), estimated from the number of fatalities, reported infections, and the death probability f.