The COVID-19 pandemic is raging around the world with an immense toll in terms of human, economic and social impact. Forecasting a pandemic dynamics and its spreading is therefore paramount in helping governments to make informed decisions on a number of social and economic measures, apt at curbing the pandemic and dealing with its aftermath.
While different empirical models already exist to describe the epidemic dynamics locally and globally, a coherent framework is missing. Using a powerful language and methodology borrowed from high energy physics, we study and forecast the spreading dynamics and containment across different regions of the world. This framework is the renormalisation group approach [1, 2] , which was successfully employed for epidemic dynamics in [3] . Here we will generalise the framework to take into account the dynamics in between different regions of the world. The approach is complementary to other methods summarised in [4] [5] [6] [7] [8] . As for the widely adopted choice to represent the data by fitting them to simple minded logistic functions we refer to [9] [10] [11] [12] [13] [14] . We will also provide a map between our framework and compartmental models such as the time-honoured SIR model [15] . Our epidemic renormalisation group (eRG) approach is based * g.cacciapaglia@ipnl.in2p3.fr † sannino@cp3.sdu.dk upon a simpler set of equations, which can be extended in a straightforward way to include interactions between multiple regions of the world, without the need for powerful numerical simulations.
As already noted in [3] , rather than the number of cases, it is convenient to discuss its logarithm, which is a more slowly varying function. We define through it an epidemic strength function
where I(t) is the total number of infected cases per million inhabitants in the region, and ln indicates its natural logarithm. The derivative of α with respect to time provides a new quantity that we interpret as the betafunction of an underlying microscopic model. In statistical and high energy physics, the latter governs the time (inverse energy) dependence of the interaction strength among fundamental particles. Here it regulates infectious interactions. More specifically, as the renormalisation group equations in high energy physics are expressed in terms of derivatives with respect to the energy µ, it is natural to identify the time as t/t 0 = − ln µ/µ 0 , where t 0 and µ 0 are respectively a reference time and energy scale. We choose t 0 to be one week so that time is measured in weeks, and will drop it in the following. Thus, the dictionary between the eRG equation for the epidemic strength α and the high-energy physics analog is It has been shown in [3] that α captures the essential information about the infected population within a sufficiently isolated region of the world. The pandemic beta function can be parametrised as
whose solution, for n = 1, is a familiar logistic-like function
The dynamics encoded in Eq. (3) is that of a system that flows from an UV fixed point at t = −∞ where α = 0 to an IR fixed point where α = a. The latter value encodes the total number of infected cases per million expected in the region under study. The coefficient γ is the diffusion slope, while b shifts the entire epidemic curve by a given amount of time. Further details, including what parameter influences the flattening of the curve and location of the inflection point and its properties can be found in [3] . Note that here we work with number of cases per million, so that our α corresponds to α − ln n m of [3] with n m the number of inhabitants per million per each sufficiently isolated region of the world. In this work, we extend the eRG formalism to include the diffusion of the epidemic between multiple nearly-isolated regions. Our work is organised as follows. In Sec. II we present the formalism in the simple case of two regions, and study how the interaction term influences the delay between the epidemic peak between the first and second regions. We also investigate the effect of closing the borders at different times between the two regions. In Sec. III we map the eRG formalism onto time-honoured compartmental models of the SIR-type. In Sec. IV we test the eRG formalism by comparing the predictions to data relative to the COVID-19 epidemic in Europe and in the United States, and generalise it to include multiple regions at the same time. We then use the model to simulate a second wave of epidemic diffusion in a sample of European countries. Finally, we present a discussion of the results in Sec. V and offer our conclusions.
Here we go beyond the state-of-the-art by considering the diffusion among multiple regions of the world, each characterised by its own α i (t), which in isolation obeys a beta function like Eq. (3), with its own γ i and a i . We exemplify the framework by first considering two regions, and we generalise to multiple ones later. To couple the two equations, we start from the following axiom: there is a constant number of travellers moving from one region to the other, and viceversa, given by ∆N trav each week. Our basic simplifying assumption is that the number of travellers is symmetric, i.e. there is no net flow of people between the two regions: this is a reasonable approximation during a short time as immigration only involves a smaller fraction of inhabitants than that involved in the epidemic. We further use the approximation that the rate of infected cases within the travelling subset of people is the same as the rate of infected cases in the total population of each region. Thus, the variation in the number of infected cases per million in region-1 is given by 1
where n m1 is the population of region-1 in millions and
For region-2, we find the analogous
where the same k applies. Physically, the parameter k measures the number of reciprocal travellers per week in units of million people. For instance, if the number of weekly travellers is ∆N trav = 1000, then k = 10 −3 .
Using the identity
the effect of this exchange can be encoded in the two beta functions, c.f. Eq. (3), as follows:
The above equations describe the evolution of the epidemic across the two regions, once a small fraction of the population travels between the two. However, for large k, they have the interesting property of forcing α 1 = α 2 , which in turn modifies the value of the fixed point for the two regions. The fact that the α's become equal in the long run indicates that the two regions have merged into one. One surprising finding is that the total number of infected cases across the two regions, for large k, may be reduced compared to the isolated case (see details in Appendix A). While mathematically intriguing, we do not consider this result physical, as having large k modifies the values of α i and γ i in the two regions compared to the values one would have in case of isolation. In other words it would violate our initial assumption that the two regions are nearly-isolated, with small k. One can [3] . We use data updated to the 4th of May (from www.worldometers.info), so the values should be considered as averaged over the whole period of pandemic diffusion.
go beyond the realistic case envisioned here by increasing k. This would require modifying the set of equations substantially and goes beyond the scope of this work.
To quantitatively estimate the interaction between two regions of the world, we consider benchmark values for the parameters in the two beta functions using the results given in [3] . We show in Table I the values of a and γ for various regions of the world for the COVID-19 pandemic [3] . a is normalised per million inhabitants and all values are adjourned to the 4th of May, 2020. 2 The values of γ and a are average values over the whole duration of the epidemic diffusion in each country/region. We observe that the value of γ tends to diminish over time as a consequence of the effect of gradual implementing of social distancing measures in each region. At the early stages of the epidemic, we observe γ ∼ 1, so that we will consider this as a benchmark value for the epidemic diffusion without any restriction.
With the exception of South Korea and China (Hubei province), the range for a is roughly [7.5, 8.6] , while for γ we find [0.4, 0.76]. Thus, we defined the following benchmark scenario for the two regions:
while we vary the values of γ 1 and γ 2 as specified in the figures. The value of b 2 = 200 is chosen such that the peak in the two regions in isolation have a relative delay of 14 weeks. The peak is here defined as the week where the maximum number of new infected cases per million is registered and corresponds to the inflection point of the total number of infected cases curve. The explicit formula for the inflection point time as function of the parameters of the theory can be found in [3] .
As a first sanity check, we computed the total number of infected cases across the two regions, per million, at the end of the pandemic, i.e. at infinite time. This is given by
as a function of k. The result allows us to determine the largest value of k that does not affect the total number, i.e. the largest value that k can have before the two regions effectively merge into one. In Fig. 1 we show the results for two different populations. The plot shows that k as large as 0.1 is allowed before our description of the coupled system breaks down. Note that the maximal value of k grows linearly with the population in region-2, as it enters as the ratio k/n m2 in the coupled differential equations.
To understand how the interaction encoded by k affects the diffusion of the epidemic in the two regions, we study the same benchmark of Eq. (11), except that we set b 2 = ∞, i.e. the region-2 remains with zero infected cases if isolated. One caveat that should be kept in mind is that the values for γ in Table I are obtained by fitting the data during the whole period of the epidemic, i.e. they take into account the effect of social distancing measures in each region. However, at the early stages of the epidemic, when social distancing measures were not yet being enforced, we expect larger values of γ. That is the reason why, in the following, we assume γ = 1 as the natural initial benchmark value. Nevertheless, we show how different social distancing measures impact the results for region-2 by showing also the results for smaller values of γ 2 . We discover that the interaction among the two regions of the world, controlled by the parameter k, is sufficient to ignite the spread of the epidemic to region-2 and it also controls the timing of the peak. This is shown in Fig. 2 , where we plot the time of the peaks in the two regions as a function of k/n m2 . The result does not depend on n m1 . Also, the time of the peak for region-1 is unaffected by the value of k (dashed curve) while it affects the timing of the peak for region-2 (solid curves). Note that the kterm in Eq. (10) sparks the epidemic diffusion in region-2 as soon as k nm2 e α1(t) becomes sizeable. After this point, the epidemic evolution follows the solution of the initial equation (3), as encoded in the first term of the region-2 beta function.
The numerical results for the peak delay show a linear dependence on ln k, with a change in slope appearing for k/n m2 ∼ 10 −3 . This value corresponds to k n m2 e a1 = 1 ,
i.e. it marks the threshold (grey line in the plots) between the regime where the interaction term is always smaller than one, and the one where strength 1 is attainable.
To test our approach we consider the COVID-19 epidemic spread from China (Hubei province) to Europe (Italy). From data it is known that the peaks in the two regions are about 7 weeks apart. A reasonable estimate of weekly travellers between the two regions is in the order of the thousands, so we consider k = 5 × 10 −3 as a benchmark. This means that
For this value, the bottom plot in Fig. 2 allows us to estimate the peak delay to be around 6 weeks for γ 2 = 1, i.e. for unrestricted diffusion within Italy. This nicely confirms our expectations while validating the model. It is useful to note that both k and b 2 lead to a temporal shift of the epidemic curve for region-2, however the underlying mechanisms are distinct. The former is due to an interaction between two different regions of the world, while the latter is a constant of integration that depends on the number of cases at the initial time t = 0 in region-2. This also means that a specific peak time for region-2 relative to region-1 can emerge as a combination of the two effects, interpolating between the two limiting cases: the peak delay is entirely due to the interaction with region-1, or it is due to the presence of cases in region-2 (which may have different origin) and the coupling to region-1 is negligible. We will discuss this interplay in more details in the next section.
We now turn out attention on the impact of closing the borders between two regions of the globe versus different degree of social distancing. In the eRG approach this is implemented by setting to zero k after the closing time t cl . We consider the benchmark values given in Eq. (11) , while the impact of social distancing is encoded in region-2 in varying the value for γ 2 . Furthermore we consider two scenarios: one in which region-2 has zero initial cases, meaning that the epidemic would not occur for k = 0 (corresponding to b 2 = ∞) and another where we fix the initial condition according to the benchmark (corresponding to b 2 = 200).
The results are shown in Fig. 3 , where we report the delay in the peak of region-2 caused by closing the borders (i.e., delay relative to the case of t cl = ∞). Such a delay depends crucially on the value of γ 2 in region-2 as shown in the top plot when the epidemic in region-2 is only driven by the interaction term. In particular the results show that a significant delay in the spreading of the epidemic can be achieved only if the closing is enacted before the peak in region-1 (which is unaffected by k).
In the bottom panel we show the case where region-2 features already some initial cases, so that t cl = 0 would correspond to isolated regions with both featuring infected cases. In this case, we also see that the effect of the interaction is more pronounced for small values of γ 2 in region-2, indicated by the red curve for γ 2 = 0.4. For this value of γ 2 , isolation would yield a delay of 4.5 weeks in the peak. For larger values of γ 2 (less social distancing) the peak delay is strongly reduced to within Week of closed frontiers, t cl
[weeks] one or two weeks. In any case, closing the borders is only relevant if done before the peak in region-1 is attained.
Our results, obtained using the simple and effective eRG approach, agree qualitatively with the ones presented in [16] obtained using a numerical analysis. The take home message is that social distancing plays the dominant role in curbing and delaying the epidemic spread in region-2 with respect to seed region-1.
Epidemic dynamics is often described in terms of simplistic compartmental models introduced long time ago in [15] . Here, the affected population is described in terms of compartmentalised sub-populations that have different roles in the dynamics. Then, differential equations are designed to describe the time evolution of the various sub-populations. For an application to the COVID-19 epidemic, see [17, 18] . The sub-populations can be chosen to represent (S)usceptible, (I)nfected and (R)ecovered individuals (SIR model), obeying the following differential equations:
where P = S + I + R is a constant, measuring the total number of individuals affected. As the equations do not depend on the normalisation of the number of individuals, we can consider them for cases per million. Due to the constant P , only two equations are independent, so that we can drop the one for S. The total number of infected, I(t), we study in our model is related to the above sub-populations as
We can therefore re-write the two independent SIR equations as
Eq. (19) has a form similar to Eq. (3), except for the following: it is written in terms of the total number I(t) instead of its log α(t); it contains a dependence on the number of recovered cases, R(t). Thus, our eRG approach would be equivalent to the SIR model if we could drop the R(t) dependence in the differential equation for I(t). It is conceivable that this is the case: in fact, in Eq. (19) we can already see that the second factor drives the number of infected cases to the fixed point I(∞) → P ≡ e a , which corresponds to the IR fixed point in the eRG approach. R(t), instead, is zero at early times and only grows slowly as long as the recovery rate is small, thus its effect should remain negligible once the dynamics of I(t) is driven towards the fixed point. As investigated in [3] , the dynamics of I(t) and α(t) can be described by the same equation, as they both are driven to flow between the two fixed points.
Once a solution for I(t) is found following the eRG approach, i.e. Eq. (3), the number of recovered cases can be calculated by solving Eq. (20), with solution
To validate this approach and calibrate , we compared the above formula to the number of recovered cases for the United States (US), where I(t) is obtained using the fit values in Table I : the results are shown in Fig. 4 , where R(t) (in red) reproduces the data for = 0.09. We checked that for other countries, a similarly good fit can be obtained for ∼ 0.1, thus we consider this description consistent. To establish a more quantitative dictionary between the eRG approach and the SIR model, we compared the numerical solutions of the SIR equations (19) and (20) to the solutions of the beta function in Eq. 3 (with R(t) given by Eq. (21)). We find that the solutions overlap as long as matching values of γ andγ are used. In Fig. 5 we show the numerical relation between the matching values of the couplings for 3 choices of the recovery rate : the result shows a linear relation between the couplings in the two models.
Being able to reproduce the number of recovered cases for one region in isolation, we can now address the issue of the effect of the recovered cases in the coupled system. In fact, the transmission of the epidemic due to travel of individuals between the two regions is only due to the presence of people actively infected, namely it depends on
Thus, it suffices to replace the expression in Eq. (5) with
and similarly for Eq. (7). We have compared the solutions of the coupled differential equations with and without taking into account the recovered cases, and found that including R i (t) only affects the epidemic diffusion in region-2 by a few days. Thus, this effect can be neglected in first approximation.
We now confront the eRG framework to data from the COVID-19 pandemic collected from www.worldometers.info and adjourned to the 4th of May, 2020. Although we are well aware of the pitfalls stemming from comparing data provided by different countries due to the inhomogeneous way infectious cases were tested and reported, it is still possible to extract from these reliable time behaviour and structure. Of course, when coupling two regions of the world, part of the initial uncertainty also affects the epidemic transmission probability without affecting the overall picture.
Nevertheless, we will now see that the eRG formalism can be used simultaneously to quantiately project the spreading dynamics across different regions of the world, or as an a-posteriori way to learn how this spreading came to be. We focus on two examples, one intra European (Italy-Denmark) and the other between Europe and the US. The values for γ and a for each country are taken from the fit in Table I , which assumes isolation.
In the right plot of Fig. 6 we show the total number of infected cases in Italy and Denmark (blue and red dots, respectively), compared to the fit in Table I (solid  curves) : the latter assumes that the epidemic occurred in the two regions while in isolation. We now wish to understand how and whether the virus spread from Italy to Denmark: thus, we used the coupled Eqs (9) and (10), while we set the number of initial cases in Denmark to be null, i.e. b 2 = ∞. All the other parameters are fixed to the values in the Table. We find that the two curves can be reasonably fit by assuming k = 0.16, as shown in the left plot of Fig. 6 : the dashed orange curve, corresponding to Italy, overlaps to the isolated fit (solid blue), while the new curve for Denmark (dashed green) is close to the isolated fit (solid red). Let us now comment on the actual value of k. If we take it literally this it would correspond to a rate of 160.000 travellers between the two regions each week. This is an unreasonably large value but it can be alternatively and conservatively interpreted in the following ways:
i) More countries contributed to the epidemic spread in Denmark;
IT and DK k=0.16 i) The original spreading dynamics in Denmark is due to few very socially active infected individuals that traveled back from Italy and/or were superspreaders;
iii) A combination of the above.
Whatever the reason, it is naturally incorporated in a larger value of k. One can also take into account the various scenarios by effectively re-instating an initial value for α 2 at t = 0 while reducing the k value.
To further test our model we consider the system consisting of Europe as region-1 and the United States as region-2. For simplicity, we modelled Europe on the European Union (with n m1 = 445) with parameters from the fit of the epidemic diffusion in Italy (c.f., Table I ). After setting to zero the initial cases in the US, we were able to reproduce the diffusion of the epidemic in region-2 (US) for k = 10, as shown in the right plot of Fig. 6 . While it is still possible that the large value for k may be interpreted as in the above case, it has the further effect of distorting the epidemic curve for region-1, the EU, thus suggesting that it may be hard to explain the diffusion of the COVID-19 epidemic in the US as originating solely from the EU.
At this stage, we cannot exclude that adjusting the epidemic parameters in the EU could improve the agreement. This exercise, nevertheless, proves the effectiveness of our simple eRG model to describe the diffusion of the epidemic among different regions of the world. A more accurate fit may be obtained if more than one region is included in the analysis, which will be considered in a future work.
We now use the eRG framework to model the impact of a new wave of epidemic spread of the COVID-19 virus (or a related one) in Europe. To do so, we simulate the effect of transmission among countries in a pool of European countries, namely Italy, Spain, France, the United Kingdom, Germany, Denmark and Switzerland. We also include an unspecified "seed region", with a population of n m0 = 50, which has some initial cases, while no case is initially present for the simulated European countries. This is achieved by setting b i = ∞, where i = 1, . . . 7 spans over the 7 sample countries mentioned above.
We generate randomly the diffusion factors γ i in the range [0.4, 0.76], based on the data of the current COVID-19 epidemic in Europe, and also generate random values of a i in the range [7.5, 8.5] . This also includes the seed region. Finally, we provide randomly generated numbers of travellers between each of the regions, including the seed one, giving coupling values k ij in the range [1, 10] × 10 −3 . We then solve the 8 coupled differential equations:
where i, j = 0, . . . 7 and α 0 corresponds to the seed region. The result is shown in the top row of Fig. 7 , where black indicates the seed region and the coloured curves correspond to the 7 sample European countries. The top-right plot, where the distribution of new cases is displayed, clearly shows that the peaks in the infected regions occur between 3 to 12 weeks after the peak in the seed region. This effect, however, is mainly due to the values of the γ's in those regions, and not on the values of the interaction couplings k ij .
To prove this, we have run the same simulation again, by fixing γ i = 1, i = 1, . . . roughly correspond to unrestricted diffusion of the virus in the target regions. The result shows that all the peaks are now occurring within 4 weeks after the peak in the seed region.
The results nicely demonstrate that our eRG framework not only is useful, simple and effective to understand the current pandemic, but can also be used to model future ones.
We extended the epidemic renormalisation group approach to analyse the dynamics of disease transmission and spreading across different regions of the world. We have shown that the eRG framework constitutes an effective way to understand the relative impact of border control versus social distancing measures on the global spread of the epidemic. The simplicity of the approach, stemming from an effective description of complex phenomena, make it a reliable alternative to the use of expensive high-performance numerical computations.
We calibrated our approach via internationally reported cases. The approach elucidates the underlying mechanism that governs the delay in the relative peaks of newly infected cases across different regions of the world. Among our results, we were able to demonstrate that social distancing measures are more efficient than border control in delaying the epidemic peak.
In order to connect with widely used time-honoured compartmental models of the SIR-like type, we established the proper map with our eRG framework. We have also shown how to generalise the eRG framework to account for the epidemic interactions across multiple regions of the world.
We foresee a number of future applications and extensions of our seed work. From a more phenomenological point of view, of immediate impact for society, we plan on embarking on a world-wide monitoring to make global projections that will help governments and industries make containment plans and strategise about reopening society and how to best implement border control. We also wish to improve on understanding the link between the eRG approach and microscopic models of population dynamics and epidemic spread including a number of granular effects that are, by construction, averaged over by effective descriptions such as the eRG approach.
Appendix A: Asymptotic behaviour for large k
The coupled system of beta functions in Eqs (9) and (10) has an interesting solution in the limit for large k. While this limit should be considered unphysical, it is intriguing from a mathematical point of view. Moreover, it may lead to some insights on the dynamics of two regions merging into a single one.
In the limit of large k → ∞, the interaction term in the two beta functions dominates. It thus forces α 1 = α 2 asymptotically. To check if there is a fixed point, we can take the sum of the two functions, eliminating the k-term, and search for zeros:
n m1 e α1 β(α 1 ) + n m2 e α2 β(α 2 ) = e α1 n m1 γ 1 α 1 1 − α 1 a 1 + e α2 n m2 γ 2 α 2 1 − α 2 a 2 = 0 .
Imposing the condition α 1 = α 2 = α, we find: e α α n m1 γ 1 1 − α a 1 + n m2 γ 2 1 − α a 2 = 0(A2) that is solved by α = 0 (UV fixed point) or α * = (n m1 γ 1 + n m2 γ 2 )a 1 a 2 n m1 γ 1 a 2 + n m2 γ 2 a 1 ,
which defines the new IR fixed point. It is now interesting to ask: for fixed a 1 and a 2 , is the total number of cases in the two regions in the k → ∞ limit larger or smaller than that in the k → 0 limit? In fact, I tot k→∞ = (n m1 + n m2 )e α * = (n m1 + n m2 )e 
Similarly, I tot k→0 = n m1 e a1 + n m2 a a2 = e A (n m1 e δ + n m2 e −δ ) .
Thus I tot k→∞ I tot k→0 = n m1 + n m2 n m1 e δ + n m2 e −δ e − δ(δ+gA) A+gδ .
(A7)
Interestingly, as the exponential in the numerator may be negative, it is possible to have a reduction in the total number of cases when the exchange of people is large. In Fig. 8 we show numerically the above ratio, for n m2 = n m1 (top plot) and n m2 = 10 n m1 (bottom plot), and for various values of g (we recall that −1 < g < 1). To make sense of the result, let's consider the case of similar slope in the two countries, i.e. γ 1 = γ 2 . In the case n m1 = n m2 , this would give g = 0, i.e. the blue curve in the top plot: this case features a reduction of the total cases, as long as a 1 = a 2 . For n m2 = 10 n m2 , we have g ≈ 1, corresponding to the orange curve in the bottom plot: again, there is a significant reduction for a 1 > a 2 , k→∞ /I tot k→0 as a function of delta for g = 0 (blue), g = 0.5 (red), g = −0.5 (purple), and the extreme values g = 1 (orange) and g = −1 (green). We fixed A = 12, but the results have a mild dependence on its value.
while the increase if minor if a 2 > a 1 . Of course, these results are just illustrative, because the values of a i and γ i should be influenced by the policies concerning the movement of people between infected regions.