A Feynman Path Integral-like Method for Deriving Reaction–Diffusion Equations

This work is devoted to deriving a more accurate reaction–diffusion equation for an A/B binary system by summing over microscopic trajectories. By noting that an originally simple physical trajectory might be much more complicated when the reactions are incorporated, we introduce diffusion–reaction–diffusion (DRD) diagrams, similar to the Feynman diagram, to derive the equation. It is found that when there is no intermolecular interaction between A and B, the newly derived equation is reduced to the classical reaction–diffusion equation. However, when there is intermolecular interaction, the newly derived equation shows that there are coupling terms between the diffusion and the reaction, which will be manifested on the mesoscopic scale. The DRD diagram method can be also applied to derive a more accurate dynamical equation for the description of chemical reactions occurred in polymeric systems, such as polymerizations, since the diffusion and the reaction may couple more deeply than that of small molecules.


Introduction
Finding correct and thermodynamically consistent reaction-diffusion equations [1][2][3][4][5] is important for the study of both small-molecule and macromolecular reacting systems. Previously, the reaction-diffusion equations were obtained by simply adding the contributions from the reaction and the diffusion together [1][2][3][4][5], and they were, thus, assumed to be independent of each other.
However, experiments [6,7] revealed that diffusions and chemical reactions may affect each other, and these coupling effects will cause observable consequences. For example, Wang et al. [6] found that chemical reactions will change the mobility of molecules or affect diffusions by kicking the molecules near the reaction sites. Accordingly, diffusions can also affect chemical reactions. Liu et al.'s recent work [7] showed that because in the polymerization, monomers have to be transported, by passing through a complicated entanglement of the long chain to the right location of the reaction site of a polymer chain in order for the chain growth reaction to actually happen, a wait-and-jump phenomenon was observed in the chain growth polymerization.
Therefore, it is necessary to propose a theoretical framework that can be used to describe the coupling between chemical reactions and diffusions. The common framework used to investigate this coupling is the diffusion-controlled reaction model. For example, based on the formulation of second quantization or the path integral, Doi and Peliti [8][9][10] successively proposed a very general theoretical framework for the dynamics of many-particle system, which can, in principle, be used to derived the dynamical equation for describing the reactiondiffusion coupling. However, as pointed out by Doi himself, the second-quantization method might not be very useful for the usual classical systems. Another method is the Smoluchowski diffusion limited reaction model [11,12], which is mainly based on the Smoluchowski equation of this system. Recently, the reaction-diffusion master equation (RDME) model has also been attracting the attention of researchers for stochastic chemical kinetics [13][14][15][16], where the particles are treated as points and evolve according to master equations.
However, these models studied diffusion-reaction coupling by investigating bimolecular reactions. One may, therefore, be interested to know whether this coupling still exists in monomolecular reactions. Furthermore, if this coupling exists in the monomolecular case, can we simply add the contributions from the reaction and diffusion together to derive dynamic equations? Additionally, to what extent can we neglect the coupling term? Therefore, this work is devoted to re-deriving the reaction-diffusion equation by summing over the contributions from all possible trajectories of reactions and diffusions. Hopefully, we can identify the coupling terms of these two processes in the newly derived reaction-diffusion equation.
Inspired by the Feynman diagram in particle physics, this work will introduce a series of diffusion-reaction-diffusion (DRD) diagrams to derive a more accurate reaction-diffusion equation. The basic idea of the DRD diagram, taking a binary reacting system (A B) as an example, is as follows. Previously, most of us might see a trajectory (or more precisely, a trajectory ensemble) of an A-type molecule moving from r to another point r at ∆τ later as just one physical process denoted as (A; r, t) → (A; r , t + ∆τ), and summing all possible r and ∆τ will lead to the diffusion equation of A. This strategy has been employed by many previous works [17][18][19] to derive field-based dynamical equations based on the particle-based microscopic dynamics. However, this work notices that when A can turn into B and B can turn into A through chemical reactions, in between the trajectory of A → A, other physical processes can also occur. For example, even though the physical process (see a 2 in Figure 1a) involves two chemical reaction processes in between r and r , it can be still seen as the above process (A; r, t) → (A; r , t + ∆τ) by us. Therefore, a simple physical process (A; r, t) → (A; r , t + ∆τ) might have many possible variants. This work is going to consider these possibilities, and we hope to derive a more accurate reaction-diffusion equation based on this idea. In particular, we are specially interested in what conditions the reaction will couple with diffusions according to the newly derived equations, which might also have potential applications for the description of the dynamics in the mesoscopic scale or the dynamics of polymerization. Figure 1. Diffusion-reaction-diffusion (DRD) diagrams of four processes related with the A/B binary reaction-diffusion system. In the diagram, the solid arrow indicates diffusion, and the dotted arrow means chemical reaction. (a,c) The diffusions of A or B component from r to r + ∆r within a given amount of time ∆τ actually consists of a bunch of possible DRD diagrams in the A/B reaction system. The reactions A → B (b) and B → A (d) also consist of numerous DRD diagrams. Note that a given diagram actually corresponds to an ensemble of trajectories.

