Multi-patch epidemic models have been used to study various infectious diseases, for example, nosocomial infection [20] , vector-borne diseases [19] , HIV/AIDS transmission [18] , SARS epidemic [22] , and so on. They are often used to capture the heterogeneity between different geographic locations, for example, a densely populated city and a less populated rural area. It also helps to study the effect of migrations or lock-down measures among different population groups or locations. In the Covid-19 pandemic, it has been observed that the infectivity in different regions may vary and is impacted by various social-distance and lock-down measures [27] .
ODE models are often used to study the dynamics of such multi-patch epidemic models. It is well known that the ODE dynamics arises from the Markovian assumptions in the stochastic multi-patch epidemic model, that is, the infection process is Poisson, the infectious (and/or exposed/latent) periods are exponentially distributed and the migration processes are also Markovian [1, 8, 21, 22, 4, 19] . Some ODE/PDE models are also used to study their dynamics when the infection rates are age-dependent (depending on how long the population has been infected), see, for example, [20, 11] . These models also assume exponentially distributed infectious periods and Markovian migration processes.
In this paper, we study multi-patch SIR and SEIR models, in which the infectious (and exposed/latent) periods have a general distribution, while the migration processes are Markovian. The infection is assumed to be local, that is, the infection rate depends on the susceptible and infectious population in each patch only. Individuals may migrate from one patch to another in each of the Susceptible, Exposed(Latent), Infected and Recovered stages. The epidemic dynamics evolves in a manner of mixed Markovian and non-Markovian components, which seems challenging to directly derive its first-order fluid equations.
We describe the evolution dynamics by tracking the time epochs of becoming exposed and/or infectious and the location of an individual at these event times. Specifically, in the multi-patch SIR model, each individual tracks the time epochs of becoming infected and recovered, and is associated with a continuous time Markov chain taking values of the patches/locations. In the multi-patch SEIR model, each individual tracks the time epochs of becoming exposed, infected and recovered, and is associated with two Markov chains that are used to track their movement starting when the individual becomes exposed / infectious. For the initially exposed and/or infected individuals, we also assume that their remaining exposing and/or infectious periods have general distributions, which may be different from those of the newly exposed/infected individuals. For these initially exposed/infectious individuals, we also track their movement while being exposed/infectious. Given the representations with these time epochs and location processes, we show a functional law of large numbers (FLLN) and a functional central limit theorem (FCLT) for the associated counting processes in the multi-patch SIR and SEIR models. The fluid limits are determined by a set of Volterra integral equations. When the infectious (and exposed/latent) periods are deterministic, we can write the fluid integral equations as a set of ODEs with delay (Remarks 2.1 and 3.1). The limit processes in the FCLT are determined by a set of stochastic Volterra integral equations, driven by a sum of independent Brownian motions and continuous Gaussian processes with a certain covariance structure. When the infectious (and exposed/latent) periods are deterministic, the limits become stochastic differential equations with linear drifts and delay (Remarks 2.3 and 3.2). We also remark how the approach and results can be extended to study multi-patch SIS and SIRS models (Remarks 2.4 and 3.3).
In the proofs of these results, we employ Poisson random measures (PRMs) that are constructed as the sums of the Dirac masses at the time epochs of becoming infectious (and/or exposed), the infectious and exposing periods and the Markov chain starting from the location of each individual at those epochs. We can then use the martingale properties and convergence theorems as critical tools in the proofs. For a single patch SIR and SEIR models with general infectious and exposing periods, an approach using PRMs that are constructed at the time epochs of becoming infectious (and/or exposed), was developed in [24] . The approach is further developed in this paper for multi-patch SIR and SEIR models, to track the locations of each individual at each event epoch.
This paper contributes to the limited literature on stochastic epidemic models with general infectious periods, see the overview in Chapter 3.4 of [8] on the common approaches to study non-Markovian epidemic models and LLN and CLT for the final sizes of the epidemic; see also the recent method using piecewise Markov deterministic processes in [9] and [15] for the SIR model. We also refer to Reinert [28] that proves an FLLN for the empirical measure of the SIR epidemic dynamics using Stein's method. In [24] , both FLLN and FCLT were proved for the epidemic dynamics in the classical models (SIR, SIS, SEIR, SIRS) where the PRM representations of the dynamics plays a fundamental role in the proofs. Although Volterra integral equations were used to describe the proportion of infectious population in the SIS, SIR or SEIR model without proving an FLLN (see [6, 7, 10, 12, 17, 29] ), no Volterra integral equations have been proposed for multi-patch epidemic models with general infectious (and/or exposing) periods. Our work shows both FLLN and FCLT for non-Markovian multi-patch models, and identify (stochastic) multidimensional Volterra integral equations as their limits.
It is also worth mentioning the multi-type epidemic models where the population splits up into multiple groups of individuals and each group may infect any other group in addition to itself (no migration, and different from the local infection assumption in our model) and proportionate mixing taking into account control measures like social distance or lockdowns is also incorporated, see Chapters 6.1 and 6.2 in [1] and [2, 3] . This is very different from our model since infection is local in our setup and individuals in any of the stages may migrate from one group to another.
We also highlight that Fodor et al. [13] argue that Volterra integral equations are more effective than the ODE models in the study of the Covid-19 pandemic. In particular, they show that the ODE models can significantly underestimate the initial basic reproduction number R 0 , while the integral equations results in more accurate epidemic dynamics after sharp changes in R 0 due to lockdown measures. The integral equations are also used to estimate the state of the Covid-19 epidemic [14] .
1.1. Notation. Throughout the paper, N denotes the set of natural numbers, and R k (R k + ) denotes the space of k-dimensional vectors with real (nonnegative) coordinates, with R(R + ) for k = 1. For x, y ∈ R, denote x ∧ y = min{x, y} and x ∨ y = max{x, y}. Let D = D([0, T ], R) denote the space of R-valued càdlàg functions defined on [0, T ]. Throughout the paper, convergence in D means convergence in the Skorohod J 1 topology, see chapter 3 of [5] . Also, D k stands for the k-fold product equipped with the product topology. Let C be the subset of D consisting of continuous functions. Let C 1 consist of all differentiable functions whose derivative is continuous. For any function x ∈ D, we use x T = sup t∈[0,T ] |x(t)|. For two functions x, y ∈ D, we use x • y(t) = x(y(t)) denote their composition. All random variables and processes are defined in a common complete probability space (Ω, F, P). The notation ⇒ means convergence in distribution. We use 1 {·} for indicator function. We use small-o notation for real-valued functions f and non-zero g: f (x) = o(g(x)) if lim sup x→∞ |f (x)/g(x)| = 0.
2. The multi-patch SIR model with general infectious periods 2.1. Model description. We consider a multi-patch epidemic model, where individuals in each patch experience the Susceptible-Infectious-Recovered (SIR) process. The patches may refer to populations in different locations, for example, a densely populated city and a less populated rural area. Individuals in each patch are infected locally, while they may migrate/travel from one to the other in any of the susceptible, infectious and recovered states. The rate of infection is different in the patches (because of the differences in the density of population or the use of public transportations), while the law of the infectious period is the same (same illness).
Let N be the total population size and L be the number of patches. For each i = 1, . . . , L, let S N i (t), I N i (t) and R N i (t) denote the numbers of individuals in patch i that are susceptible, infectious and recovered at time t, respectively. Then we have the balance equation:
Assume that S N i (0) > 0, L i=1 I N i (0) > 0 and R N i (0) = 0, i = 1, . . . , L. Let λ i be the infection rate of patch i, i = 1, . . . , L. Define the following processes
, i = 1, . . . , L .
Let A N i (t) be the cumulative counting process of individuals that become infectious during (0, t]. Then we can give a representation of the process A N i (t) via the standard Poisson random measure Q i,inf :
Equivalently, we can write
where P A,i is a unit-rate Poisson process, and independent from each other for i = 1, . . . , L. We let Susceptible (resp. infectious, resp. removed) individuals migrate from patch i to patch j at rate ν S,i,j (resp. at rate ν I,i,j , resp. at rate ν R,i,j )). Let X(t) denote the location (i.e., the patch) at time t ≥ 0 of an infected individual. X(t) is a Markov process which alternates between states 1, . . . , L. Define p i,j (t) = P(X(r + t) = j|X(r) = i) for i, j = 1, . . . , L and r, t ≥ 0. Then these probabilities can be computed, see, for example, Corollary 2.2 in Chapter 7 of [26] .
Example 2.1. Suppose that there are two patches. While in state 1 (resp. 2) it stays there for a duration of time which follows the exponential distribution with parameter ν I (resp. aν I ). We define p 1,1 (t) = P(X(r + t) = 1|X(r) = 1) , p 1,2 (t) = P(X(r + t) = 2|X(r) = 1) , p 2,1 (t) = P(X(r + t) = 1|X(r) = 1) , p 2,2 (t) = P(X(r + t) = 2|X(r) = 2) .
A standard computation (based upon, e.g., Corollary 2.2 in Chapter 7 of [26] ) yields the following explicit formulas:
We will use X 0,k i and X k i to indicate the associated process for individual k in patch i, for the initially and newly infected ones, respectively. Note that they are all mutually independent and have the same law as described for the process X(t) above. Also note that the processes X k i start from the time becoming infected τ N k,i while the processes X 0,k i start from time 0. We now provide a representation of the epidemic evolution dynamics:
where P S,i,ℓ , P R,i,ℓ , i, ℓ = 1, . . . , L, are all unit-rate Poisson processes, mutually independent, and also independent of P A,i . Here, the first term in I N i (t) represents the number of initially infected individuals from patch ℓ = 1, . . . , L that remain infected and are in patch i at time t, and the second term represents the number of newly infected individuals from patch ℓ = 1, . . . , L that remain infected and are in patch i at time t. The first term in R N i (t) represents the number of initially infected individuals from patch ℓ = 1, . . . , L that have recovered by time t and were in patch i at the time of recovery, and the second term represents the number of newly infected individuals from patch ℓ = 1, . . . , L that have recovered by time t, and were in patch i at the time of recovery.
Define a PRM Q ℓ,inf (ds, du, dv, dθ) on R 3 + × {1, . . . , L}, which is the sum of the Dirac masses at the points (τ N j,ℓ , U N j,ℓ , η j,i , X j ℓ (η j,ℓ )) with mean measure ds
, and an infection occurs at time τ N j,ℓ in case U N j,ℓ ≤ λ ℓ Φ N ℓ (τ N j,ℓ ). We can then write for ℓ, ℓ ′ = 1, . . . , L,
We denote the corresponding compensated PRM Q ℓ,inf (ds, du, dv, dθ 
Define the PRMQ ℓ,inf (ds, du, dθ) on R 2 + × {1, . . . , L}, with mean measure ds du µ ℓ (dθ), where µ ℓ ({ℓ ′ }) = p ℓ,ℓ ′ (t o ) for ℓ, ℓ ′ = 1, . . . , L Then we can represent
FLLN. For any process Z N , letZ N := N −1 Z N be the fluid-scaled process. We make the following assumptions on the c.d.f. F .
for a finite or countable number of positive numbers a i and the corresponding t i such that i a i ≤ 1 and t 0 < t 1 < · · · < t k < . . . , and F 2 is Hölder continuous with exponent 1 2 + θ for some θ > 0, that is, F 2 (t + δ) − F 2 (t) ≤ cδ 1/2+θ for some c > 0. In addition, assume F 0 is continuous. 8) in probability, locally uniformly on [0, T ], where (S i (t),Ī i (t),R i (t), i = 1, . . . , L) ∈ C 3L is the unique solution to the following set of deterministic integral equations:
.
(2.12) 
Remark 2.2. Let the infectious periods be exponentially distributed with rates γ, η j,ℓ ∼ exp(γ) for all i = 1, . . . , L and j ∈ N. Due to the memoryless property, the remaining infectious periods of the initially infected individuals have the same exponential distribution. It is well known that the epidemic evolution dynamics can be described by the following set of ODEs:
It can be checked that the set of integral equations in Theorem 2.1 reduces to the ODEs in this special case.
FCLT. For any process Z N , letẐ N := √ N (Z N −Z) be the diffusion-scaled process wherē Z N is the fluid-scaled process andZ is its limit. 
where the limits are the unique solution to the following set of stochastic Volterra integral equations driven by Gaussian processes:
(2.14)
with B A,i , B S,i,j , B I,i,j , B R,i,j being mutually independent standard Brownian motions, and with the deterministic functionsS i ,Ī i ,R i being the limits in Theorem 2.1. The processesÎ 0 i,j andÎ i,j are continuous Gaussian processes with mean zero and covariance functions:
otherwise.
In addition,Î 0 i,j andÎ i,j are independent, and also independent of the Brownian terms. Remark 2.3. As discussed in Remark 2.1, when the infectious periods are deterministic and equal to t o , the stochastic integral equations become linear stochastic differential equations with delay. In particular, the term
Remark 2.4. The analysis can be easily extended to the multi-patch SIS model, where the population in each patch has susceptible and infectious groups, and when infectious individuals recover, they become susceptible immediately. The epidemic evolution dynamics is described as
Thus, in the FLLN, we obtain the same limitĪ i in (2.10) as in the multi-patch SIR model, and the limitS i (t):
Similarly in the FCLT, we obtain the same limitÎ i as in (2.15) for the multi-patch SIR model, and the limitŜ i (t):
3. The multi-patch SEIR model with general exposing and infectious periods 3.1. Model description. In the multi-patch SEIR model, individuals in each patch experience the Susceptible-Exposing (latent)-Infectious-Recovered (SEIR) process, and they may migrate from one patch to another in either of the four stages. As in the SIR model, the infection process is local, the rates of infection vary among patches, and the exposing and recovery rates are the same among different patches. Let N be the total population size and L be the number of patches. For each patch i, let S N i (t), E N i (t), I N i (t), R N i (t) count the numbers of individuals that are susceptible, exposed (latent), infectious and recovered in patch i at time t, respectively. We have the balance equation:
be the cumulative counting process of individuals that become infectious during (0, t], which can be represented as (2.1) and (2.2) in the SIR model with infection rate λ i and
The E N i (0) initially exposed individuals experience the exposing and infectious periods before recovery. Let {ζ 0 k,i : k = 1, . . . , E N i } be the remaining exposing/latent periods of the initially exposed individuals in patch i. After the exposing period, let {η −k,i : k = 1, . . . , E N i } be their infectious times. The I N i (0) initially infected individuals experience a remaining infectious period before recovery. We have the same notation {η 0 k,i } as in the SIR model for them. The A n i (t) newly infected individuals experience the exposing and infectious periods. Let {ζ j,i : j ∈ N} and {η j,i : j ∈ N} be the associated exposing and infectious periods.
Assume that {ζ 0 k,i }, {ζ j,i }, {η 0 k,i } and {η j,i } are all i.i.d. sequences of random variables having c.d.f.'s G 0 , G, F 0 , F , respectively, and they are also mutually independent. Note that η j,i is defined for j ∈ Z and 1 ≤ i ≤ L (those with j < 0 code the infectious periods of the initially exposed individuals, while those with j > 0 code the infectious periods of the newly exposed individuals). Let G c 0 = 1 − G 0 , G c = 1 − G, F c 0 = 1 − F 0 and F c = 1 − F . Individuals may migrate from patch i to j in any of the four epidemic stages, with rates ν S,i,j , ν E,i,j , ν I,i,j and ν R,i,j for the susceptible, exposing, infectious and recovered ones, respectively. Times between migration in each of the stages are exponentially distributed. Let Y (t) be the Markov process denoting the location of an individual k that is exposed from patch i at time t, and let q i,j (t) = P(Y (r + t) = j|Y (r) = i) for i, j = 1, . . . , L and r, t > 0. We use Y 0,k i (t) and Y j i (t) for the initially and newly exposed individuals in patch i, and note that they have the same distribution as Y (t). These processes are only used before an individual moves from the exposing stage to the infectious stage. We use the same Markov processes X 0,k i (t) and X j i (t) for the migration of infectious patients as in the SIR model.
The multi-patch SEIR epidemic evolution dynamics can be described as follows:
and
where P S,i,j , P R,i,j , i, j = 1, . . . , L, are all unit-rate Poisson processes, mutually independent, and also independent of P A,i . It is clear that the dynamics of E N i (t) is the same as that of I N i (t) in the SIR model. We remark that in the expression for the I N i (t) in the SEIR model, the first term counts the number of initially exposed individuals from all the patches that remain exposed and are in patch i at time t, and the second term counts the numbers of initially exposed individuals from all the patches that have become infectious in patch i at time t (for tracking purposes, the location at the epochs of becoming infection is recorded, since we use η k,ℓ and X 0,k ℓ for different patches ℓ). Also note that we use the Markov process X 0,k ℓ to indicate that these are for the initially exposed individuals. The third term counts the number of newly exposed individuals at all patches that have become infectious and in patch i at time t, and we also track the patch in which each individual becomes infectious.
Note that some of the key components in the dynamics above can be represented via PRMs. The infection process A N ℓ has the same representation in (2.1) using the PRM Q ℓ,inf . Define a PRMQ ℓ,inf (ds, du, dv, dϑ) on R 3 + × {1, . . . , L}, which is the sum of the Dirac masses at the points
. We can then write for ℓ, ℓ ′ = 1, . . . , L,
. . , L}, which is the sum of the Dirac masses at the points
, and for each z > 0, µ X ℓ (z, {ℓ ′ }) = p ℓ,ℓ ′ (z), and again an infection occurs at time τ N j,ℓ if and only if
. We can then write for ℓ, i = 1, . . . , L,
We denote the corresponding compensated PRM Q ℓ,inf (ds, du, dy, dϑ, dz, dθ) = Q ℓ,inf (ds, du, dy, dϑ,
Assumption 3.1. The c.d.f.'s G and F satisfy the conditions of F in Assumption 2.1. In addition, assume that G 0 and F 0 are continuous. 
is the unique solution to the following set of deterministic integral equations:
and
(3.14)
Remark 3.1. Suppose the exposing and infectious periods are deterministic, taking values of t e > 0 and t o > 0. Also, assume that the remaining exposing and infectious periods of the initially exposed and infectious are uniformly distributed over (0, t e ) and (0, t o ), respectively. Then the fluid equations
It is easy to see that we obtain a set of ODEs with delay after taking derivative. 
where the limits are the unique solution to the following set of stochastic Volterra integral equations driven by Gaussian processes:
HereΦ 
are continuous Gaussian processes, independent of the Brownian motions above, with mean zero and covariance functions:
are independent from each other, and
Then the limits in Theorem 3.2 become stochastic differential equations with linear drifts and delay. In particular,
Remark 3.3. The analysis for the multi-patch SEIR model can be easily extended to multi-patch SIRS model, where in each patch, the population is grouped into susceptible, infectious, and recovered individuals and individuals become susceptible after experiencing a recovery period. In this model, the infectious and recovered processes I N , R N correspond to the exposed and infectious processes E N , I N in the SEIR model. In the description of the epidemic dynamics, we need to change the dynamics of S N i in (3.1) by adding the individuals that have become susceptible after recovery, i.e., the first three terms in R N i in (3.4) . This is similar to the susceptible process S N i in (2.18) for the SIS model. Then it is straightforward to write down the limit processes in the FLLN and FCLT for the processes (S N , I N , R N ) (corresponding to (S N , E N , I N ) in the SEIR model).
In this section we prove Theorem 2.1. We first provide a representation of the processes
Then for each i = 1, . . . , L,
It has the predictable quadratic variation: 
are stochastically bounded in D by Lemma 5.8 in [25] , and thus,M 
Proof. We can rewrite the processesS N i as
The processesM N S,ℓ,i are square integrable martingales with respect to the filtration F N S (t) : t ≥ 0}, defined by
They have the predictable quadratic variation:
Thus, the processes
Under Assumption 2.2, by Lemma 4.1 and Lemma 8.1 (with y j ≡ 0 in the map Υ), applying the continuous mapping theorem, we conclude the convergence of {(S N 1 , . . . ,S N L )}. We next study the processes I N i (t), for which we give the following representation. Lemma 4.3.
where P I,i,j , i, j = 1, . . . , L, are all unit-rate Poisson processes, mutually independent, and also independent of P A,i , P S,i,j and P R,i,j .
Proof. In the representation of I N i (t), we observe that
and for ℓ = i,
is the number of initially infected individuals from patch i that are in patch ℓ at the time t ∧ η 0 k,ℓ for k = 1, . . . , I N ℓ (0).We also observe that
and for ℓ = i,
Thus, using the above identities, we obtain the expression in (4.5).
For convenience, we let
for ℓ, i = 1, . . . , L and t ≥ 0. We next treat these processes in the following lemmas. 
Proof. We first focus onĪ N,0 1,1 (t). Note that, since η 0 k,1 and X 0,k 1 are independent,
where the expectation is taken under the condition that X 0,k 1 (0) = 1. Note that the pairs (η 0 k,1 , X 0,k 1 (·)) are independent over k, and have the same distributions. Thus, by the LLN of i.i.d. variables, we obtain that for each t ≥ 0,
Note that the convergence holds in fact in probability, since the limit is deterministic. This can be extended to convergence of finite dimensional distributions. In order to establish tightness in D, we will show that lim sup
The fact that this implies tightness in D follows from Corollary page 83 of [5] .By the independence of the pairs
The convergence of the other processesĪ N,0 ℓ,i follows similarly. It remains to prove the joint convergence. Since the initial variables of the patches and their migration processes are independent, it suffices to prove the joint convergence of Ī N,0 ℓ,1 , . . . ,Ī N,0 ℓ,L for different ℓ's separately. On the other hand, the joint convergence of Ī N,0 ℓ,1 , . . . ,Ī N,0 ℓ,L is straightforward since they count an exclusive partition of individuals and the associated pairs η 0 k,ℓ , X 0,k ℓ (·) are independent. This completes the proof.
ForĪ N ℓ,i , we first consider their conditional expectations
Lemma 4.5. With the limit Ā 1 , . . . ,Ā L of the convergent subsequence of
Proof. We first focus onȊ N 1,1 . Observe that
Then we can writȇ
where dG 1,1 (t − u) is the differential of the map u → G 1,1 (t − u). Thus, by the continuous mapping theorem, we obtainȊ N 1,1 ⇒Ī 1,1 in D as n → ∞, wherē
For the joint convergence, it can be shown that the mapping from (x 1 , . . . , x K ) ∈ D L to
in D L 2 , is continuous in the Skorohod J 1 topology. Then we can apply the continuous mapping theorem to prove the joint convergence of the processes. This completes the proof.
We then show that the processesĪ N ℓ,i andȊ N ℓ,i are asymptotically negligible for each ℓ, i = 1, . . . , L. 
Proof. We focus on the case ℓ, i = 1. We havē
Recall G 1,1 (t) in (4.10). Then it is clear that for each j, E χ N j,1 (t)|τ N j,1 = 0, and
And by the independence of the pairs η j,1 , X j 1 (·) and η j ′ ,1 ,
which implies that for any ǫ > 0,
Observe that the three terms on the right hand side are nondecreasing in u. Thus we obtain P sup
Using the PRM Q 1,inf (ds, du, dv, dθ) and its compensated PRM Q 1,inf (ds, du, dv, dθ), we have
The first term converges to zero as N → ∞, and the second term satisfies
which follows from Assumption 2.1, as shown in Lemma 4.3 in [24] . The second term on the right hand side of (4.13) can be treated similarly as the second term right above. Now for the third term, using (4.1),
By (4.4), the first term converges to zero as N → ∞. The second term is bounded by 2λ 2 1 δ 2 by (4.3).
Thus, combining the above, we obtain
The result now follows from the next Lemma. 
Proof. The Lemma is a direct consequence of the inequality (4.21) in [24] , which we repeat here for the reader's convenience:
As a consequence of Lemmas 4.5 and 4.6, we obtain the following. 
We are now ready to prove the convergence of (Ī N 1 , . . . ,Ī N L ). 
withĪ 0 ℓ,i andĪ ℓ,i being defined in (4.6) and (4.8), respectively.
Proof. By the representations of I N i (t) in (4.5), we havē
where 
withĪ 0 i,j andĪ i,j being defined in (4.6) and (4.8), respectively.
Proof. We can represent the processesR N i (t) bȳ 18) where ∆ N I,ℓ,i (t) is given (4.15), and for ℓ = i, Thus, again, by Lemmas 4.1, 4.4, 4.5, 4.6 and 8.1, we apply the continuous mapping theorem to the mappingΥ to conclude the convergence of (R N 1 , . . . ,R N L ). From the argument above, since we have the joint convergence in Lemmas 4.4 and 4.6, we can conclude the joint convergence of (S N i ,Ī N i ,R N i , i = 1, . . . , L). In addition, the mapping (x, y, z) ∈ D 3 → xy x+y+z ∈ D is Lipschitz and continuous in the Skorohod topology. Thus, we obtain the convergenceΦ
and thus,
Therefore, all the limits satisfy the integral equations given in Theorem 2.1. Finally, the uniqueness of solutions to the set of integral equations in Theorem 2.1 follows from the Lipschitz continuity of the mapping (x, y, z) ∈ D 3 → xy x+y+z ∈ D and applying Gronwall's inequality. This completes the proof of Theorem 2.1.
In this section we prove Theorem 2.2. We first provide the following representations of the diffusion-scaled processes. The processÂ N i (t) can be decomposed as:
,
For the processesŜ N i (t), we havê For the processesÎ N i (t), we havê
where for ℓ, i = 1, . . . , L,
For the processesR N i (t), we havê
where for ℓ, i = 1, . . . , L, and ℓ = i,
We establish the convergence of some key components in these representations. Proof. This follows from a standard martingale convergence argument, see, e.g., [25] . The main steps include proving that the quadratic variations converge (involving the convergence of fluidscaled processes) and then applying the FCLT for martingales. We omit the details for brevity. where the limits are as given in Theorem 2.2.
Proof. We first focus on the convergence ofÎ N,0 1,1 . Recall that
Observe that the pairs η k,1 , X 0,k 1 (·) and η k ′ ,1 , X 0,k ′ 1 (·) are independent and have the same law. Thus, its proof follows in a similar approach for empirical processes, see, e.g., Theorem 14.3 in [5] . There are some differences due to the process X 0,k 1 , which we highlight below. So, we apply Theorem 13.5 in [5] .
For each t > 0 and ϑ > 0, we have
Similarly, it can be also shown that for any 0 < s < t,
as N → ∞. Thus, for the convergence of finite dimensional distributions, with t 1 < t 2 < · · · < t k and ϑ ℓ , ℓ = 1, . . . , k, we can write k ℓ=1 iϑ ℓÎ N,0 1,1 (t ℓ ) using the incrementsÎ N,0 1,1 (t ℓ )−Î N,0 1,1 (t ℓ−1 ), which have covariances equal to zero over disjoint intervals.
Next, to prove tightness, we employ Theorem13.5 and verify condition (13.14) in [5] . We show that for r ≤ s ≤ t and for N ≥ 1,
for some constant C and φ(t) = t 0 p 1,1 (u)dF 0 (u) which is a nonnegative, nondecreasing and continuous function. Recall F 0 is continuous. This will enforce condition (13.14) in [5] , which according to Theorem13.5 implies tightness in D. Let By direct calculations, following similar steps in the proof of (14.9) in [5] , we obtain
Thus we have shown the convergence of the processÎ N,0 1,1 ⇒Î 0 1,1 . For the joint convergence Î N,0 ℓ,i , ℓ, i = 1, . . . L , since the variables and processes associated with patch ℓ and patch ℓ ′ are independent, it suffices to show the joint convergences Î N,0 ℓ,i , i = 1, . . . , L for different ℓ's separately. For the joint convergence Î N,0 ℓ,i , i = 1, . . . , L , we obtain tightness from that of each process as established above, so it suffices to show the joint convergence of their finite dimensional distributions. Take ℓ = 1, i = 1, 2 as an example. For 0 < t 1 < t 2 and ϑ 1 , ϑ 2 > 0,
This calculation can be extended to the computation of final dimensional distributions of Î N,0 1,1 ,Î N,0 1,2 . Therefore, we can conclude the joint convergence.
Finally, to prove that the limit processes are continuous when F 0 is continuous, since they are Gaussian, it suffices to show continuity in the quadratic mean [16] , that is, for all t > 0, lim s→t E Î N,0 ℓ,i (t) −Î N,0 ℓ,i (s) 2 = 0. This is easily checked from the continuity of the covariance functions. Therefore, the proof of this lemma is complete. where the limits are as given in Theorem 2.2.
Proof. Recall the PRM Q i,inf (ds, du, dv, dθ) and its compensated PRM Q i,inf (ds, du, dv, dθ), and the representation in (2.7). We write for ℓ, i = 1, . . . , L,
To prove the convergence, we define the auxiliary processes:
where Φ ℓ (t) is given in (2.12) and is a deterministic function. It can then be shown that the processesĨ N ℓ,i (t) are square-integrable martingales with respect to the filtration {F N I (t) : t ≥ 0} defined by
. . , L ∨ σ (η j,ℓ , X j ℓ (s)) : ℓ = 1, . . . , L, j = 1, . . . , A N ℓ (s), 0 ≤ s ≤ t .
They have quadratic variations
In addition, their cross quadratic variations: for i = i ′ , Ĩ N ℓ,i ,Ĩ N ℓ,i ′ (t) = 0, and for ℓ = ℓ ′ and any i, i ′ , Ĩ N ℓ,i ,Ĩ N ℓ ′ ,i ′ (t) = 0. Thus, by the FCLT for martingales (see, e.g., [30] ), we obtain Ĩ N ℓ,i , ℓ, i = 1, . . . , L ⇒ Î ℓ,i , ℓ, i = 1, . . . , L in D L 2 as N → ∞. for each ℓ, i = 1, . . . , L. We focus on the processÎ N 1,1 . It is clear that
where the convergence holds by Theorem 2.1 and the dominated convergence theorem. We next show that the sequence {Î N 1,1 −Ĩ N 1,1 } N is tight. Observe that
We can decompose sign
such that each of these will induce a process that is nondecreasing in t. It is also clear that tightness of these processes will be implied by the tightness of the following processes:
Since these two processes are nondecreasing in t, it suffices to show (see the Corollary on page 83 in [5] or Lemma 4.7) that for any ǫ > 0, and j = 1, 2, 
for each s ≥ 0, and i = 1, . . . , L. In the representations ofŜ N i (t),Î N i (t) andR N i (t) in (5.4)-(5.7), respectively, the following hold: there exists a constant C > 0 such that for all N ,
Thus, by taking squares of the processesŜ N i (t),Î N i (t) andR N i (t) in (5.4)-(5.7), we can apply Cauchy-Schwartz inequality and Gronwall's inequality to conclude that 
Consider the two cases of F in Assumption 2.1. In the case F = F 1 ,
Thus, in both cases, we obtain (5.18). For the process Ξ N 2 (t), we have
The argument for these two terms follow from that for the second and fourth terms above for Ξ N 1 (t). This completes the proof.
Completing the proof of Theorem 2.2. Note that the processes associated with the initial variables in Lemma 5.2, those associated with the newly infected individuals in Lemma 5.3, and those associated with the migrations in Lemma 5.1 are mutually independent. Thus, all these processes converge jointly. Given the convergence of the initial variables in Assumption 2.3, by the convergence results in Lemmas 5.1, 5.2 and 5.3, we can apply the continuous mapping theorem to the map ̥ defined in Lemma 8.2 below (with m = L). The proof is complete. ✷
In this section we prove Theorem 3.1. Note that for the SEIR model, the tightness result in Lemma 4.1 ofĀ N i , and the convergence ofS N i in Lemma 4.2 hold with the same argument. The analysis of the processĒ N i follows the same as that of the processĪ N i in the SIR model, so we will omit its proof. We give the following representation ofĒ N i as in Lemma 4.3:
Similarly, we obtain the following representations for the process I N i (t):
We focus on the convergence of I N i (t). For ℓ, i = 1, . . . , L, let
We first treat the components associated with the initial quantities. 
which implies by the LLN of i.i.d. variables that for each t ≥ 0,
as n → ∞. The convergence of finite dimensional distribution is a straightforward extension. For tightness we use the same argument as in Lemma 4.4.
We start with Ī N,0,2
We newt note that each of the two terms on the right hand side are increasing in s, so that
as N → ∞. And the two terms in the limit both vanishes as δ → 0. Thus, for any ǫ > 0,
Thus, we can conclude thatĪ N,0,2 ℓ,i ⇒Ī 0,2 ℓ,i in D as N → ∞. Then, since G 0 is continuous, we can verify the continuity of the covariance function, and thus the continuity of the limit processesĪ 0,2 ℓ,i . ), which is straightforward. This completes the proof.
As in the proofs of Lemmas 4.5 and 4.6 for the SIR model, we obtain the convergence of the processes E N ℓ,i (t) and I N ℓ,i (t). where for ℓ, i = 1, . . . , L and t ≥ 0,
with H ℓ,i defined in (3.14).
Proof. The proof of the convergence ofĒ N ℓ,i follows from the exact same argument as in Lemmas 4.5 and 4.6. So we focus on the convergence ofĪ N ℓ,i . The main steps are similar, so we highlight the differences.
Define the auxiliary processes: for ℓ, i = 1, . . . , L,
We first show that Ȋ N ℓ,i , ℓ, i = 1, . . . , L ⇒ Ī ℓ,i , ℓ, i = 1, . . . , L in D L 2 as N → ∞. (6.10)
Observe that
Recall H ℓ,i (t) in (3.14). Then we can writȇ
and apply the continuous mapping theorem (as in the proof of Lemma 4.5) to obtain the convergencȇ I N ℓ,i ⇒Ī N ℓ,i in D as n → ∞, and then the joint convergence in (6.10) using the mapping (4.11). We next show that for any ǫ > 0, and for ℓ, i = 1, . . . , L, P sup
This follows from similar steps as in the proof of Lemma 4.6. We highlight the main differences below.
where
It is clear that E χ N j,ℓ,i (t)|τ N j,ℓ = 0 and E χ N j,ℓ,i (t) 2 |τ N j,ℓ = H ℓ,i (t−τ N j,ℓ )(1−H ℓ,i (t−τ N j,ℓ )) where H ℓ,i (t) is defined in (3.14) . Moreover E χ N j,ℓ,i (t)χ N j ′ ,ℓ,i (t)|F N A,ℓ (t) = 0 due to the independence of the pairs ζ j,ℓ , η j,ℓ ′ , Y j ℓ (·), X j ℓ ′ (·) and ζ j ′ ,ℓ , η j ′ ,ℓ ′ , Y j ′ ℓ (·), X j ′ ℓ ′ (·) . Thus, we obtain
. Observe that each of the five terms on the right hand side is increasing in s. Thus, we have P sup
The last term is treated in the same way as in the proof of Lemma 4.6 for the SIR model. For the first two terms, we use the PRM representation in (3.6). We have
It is clear that the first term converges to zero as N → ∞, and the second term satisfies
The convergence follows the same argument as for the proof of (6.11) in [24] under Assumption 3.1.
Similarly, for the second term on the right hand side of (6.13), we have
Here it is clear that the first term converges to zero as N → ∞, and the second term satisfies
The convergence follows a similar argument as for the proof of (6.9) in [24] under Assumption 3.1. Now for the third and fourth terms on the right hand side of (6.13), we have
These two terms can be treated similarly as above. Thus we have shown that
and we can apply Lemma 4.7 to conclude the tightness.
For the joint convergence of the processes Ē N ℓ,i , ℓ, i = 1, . . . , L and Ī N ℓ,i , ℓ, i = 1, . . . , L , it suffices to prove the joint convergence of the processes Ȇ N ℓ,i , ℓ, i = 1, . . . , L and Ȋ N ℓ,i , ℓ, i = 1, . . . , L , whereȊ N ℓ,i(t) is given in (6.11) , and
14)
It is easy to prove that the mapping from (x 1 , . . . ,
is continuous in the Skorohod J 1 topology. Thus, applying the continuous mapping theorem, we obtain the convergence Ȇ N ℓ,i ,Ȋ N ℓ,i , ℓ, i = 1, . . . , L ⇒ Ē ℓ,i (t),Ī ℓ,i , ℓ, i = 1, . . . , L in D 2L 2 as N → ∞. Then the asymptotic equivalence betweenĒ N ℓ,i andȆ N ℓ,i , and betweenĪ N ℓ,i andȊ N ℓ,i results in the joint convergence in (6.7). This completes the proof.
We are now ready to prove the convergence of (Ī N 1 , . . . ,Ī N L ). 
withĒ 0 ℓ,i ,Ī 0,1 ℓ,i andĪ 0,2 ℓ,i being given in (6.4), (6.5) and (6.6) in Lemma 6.1, and withĒ ℓ,i andĪ ℓ,i being defined in (6.8) and (6.9) in Lemma 6.2, respectively.
Proof. By the representations of I N i (t) in (6.2), we havē (6.16) andM N I,ℓ,i (t) is given in (4.16) , and its convergence in (4.17) holds. Recall the representations of E N ℓ,i (t) andȊ N ℓ,i (t) as integrals with respect toĀ N ℓ in (6.14) and (6.11), respectively. Thus, by Lemmas 4.1, 6.1, and 6.2 and equation (6.12) , and by Lemma 8.1, we apply the continuous mapping theorem to the mapping Υ to conclude the convergence of Ī N 1 , . . . ,Ī N L . We next prove the convergence of (R N 1 , . . . ,R N L ). Similar to (6.2), we obtain the following representations for the process R N i (t): 
withĪ 0,1 ℓ,i ,Ī 0,2 ℓ,i andĪ ℓ,i being given in (6.5), (6.6) and (6.9), respectively. Proof. We can represent the processesR N i (t) bȳ 18) where ∆ N I,ℓ,i (t) is given (4.15) , andM N R,ℓ,i (t) is given in (4.19) , and its convergence in (4.20) holds. Recall the representation ofȊ N ℓ,i (t) as an integral with respect toĀ N ℓ in (6.11). Thus, again, by Lemmas 4.1, 6.1, and 6.2 and equation (6.12) , and by Lemma 8.1, we apply the continuous mapping theorem to the mappingΥ to conclude the convergence of (R N 1 , . . . ,R N L ).
The rest of the proof follows from the same argument as for the SIR model. This completes the proof of Theorem 3.1.
In this section, we prove Theorem 3.2. We first give the following representations of the diffusionscaled processes. The processÂ N i (t) has the same decomposition as in (5.1), but witĥ
We have the same representation of the processŜ N i (t) in (5.4 ). For the processÊ N i (t), by the representation in (6.1), using the definitions of E N,0 ℓ,i (t) and E N ℓ,i (t) in (6.3), we obtain
where for ℓ, i = 1, . . . , L,
and
For the processÎ N i (t), we obtain 5) whereÊ N,0 ℓ,i (t) andÊ N ℓ,i (t) are defined in (7.3) and (7.4), respectively, and
andM N I,ℓ,i (t) is as defined in (5.6) for the SIR model. For the processR N i (t), we havê
whereM N R,ℓ,i (t) is as defined in (5.8) for the SIR model. The convergence of the processesÊ N i follows the same argument for that ofÎ N i in the SIR model, so we focus on the convergence ofÎ N i in the SEIR model which is given in (7.5) . First, as in Lemma 5.1, we have the joint convergence of the following martingales, whose proof is omitted. We now prove the convergence of the components associated with the initially exposed and infected individuals. . We again apply Theorem 13.5 in [5] . By direct calculations, we obtain for t ≥ 0,
and for t > t ′ ≥ 0,
and then establish the convergence of finite dimensional distributions ofÎ N,0,2 ℓ,i . Also, we obtain for t ′ ≤ t ≤ t ′′ and for N ≥ 1, , ℓ, i = 1, . . . , L , which also follows similarly as in Lemma 5.2. Finally for the continuity of the limit processes, it suffices to show the continuity in the quadratic mean [16] , which follows from the continuity of the covariance functions.
We next prove the convergence of the components associated with the newly exposed individuals. where the limits are as given in Theorem 3.2.
Proof. The convergence of Ê N ℓ,i ,Î N ℓ,i , ℓ, i = 1, . . . , L follows from a similar argument as that of Î N ℓ,i ,Î N ℓ,i , ℓ, i = 1, . . . , L for the SIR model in Lemma 7.3. It uses the PRMQ ℓ,inf (ds, du, dv, dϑ) and the representation in (3.5) . In particular, we can writê
. Define the auxiliary process
where Φ ℓ (t) is given in (3.12) . It can be proved in the convergence ofÊ N ℓ,i , we also havê E N ℓ,i −Ẽ N ℓ,i ⇒ 0 in D as N → ∞ for each ℓ, i = 1, . . . , L.
We focus on the convergence ofÎ N ℓ,i (t). Recall the PRM Q ℓ,inf (ds, du, dy, dϑ, dz, dθ) and the representation in (3.6). We can then writê
We define the auxiliary process
where Φ ℓ (t) is given in (3.12) . Since Φ ℓ (t) is a deterministic function, thenĨ N ℓ,i (t) is a squareintegrable martingale with respect to the filtration {F N I (t) : t ≥ 0} defined by
It has the quadratic variation
where H ℓ,i (t) is defined in (3.14), and the cross quadratic variations of the processesĨ N ℓ,i (t) and I N ℓ ′ ,i ′ (t) satisfy: for i = i ′ , Ĩ N ℓ,i ,Ĩ N ℓ,i ′ (t) = 0, and for ℓ ′ = ℓ and any i, i ′ = 1, . . . , L, Ĩ N ℓ,i ,Ĩ N ℓ ′ ,i ′ (t) = 0. Thus, by the FCLT for martingales (see, e.g., [30] ), we obtain Ĩ N ℓ,i , ℓ, i = 1, . . . , L ⇒ Î ℓ,i , ℓ, i = 1, . . . , L in D L 2 as N → ∞. Now it remains to show thatÎ N ℓ,i −Ĩ N ℓ,i ⇒ 0 in D as N → ∞ for each ℓ, i = 1, . . . , L. It is clear that
where the convergence holds by Theorem 3.1 and the dominated convergence theorem. We next show that the sequence {Î N ℓ,i −Ĩ N ℓ,i } is tight. We havê
As in the proof of Lemma 5.3, it suffices to prove tightness of the following processes
By the monotone property of these two processes in t, we then show that for any ǫ > 0, and κ = 1, 2,
Similar to the derivation in (5.13) of Lemma 5.3, we obtain
(7.11)
The first and third terms converge to zero as N → ∞ by the convergence E Φ N ℓ (s) − Φ ℓ (s) → 0 and applying the dominated convergence theorem. By a similar argument as in Lemma 5.3, we also obtain lim sup N sup and applying Cauchy-Schwartz inequality and Gronwall's inequality, we obtain (5.17) together with the bound forÊ N i (t), and thus the claim follows. For the last term in (7.11), we have
We can bound the sum of those two terms by
By the Lipschitz continuity condition in the case of G 2 and F 2 in Assumption 3.1, say, with Lipschitz coefficient θ ′ and θ and constants c ′ and c, respectively, we obtain the bounds for the right hand side:
Now suppose that both functions G and F are discrete, say, G(t) = G 1 (t) = k b k 1 t≥s k for a finite or countable number of b k and the corresponding s k such that k b k ≤ 1 and s 0 < s 1 < · · · < s k < . . . , and F (t) = F 1 (t) = i a i 1(t ≥ t i ) for a finite or countable number of positive numbers a i and the corresponding t i such that i a i ≤ 1 and t 0 < t 1 < · · · < t k < . . . . The first term on the right hand side of (7.12) is equal to
Similarly, the second term on the right hand side of (7.12) is bounded by
Thus, in both cases, we obtain lim δ→0 lim sup
The other cases where G (resp. F ) is discrete and F (resp. G) is Lipschitz continuous can be treated similarly without difficulty. Thus, (7.10) holds for Ξ N 1 . Consider next
We have
Here the first term is bounded by 2λ 2 ℓ δ 2 sup s∈[0,T ] E Φ N ℓ (s) 2 and the second term is treated as above in (7.13) . This completes the proof of (7.10), and thus the convergence Î N ℓ,i , ℓ, i = 1, . . . , L ⇒ Î ℓ,i , ℓ, i = 1, . . . , L in D L 2 as N → ∞. It remains to prove the joint convergence of (Ê N ℓ,i ,Î N ℓ,i ), for which it suffices to prove the joint convergence of (Ẽ N ℓ,i ,Ĩ N ℓ,i ). Recall the representations of them in (7.8) and (7.9) using the PRMš Q ℓ,inf (ds, du, dv, dϑ) and Q ℓ,inf (ds, du, dy, dϑ, dz, dθ), respectively. By their definitions of the two PRMs, we can regardQ ℓ,inf as projections of Q ℓ,inf from R 3 + × {1, . . . , L} × R + × {1, . . . , L} onto R 3 + × {1, . . . , L}. Thus it is straighrforward to prove the convergence of the finite-dimensional distributions of (Ẽ N ℓ,i ,Ĩ N ℓ,i ) by computing
and their extensions to multiple time points. Since tightness of (Ẽ N ℓ,i ,Ĩ N ℓ,i ) follows from that of each individual process, we can conclude their joint convergence in D L 2 . This completes the proof.
Completing the proof of Theorem 2.2. We observe the correspondence of the processeŝ Φ N i in (7.1), andŜ N i ,Ê N i ,Î N i ,R N i in (5.4), (7.2), (7.5), and (7.6) with the mapping̥ defined below in Lemma 8.3. We have shown the convergence of the components associated with the initial variables, newly exposed individuals, and migration processes in Lemmas 7.2, 7.3, 7.1, which are all mutually independent. Thus, given the convergence of the initial variables in Assumption 3.3, we can apply the continuous mapping theorem to the map̥ with m = L, and conclude the joint convergence of the processes (Ŝ N i ,Ê N i ,Î N i ,R N i ). ✷
In this section we define several integral mappings and prove their continuity in the Skorohod J 1 topology.
Define the m-dimensional integral mapping Υ : (y 1 , . . . , y m , z 1 , . . . , z m ) ∈ D 2m → (x 1 , . . . , x m ) ∈ D m , where the y j 's are nondecreasing:
x i (t) = x i (0) + y i (t) + z i (t) + for given constants a i,j , b i,j ∈ R and fixed functions φ i,j in D, i, j = 1, . . . , m. LetΥ be the mapping without y i (t) terms but with the integral with respect to y i (t) only. This is an extension of the mapping in Section 8.1 of [24] , whose proof can be easily extended to this mapping. Thus we omit the proof for brevity.
Lemma 8.1. Given (y 1 , . . . , y m , z 1 , . . . , z m ) ∈ D 2m with (y 1 , . . . , y m ) being nondecreasing, a i,j , b i,j ∈ R, φ i,j ∈ D in D, i, j = 1, . . . , m, and (x 1 (0), . . . , x m (0)) ∈ R m , there exists a unique solution (x 1 , . . . , x m ) ∈ D m to the integral mapping Υ. If (y 1 , . . . , y m , z 1 , . . . , z m ) ∈ C 2m , then (x 1 , . . . , x m ) ∈ C m . The map Υ is continuous in the Skorohod J 1 topology, that is, if (y n 1 , . . . , y n m , z n 1 , . . . , z n m ) → (y 1 , . . . , y m , z 1 , . . . , z m ) in D 2m , then (x n 1 , . . . , x n m ) = Υ(y n 1 , . . . , y n m , z n 1 , . . . , z n m ) → (x 1 , . . . , x m ) in D m , and if, in addition, (y 1 , . . . , y m , z 1 , . . . , z m ) ∈ C 2m , then (x 1 , . . . , x m ) ∈ C m , and the convergence holds uniformly on compacts. The same conclusions hold for the mappingΥ.
We next define a 3m-dimensional integral mapping ̥: given a i , b i , c i , φ i , ψ i , ϕ i ∈ D, some constants α i , β i , γ i > 0 and functions F ℓ,i for ℓ, i = 1, . . . , m, let x i , y i , z i be the solutions to the following set of integral equations: In the next lemma, we study the existence and uniqueness of its solution and the continuity property.
Lemma 8.2. Assume that F ℓ,i , ℓ, i = 1, . . . , m are bounded and continuous functions satisfying F ℓ,i (0) = 0, and let the constants α i , β i , γ i > 0 and the functions φ i , ψ i , ϕ i be given. There exists a unique solution (x i , y i , z i , i = 1, . . . , m) ∈ D 3m to the set of integrable equations defining the mapping ̥. The mapping is continuous in the Skorohod J 1 topology, that is, if (a n i , b n i , c n i , φ n i , ψ n i , ϕ n i , i = 1, . . . , m) → (a i , b i , c i , φ i , ψ i , ϕ i , i = 1, . . . , m) in D([0, T ], R 6m ) as n → ∞ and (x n i (0), y n i (0), z n i (0), i = 1, . . . , m) → (x i (0), y i (0), z i (0), i = 1, . . . , m), then (x n i , y n i , z n i , i = 1, . . . , m) → (x i , y i , z i , i = 1, . . . , m) in D([0, T ], R 3m ) as n → ∞. In addition, if φ i , ψ i , ϕ i are continuous, then (x i , y i , z i , i = 1, . . . , m) ∈ C 3m and the mapping ̥ is continuous uniformly on compact sets in [0, T ].
Proof. For the existence and uniqueness of solutions, we can apply the Schauder-Tychonoff fixed point theorem, and modify the proofs in Theorems 1.2 and 2.3 in Chapter II of [23] (where these results are shown for Volterra integral equations with continuous functions).
We now prove the the continuity of the mapping in the Skorohod J 1 topology. Note that functions in D are necessarily bounded. For the given functions φ i , ψ i , ϕ i ∈ D, let T be a common continuity point. Then there exists increasing homeomorphisms λ n on [0, T ] such that λ n − e T → 0, φ n i −φ i •λ n T → 0, ψ n i −ψ i •λ n T → 0, ϕ n i −ϕ i •λ n T → 0, a n i −a i •λ n T → 0, b n i −b i •λ n T → 0, and c n i − c i • λ n T → 0, as n → ∞. Here e(t) := t for all t ≥ 0. It also suffices to consider homeomorphisms λ n that are absolutely continuous with respect to the Lebesgue measure on [0, T ] having derivativesλ n satisfying λ n − 1 T → 0 as n → ∞. Let C 0 be some constant such that sup ℓ,i sup t∈[0,T ] |F ℓ,i (t)| ≤ C 0 .
For brevity, we consider one of the functions, y 1 , since the others can be treated similarly. We have |y n 1 (t) − y 1 (λ n (t))| ≤ |y n 1 (0) − y 1 (0)| + ψ n 1 − ψ 1 • λ n T + t 0 (a n 1 (s)x n 1 (s) + b n 1 (s)y n 1 (s) + c n 1 (s)z n 1 (s))ds − λ n (t) 0 (a 1 (s)x 1 (s) + b 1 (s)y 1 (s) + c 1 (s)z 1 (s))ds + m ℓ=1 t 0 F ℓ,1 (t − s)(a n ℓ (s)x n ℓ (s) + b n ℓ (s)y n ℓ (s) + c n ℓ (s)z n ℓ (s))ds − λ n (t) 0 F ℓ,1 (λ n (t) − s)(a ℓ (s)x ℓ (s) + b ℓ (s)y ℓ (s) + c ℓ (s)z ℓ (s))ds + m ℓ=2 t 0 (β ℓ,1 y n ℓ (s) − β 1,ℓ y n 1 (s))ds − λ n (t) 0 (β ℓ,1 y ℓ (s) − β 1,ℓ y 1 (s))ds .
By the change of variables for the second integral, and by the boundedness of the functions x n 1 , y n 1 , z n 1 and a 1 (λ n (s)), b 1 (λ n (s)), c 1 (λ n (s)), the third term on the right hand side of (8.2) is bounded by λ n − 1 T T 0 |a 1 (s)x 1 (s) + b 1 (s)y 1 (s) + c 1 (s)z 1 (s)|ds + t 0 |a n 1 (s) − a 1 (λ n (s))||x n 1 (s)| + |b n 1 (s) − b 1 (λ n (s))||y n 1 (s)||c n 1 (s) − c 1 (λ n (s))||z n 1 (s)| ds + t 0 |a 1 (λ n (s))||x 1 (λ n (s)) − x 1 (s)| + |b 1 (λ n (s))||y 1 (λ n (s)) − y 1 (s)| + |c 1 (λ n (s))||z 1 (λ n (s)) − z 1 (s)| ds ≤ λ n − 1 T T 0 |a 1 (s)x 1 (s) + b 1 (s)y 1 (s) + c 1 (s)z 1 (s)|ds + C 1 a n 1 − a 1 • λ n T + b n 1 − b 1 • λ n T + c n 1 − c 1 • λ n T + t 0 |x 1 (λ n (s)) − x 1 (s)| + |y 1 (λ n (s)) − y 1 (s)| + |z 1 (λ n (s)) − z 1 (s)| ds (8.3)
for some constant C 1 , where the first term converges to zero since λ n − 1 T → 0 and the next three terms converge to zero because of the convergences of a n 1 → a 1 , b n 1 → b 1 and c n 1 → c 1 in the Skorohod topology.
For the fourth term on the right hand side of (8.2), by the change of variables for the second integrals, we bound the term with ℓ = 1 by λ n − 1 T T 0 |F 1,1 (λ n (t) − λ n (s))||a 1 (s)x 1 (s) + b 1 (s)y 1 (s) + c 1 (s)z 1 (s)|ds + t 0 F 1,1 (t − s) − F 1,1 (λ n (t) − λ n (s)) a n 1 (s)x n 1 (s) + b n 1 (s)y n 1 (s) + c n 1 (s)z n 1 (s) ds + t 0 F 1,1 (λ n (t) − λ n (s)) a n 1 (s)x n 1 (s) + b n 1 (s)y n 1 (s) + c n 1 (s)z n 1 (s) − a 1 (λ n (s))x 1 (λ n (s)) − b 1 (λ n (s))y 1 (λ n (s)) − c 1 (λ n (s))z 1 (λ n (s)) ds.
Here the first two terms converge to zero since λ n − 1 T → 0 and F 1,1 is continuous, respectively, and also because the functions in D are bounded over [0, T ]. For the third term, similar to the second and third terms in (8.3) , it can be bounded by C 0 C 1 a n 1 − a 1 • λ n T + b n 1 − b 1 • λ n T + c n 1 − c 1 • λ n T + t 0 |x 1 (λ n (s)) − x 1 (s)| + |y 1 (λ n (s)) − y 1 (s)| + |z 1 (λ n (s)) − z 1 (s)| ds .
For the last term on the right hand side of (8.2), again by the change of variables, we can bound the ℓ-component by λ n − 1 T T 0 β ℓ,1 y ℓ (λ n (s)) − β 1,ℓ y 1 (λ n (s)) ds + β ℓ,1 t 0 y n ℓ (s) − y ℓ (λ n (s)|ds + β 1,ℓ t 0 y n 1 (s) − y 1 (λ n (s)) ds.
Thus, the continuity property in the Skorohod J 1 topology holds by applying Gronwall's inequality to the set of integral inequalities for x i , y i , z i . The uniform continuity property on compacts when the functions are continuous is straightforward. This completes the proof.
Finally, we define a 4m-dimensional integral mapping̥: given a i , b i , c i , d i , φ i , ψ i , ϕ i , χ i ∈ D, some constants α i , β i , γ i , κ i > 0 and functions F ℓ,i , G ℓ,i for ℓ, i = 1, . . . , m, let x i , y i , z i , w i be the solutions to the following integral mapping: 