{'fraction_non_alphanumeric': 0.05548143568383717, 'fraction_numerical': 0.015261462366990986, 'mean_word_length': 3.8853827199745647, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fluctuation analysis for particle-based stochastic reaction-diffusion models

Recent works have derived and proven the large-population mean-field limit for several classes of particle-based stochastic reaction-diffusion (PBSRD) models. These limits correspond to systems of partial integral-differential equations (PIDEs) that generalize standard mass-action reaction-diffusion PDE models. In this work we derive and prove the next order fluctuation corrections to such limits, which we show satisfy systems of stochastic PIDEs with Gaussian noise. Numerical examples are presented to illustrate how including the fluctuation corrections can enable the accurate estimation of higher order statistics of the underlying PBSRD model.


Introduction
In this work we consider the large-population limit of particle-based stochastic reaction-diffusion (PB-SRD) models. Such models are used to study the interplay between spatial transport via diffusion and chemical reactions for systems of interacting particles/molecules, and have been used in a variety of biological and population process studies. Several recent works [LLN19,IMS22,PV22] have investigated the large-population (thermodynamic) limit of such models, first studied in [O89], where the number of particles and a system size parameter (typically domain volume and/or Avogadro's number) are allowed to become infinite such that the initial concentration of different chemical species (i.e. different types of particles) is held constant. In such limits one obtains systems of partial-integral differential equations (PIDEs) for the associated macroscopic concentration field of each species.
Here we build on our previous studies in [IMS22] and [IMS21], in which we investigated the largepopulation mean-field limit for the volume reactivity (VR) PBSRD model, which is a generalization of the classical Doi model [TS67,D76a,D76b]. In the VR model, molecules are approximated as point particles that move by Brownian Motion, unimolecular reactions such as A → B occur as an internal Poisson process with a fixed probability per time for each substrate molecule (i.e. each A molecule), and bimolecular reactions between two individual reactant particles, like A + B → . . . , occur with a probability per time based on the positions of the two substrate molecules. Products of a reaction are placed randomly based on specified particle placement densities. In [IMS22] we derived a set of reaction-diffusion PIDEs that we then proved corresponded to the rigorous large-population limit of the VR PBSRD model, with our analysis being valid for general reaction networks involving both first and second order reactions in which particle numbers can change. Our method of proof was based on formulating a weak representation for the dynamics of measure-valued stochastic processes (MVSPs), which correspond to concentration fields of different chemical species in the VR PBSRD model. We then adapted the classical Stroock-Varadhan Martingale approach to calculating mean-field limits [EK86,SV06] to prove the large-population limit for the MVSPs. A similar approach was used in [O89] to derive the analogous mean-field PIDEs for systems that involved only linear zero and first-order reactions (and therefore no interactions between particles). In contrast, a different approach was considered in [LLN19], which derived the analogous mean-field limit using relative entropy methods, but for systems involving only reactions that conserve particle numbers. Using entropy methods allowed for more quantitative estimates, but as pointed out by the authors, appears restricted to systems in which particle numbers are unchanged by each reaction. Finally, we note that [PV22] considers a novel underlying PBSRD model that is distinct from the VR model. This enables the consideration of higher than second order reactions, and is used in studying a large-population limit in which some species remain discrete with a stochastic evolution while others converge to continuous concentration fields.
We next investigated in [IMS21] how the derived mean-field PIDEs relate to commonly used, but phenomenological, mass action reaction-diffusion PDE models. We demonstrated that the latter could be interpreted as an approximation to the rigorous mean-field PIDEs in the limit that bimolecular interactions become short-ranged when bimolecular reaction kernels have an averaging form. Collectively, the works [O89, LLN19, IMS22, IMS21] establish the rigorous large population limit for VR PBSRD models, and then rigorously demonstrate how they relate to commonly used macroscopic reaction-diffusion PDE models.
The derived mean-field model has several clear limitations. In particular, as the resulting PIDEs are deterministic it is unable to give any approximation to higher order statistics of the underlying VR PBSRD models. To address this issue, in this work we derive the next order fluctuation corrections to the mean-field PIDEs we derived in [IMS22], giving a proof that the scaled difference between the PBSRD MVSPs and the mean-field PIDE solutions, see (1.1), converges to the derived fluctuation corrections in the large-population limit.
Our analysis allows us to obtain an approximation to the distribution of the empirical measure of the system to leading order in the system size parameter. Denote by µ γ,j t the marginal distribution (molar concentration) of species j = 1 · · · , J, where γ represents the system size parameter (e.g. Avogadro's number in R d ). Lettingμ t be the limit of µ γ,j t as γ → ∞, we study the behavior of the fluctuation corrections Ξ γ,j t = √ γ µ γ,j t −μ j t , j = 1, · · · J, (1.1) as γ → ∞. We prove that, in the appropriate sense, Ξ γ,j t →Ξ j t as γ → ∞ withΞ j · being Gaussian. We therefore obtain the approximation which expands on earlier mean-field limit studies by allowing the approximation of the distribution of µ γ,j t , as well as associated statistics.
The proof of convergence is based on weak convergence analysis of interacting particle systems. One of the main challenges is that the fluctuations of the empirical distribution correspond to a signed measure-valued process. As discussed in [BDG07,SSG14], the space of signed measures endowed with the weak topology is in general not metrizable. For this reason, we study convergence of the distribution-valued fluctuation process in the dual space of an appropriately weighted Sobolev space (weights are introduced to control the behavior at infinity). We use the weights introduced in [FM97] in studying the McKean-Vlasov equation, and show fluctuations arising from our VR PBSRD model are elements of these spaces, see Section 6 for a short exposition. The limit,Ξ j t , is then obtained as a distribution-valued process in an appropriate weighted negative Sobolev space. In the process of studying the limit of Ξ γ,j as γ → ∞, we need to linearize the dynamics of µ γ,j around its limiting behaviorμ j . Doing so produces a number of remainder terms that need to be appropriately controlled; details of these bounds are in Section 8.
As stating our main result, Theorem 5.6, requires significant notation and a number of assumptions, we begin in Section 2 by summarizing our results for the special case of the reversible A + B C reaction. In Section 3 we then summarize the main notation we will use in describing the weak MVSP representation for the VR PBSRD model. Section 4 summarizes our earlier work on the large-population mean-field limit from [IMS22], along with recalling the main assumptions under which we proved the large-population mean-field limit. In Section 5 we present the main result of this work, Theorem 5.6, along with stating the additional assumptions we introduce to prove this theorem. We summarize in Section 6 the weighted Sobolev spaces used to control behavior at infinity, and in which we prove convergence. In Section 7 we illustrate how including fluctuation corrections can enable the accurate approximation of higher-order (i.e. non-mean) statistics of an underlying VR PBSRD model. This is demonstrated via comparison of numerical solutions to the SPIDEs of the fluctuation process model and a jump process approximation to the VR PBSRD model [I13,IZ18]. Finally, in Section 8 we provide the proof of Theorem 5.6 via proving uniform spatial moment bounds of the PBSRD model in Section 8.2, proving tightness in Section 8.3, identifying the limiting SPIDEs satisfied by the fluctuation corrections in Section 8.4, and proving uniqueness of solutions to the limiting equations in Section 8.5. An appendix proves several results that are used in Section 8.

Summary of fluctuation corrections for the A + B C reaction
In order to introduce the main result of this paper, let us first present the result in the special case of the reversible A + B C reaction. We first summarize the large population limit derived in [IMS22], and then we show the fluctuation corrections we obtain in this work.
Let γ denote a system size parameter (i.e. Avogadro's number, or in bounded domains the product of Avogadro's number and the domain volume). We assume all molecules move by Brownian Motion in R d , with species-dependent diffusivities, D A , D B and D C respectively. Let K γ 1 (x, y) = K 1 (x, y)/γ denote the probability per time an individual A molecule at x and B molecule at y can react, with m 1 (z|x, y) giving the probability density that when the A and B molecules react they produce a C molecule at z. We define K γ 2 (z) = K 2 (z) and m 2 (x, y|z) similarly for the reverse reaction. Finally, denote by A(t) the stochastic process for the number of species A molecules at time t, and label the position of ith molecule of species A at time t by the stochastic process Q corresponds to the molar concentration of species A at x at time t. We can similarly define B γ (x, t) and C γ (x, t). In [IMS22] we derived the large population (thermodynamic) limit where γ → ∞ and A γ (x, 0) converges to a well defined limiting molar concentration field (with similar limits for the molar concentrations of species B and C). We proved, in a weak sense, that as γ → ∞, K 1 (y, z)m 1 (x|y, z)A(y, t)B(z, t) dy dz.
Here (A(x, t), B(x, t), C(x, t)) corresponds to the deterministic mean-field large-population limit.
In this work we derive the next order correction to these equations, obtaining an approximation as γ → ∞ of with similar equations for B(x, t) and C(x, t). We prove that in an appropriate (weak) sense, Here the fluctuation correction processes, (Ā(x, t),B(x, t),C(x, t)), satisfy the stochastic PIDEs where the noise processes, (ξ A t , ξ B t , ξ C t ), are mean-zero Gaussian processes. For test functions f and g they have the covariance structure for ∈ {1, . . . ,L}. We assume the remaining L −L reactions have one or more product particles. Let D i label the diffusion coefficient for the ith molecule, taking values in {D 1 , . . . , D J }, where D j is the diffusion coefficient for species S j , j = 1, · · · , J. We denote by Q i t ∈ R d the position of the ith molecule, i ∈ N + , at time t. A particle's state can be represented as a vector inP = R d × S, the combined space encoding particle position and type. This state vector is subsequently denoted byQ i t def = (Q i t , p i ). We now formulate our representation for the (number) concentration, equivalently number density, fields of each species. Let E be a complete metric space and M (E) the collection of measures on E. For f : We will frequently have E = R d . In this case we omit the subscript E and simply write f, µ . For each t ≥ 0, we define the concentration of particles in the system at time t by the distribution where, borrowing notation from [BM15], N (t) = 1, ν t P represents the stochastic process for the total number of particles at time t. To investigate the behavior of different types of particles, we denote the marginal distribution on the jth type, i.e. the concentration field for species j, by will similarly label the total number of particles of type S j at time t.
For ν any fixed particle distribution of the form (3.1), we will also use an alternative representation in terms of the marginal distributions ν j ∈ M (R d ) for particles of type j, In addition to having notations for representing particle concentration fields, we will also often make use of state vectors for the positions of particles of a given type. Define the particle index maps {σ j (k)} Nj (t) k=1 , which encode a fixed ordering for particles of species j, Q σj (1) t · · · Q σj (Nj (t)) t , arising from an (assumed) fixed underlying ordering on R d . We denote by the position state vector for type j particles, using the same ordering on R d , with H i (ν j t ) ∈ R d labeling the ith entry in H(ν j t ). Note, as particles of the same type are assumed indistinguishable, there is no ambiguity in the value of H(ν j t ) in the case that two particles of type j have the same position. With the preceding definitions, we last introduce a system of notation to encode reactant and particle positions and configurations that are needed to later specify reaction processes.
Definition 3.1. In describing the dynamics of ν t , we will sample vectors containing the indices of the specific reactant particles participating in a single -type reaction from the reactant index space For the allowable reactions considered in this work we label the elements of I ( ) in a manner that shows which species they belong to: (1) For R of the form ∅ → · · · I ( ) = ∅.
(2) For R of the form S j → · · · (3) For R of the form S j + S k → · · · with j < k A generic element i ∈ I ( ) can then be written as i = (i Definition 3.2. We define the product index space J ( ) analogously to I ( ) , with j ∈ J ( ) given by j = (j Definition 3.3. We define the reactant particle position space analogously to I ( ) , . For x ∈ X ( ) a sampled reactant position configuration for one individual R reaction, x (j) r then labels the sampled position for the rth reactant particle of species j involved in the reaction.
r ) be the corresponding volume form on X ( ) , which also naturally defines an associated Lebesgue measure.
Definition 3.4. For reaction R withL + 1 ≤ ≤ L, i.e. having at least one product particle, define the product position space analogously to X ( ) , with an element y ∈ Y ( ) written as y = (y be the corresponding volume form on Y ( ) , which also naturally defines an associated Lebesgue measure.
Definition 3.5. Consider a fixed reaction R , with i ∈ I ( ) and ν corresponding to a fixed particle distribution given by (3.1) with representation (3.2). We define the th projection mapping P ( ) : M (P ) × I ( ) → X ( ) as When reactants with indices i in particle distribution ν are chosen to undergo a reaction of type , P ( ) (ν, i) then gives the vector of the corresponding reactant particles' positions. For simplicity of notation, in the remainder we will sometimes evaluate P ( ) with inconsistent particle distributions and index vectors. In all of these cases the inconsistency will occur in terms that are zero, and hence not matter in any practical way.
Definition 3.6. Consider a fixed reaction R , with ν a fixed particle distribution given by (3.1) with representation (3.2). Using the notation of Def. 3.1, we define the allowable reactant index sampling space Note that in the calculations that follow Ω ( ) (ν) will change over time due to the fact that ν = ν t changes over time, but this will not be explicitly denoted for notational convenience.
Definition 3.7. Consider a fixed reaction R , with ν any element of M (P ) having the representation (3.2). We define the th reactant measure mapping λ ( ) Definition 3.8. For reaction R , define a subspaceX ( ) ⊂ X ( ) by removing all particle reactant position vectors in X ( ) for which two particles of the same species have the same position. That is

Previous Results and Assumptions for the Mean Field Limit
In this section we review our previous results on the mean field large-population limit from [IMS22] that are relevant to the current work.
Let us consider the time evolution of the processes ν ζ,j Sj which give the spatial distribution of particles of type j (i.e. number density or concentration). Here N ζ (t) = J j=1 N ζ,j (t) denotes the total number of particles in the system at time t, with N ζ,j (t) = 1, ν ζ,j t P the number of particles of type j at time t. ζ = ( 1 γ , η) ∈ (0, 1) 2 is a two-vector consisting of a scaling parameter, γ, and a displacement range parameter, η (which will be explained in Section 4.2).
In the large population limit we consider γ plays the role of a system size, and is considered to be large (e.g. Avogadro's number, or in bounded domains the product of Avogadro's number and the domain volume) [DK15]. On the other hand, η is a regularizing parameter allowing us to be able to consider and rigorously handle delta-function placement measures for reaction products (a common choice in many PBSRD simulation methods). We will further clarify these parameters later on, focusing on the (large-population) limit that γ → ∞ and η → 0 jointly, denoted as ζ → 0. 4.1. Generator and process level description. To formulate the process-level model, it is necessary to specify more concretely the reaction process between individual particles. For reaction R , denote by K γ (x) the rate (i.e. probability per time) that reactant particles with positions x ∈ X ( ) react. As described in the next section, we assume this rate function has a specific scaling dependence on γ. Let m η (y | x) be the placement measure when the reactants at positions x ∈ X ( ) react and generate products at positions y ∈ Y ( ) . We assume this placement measure depends on the displacement range parameter η.
With the preceding definitions, we can formulate a weak representation for the dynamics of ν ζ t and the concentration (i.e. number density) fields for each species, ν ζ,j t , see [IMS22] for the resulting equations. In this work we only require a weak representation for the time evolution of the scaled empirical measures, i.e. the molar concentration field for species j, µ ζ,j t = 1 γ ν ζ,j t with j = 1, · · · , J. Note, with this definition t } n∈N+ as a countable collection of standard independent Brownian motions in R d for species j, j = 1, 2, · · · , J. We can write the marginal distribution (molar concentration) of species j as which satisfies the coupled system (4.1) In this formulation, one important fact is that for fixed ζ, γ 1, µ ζ,j s− is assumed finite, see Assumption 4.9 in the next section, which provides exchangeability of the sum and Lebesgue integral. It also implies that the stochastic integrals with respect to Brownian motions in (4.1) are martingales (for a fixed ζ).

4.2.
Assumptions for the Mean Field Limit. In this section we summarize the assumptions we previously used in proving the mean field large population limit [IMS22]. We will assume they hold in the remainder in studying the fluctuation corrections.
Assumption 4.1. We assume that for all 1 ≤ ≤ L, the reaction rate kernel K (x) is uniformly bounded for all x ∈ X ( ) . We denote generic constants that depend on this bound by C(K).
Assumption 4.2. We assume that for any η ≥ 0,L + 1 ≤ ≤ L, y ∈ Y ( ) and x ∈ X ( ) , the placement density m η (y | x) is a bounded probability density, i.e. Y ( ) m η (y | x) dy = 1. We want to allow for placement densities involving delta-functions. To do so in a mathematically rigorous way we introduced the displacement (i.e. smoothing) parameter η, through which we can define a corresponding mollifier in a standard way, as given by Definition (4.1) below. This is needed for (4.1) to be well-defined, since expressions like {θ 2 ≤ m η y | P ( ) (ν ζ s− , i) } are nonsensical when η = 0 and the placement density is a Dirac delta function.
Definition 4.1. For x ∈ R d , let G(x) denote a standard positive mollifier and G η (x) = η −d G(x/η). That is, G(x) is a smooth function on R d satisfying the following four requirements (1) G(x) ≥ 0, (2) G(x) is compactly supported in B(0, 1), the unit ball in R d , is the Dirac delta function and the limit is taken in the space of Schwartz distributions.
The allowable forms of the placement density for each possible reaction are given by Assumptions 4.3-4.6: Assumption 4.3. If R is a first order reaction of the form S i → S j , we assume that the placement density m η (y | x) takes the mollified form of m η (y | x) = G η (y − x), with the distributional limit as η → 0 given by This describes that the newly created S j particle is placed at the position of the reactant S i particle.
Assumption 4.4. If R is a second order reaction of the form S i + S k → S j , we assume that the binding placement density m (z | x, y) takes the mollified form of with the distributional limit as η → 0 given by where I is a fixed finite integer and i p i = 1. This describes that the creation of particle S j is always on the segment connecting the reactant S i and reactant S k particles, but allows some random choice of position. A special case would be I = 2, p i = 1 2 , α 1 = 0 and α 2 = 1, which corresponds to placing the particle randomly at the position of one of the two reactants. One common choice is taking I = 1, p 1 = 1 and choosing α 1 to be the diffusion weighted center of mass [IZ18].
Assumption 4.5. If R is a second order reaction of the form S i + S k → S j + S r , we assume that the placement density m (z, w | x, y) takes the mollified form of with the distributional limit as η → 0 given by This describes that newly created product S j and S r particles are always at the positions of the reactant S i and S k particles. p is typically either 0 or 1, depending on the underlying physics of the reaction.
Assumption 4.6. If R is a first order reaction of the form S i → S j + S k , we assume the unbinding displacement density is in the mollified form of with the distributional limit as η → 0 given by with i p i = 1. Here we assume the relative separation of the product S j and S k particles, |x − y|, is sampled from the probability density ρ(|x − y|). Their (weighted) center of mass is sampled from the density encoded by the sum of δ functions. Such forms are common for detailed balance preserving reversible bimolecular reactions [IZ18].
We further assume some regularity of the separation placement density, ρ(r), introduced in Assumption 4.6: Assumption 4.7. For Assumption 4.2 to be true, we'll need that the probability density ρ is normalized, i.e. In the remainder we abuse notation, and write ρ ∈ L 1 (R d ) to mean ρ(|·|) ∈ L 1 (R d ).
Finally, to study the large-population limit of the population density measures, we must specify how the reaction kernels depend on the scaling parameter (i.e. system size parameter) γ. Motivated by the classical spatially homogeneous reaction network large-population limit [DK15], we choose Assumption 4.8. The reaction kernel is assumed to have the explicit γ dependence that for any x ∈ X ( ) , 1 ≤ ≤ L.
Recall that |α ( ) | represents the number of reactant particles needed for the -th reaction. As we assume |α| ≤ 2, we obtain three scalings for the three allowable reaction orders: • |α ( ) | = 0 corresponds to a pure birth reaction. By Assumption 4.8, the scaling is γ; i.e. a larger system size implies more births. In a well-mixed model this would imply that as γ and the initial number of molecules are increased, we maintain a fixed rate with units of molar concentration per time for the birth reaction to occur.
• |α ( ) | = 1 corresponds to a unimolecular reaction. By Assumption 4.8, there's no rescaling as it's linear. We assume the rates of first order reactions are internal processes to particles, and as such independent of the system size. • |α ( ) | = 2 corresponds to a bimolecular reaction. By Assumption 4.8, the scaling of reaction kernel is γ −1 . As the system size increases it is harder for two individual reactant particles to encounter each other and react. See [IMS22] for a more detailed discussion on how such scalings could arise.
Finally, we made two assumptions on the molar concentration fields, µ j t , in [IMS22], requiring boundedness of the amount of particles in the system and convergence of the initial conditions in the large-population limit.
Assumption 4.9. We assume that the total (molar) population concentration satisfies J j=1 1, µ ζ,j t ≤ C(µ) for all t < ∞, i.e. are uniformly bounded in time by some constant C(µ). In the remainder we abuse notation and also denote generic constants that depend on this bound by C(µ).
Here L j := D j ∆ x denotes the diffusion operator associated with a particle of type j.

Main result: fluctuation theorem
In this work, we aim to prove a central limit theorem result for the particle-based stochastic reaction diffusion (PBSRD) model. We define the fluctuation process by In Theorem 5.6 we show that the signed measure-valued process T ]} as ζ → 0 in the appropriate space. As is commonly found for problems involving fluctuation analysis of interacting particle systems, weighted Sobolev spaces are needed to control behavior of the fluctuation process as x → ∞. Such spaces have been previously used to study central limit theorems of mean field systems in a variety of studies, include [FM97], [KX04] and [SSG14]. The weighted spaces introduced in [ are what is needed for our analysis. We subsequently denote by W −Γ,a = W −Γ,a (R d ) the dual space to W Γ,a 0 . In Section 6 we review the precise definitions and key properties of such spaces.
In order to prove the fluctuations theorem, we make a few additional assumptions: Assumption 5.1. In the remainder we will make use of the non-negative integer parameters {D, Γ, Γ 1 , a, b} in specifying weighted Sobolev spaces. These parameters are chosen to satisfy the following constraints: Assumption 5.2. We assume the large population limit γ → ∞ is taken such that √ γη → 0. As such η → 0 simultaneously and we write ζ → 0 for the dual limit.
Assumption 5.3. For the density ρ, in addition to Assumption 4.7 we assume Assumption 5.4. For Γ = 2D + 2, we assume that K ∈ C Γ (X ( ) ), that max j≤Γ ∂ (j) x K ∞ < ∞, and that K (x) is symmetric with respect to permutations in the ordering of components of x that correspond to the same species (e.g. to interchanging x Note that this assumption means our results do not rigorously encompass the commonly used Doi bimolecular reaction model for an A + B → C reaction, in which an A at x and B at y would react with probability per time λ when separated by ε or less, i.e. K(x, y) = λ1 [0,ε] (|x − y|). Similar to how we have introduced η as a regularization of the δ-functions that arise in the particle placement densities and then incorporated the limit η → 0 in our analysis, we could have introduced a second regularization parameter into the Doi rate function K(x, y). For simplicity of exposition we have chosen to ignore this additional complication, but note that our results would hold for any mollification of such a Doi kernel. We expect these smoothness requirements are an artifact of our method of proof, and not intrinsic to the existence of a limiting fluctuation process correction.
Assumption 5.5. Let Γ = 2D + 2 and a = D. We assume that sup For reasons that will become clearer in Subsection 8.1 we have the following definition.
Definition 5.1. For any -th reaction, particle distribution ν = J j=1 ν j δ Sj ∈ M (P ) and particle fluctuation if the -th reaction is of the form S k + S r → · · · Theorem 5.6. Assume T < ∞. Given Assumptions 4.1-4.10 as well as Assumptions 5.1-5.5, the sequence . For any subsequence of this sequence, there is a further subsubsequence that converges in law to the process is a distribution-valued mean-zero Gaussian martingale with marginal variance covariance structure for 0 ≤ s, t ≤ T and f, g ∈ W 2+Γ,a 0 (R d ) Finally, the limiting stochastic evolution equation has a unique solution in W −(2+Γ),a ⊗J for t ∈ [0, T ], and thus the limit accumulation point

On the appropriate Sobolev space
We study convergence in Sobolev spaces, see for example [A78] for a general exposition. The weighted Sobolev spaces introduced in [FM97] are designed to control growth as x → ∞, and are just what is needed for the problem at hand. In this section, we recall their definitions and main properties.
For any integers Γ, a ∈ N, consider the space of real valued functions f with partial derivatives up to order Γ which satisfy Define the space W Γ,a 0 as the closure of functions of class C ∞ 0 in the norm defined above, where C ∞ 0 is the space of all functions in C ∞ with compact support. W Γ,a 0 is a Hilbert space (see Theorem 3.5 and Remark 3.33 in [A78] and also [FM97] for the weighted version) and has the inner product When Γ, a = 0, we write f, g 0 = f, g . W −Γ,a denotes the dual space of W Γ,a 0 that is equipped with the norm Let C Γ,a denote the space of continuous functions f that have continuous partial derivatives up to order Γ such that We refer the interested to [FM97] for details on this class of weighted spaces.
A few properties that we will use in this paper include the embeddings that A key property we will subsequently make use of is that the latter is Hilbert-Schmidt, and implies that the embedding is also Hilbert-Schmidt.

Numerical results
In this section, we numerically compare the mean field with fluctuation corrections SPIDE system to the underlying PBSRD model for a 1D (d = 1) model of the three species reversible reaction. The model we consider corresponds to that presented in Section 2, with the specific choice of reaction kernels and placement densities summarized in Section 7.1.
We start by describing the model problem and discretization schemes for the SPIDE and PBSRD models in detail before giving numerical results. Note that our notation in this section differs slightly from the more abstract notation used previously; for example, we parametrize the particle types by uppercase letters A, B, C instead of natural numbers S 1 , S 2 , S 3 .
We demonstrate that for the total molar mass (i.e., the integral or L 1 -norm of the molar concentration field) of the type C particles, the fluctuation process gives an increasingly accurately approximation of the PBSRD variance as γ increases. For γ = 8000, the largest value of γ that we consider, the variances for the two models agree to statistical error. Since the concentration field is a Gaussian process, the distribution of the molar mass is Gaussian at each time t. We observe empirically that the PBSRD molar mass is also approximately Gaussian for large γ, and so, given the close agreement of the means and variances, we also see close agreement between the entire statistical distribution for the molar mass for large enough γ.
Throughout the following section, we use the terms PBSRD and particle model interchangeably. We refer to the jump process discretization of the particle model that we actually simulate as the CRDME (see discretization section). Finally, we subsequently refer to the solution of the combined mean field with fluctuation correction SPIDE model as the fluctuation process.
7.1. Description of model problem. We restrict the reaction system to the periodic domain Ω = (0, 2π). The function λK(x, y) denotes the forward reaction kernel, i.e., the probability per unit time that a forward reaction A + B → C occurs given an A at x and a B at y, where K(x, y) is a normalized function giving the spatial distribution and λ > 0 is the total reaction rate. Here we use the Gaussian with ε a parameter determining the kernel width, |x − y| the periodic distance between x, y ∈ Ω, i.e.
and Z a normalization constant, For the placement density m(z|x, y), the probability of placing a C at z given a reaction between an A at x and B at y, we use the combination of δ-functions meaning that in the event of an A + B → C reaction the product C is placed at the location of the A or the B each with probability 1 2 . Reversibility (equivalently detailed balance) implies that the reverse reaction in which a C at z unbinds into an A at x and a B at y occurs with rate µK(x, y)m(z|x, y), where µ > 0 is a constant [ZI22]. Integrating over x and y, we find that the probability per time of an unbinding reaction for a C at z is given by Similarly, we can derive that given an unbinding reaction of a C at z, either the A or B is placed at z with probability 1 2 and the other molecule is placed with the Gaussian distribution K(x, z) (respectively K(z, y)) centered at z.
In the results that will follow, we have used the parameters D A = 1, D B = .5, and D C = .1. We also set λ = 1, µ = .05, and ε = 2 −7 . We will also need an initial condition, which we set to be proportional to We use this initial condition for the PIDEs satisfied by the mean-field concentration fields, and to generate initial conditions for the particle model (see the next section). The initial conditions for the SPIDEs that correspond to the fluctuation corrections, a scaled difference between the particle and mean-field initial conditions, were therefore zero.
7.2. Discretization of particle model. We numerically study the particle model using the convergent reaction-diffusion master equation (CRDME), a convergent spatial discretization of the forward Kolmogorov equation associated with the PBSRD model [I13,IZ18]. The CRDME corresponds to the forward equation for a system of continuous-time jump processes on a mesh, and can therefore be simulated using optimized versions of the stochastic simulation algorithm (SSA), also known as Gillespie's method or Kinetic Monte Carlo [G76]. The PBSRD particles' Brownian motions are approximated by continuous-time random walks on the grid, and their reaction interactions by jump processes that depend on the relative positions of reactants on the mesh. For these simulations we use a uniform mesh with nodes N . We denote the compartments, or voxels, that particles hop between by . As we require integer numbers of particles for each simulation, we converted from the molar concentration field at x i to the number of particles, a i , in voxel V i as This choice ensured the average number of particles of species A in the ith voxel is γhA(x i , 0), consistent with the mean-field model's initial condition. The initial number of particles of species B and C were chosen similarly.
The diffusion rate for each species s ∈ {A, B, C} from one voxel V i to either of its neighbors is given by Ds h 2 , where h is the distance between the centers of neighboring voxels. Pair reaction rates between a molecule of type A in voxel V i and a molecule of type B in voxel V j are computed as λK γ (x i , y j ) = K(xi,yj ) γ . The resulting C is placed in voxel V i or V j with probability 1 2 each. When a C at location V k unbinds, an A or B is placed in V k with probability 1 2 . The location of the other product particle is determined by the over the voxel sites. The total rate of unbinding in the CRDME model then is given by the discrete marginal integral of the reaction kernel, which for our chosen parameters is essentially µ (i.e. to numerical precision). Specifying the rates with the midpoint rule quadrature scheme as we have just described preserves detailed balance and also O(h 2 ) convergence for statistics as compared to the version of the CRDME presented in [IZ18].
The CRDME is simulated as a continuous-time jump process using a variant of the direct Gillespie stochastic simulation algorithm (SSA) in which we view the simulation mesh as the top level (also referred to as the fine grid) of a nested hierarchy of uniform meshes so that |Ω i | = 2 i . We keep track of the total reaction rate for each subvolume on each mesh throughout the simulation, which allows us to identify the next fine grid subvolume where a reaction occurs in logarithmic time using a binary search. When a reaction or a diffusion event occurs, we update the total reaction rate for each subvolume on each mesh that contains a fine grid subvolume involved in the reaction. The complexity for these updates is also logarithmic, although the total cost of rate updates after a B molecule reacts or diffuses will be O(N log N ) (with a small constant) due to the cost of updating the reaction rate for all A in the range of the reaction kernel.
In our CRDME simulations, we chose I = 10 so that the simulations used N = 2 10 voxels. This choice was informed by looking at successive differences of statistics for simulations at different mesh resolutions to check for convergence, and also by the requirement for the mesh spacing h = 2π N to be below the kernel spacing ε in order to achieve asymptotic O(h 2 ) convergence rates. 7.3. Fluctuation Process Model and Discretization. The mean field model for the reversible reaction is given by the system of PIDEs where A, B, and C are the mean field molar concentrations for the corresponding particle types, and we define the integral operator This model has been previously explored and compared to CRDME simulations of the PBSRD model in [IMS21]. We denote byĀ,B, andC the fluctuation corrections that correspond to the limiting process in Theorem 5.6. By augmenting the mean field solution with the fluctuation corrections, we obtain the fluctuation processes which we expect to provide an improved approximation of the particle model for finite γ compared to the mean field model. We note that the fluctuation processes are Gaussian processes. From the rigorous weak MVSP equations (5.2) and (5.3) we obtain the following (formal) strong form representation the fluctuation corrections should satisfy Here the noise terms, ξ A t , ξ B t , ξ C t , are mean zero Gaussian processes with covariance structure with the remaining covariances following by symmetry. Note that in the preceding formula we have suppressed writing the spatial dependence of the mean field solutions and test functions, and we have denoted by ., . the L 2 inner-product on Ω.
To compute the SPIDE solution first requires the mean field solution. We use Fourier spectral discretizations for both problems. For the mean field solver, we use a total of 2 9 basis functions {e inx } 2 9 −1 n=0 to represent the concentration fields for each of the three species. We use fast Fourier transforms provided by Scipy to convert between Fourier representations and values at collocation points and compute all integrals using the midpoint rule on intervals centered at collocation points.
Solving for the fluctuation process is more computationally expensive because of the covariance matrix computation and factorization, and the large number of trials required to resolve statistics to sufficient accuracy for convergence testing. We therefore used 61 basis functions {sin(nx), cos(nx), 1} 30 n=1 for approximating the fluctuation corrections. Quadratures were computed at the original 2 9 collocation points for integrals involving the mean field solution using the midpoint rule.
For the time discretization we used a one step IMEX Euler method for both the mean field and fluctuation processes, treating the reaction and noise terms explicitly and the diffusion terms implicitly. In both solvers we used a time step of .001. 7.4. Numerical results. In this section we present numerical results for the reversible A + B C reaction using the CRDME approximation to the particle model, and Fourier spectral approximation to the fluctuation process SPIDEs, as previously described. We examine several statistics that are not available from the (deterministic) mean-field approximation, to illustrate the additional information gained by having the fluctuation corrections. To generate the empirical variances and distributions, we used 280,000 simulations for the CRDME and 1,600,000 simulations for the fluctuation process for each individual set of parameter values.
We first compare the variance of the ( √ γ-scaled) molar mass within the interval for the fluctuation process at various times t. Denote by V SP IDE (t) the variance of the total (scaled) molar mass of C molecules from the fluctuation process, i.e.
We denote by V γ CRDM E (t) the raw (unscaled) variance of the CRDME molar mass for a given γ at time t, i.e.
where C γ CRDM E (x, t) is the molar concentration field process in the CRDME. Then γV γ CRDM E (t) denotes the rescaled variance of the CRDME molar mass, which is directly comparable to V SP IDE (t). That is, we expect γV γ CRDM E (t) → V SP IDE (t) as γ → ∞. Directly comparing the unscaled variances would be more difficult as we expect the unscaled CRDME variance V γ CRDM E (t) → 0 and unscaled SPIDE variance V SP IDE /γ → 0 as γ → 0.
The top panels in Figure 1 show the rescaled CRDME variances over the whole time interval over the entire interval [0, 1] the difference between the variances decreases as γ increases. In the right bottom panel we observe a noisy but roughly consistent rate of convergence as we increase γ. The orange line represents an O(γ − 3 4 ) convergence rate, which is similar to the best fit rate of convergence for the first five values of γ of γ −.71 .
As a final demonstration of the accuracy of the SPIDE solution, in Figure 2 we show the entire distribution of the CRDME molar mass, i.e. the distribution of compared to the distribution of the fluctuation process SPIDE molar mass, i.e. the distribution of In the figure we have normalized the histograms such that both integrate to one, and use a fixed histogram bin-width corresponding to the molar concentration of one-particle for the CRDME curves. We see that, as expected, the SPIDE molar mass has an approximately Gaussian distribution. For large γ, the SPIDE and CRDME distributions have very little visible difference. In our tests, the SPIDE solver is about 100 times faster than the CRDME simulation for large enough γ; it took approximately 21 hours on 28 threads for 1,600,000 SPIDE simulations, or 588 compute hours, for an approximate time of 1.3s/simulation. On the other hand, the γ = 8000 CRDME simulations took around 132s/trial. Hence, the SPIDE provides an intermediate fidelity model, more accurate and more In this section we prove Theorem 5.6. We begin in the next section by formulating an evolution equation for the fluctuation process, Ξ ζ,j t . In Section 8.2 we prove a lemma giving moment bounds for the pre-limit measures, µ ζ,j t , and the post-limit mean-field measures,μ ζ,j t . We then prove tightness for the fluctuation process in Section 8.3, identify the limiting equation it satisfies in Section 8.4, and finally prove uniqueness of the solution to this equation in Section 8.5. 8.1. Preliminary calculations for the prelimit fluctuation process. From Eqs. (4.1) and (4.2), we directly obtain Note that in the last equality of Eq. (8.1), we are using the notation as the martingale part with quadratic variation , denotes the discrete part, i.e. Poisson process, of the martingale.
The quadratic variation of C ζ,j t (f ) is and the quadratic variation of D ζ,j t (f ) is In particular, we have the following cases. If the -th reaction is a first order reaction, If the -th reaction is a second order reaction, When ζ → 0, we would expect µ ζ,j s →μ j s for all j = 1, · · · , J. These considerations led to the notation used in Definition 5.1.
For the remainder of the proof of Theorem 5.6, and without loss of generality, we assume thatL = 0. The case whenL > 0 follows by similar arguments. 8.2. Preliminary moment bounds. The purpose of this section is to obtain moment bounds on µ ζ,j t that are uniform with respect to ζ, and on the limit measureμ j t , for all t ∈ [0, T ] for any T < ∞. We are interested in obtaining a uniform bound with respect to ζ ∈ (0, 1) 2 of the quantity For this purpose, we have the following lemma.
Lemma 8.1. Let k ∈ N be given. Let us assume ρ ∈ L 1 (R d ) and that Assumptions 4.9, 5.3, 5.4, 5.5 hold. Then, we have that there exists a finite constant C(µ, K, ρ, T, D) < ∞ that depends on the upper bounds of the quantities that appear in those assumptions, such that the following moment bound holds Proof of Lemma 8.1. We will only prove the statement for the moments associated with the prelimit measure µ ζ,j · . The statement for the moments of the limit measureμ j · follows by a similar argument. Without loss of generality we assumeL = 0. We start by defining We will first prove that J j=1 h ζ,j,8D T < C for some finite constant C < ∞. Beginning with the equivalent Fock space forward equation representation of (4.1), we show in Appendix A that Now we derive an a-priori bound. With f (x) = |x| 8D in (8.8) we have that We first estimate by analyzing all possible cases.
(1) If the -th reaction is of the form S i → S j , so that α ,i = β ,j = 1, then (2) If the -th reaction is of the form S i → S j + S n , (3) If the -th reaction is of the form S i + S n → S j , (4) If the -th reaction is of the form S i + S n → S j + S r , Thus, we obtain that It remains to bound the second term in (8.9), The analysis of this is done as in Lemma B.2 of [IMS22] using the decomposition in terms of the different types of reactions as was done above. Following the same process, one gets that there exists a constant C = C(µ, K) that depends on t∈[0,T ] J j=1 1, µ ζ,j t and on the upper uniform bound for the reaction Combining the last bound with (8.10) implies the following bound for (8.9) We now prove the uniform bound for θ ζ,j,4D t . Recall that we can assume, without loss of generality, that L = 0. Due to the preceding bound, we can now write (4.1) for f (x) = |x| 4D (8.12) where M ζ,j t (f ) = C ζ,j t (f ) + D ζ,j t (f ) is the martingale given by (8.2). By the Burkoholder-Davis-Gundy inequality with f (x) = |x| 4D , we have that there exists a finite constant C that depends on D j and k, such Similarly, using Doob's maximal inequality and calculations for the different reaction possibilities as before, we obtain that there exists a constant C < ∞ that may depend on t∈[0,T ] J j=1 1, µ ζ,j t , on the upper uniform bound for the reaction kernel K (x) , on Using then the bound (8.11), we obtain that for some constant C < ∞ that may depend on t∈[0,T ] J j=1 1, µ ζ,j t , on the upper uniform bound for the , and on T < ∞, but not on ζ. Taking now supremum over t ∈ [0, T ] and then the expectation operator in (8.12), using (8.14) for the martingale terms and (8.11) to bound the supremum of the Riemann integral terms, we obtain the desired bound 8.3. Tightness. In this section, we aim to rigorously show the relative compactness of the fluctuation process Ξ ζ,1 t , · · · , Ξ ζ,J t , t ∈ [0, T ] ζ∈(0,1) 2 and of the martingale vector Let κ > 0, q ≥ 0 and define the stopping time The estimate in Lemma 8.1 implies that for T < ∞ and q = 2k ≥ 0, we have that . Before we prove our main Lemma 8.9, let us first present the technical Lemmas 8.2 to 8.5.
For a fixed γ, and assuming from Assumption 5.1 that b > a + d/2, Γ 1 ≥ 2D + 1, and a ≥ D, we have that p≥1 t∧θ 2b is a square integrable martingale with the expected value of its quadratic variation bounded by for a finite constant C < ∞.
Proof. For a fixed γ, where we used the moment bound from Lemma 8.1 and Sobolev embedding since Γ 1 > d/2. For a fixed γ, b > a + d/2, by Eq. (8.6), we have that p≥1 The derivation of the last line used the Sobolev embedding theorem, and that since Γ 1 − (1 + D) > d/2 and b − a > d/2 the embedding W Γ1,a 0 → W 1+D,b 0 is of Hilbert-Schmidt type, so p≥1 ||f p || 2 1+D,b < ∞. This concludes the proof of the lemma.
We now develop several useful bounds involving the function is uniformly bounded in ζ ∈ (0, 1) 2 . The constant C depends on L, C(K), C(µ), κ and the upper bound from the moment bound of Lemma 8.1.
< CT ||f p || 2 D,b < ∞ which follows using the Sobolev embedding theorem as D > d/2, using Assumption 4.8 on the γ-scaling of K γ (x), and using the uniform bound on the moments of µ ζ s from Lemma 8.1. Note, in the preceding estimate the final bound has no γ dependence.
Lemma 8.4. For a fixed γ, and any is a square integrable martingale with quadratic variation bounded by The constant C depends on L, C(K), C(µ), κ, and the upper bound from the moment bound of Lemma 8.1.
Proof. As shown in Eq. (8.17), for a fixed γ, ||Ξ ζ,1 s || −Γ1,a < ∞, so that the quadratic variation (8.20) Lemma 8.5. For Ξ ∈ W −Γ1,a and a ≥ 1, recalling that L j = D j ∆ x , we have that Proof. By the Riesz Representation Theorem for Hilbert spaces, there exists a unique Ψ ∈ W Γ1,a 0 such that Let us denote F (Ξ) = Ψ, and note that S := {Ξ ∈ W −Γ1,a : F (Ξ) ∈ W Γ1+2,a 0 } is dense in W −Γ1,a . We will first focus on such Ξ ∈ S. Then L j F (Ξ) ∈ W Γ1,a 0 and Ξ, Then, by Lemma B.1 we obtain that for some constant C < ∞ By a density argument, the above result holds for any Ξ ∈ W −Γ1,a .
Lemma 8.6. Let {f p } p≥1 be a complete orthonormal system in W Γ1,a 0 (R d ) of class C ∞ 0 (R d ). Under Assumptions 5.1 and 5.3 we have that there exists a finite constant C < ∞ such that Proof. We begin by considering the term (1) If the -th reaction is a unimolecular reaction S j → · · · , (2) If the -th reaction is a bimolecular reaction 2S j → · · · , by expanding out all the allowable reaction types: (1) If the -th reaction is of the form S i → S j , (2) If the -th reaction is of the form S i → 2S j , we set the operator and we have (3) If the -th reaction is of the form S i → S j + S k for j = n, mimicking the previous case we find (4) If the -th reaction is of the form S i + S k → S j , (5) If the -th reaction is of the form S i + S k → 2S j , If the -th reaction is of the form S i + S k → S j + S r for r = j, mimicking the previous case we find In each case studied above, the upper bound is of the form Df p , Ξ ζ s where the operator D is either a local operator as Df p (x) = K (x)f p (x), Df p (x, s) = K (x), µ ζ,i Examining the left hand side of (8.21), for such operators D we see that in all of these cases estimates of the following form hold: Let us prove this for one possible form of the operators D, say for Df p (x, s; ζ, i) = K (x, y)f p (y), µ ζ,i s (dy) . The rest of the cases can being treated similarly 1 . Let us first show that D : W Γ1,a 0 → W Γ1,a 0 and for f p ∈ W Γ1,a 0 the bound Df p 2 Γ1,a ≤ C f p 2 Γ1,a holds, for some unimportant constant C < ∞. Indeed, the following computation holds where we have used the Sobolev embedding theorem, the boundedness of max k≤Γ1 ∂ (k) x K 2 ∞ , and that 2a > d by assumption. This estimate also demonstrates why weighted spaces are needed. Then, as in Lemma 8.5, the Riesz representation theorem for Hilbert spaces allows us to show that there exists a unique Ψ j , Ψ i ∈ W Γ1,a 0 such that (omitting the ζ for notational simplicity) Hence, using Young's inequality, pulling the supremum of ∂ (k) x K outside the integration and omitting the time dependence for now, we shall have (following a similar computation as in (8.23)) 1 Even though we do not show this here, showing the desired bound for the choice Dfp(x) = D i (1−αn) fp(x) leads to the requirements ρ L 1 < ∞ and R d |w| 2a ρ(w)dw < ∞.
where to go from the first inequality to the second we did the same computation as in (8.23). From the preceding estimate we find We conclude that there is a constant C < ∞ that depends on κ, C(K), C(µ), ρ L 1 , R d (1 + |x| 2a )ρ(x)dx, and the moment bound from Lemma 8.1 such that (8.25) The constant C < ∞ depends on κ, C(K), C(µ), ρ L 1 , R d (1 + |x| 2a )ρ(x)dx and the moment bound from Lemma 8.1.
Proof. Let {f p } p≥1 be a complete orthonormal system in W Γ1,a 0 (R d ) of class C ∞ 0 (R d ). By applying Ito's formula on Eq. (8.1), we can obtain We will sum over all p ≥ 1, take the supremum over t ∈ [0, T ] and then take expectations on both sides of Eq. (8.27). Using Parseval's identity, the first line of Eq. (8.27)'s right side becomes For the third line of Eq. (8.27), using Lemma 8.6 we obtain the bound 2E sup For the fourth line, using Lemma 8.7 we obtain the bound 2E sup The fifth line in Eq. (8.27) is a martingale by Lemma 8.4 and thus by applying Jensen's inequality, Doob's martingale inequality, and Young's inequality, similarly to the derivation of Eq. (8.28), it is bounded by The last term in Eq. (8.27) by Lemma 8.3 is bounded by Therefore, for some finite constant C that depends on κ, C(K), C(µ), T , ρ L 1 , b, D, R d |x| 2a ρ(x)dx and on the upper bound of Lemma 8.1 for the moments, we find that As Γ 1 > 1 + d/2 + D and b − a > d/2, the embedding W Γ1,a 0 → W 1+D,b 0 is of Hilbert-Schmidt type, so p≥1 ||f p || 2 1+D,b < ∞. Summing over j = 1, · · · , J, using Parseval's identity, and using Assumption 5.2 to make √ γη sufficiently small, we find for j = 1, · · · , J.
Proof. Notice that the quadratic variation of M ζ,j t∧θ 2b with D > d/2, for some C 1 depending on D j s, C(µ), C(K), T, L, κ but independent of ζ. Above the last inequality is due to the Sobolev embedding theorem since m > d/2. Doob's martingale inequality gives concluding the proof of the lemma.
8.3.3. Control of the increments. Since we showed Lemmas 8.8 and 8.9, to show tightness (see [FM97] pg. 46), we only need to verify the following Aldous conditions.
Lemma 8.10 (Aldous condition for martingale). For every ε 1 , ε 2 > 0, there exists δ > 0 and an integer n 0 such that for every F n t -stopping time τ n < T , j = 1, · · · , J, (8.32) sup Proof. Let {f p } p≥1 be a complete orthonormal system in W Γ1,a 0 (R d ) of class C ∞ 0 (R d ). Then with D > d/2, for some C depending on the D j 's, C(µ), C(K), T , L, and κ, but independent of ζ. Here the last inequality is due to Sobolev embedding theorem since D > d/2.
Thus, by Parseval's identity, we obtain where C is a constant depending on C(K), C(µ), and the D j 's. By the Markov inequality, we have Cδ ≤ ε 2 , j = 1, · · · , J, for δ sufficiently small.
Lemma 8.11 (Aldous condition). For every ε 1 , ε 2 > 0, there exists δ > 0 and an integer n 0 such that for every F n t -stopping time τ n < T , j = 1, · · · , J, (8.34) sup Proof. By the Markov inequality, we have as long as we choose δ sufficiently small, using the following argument to derive the bound on E||Ξ ζ,j . By definition (let us ignore the stopping time θ 2b ζ,κ for the moment), Applying Ito's formula, we obtain In the second line of Eq. (8.35), by Lemma 8.2, the stochastic integral 2 is a martingale and therefore its expectation is zero. By Eq. (8.6), the term C ζ,j The third line of Eq. (8.35), according to Lemma 8.6, is bounded by The fourth line of Eq. (8.35), by Lemma 8.7, is bounded by Therefore, we obtain for some generic finite constant C < ∞ and b > a + d/2. Since Γ ≥ 2D > D + d/2 and b − a > d/2, the embedding W Γ,a 0 → W D,b 0 is of Hilbert-Schmidt type, so p≥1 ||f p || 2 D,b < ∞. Thus, by Parseval's identity, we obtain, and therefore Lemma 8.14. For every f, g ∈ W Γ,a 0 , the process M ζ,1 t , · · · , M ζ,J t defined by Eq. (8.2) is a martingale in Proof. Let f k ∈ W Γ,a 0 , k = 1, · · · , J. Consider the martingale J k=1 σ k M ζ,k t (f k ), for some σ k ∈ R, k = 1, · · · , J. By definition of M ζ,k t (f k ) in Eq. (8.2), we see that For two different species, S i and S j , their Brownian motions are independent, and therefore there are no cross terms in the quadratic variation arising from the integrals with respect to Brownian motions that involve S i and S j (see Chapter 3 Proposition 2.17 in [KS98]). Note, however, species can share the same Poisson "random generators" via reactions in which they both participate.
The second term, again depends linearly on {µ ζ,j 8.5. Uniqueness of Limiting Solution.

A. Forward Equation Fock Space Representation
Our goal in this appendix is to derive (8.8). Our basic approach is to begin with the forward equation Fock space representation for the dynamics, derive the analogous equation in the Fock space representation, and then show it is equivalent under a change of notation to (8.8). This approach allows us to bypass estimates that would be needed to rigorously derive (8.8) directly from (4.1) in the case that the test function can be unbounded (as needed for moment estimates).
To describe the Fock space representation, we first define some alternative notation for representing the numbers of particles of each species and their positions. Let N (t) := (N 1 (t), . . . , N J (t)) = 1, ν 1 t , . . . , 1, ν J denote the vector stochastic process for the number of each species at time t, with n = (n 1 , . . . , n J ) a possible value for this process. We treat N (t) as a multi-index, so that N (t) = J j=1 N j (t) = |N (t)| gives the total number of particles at time t. j n j = |n| is defined analogously. Note, we previously used N to denote a Poisson random measure, but within this appendix N j (t) will always denote the stochastic process for the number of type j molecules at time t.
We denote by Q Nj (t) i (t) = H i (ν j t ) ⊂ R d the stochastic process for the position of the ith particle of species j at time t, with Q Nj (t) (t) = Q Nj (t) 1 (t), . . . , Q Nj (t) Nj (t) (t) ⊂ R dNj (t) the stochastic process for the state vector of species J and Q N (t) (t) = Q N1(t) (t), . . . , Q N J (t) (t) ⊂ R dN (t) the stochastic process for the state vector of all particles. Possible values for each of these are given by q nj i , q nj and q n respectively.
With these definitions we denote by P n (q n , t) the probability density that N (t) = n with Q N (t) (t) = q n , and define the Fock space probability vector P (t) = {P n (q n , t)} n . Note that because particles of the same species are assumed to be identical, P n (q n , t) is symmetric with respect to reorderings of the particle positions of type j within q nj . As in [IMS22], we consider P (t) for t fixed to be an element of an L 2 Fock space, F = L 2 (X), where X = n R d|n| .
Note that Assumption 4.9 implies N (t) is bounded from above, and hence X is finite dimensional. Letting G = {G n } n ∈ F andG = {G n } n ∈ F , we construct an inner product on F by where n! is defined in the usual multi-index sense. The forward equation for P (t) is then given in strong form by (A.1) P (t) = (L + R)P (t), where L denotes a diffusion operator and R a reaction operator. The former is defined by (LP ) n (q n , t) = J j=1 D j ∆ q n j P n (q n , t), where ∆ q n j denotes the Laplacian in q nj . To define the reaction operator we need to introduce notation for adding and removing particles from a state q n , which requires notation for configurations of particles that could have been substrates or products of one reaction. Abusing notation, we let Ω ( ) (n) := Ω ( ) (ν), see Definition 3.6, label the collection of allowable reactant particle indices when there are n particles in the system, and note that the number of elements in this set is given by We define the allowable product index sampling spaceΩ ( ) (n) ⊂ J ( ) bỹ β ( ) = 2, β j = β k = 1, j < k.
Note, this is a slightly different definition than that of Ω ( ) (n), needed due to our choice of normalization for m (y | x). The number of elements inΩ ( ) (n) is Ω ( ) (n) = n! (n − β ( ) )! .
Given the current state of a system, q n , we will write q n i for i ∈ Ω ( ) (n) to label the position of the substrates for the th reaction determined by i. Similarly, for j ∈Ω ( ) (n), q n j will denote one product collection of particles in q n that could have been produced through one occurrence of reaction . To represent the state vector with particles q n i removed from q n we use the notation q n \ q n i . We analogously use the notation q n ∪ x to denote adding the particles in x into the state vector q n . Finally, we define the overall reaction interaction function, encoding both the rate of a reaction and the placement of products, by K γ (y | x) = K γ (x)m η (y | x) , and note that by assumption. With these notations, the reaction operator is given by (RP ) n (q n , t) = − i∈Ω ( ) (n) K γ (q n i )P n (q n , t) + 1 α ( ) ! i∈Ω ( ) (n) where υ ( ) = β ( ) − α ( ) is the net stoichiometry vector for the th reaction.
In the special case of the A + B C reaction, we proved in [IMS22] that P (t) is non-negative, normalized in time, i.e. . We do not show here, but subsequently assume, that these properties still hold for the general class of reaction systems considered in this work under the assumptions of our main result, Theorem 5.6. Recall this includes an assumption that the total number of molecules is uniformly bounded from above, which rules out finite-time blowup.
To derive (8.8), we begin by noting that for a given scalar test function g(x) ∈ C(R d ), we have that

B. Auxiliary lemmas
Lemma B.1. Let a ≥ 1 and J ≥ 0. Then, for any Ψ ∈ W J+2,a 0 , there is constant C < ∞, that may depend on D j and a but not on Ψ, such that L j Ψ, Ψ J,a ≤ C||Ψ|| 2 J,a Proof of Lemma B.1. For notational convenience we set, without any loss of generality, D j = 1. Let us first consider the case that J = 0. By definition of the norm in question and integration by parts we have Integration by parts of the last expression gives Analyzing now the latter integral separately for the regions |x| < 1 and |x| ≥ 1, and using the assumption a ≥ 1, we obtain that there is a constant C < ∞ that depends on a such that R d 4a 2 |x| 4a−2 (1 + |x| 2a ) 2 1 1 + |x| 2a |Ψ(x)| 2 dx ≤ C||Ψ|| 2 0,a .
The latter then readily gives that 0,a , completing the proof for J = 0. The computation for J > 0 follows along exactly the same lines.