Diffusion-Reaction-Diffusion Diagram
In this section, we will give a brief introduction about the overall strategy on how to derive the dynamical equation of a binary system undergoing chemical reactions and diffusions. In particular, the diffusion-reaction-diffusion (DRD) diagram, which is similar to the Feynman diagram in particle physics, will be specially introduced.
The overall derivation strategy is simple and direct. First, calculate the probabilities of all possible physical processes or trajectories. Second, evaluate the contributions of these trajectories to the number density variation based on their probabilities. Third, obtain the number density variation (or dynamical equation) of a given component by summing over the contributions from all possible trajectories.
In the present binary system undergoing chemical reactions, there are mainly four physical processes as shown in Figure 1, (a) diffusion of an A molecule from r to r + ∆r within a small time span ∆τ, (b) an A-type molecule turns into a B-type molecule through chemical reaction, (c) diffusion of a B molecule, and (d) a B-type molecule turns into an A-type molecule. Just as pointed out in the introduction section, each of these physical processes is not just a single process but actually represents a group of physical processes.
For example, for the first group physical processes (a), one might observe an A molecule disappearing at r at t = 0 and re-appearing at r + ∆r at t = ∆τ. This phenomenon may correspond to many possible physical processes as shown in Figure 1, where the A molecule can turn into B molecule through chemical reaction. Take the a 2 DRD diagram for example ( Figure 1a); the A molecule first diffuses into another position and turns into B at that position, which will further diffuse into somewhere, turn back into A, and finally diffuse to the position r + ∆r with the whole process taking a time of ∆τ. Obviously, in between A → A diffusion process, one can add as many (A B) + (B A) chemical reaction pairs as possible as long as the whole process takes the time of ∆τ and the final position is r + ∆r.
More information about these DRD diagrams can be seen in Figure 2.  Figure 1b, which indicates that in this process, an A molecule will firstly diffuse from r to r + ∆r within ∆τ 1 , then turn into the B molecule through chemical reaction (green dot) at the position r + ∆r , and finally B diffuses to the position r + ∆r within ∆τ − ∆τ 1 . The probability of this process is denoted as P b 1 (r, r + ∆r; ∆τ), which is obtained by path integrating over all possible trajectories r → r + ∆r and r + ∆r → r + ∆r and by integrating over the intermediate variables ∆r and ∆τ 1 .
If the probabilities P x i (r, r + ∆r; ∆τ) (x = {a, b, c, d} and i = 1, 2, 3 . . .) of these DRD diagrams can be calculated, then the probabilities of these four groups of processes can be obtained as P x (r, r + ∆r; ∆τ) = ∑ i P x i (r, r + ∆r; ∆τ). Luckily, the order of P x i in ∆τ is completely determined by how many chemical reactions occur in between the process (e.g., P a 1 ∼ ∆τ 0 and P b 2 ∼ ∆τ 3 ), then most of the DRD diagrams of a higher order can be neglected if ∆τ is small.
With P x i 's, the dynamical equation of A can be derived as follows: where φ A/B are the number densities of A and B, respectively. Note that this equation is similar to the master equation in Doi's work [9]. The equation for B is similar. We shall explain the four terms on the right hand side of the above equation one by one as follow. The first two terms account for the diffusions of A from r to r + ∆r and from r + ∆r to r, respectively. The third term φ B (r + ∆r)P d i (r + ∆r, r; ∆τ) accounts for the density increase due to the chemical reaction B A and it is obviously proportional to number of B molecules found at r + ∆r. The last term is related to the reaction A B. Evaluations of these four terms are quite involved and will be presented in Sections 3-5.

DRD Diagram Rules
If A/B molecules are assumed to obey Langevin dynamics, it can be proved that the probability of a pure diffusion DRD diagram (a 1 and c 1 in Figure 1) can be obtained as follows, by path integrating all possible trajectories connecting r and r + ∆r (see Equation (A7) in Appendix B): where D is the diffusion coefficient. U A (r) is the potential field felt by a single molecule A at r. For example, when there is interaction χφ A φ B between A and B, then U A (r) = χφ 2 B . Note that chemical potential of A can be expressed as µ A = ln φ A + U A + A with µ 0 A = A the inner energy of A. ∆τ is chosen such that it should be larger than the minimum time scale that makes Langevin dynamics valid for A and B molecules. Note that, for the sake of simplicity, we set k B T = 1 throughout this work.
Similarly, for the pure diffusion of B, we have With Equations (2) and (3), we have the following diagram rules (similar to the Feynman rules): (i) The two end points of the DRD diagram can be either at r or r + ∆r.
(ii) Each internal vertex (dotted circle) represents a chemical reaction, either A B or B A, but it should be consistent with its incoming or outgoing molecule. For example, if the incoming molecule is A, then the reaction should be A B (we use the dash arrow to denote the reaction). (iii) A B vertex contributes a factor k + exp{U A (r + ∆r )} with r + ∆r the spatial position of the vertex, while B A vertex contributes a factor k − exp{U B (r + ∆r )} with k + and k − , the reaction rate constants of the forward and the reverse reactions, respectively. (iv) A solid-line arrow represents a pure diffusion process, and it contributes a factor P A/B (r i , r i+1 ; ∆τ i ) as given in Equations (2) and (3), where r i and r i+1 denotes the spatial positions of the two end points of the arrow, and ∆τ i denotes the time span of the diffusion. Note that ∑ i ∆τ i = ∆τ and ∆τ i > 0. (v) Multiplying all the factors of arrows (P A/B ) and vertices (k ± exp{U A/B (r + ∆r )}) and integrating over all possible spatial positions of vertices and time spans of the arrows, the probability of the diagram can be finally evaluated.
Comments about (iii). Note that even though we have not explicitly let the probability of reaction depend on the time span, when being incorporated with the diffusion process, the reaction DRD diagram will be found to depend on the time span (see Equation (6)). The forward and backward reaction rates are related by k + = k 0 e A and k − = k 0 e B with A = µ A0 and B = µ B0 , the inner energies of A and B, respectively.
The reason we choose the factor (the local reaction probability) of the reaction vertex of A B as k + exp{U A (r + ∆r )} is as follows. First, the probability of the forward reaction of A B should locally depend on only the state of the reactant A other than that of the product B. Physically, A's state is totally determined by its inner energy A and the external field felt U A by A, so the probability function of A B should take the form of f ( A , U A ). The expression of k + = k 0 e A obviously indicates f ( A , U A ) = e A g(U A ), and it is natural to conjecture g(U A ) = e U A . Second, as can be seen in Sections 6 and 7, the derived reaction-diffusion equation using this factor is more reasonable than those obtained using other methods, especially for the small density region.

Probability of DRD Diagram
For example, based on the above five rules, we can evaluate the DRD diagram of b 1 ( Figure 2) as follows: Plugging Equations (2) and (3) into the above equation and re-arranging the relevant terms, the above equation can be rewritten as where By Equation (A5) in Appendix A, P b 1 can be evaluated as where only terms up to the second order of ∆τ are kept, and the meaning of ∂ α ∂ β g∆x α ∆x β is referred to the context around Equation (A5). The probability of the d 1 diagram (Figure 1d) can be easily obtained by simply switching A and B and replacing k + with k − in Equation (6),

Probability of More Complex Diagrams
After the b 1 and d 1 diagrams' probabilities are obtained, other diagrams' probabilities can be derived one by one.
For example, a 2 diagram can be obtained by adding a B A reaction and an A-type diffusion to the b 1 diagram (see Figure 1a,b). Accordingly, by the DRD diagram's rule, it can be realized by multiplying b 1 diagram's probabilityP b 1 (r, r + ∆r ; ∆τ ) with k − e U B (r+∆r ) (B A reaction) and P A (r + ∆r , r + ∆r; ∆τ − ∆τ ) (A-type diffusion) and integrating over ∆r and ∆τ . Closely following Equation (4), the probability of a 2 diagram can be written as Even without working out the above integration, one can conclude that the leading order of the a 2 diagram is of the second order of ∆τ.
Similarly, multiplying one more forward reaction and a pure B-type diffusion to a 2 diagram and integrating the corresponding positions and time intervals, one obtains the b 2 diagram's probability. All other diagrams can be derived step by step in this way.

Order, Symmetry and Duality of the Diagram
Here, the order of the diagram is defined as the leading order in ∆τ of the diagram's probability. By the above subsection, one can easily prove that the order of the diagram amounts to the number of reactions or the number of dotted circles in the diagram. Accordingly, a i and c i diagrams are of ∼∆τ 2i−2 , while b i and d i diagrams are of ∼∆τ 2i−1 .
For diagrams a i and c i (Figure 1a,c), reversing the directions of all the sub processes results in the same diagram. For example, we will still obtain the a 2 diagram after reversing the directions of the a 2 diagram. Therefore, a i and c i are symmetric in this respect.
On the other hand, b i is not symmetric, but is dual to d i , i.e., reversing b i will result in d i .
The order, symmetry and duality of the diagram are important in the following derivation of the dynamical equation.

Contribution of DRD Diagram to the Density Variation
We first present the following two general formulae, which can be used to evaluate the contribution of a DRD diagram to the number density variation, as follows: ∆φ Y (r; P X→Y , ∆τ) = d∆rφ X (r + ∆r)P X→Y (r + ∆r, r; ∆τ), where the subscript of P X→Y indicates that the two chemicals at the two ends of the diagram are X and Y. Note that X and Y can be the same species. These two formulae can be applied to all diagrams in Figure 1.

Asymmetric Diagrams' Contributions
The contribution of the b 1 diagram to the density variation of A in Equation (1) can be obtained by multiplying the density φ A (r) and integrating over ∆r, as follows: where the minus sign indicates that the b 1 diagram causes φ A to decrease, The detailed derivation of Equation (11) is referred to in Appendix C.
The b 1 diagram will also contribute to the variation of φ B as follows: where ∆U ≡ U A − U B . As the dual diagram of b 1 , the d 1 diagram's contribution to the concentration variation B will be exactly the same as b 1 's contribution to A (Equation (11)), except that A needs to be replaced with B, and it can be expressed as follows: Similarly, d 1 diagram's contribution to the concentration variation A can be expressed as The net contribution of b 1 and d 1 diagrams to A's variation can be, thus, obtained as Similarly, the net contribution of the b 1 and d 1 diagrams to B's variation is The contributions of other higher-order asymmetric diagrams (b i and d i ) can be calculated accordingly.

Symmetric Diagrams' Contributions
Obviously, the probabilities of a 1 and c 1 diagrams in Figure 1a,c are just P A (Equation (2)) and P B (Equation (3)), respectively. The contribution of a 1 to density variation is as follows: ∆φ A (r; P B a 1 , ∆τ) = d∆rφ A (r + ∆r)P A (r + ∆r, r; ∆τ) where the superscript F in P F a 1 indicates the forward trajectory, while B indicates the backward trajectory. It can be seen that even though the leading-order contribution is of ∆τ 0 , the forward and backward zeroth-order contributions will exactly cancel each other, leaving terms of the first, second and other higher-order terms. These two equations are easily derived by applying Equation (A3).
By adding up Equations (17) and (18), we obtain the total contribution of the forward and backward trajectories of the a 1 diagram (Figure 1a) as where µ A = U A + A + ln φ A with A being independent of r. Detailed derivation of the first-order term is referred to in Appendix D. Note that the first-order term will obviously reproduce the Smoluchowski equation [20,21] and when U A = 0, it is reduced to the familiar expression D∆τ∇ 2 φ A encoding Fick's law. Similarly, the expressions for ∆φ B (r; P F c 1 , ∆τ), ∆φ B (r; P B c 1 , ∆τ) and ∆φ B (r; P c 1 , ∆τ) can be also obtained as The contributions of other higher-order symmetric diagrams (a i and c i ) can be calculated accordingly.

Summary of DRD Diagram in a Table
Here we summarize the properties of the DRD diagrams of Figure 1 in Table 1.

Dynamical Equation of Chemical Reactions Coupling with Diffusions
With the results of the above section, dynamical equations can now be derived as follows: where µ A/B = A/B + U A/B + ln φ A/B , ∆U = U A − U B , and ∆τ 2 terms and other higherorder terms are neglected. Note that we set k B T = 1 through this work for simplicity. Similarly, for B, we have The above two equations should be the most accurate equations to date that can be derived to describe the dynamics of chemical reaction coupling with diffusions.
First, we would like to analyze the origins of the terms on the right-hand side of the above equations. (i) Even though there are no gradient terms in the leading-order term of the diffusion diagram (e.g., Equation (17) or Figure 1a), the leading-order term of the diffusion diagram and that of its corresponding backward diagram will always cancel each other, leaving terms that have gradient terms. (ii) At the same time, the leading-order term of the reaction diagram and that of its dual diagram will not cancel each other. So the main reaction term in the final dynamical equation has no gradient term. (iii) It is interesting to find that there is a first-order pure diffusion term proportional to D 2 ∆τ, which is actually produced by diagrams a 1 and c 1 , while there is no such first-order pure reaction term proportional to k 2 0 ∆τ, which can be only produced in the leading-order term in diagrams a 2 and c 2 but will be canceled by their backward diagrams. (iv) The left first-order terms in the dynamical equations are the true coupling terms between the diffusion and the reaction that are proportional to Dk 0 ∆τ, which are produced by the diagrams b 1 and d 1 . Therefore, we have shown that even for monomolecular reactions, there will be reaction-diffusion coupling terms. The reason for this coupling might be that the reaction rate depends on the state of the reactant (e.g., ∝ e U A (r) ). However, on the mesoscopic scale, this dependence should be replaced with a more precise one that is obtained by averaging over a small volume, and it is like The terms a(∇U A ) 2 D∆τ and b∇ 2 U A D∆τ will finally lead to the reaction-diffusion coupling.
There are two circumstances for which the first-order terms of ∆τ cannot be neglected. Here, ∆τ should be at least several τ f 's, with τ f being the mean free time of A/B molecules such that A/B molecules obey Langevin dynamics during ∆τ.
First, in the mesoscopic world, the time scale of interest will be close to ∆τ, then the first-order terms cannot be neglected, and it would be interesting to investigate the effect of the coupling terms ∼ Dk 0 ∆τ on the dynamics. We anticipate that these two equations will be useful in the future when researchers are simulating the chemical reactions coupling with diffusions inside the cell, which can be seen as a mesoscopic system.
Second, during the polymerization, polymers or reacting monomers have to diffuse a longer distance than the small-molecule reactions in order to initiate another chemical reaction, which will, in turn, make ∆τ larger such that the first-order terms are not negligible, even when viewed from the macroscopic world.
On the contrary, for a macroscopic system composed of small molecules, ∆τ is extremely small compared to the time scale of interest. Thus, the first-order terms can be neglected, and the above two equations are reduced to Now, the diffusion and the reaction are actually decoupled since there are no terms ∼ Dk 0 .
Here are some comments about these two equations: (i) When there are no chemical reactions, the above equations are actually Smoluchowski equations [20], and these two equations cannot guarantee the local incompressibility of the total concentration since we have not considered this constraint in this work; this will be specially discussed in our future work. (ii) When the intermolecular interaction is neglected or U A/B = 0, the above two equations can be further reduced to the traditional reaction-diffusion equations [1][2][3][4][5], which validates the derivation of this work in some sense. Here, we present the traditional reaction-diffusion equations as follows: In particular, the diffusion terms . At the same time, the chemical reaction terms are also reduced to the familiar ones, i.e., is strongly recommended by Carati and Lefever [22], behaves well when φ A or φ B or both are close to zero. For example, if φ A = 10 −1000000 and φ B = 10 −1000 , the reaction rate given by our dynamical equation is obviously very close to zero, but the reaction rate [23,24] given by the term, like k 0 (µ B − µ A ), can be extremely big. The term k 0 (µ B − µ A ) has been used by some previous works [23,24].

Comparison with Previous Work
To show the difference between our newly derived dynamical equations (Equations (23) and (24)) and the other theory, we give a simulation example. Consider an A/B binary system. The two components can transform into each other through chemical reactions A B, and the interaction between the two components can be expressed as with χ the dimensionless interaction parameter and b A/B the interfacial energy coefficients and U A/B = δF int /δφ A/B .
By numerically solving the corresponding dynamical equations with the reactiondiffusion (RD) coupling terms (Equations (23) and (24)), the equations without RD terms (Equations (25) and (26)) and the equations proposed by previous works [23,24], respectively, we see the difference. First, we compared the newly derived dynamic equations in this work with the k 0 (µ A − µ B )-type equations proposed in previous works [23,24] and the difference is significant. The k 0 (µ A − µ B )-type equation is much slower than the newly derived dynamic equations (see Figures 3 and 4) to reach the equilibrium state. We found that k 0 (µ A − µ B )-type equations are numerically instable when the density is extremely small, as mentioned in the preceding section. There is no such problem for our theory. Second, we examined the effect of RD coupling terms (or ∆τ terms) in the dynamic equations. By comparing the light scattering results of with (Equations (23) and (24)) and without (Equations (25) and (26)) coupling terms, we can easily find that the dynamics with coupling terms will reach the equilibrium state slightly faster (see Figure 5). In our simulation, we set ∆τ = 0.05, which is not a big number, but the difference is apparent. Therefore, on the mesoscopic scale, the coupling effect between the reaction and the diffusion might be observable.  (23) and (24) of this work at (a) t = 0τ, (b) t = 5τ, (c) t = 20τ, (d) t = 100τ and (e) t = 1000τ, while the bottom snapshots are obtained by the dynamical equations of previous works [23,24] at (f) t = 0τ, (g) t = 500τ, (h) t = 1000τ, (i) t = 5000τ and (j) t = 10000τ. The numerical simulations are performed on a 256 × 256 lattice.

Discussion
Compared with previous works, the newly proposed method and the derived reactiondiffusion equations of this work have the following important features.
(i) The derived dynamical equations in this work contain higher-order terms in ∆τ, which show that reactions and diffusions are no longer independent but are coupled. (ii) The newly derived reaction-diffusion equations are more suitable for the description of reaction-diffusion processes on the mesoscopic scale because the time scale of interest on the mesoscopic scale is comparable to the time scale ∆τ. It is also suitable for the description of polymerizations since polymers diffuse much more slowly in polymerization than the small-molecule diffusion, which might make the reactiondiffusion coupling effect observable on the macroscopic scale. Nevertheless, the above anticipations still need to be validated by more simulations or experiments in the future. (iii) The newly proposed method helps us derive a more consistent reaction-diffusion equation directly from the microscopic dynamics. The reaction term in the newly derived equation is more reasonable than the previous one, especially for the extremely small density region. (iv) Physical origins of every term in the derived dynamical equations can be clearly seen in the DRD diagrams. For example, one can clearly see that the reason for the diffusion terms containing no pure local terms but gradients is because even though the leading-order term of a pure diffusion diagram (e.g., a 1 in Figure 1a) is purely local, it will be exactly canceled by the leading-order term of its backward diagram, leaving next-and higher-order terms, which are all non-pure local terms. On the contrary, the reaction term can be purely local, which is because the leading order term of the reaction diagram cannot be canceled by that of its dual diagram. Actually, the above results are also in accord with Curie's principle about the symmetry of cause and effect.
However, improvements can be also made on the present work in the future.
(i) In the derivation, the definition of density is actually not strict and accurate. A more accurate definition of the density should depend on the way of how we define the density or depend on the box size which we use to define the local density. In our future work, we will try to introduce the more accurate definition of density that has been employed by other works [17,18]. (ii) In this work, the binary reversible reaction is actually very special in the sense that it involves no explicit intermolecular collisions, which is normally impossible for other reactions. In our next work, we will consider more complicated reactions involving more than two types of reacting molecules and explicitly consider intermolecular collisions in theory. Note that to include collisions, the new definition of the density mentioned in (i) will play an important role. (iii) In this work, we have not considered the local incompressibility of the total density.
In our future work, we will specially discuss how the local incompressibility spontaneously arises from the microscopic world without using any penalty energy to impose the incompressibility. (iv) In the future, particle-based simulations are still required to be done to validate the higher-order terms in the reaction-diffusion equations derived in this work.

Conclusions
In this work, we proposed a Feynman path integral-like method to derive the reactiondiffusion equations that are suitable for the description of the reaction-diffusion dynamics on the mesoscopic scale, where, different from those on the macroscopic scale, reactions and diffusions might deeply couple. By noting that a previous simple physical process might consist of a series of diffusion-reaction-diffusion trajectories or diagrams, the new method was able to derive the reaction-diffusion equation for the A/B reacting system of very high accuracy. The higher-order terms in the newly derived equations show that reactions and diffusions might couple, even for monomolecular reactions, and become manifested on the mesoscopic scale. The newly derived equations are reduced to the traditional reaction-diffusion equations when the intermolecular interactions and higherorder terms are both neglected. The proposed Feynman path integral-like method can be also extended to polymeric systems to derive more accurate equations for the description of polymerization dynamics.

Appendix B. Probability of a Trajectory
To obtain the probability of a trajectory from r to r + ∆r, it is sufficient to only consider the 1D case. The extension to 3D is straightforward.
If F is the gradient of some potential U or conservative, then the probability P(∆x; ∆τ) can be obtained as, by path integrating P[x(t)] over x(t), P(∆x; ∆τ) = P[x(t)]Dx(t) ]∆t] ≈ e Note that we have changed P[r(t)] to P(r, r + ∆r; ∆τ) since now the probability only depends on the two ends and the time span ∆τ. The above equation evaluates the probability of a particle diffusing from r to r + ∆r within ∆τ (it is assumed that there is already a particle in the position r), while the probability of observing a particle diffusing from r to r + ∆r within ∆τ should be multiplied by the number density φ (or e ln φ ) at r, which can be expressed as follows: φ(r)P(r, r + ∆r; ∆τ) = ( This expression is in accord with the result obtained in the fluctuation theory [26], where the exponent − ln φ(r) − 1 2 [U(r) − U(r + ∆r)] is defined as the second entropy of the transition. (11) According to Equation (7), we have ∆φ A (r; P b 1 , ∆τ) = − d∆rφ A (r)P b 1 (r, r + ∆r; ∆τ)