Stochastic approach to Fisher and Kolmogorov, Petrovskii, and Piskunov wave fronts for species with different diffusivities in dilute and concentrated solutions
aa r X i v : . [ n li n . PS ] J u l Stochastic approach to Fisher and Kolmogorov,Petrovskii, and Piskunov wave fronts for species withdifferent diffusivities in dilute and concentratedsolutions
Gabriel Morgado , Bogdan Nowakowski , and Annie Lemarchand* Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka44/52, 01-224 Warsaw, Poland Laboratoire de Physique Théorique de la Matière Condensée, SorbonneUniversité, CNRS UMR 7600, 4 place Jussieu, case courrier 121, 75252Paris CEDEX 05, FranceJuly 22, 2020 * Corresponding author: Annie Lemarchand, E-mail: [email protected]: Wave front, stochastic description, master equation, cross-diffusion
Abstract
A wave front of Fisher and Kolmogorov, Petrovskii, and Piskunov type involvingtwo species A and B with different diffusion coefficients D A and D B is studied usinga master equation approach in dilute and concentrated solutions. Species A and Bare supposed to be engaged in the autocatalytic reaction A+B → D B > D A . In the caseof a concentrated solution, the transition rates associated with cross-diffusion arederived from the corresponding diffusion fluxes. The properties of the wave frontobtained in the dilute case remain valid but are mitigated by cross-diffusion whichreduces the impact of different diffusion coefficients. Introduction
Wave fronts propagating into an unstable state according to the model of Fisher andKolmogorov, Petrovskii, and Piskunov (FKPP) [1, 2] are encountered in many fields [3],in particular biology [4] and ecology [5]. Phenotype selection through the propagation ofthe fittest trait [6] and cultural transmission in neolithic transitions [7] are a few examplesof applications of FKPP fronts. The model introduces a partial differential equation witha logistic growth term and a diffusion term.The effect of non standard diffusion on the speed of FKPP front is currently inves-tigated [8, 9, 10, 11] and we recently considered the propagation of a wave front in aconcentrated solution in which cross-diffusion cannot be neglected [12]. Experimentalevidence of cross-diffusion has been given in systems involving ions, micelles, surface,or polymer reactions and its implication in hydrodynamic instabilities has been demon-strated [13, 14, 15, 16, 17, 18]. In parallel, cross-diffusion is becoming an active field ofresearch in applied mathematics [19, 20, 21, 22, 23, 24].The sensitivity of FKPP fronts to fluctuations has been first numerically observed [25,26]. An interpretation has been then proposed in the framework of a deterministic ap-proach introducing a cutoff in the logistic term [27]. In mesoscopic or microscopic de-scriptions of the invasion front of A particles engaged in the reaction A + B → We consider two chemical species A and B engaged in the reactionA + B k −→ , (1)where k is the rate constant. The diffusion coefficient, D A , of species A may differ fromthe diffusion coefficient, D B , of species B.In a deterministic approach, the reaction-diffusion equations are ∂ t A = D A ∂ x A + kAB (2) ∂ t B = D B ∂ x B − kAB (3)where the concentrations of species A and B are denoted by A and B . The system admitswave front solutions propagating without deformation at constant speed. For sufficientlysteep initial conditions and in particular step functions ( A ( x, t = 0) = C H ( − x ) and B ( x, t = 0) = C H ( x )), where C is constant and H ( x ) is the Heaviside function, theminimum velocity v ∗ = 2 q kC D A (4)is selected [3, 4, 27]. The parameter C = A ( x,
0) + B ( x,
0) is the sum of the initialconcentrations of species A and B. Discrete variables of space, i = x/ ∆ x , and time, s = t/ ∆ t , where ∆ x is the cell length and ∆ t is the time step, are introduced in orderto numerically solve Eqs. (2) and (3) in a wide range of diffusion coefficients D B . Weconsider a system of ℓ = 2000 spatial cells. The initial condition is a step function locatedin the cell i = ℓ / A ( i,
0) = C H ( i − i ) , (5) B ( i,
0) = C H ( i − i ) , (6)3here H ( i ) is the Heaviside function. In order to simulate a moving frame and to coun-terbalance the autocatalytic production of species A in a finite system, the followingprocedure is applied. At the time steps s such that P ℓ i =1 A ( i, s ) > P ℓ i =1 A ( i, A ( ℓ , s ) = 0 and B ( ℓ , s ) = C is created. Hence,the inflection point of the front profile remains close to the initial step of the Heavisidefunction.In small systems with typically hundreds of particles per spatial cell, the deterministicdescription may fail and a stochastic approach is required. We consider the chemicalmaster equation associated with Eq. (1) [32, 33]. The master equation is divided into twoparts ∂ P ( φ ) ∂ t = ∂ P ( φ ) ∂ t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) reaction + ∂ P ( φ ) ∂ t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) diffusion (7)where the first part corresponds to the reactive terms ∂ P ( φ ) ∂ t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) reac = P i k Ω N " ( N A ( i ) − N B ( i ) + 1) P ( { N A ( i ) − , N B ( i ) + 1 } ) − N A ( i ) N B ( i ) P ( φ ) (8)and the second part corresponds to the diffusion terms ∂ P ( φ ) ∂ t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) diff = P i " D A ∆ x ( N A ( i ) + 1) h P ( { N A ( i − − , N A ( i ) + 1 } )+ P ( { N A ( i ) + 1 , N A ( i + 1) − } ) i + D B ∆ x ( N B ( i ) + 1) h P ( { N B ( i − − , N B ( i ) + 1 } )+ P ( { N B ( i ) + 1 , N B ( i + 1) − } ) i − ∆ x (cid:16) D A N A ( i ) + D B N B ( i ) (cid:17) P ( φ ) (9)where φ = { N A ( i ) , N B ( i ) } denotes the default state, Ω , the typical size of the system, N = Ω C , the initial total number of particles in a cell, and N A ( i ) = Ω A ( i ) and N B ( i ) = Ω B ( i )are the numbers of particles A and B in cell i . We consider parameter values leading tothe macroscopic values used in the deterministic approach. The initial condition is givenby ( N A ( i ) = N , N B ( i ) = 0) for 1 ≤ i < ℓ / N A ( i ) = 0 , N B ( i ) = N ) for ℓ / ≤ i ≤ ℓ with N = 100, Ω = 10 ( C = 10). 4 . w a v e f r o n t s p ee d D B D A v ∗ v B ε v d, cut v d, ME Figure 1: Dilute system. Wave front speeds v d, ME , v d, cut , v B ε , and v d = v ∗ versus ratioof diffusion coefficients D B /D A in the dilute case. The values of v d, ME (red circles) arededuced from the direct simulation of the master equation (Eqs. (7-9)) for k = 10, Ω = 10, N = 100, D A = 1, ℓ = 2000, and ∆ x = 0 . v d, cut (blackopen triangles) are deduced from the numerical integration of the deterministic equations(Eqs. (14) and (15)) in the presence of a cutoff ε = 10 − for k = 10, C = 10, D A = 1, ℓ = 2000, ∆ x = 0 . ∆ t = 6 . × − . The values of v B ε (green crosses) are deducedfrom Eq. (16) in which the value B ε has been deduced from the numerical integration ofEqs. (14) and (15). The horizontal line gives the minimum velocity v d = v ∗ (Eq. (4)) ofan FKPP front in the absence of a cutoff.The kinetic Monte Carlo algorithm developed by Gillespie is used to directly simulatethe reaction and diffusion processes and numerically solve the master equation [34]. Theprocedure used in the deterministic approach to evaluate the front speed is straightfor-wardly extended to the fluctuating system. For sufficiently small spatial lengths ∆ x and time steps ∆ t , the numerical solution of thedeterministic equations given in Eqs. (2) and (3) leads to the same propagation speed v d , where the index d stands for dilute, in the entire range of D B /D A values [12]. Thenumber of cells created during 10 time steps once a stationary propagation is reached isused to evaluate the front speed. For the chosen parameter values, we find a propagationspeed obeying v d = v ∗ = 20 with an accuracy of 0 . v d does not depend on the diffusion coefficient D B . The front speed deduced fromthe direct simulation of Eqs. (7-9) is denoted v d, ME where the index d stands for dilute5nd the index ME for master equation. As shown in Fig. 1, the velocity v d, ME is smallerthan the deterministic prediction v ∗ given in Eq. (4).As long as D B remains smaller than or equal to D A , the velocity v d, ME is constant. Themain result of the master equation approach is that the front speed drops as D B increasesabove D A . Typically, for D B /D A = 16, the velocity v d, ME is reduced by 22% with respectto v d = v ∗ . Due to computational costs, larger D B /D A values were not investigated.In the case of identical diffusion coefficients for the two species, the decrease of thefront speed observed in a stochastic description is interpreted in the framework of thecutoff approach introduced by Brunet and Derrida [27]. For D A = D B , the dynamics ofthe system is described by a single equation. When a cutoff ε is introduced in the reactiveterm according to ∂ t A = ∂ x A + kA ( C − A ) H ( A − ε ) , (10)the velocity is given by v ε = v ∗ − π ε ) ! (11)In a particle description, the cutoff is interpreted as the inverse of the total number ofparticles in the reactive interface [28]: ε = ∆ xN W ∗ (12)where the width of the interface is roughly evaluated at [4, 12] W ∗ = 8 s D A kC (13)For the chosen parameter values, the cutoff equals ε = 10 − leading to the corrected speed v ε = 18 .
84. According to Fig. 1, the velocity v d, ME deduced from the master equationfor D A = D B agree with the velocity v ε deduced from the cutoff approach. The resultsare unchanged for D B < D A and Eq. (11) correctly predicts the velocity in a fluctuatingsystem. For D B > D A , Eq. (11) is not valid. Nevertheless, the relevance of the cutoffapproach can be checked by numerically integrating the two following equations ∂ t A = D A ∂ x A + kABH ( A − ε ) (14) ∂ t B = D B ∂ x B − kABH ( A − ε ) (15)6he values of the front speed v d, cut deduced from the numerical integration of Eqs. (14)and (15) are given in Fig. 1 and satisfactorily agree with the results v d, ME of the masterequation, including for large D B /D A values.According to Fig. 2a, the A profile is steeper than the B profile for D B > D A . Themean number of B particles in the leading edge smoothly converges to N . In average, therightmost A particle sees a number of B particles smaller than N . The significant decreaseof the front velocity v d, cut for D B > D A is qualitatively interpreted by the apparent number N ε of B particles seen by the rightmost A particle in the leading edge. The linear analysisof Eqs. (14) and (15) according to the cutoff approach [27] leads to Eq. (11) which does notaccount for the behavior at large D B . A nonlinear analysis would be necessary. Using theperturbative approach that we developed in the case of the deterministic description [4,12], applying the Hamilton-Jacobi technique [35, 36], or deducing the variance h AB i froma Langevin approach [37], we unsuccessfully tried to find an analytical estimation of thefront speed. Instead, we suggest the following empirical expression of the velocity of anFKPP front for two species with different diffusion coefficients v B ε = 2 q kB ε D A − π ε ) ! (16)where B ε denotes the concentration of B species at the abscissa x ε at which the scaledconcentration A ( x ε ) /C is equal to the cutoff ε (see Fig. 2b). The variation of B ε versus D B /D A is numerically evaluated using Eqs. (14) and (15). The result is given in Fig. 3.As shown in Fig. 1, the variation of the front speed v B ε with D B /D A deduced fromEq. (16) slightly underestimates the results v d, cut deduced from the numerical integrationof the deterministic equations (Eqs. (14) and (15)) with a cutoff. We focus on two steady properties of the wave front, the shift between the profiles ofspecies A and B and the width of species A profile [12].For a wave front propagating at speed v and using the coordinate z = x − vt in themoving frame, the shift between the profiles of the two species is defined as the difference A ( z = 0) − B ( z = 0) of concentrations between species A and B at the origin z = 07 N A N B x a) A BB ǫ x x ǫ b) Figure 2: Dilute system. (a) Numbers N A of particles A (red dashed line) and N B ofparticles B (black solid line) versus spatial coordinate x deduced from direct simulationof the master equation (Eqs. (7-9)) using Gillespie method. The snapshot is given attime t = 9 for k = 10, Ω = 10, N = 100, D A = 1, D B = 16, ℓ = 2000, and ∆ x =0 . A of species A (red dashed line) and B of species B (black solid line) versusspatial coordinate x deduced from numerical integration of the deterministic equations(Eqs. (14) and (15)) in the presence of a cutoff ε = 10 − . The snapshot is given at time t = 640 for the same other parameters as in the master equation approach. The verticaldashed line indicates the abscissa x ε for which the scaled A concentration A ( x ε ) /C reachesthe cutoff value. The horizontal line indicates the value B ε of B concentration at theabscissa x ε . 8 . . . . . B ε D B D A Figure 3: Dilute system. The green crosses give the value B ε deduced from the numericalintegration of the deterministic equations (Eqs. (14) and (15)) with a cutoff ε = 10 − versus the ratio of the diffusion coefficients D B /D A . The horizontal line indicates theconcentration C . The parameters are given in the caption of Fig. 1.chosen such that A ( z = 0) = C /
2. The shift is denoted by h d , where the index d standsfor dilute, when the concentrations are solutions of the deterministic equations withoutcutoff given in Eqs. (2) and (3). As shown in Fig. 4, the shift h d significantly varies withthe ratio D B /D A , in particular when D B is larger than D A [12]. The shift vanishes for D A = D B , is positive for D B < D A and negative for D B > D A .The direct simulation of the master equation leads to highly fluctuating profiles. Weuse the following strategy to compute the shift h d, ME . First, starting from the leftmostcell, we scan to the right to determine the label i l of the first cell in which the numberof A particles drops under N / N B ( i l , s ) for a large discrete time s at whichthe profile has reached a steady shape. Then, starting from the rightmost cell labeled ℓ , we follow a similar procedure and determine the label i r of the first cell in which thenumber of A particles overcomes N / N B ( i r , s ) for the same discrete time s .The instantaneous value of the shift deduced from the master equation at discrete time s is then given by ( N − N B ( i l , s ) − N B ( i r , s )) / Ω . The values of the shift h d, ME usedto draw Fig. 4 are obtained after a time average between the times t = 1 and t = 10 inarbitrary units, i.e. between s = 1 . × and s = 1 . × in number of time steps.The shift h d, ME between the profiles of A and B is sensitive to the fluctuations of thenumber of particles described by the master equation. Introducing an appropriate cutoffsatisfying Eq. (12) in the reactive term of the deterministic equations given in Eqs. (14)9 . − . − . − . . . s c a l e d s h i f t D B D A h d /C h d, cut /C h d, ME /C Figure 4: Dilute system. Scaled shifts h d, ME /C , h d, cut /C , and h d /C between the profilesof species A and B versus ratio of diffusion coefficients D B /D A . The values of h d, ME /C (red circles) are deduced from the master equation (Eqs. (7-9)). The values of h d, cut /C (black open triangles) are deduced from the deterministic equations (Eqs. (14 and 15))with a cutoff ε = 10 − . The values of h d /C (blue open squares) are deduced from thedeterministic equations (Eqs. (2) and (3)) without cutoff. The line gives the results for D A = D B . The parameters are given in the caption of Fig. 1.and (15) leads to values of the shift h d, cut in very good agreement with the results h d, ME of the master equation.Considering the deterministic equations, we deduce the width of A profile from thesteepness A ′ (0) in the moving frame at the origin z = 0 and find W d = C / | A ′ (0) | (17)where A is solution of Eqs. (2) and (3) without cutoff. The same definition is applied toEqs. (14) and (15) to obtain the width W d, cut in the presence of a cutoff. The definitionhas to be adapted to take into account the fluctuations of the profile deduced from themaster equation. Using the cell labels i l and i r determined for the shift between thefluctuating A and B profiles solutions of Eqs. (7-9), we define the mean cell label i m asthe nearest integer to the average ( i l + i r ) /
2. We use Eq. (17) with | A ′ (0) | ≃ ( N A ( i m − − N A ( i m + 40)) / (81 ∆ x Ω ) to compute the instantaneous width. As in the case of theshift h d, ME between the fluctuating profiles of A and B, the values W d, ME of the widthused to draw Fig. 5 are obtained after a time average between the times t = 1 and t = 10.As shown in Fig. 5, the width W d deduced from the deterministic equations withoutcutoff is smaller (resp. larger) for D B < D A (resp. D B > D A ) than the width evaluated10 . . . . . . . . . . . . p r o fi l e w i d t h D B D A W d W d, cut W d, ME Figure 5: Dilute system. Profile widths deduced from different approaches versus ratioof diffusion coefficients D B /D A . The values of W d, ME (red circles) are deduced from themaster equation (Eqs. (7-9)). The values of W d, cut (black open triangles) are deducedfrom the numerical integration of the deterministic equations (Eqs. (14) and (15)) witha cutoff ε = 10 − . The values of W d (blue open squares) are deduced from the numer-ical integration of the deterministic equations (Eqs. (2) and (3)) without cutoff. Theparameters are given in the caption of Fig. 1.at W ∗ in the case of identical diffusion coefficients D B = D A [12]. The width W d, ME deduced from the master equation (Eqs. (7-9)) and the width W d, cut deduced from thedeterministic equations (Eqs. (14) and (15)) with a cutoff obeying Eq. (12) agree and areboth smaller than the width W d of the wave front, solution of the deterministic equationswithout cutoff.According to the good agreement between the results of the master equation and thedeterministic equations with a cutoff, it is more relevant to describe the effect of thefluctuations on the wave front as the effect of the discretization of the variables than apure noise effect.Figure 6 summarizes the effect of the fluctuations on the three quantities q for q = v, h, W in the whole range of considered values of the ratio D B /D A for the dilute system.The relative differences ( q d, ME − q d ) /q d between the results deduced from the masterequation and the deterministic equations without cutoff are given in Fig. 6 for the velocity,the shift, and the width. In the whole range of D B /D A , the discrete nature of the numberof particles in the master equation induces a small decrease of 5% of the profile widthwith respect to the deterministic description without cutoff. A significant increase of 14%of the shift between the A and B profiles is observed in the presence of fluctuations in the11 . − . . . . D B D A ( h d, ME − h d ) /h d ( W d, ME − W d ) /W d ( v d, ME − v d ) /v d Figure 6: Dilute system. Relative differences between the front properties deduced fromthe master equation (Eqs. (7-9)) and the analogous properties deduced from the deter-ministic equations without cutoff (Eqs. (2 and (3)) versus D B /D A . The large red circlesgive the relative difference ( v d, ME − v d ) /v d for the front speed, the blue circles of inter-mediate size give the relative difference ( h d, ME − h d ) /h d for the shift between A and Bprofiles, and the small black circles give the relative difference ( W d, ME − W d ) /W d for thewidth of A profile. The parameters are given in the caption of Fig. 1.entire interval of ratios of diffusion coefficients. As for the width, the relative differenceof velocity ( v d, ME − v d ) /v d , with v d = v ∗ , is negative and takes the same value of −
5% for D B /D A ≤
1. However, the relative difference of velocity is not constant for D B /D A > −
22% for D B /D A = 16. Hence, a significant speed decrease is observedwhereas the shift and the width, far behind the leading edge of the front, are not affectedby large diffusion coefficients of species B with respect to the diffusion coefficient of speciesA. In a dilute system, the solvent S is in great excess with respect to the reactive speciesA and B. The concentration of the solvent is then supposed to remain homogeneousregardless of the variation of concentrations A and B . In a concentrated solution, thevariation of the concentration of the solvent cannot be ignored. In the linear domain ofirreversible thermodynamics, the diffusion fluxes are linear combinations of the concen-tration gradients of the different species. The flux j X of species X=A, B, S depends on theconcentration gradients and the diffusion coefficients of all species A, B, and S [38, 39].Using the conservation relations C tot = A + B + S , where C tot is a constant, we eliminate12he explicit dependence of the fluxes on the concentration S of the solvent and find j A = − (cid:18) − AC tot (cid:19) D A ∂ x A + AC tot D B ∂ x B (18) j B = BC tot D A ∂ x A − (cid:18) − BC tot (cid:19) D B ∂ x B (19)According to the expression of the diffusion fluxes in a concentrated system, the reaction-diffusion equations associated with the chemical mechanism given in Eq. (1) read [39] ∂ t A = D A ∂ x (cid:20)(cid:18) − AC tot (cid:19) ∂ x A (cid:21) − D B ∂ x (cid:18) AC tot ∂ x B (cid:19) + kAB (20) ∂ t B = D B ∂ x (cid:20)(cid:18) − BC tot (cid:19) ∂ x B (cid:21) − D A ∂ x (cid:18) BC tot ∂ x A (cid:19) − kAB (21)The discrete expression of the flux at the interface between cells i and i + 1 is related tothe difference of the transition rates in the master equation according to j X ( i + / ) = − ∆ x (cid:16) T − N X ( i +1) − T + N X ( i ) (cid:17) (22)where X = A, B , the transition rate T − N X ( i +1) is associated with the jump of a particle Xto the left from cell i + 1 to cell i , and T + N X ( i ) is associated with the jump of a particleX to the right from cell i to cell i + 1. Using Eqs. (18) and (19) and replacing ∂ x X by( N X ( i + 1) − N X ( i )) / Ω∆ x for X = A, B , we assign well-chosen terms of the flux j X ( i + / )to the transition rates to the left and to the right T ± N A ( i ) = D A ∆ x N A ( i ) − N A ( i ± / ) Ω C tot ∆ x [ D A N A ( i ) − D B N B ( i ± T ± N B ( i ) = D B ∆ x N B ( i ) − N B ( i ± / ) Ω C tot ∆ x [ D B N B ( i ) − D A N A ( i ± N X ( i ± / ) of particles X = A, B in the virtual cell i ± / cannot be used since it may lead to a non-zero transition rate when the departurecell is empty. Instead, we choose the harmonic mean between the number of particles incells i and i ± N X ( i ± / ) = N X ( i ) N X ( i ± N X ( i ) + N X ( i ±
1) (25)which ensures that no jump of X from cell i to cell i ± N X vanishes in cell i . We checked different definitions of the mean obeying13he latter condition and found that the results are not significantly affected when choos-ing for N X ( i ± / ) a modified arithmetic mean which vanishes if N X ( i ) = 0 and equals( N X ( i ) + N X ( i ± / q N X ( i ) N X ( i ± ∂ P ( φ ) ∂ t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) diff = P i (cid:20) T − N A ( i )+1 P ( { N A ( i − − , N A ( i ) + 1 } )+ T + N A ( i )+1 P ( { N A ( i ) + 1 , N A ( i + 1) − } )+ T − N B ( i )+1 P ( { N B ( i − − , N B ( i ) + 1 } )+ T + N B ( i )+1 P ( { N B ( i ) + 1 , N B ( i + 1) − } ) − (cid:16) T − N A ( i ) + T + N A ( i ) + T − N B ( i ) + T + N B ( i ) (cid:17) P ( φ ) (cid:21) (26)The reaction term ∂ P ( φ ) ∂ t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) reac of the master equation given in Eq. (8) for the dilute systemis unchanged in the case of a concentrated system. The kinetic Monte Carlo algorithmand the initial and boundary conditions used for the dilute system are straightforwardlyextended to the concentrated system.The front speeds v c, ME and v d, ME deduced from the master equation in concentratedand dilute cases, respectively, are compared in Fig. 7. The correction to the wave frontspeed induced by an increase of the ratio of diffusion coefficients D B /D A is smaller for aconcentrated system than for a dilute system. Indeed, in the concentrated case, the dif-fusion of a species depends on the diffusion coefficients of both species. Hence, increasing D B at constant D A has a smaller impact on the velocity since the contribution dependingon D B is partly compensated by the unchanged terms depending on D A .The effect of the departure from the dilution limit on the wave front speed v c, ME deduced from the master equation given in Eqs. (7), (8), and (26) is shown in Fig. 8.The dilution limit v d, ME ( D B /D A = 8) = 17 .
20 is recovered for C /C tot →
0. As C /C tot . w a v e f r o n t s p ee d D B D A v ∗ v ε v d, ME v c, ME Figure 7: Concentrated system. Wave front speed v c, ME deduced from the master equation(Eqs. (7), (8), and (26)) in a concentrated system (red solid disks) for C tot = 50 and speed v d, ME deduced from the direct simulation of the master equation (Eqs. (7-9)) associatedwith the dilute system (red circles) versus ratio of diffusion coefficients D B /D A . Thehorizontal solid line gives the minimum velocity v ∗ (Eq. (4)) of an FKPP front in theabsence of a cutoff. The horizontal dashed line gives the velocity v ε = 18 .
84 given in Eq.(11) for a cutoff ε = 10 − and D A = D B . The parameters are given in the caption of Fig.1. .
05 0 . .
15 0 . .
25 0 . .
35 0 . w a v e f r o n t s p ee d C C tot Figure 8: Concentrated system. Wave front speeds versus the deviation from the dilutionlimit C /C tot . The values of v c, ME (red disks) are deduced from the direct simulation ofthe master equation (Eqs. (7), (8), and (26)) for k = 10, Ω = 10, N = 100, D A = 1, D B = 8, ℓ = 2000, and ∆ x = 0 .
008 ( C = N / Ω ). The horizontal solid line gives theminimum velocity v ∗ = 20 (Eq. (4) ) of an FKPP front, solution of the deterministicequations (Eqs. (2) and (3)) without cutoff. The horizontal dashed line gives the velocity v ε = 18 .
84 given in Eq. (11) for a cutoff ε = 10 − and D A = D B .15 . − . − . − . . . s c a l e d s h i f t D B D A h c /C h c, cut /C h c, ME /C Figure 9: Concentrated system. Scaled shifts h c, ME /C , h c, cut /C , and h c /C between theprofiles of species A and B versus ratio of diffusion coefficients D B /D A . The values of h c, ME /C (red disks) are deduced from the master equation (Eqs. (7), (8), and (26)). Thevalues of h c, cut /C (black solid triangles) are deduced from the deterministic equations(Eqs. (20) and (21)) with a reactive term multiplied by the cutoff H ( A − ε ) for ε = 10 − .The values of h c /C (blue solid squares) are deduced from the deterministic equations(Eqs. (20) and (21)) without cutoff. The other parameters are given in the caption ofFig. 7. The line gives the results for D A = D B .increases, the solution is more concentrated and the cross-diffusion terms become moreimportant, so that the system is less sensitive to the difference between the diffusioncoefficients D A and D B : The wave front speed v c, ME increases and tends to the value v ε = 18 .
84 predicted by Eq. (11) for the cutoff ε = 10 − and D A = D B .The variation of the shifts h c, ME , h c, cut , and h c between the two profiles with respect tothe ratio of the diffusion coefficients D B /D A is shown in Fig. 9 in a concentrated systemfor the three approaches, the master equation and the deterministic descriptions withand without cutoff. As revealed when comparing the results given in Figs. 4 and 9, theeffect of the departure from the dilution limit on the shift is too small for us to evaluatethe difference ( h c, ME − h d, ME ) /h d, ME with a sufficient precision for the fluctuating resultsdeduced from the master equations.The effects of the departure from the dilution limit on the widths W c, ME , W c, cut , and W c of the profile are given in Fig. 10 for the three approaches. The agreement betweenthe results W c, ME and W c, cut deduced from the master equation (Eqs. (7), (8), and (26))and the deterministic equations (Eqs. (14 and 15)) with a cutoff, respectively, is satisfying16 . − . − . − . − . − . . . .
06 0 . D B D A ( W c − W d ) /W d ( W c, cut − W d, cut ) /W d, cut ( W c, ME − W d, ME ) /W d, ME Figure 10: Relative differences ( W c, ME − W d, ME ) /W d, ME , ( W c, cut − W d, cut ) /W d, cut , and( W c − W d ) /W d between the widths in a concentrated system and a dilute system fordifferent approaches versus D B /D A . The values of W c, ME and W d, ME (red disks) arededuced from the master equation (Eqs. (7), (8), and (26) and Eqs. (7-9), respectively).The values of W c, cut and W d, cut (black solid triangles) are deduced from the deterministicequations (Eqs. (20) and (21) and Eqs. (14) and (15), respectively) with a reactive termmultiplied by the cutoff H ( A − ε ) for ε = 10 − . The values of W c and W d (blue solidsquares) are deduced from the deterministic equations (Eqs. (20) and (21) and Eqs. (2)and (3), respectively) without cutoff.considering the high level of noise on the evaluation of the width W c, ME . According toFig. 5, the width in a dilute system is smaller than the width obtained for identicaldiffusion coefficients if D B < D A and larger if D B > D A . The results displayed in Fig.10 prove that, for each description method, the width in a concentrated system is largerthan the width in a dilute system if D B < D A and smaller if D B > D A . Hence, in theentire range of ratios of diffusion coefficients and for deterministic as well as stochasticmethods, the width in a concentrated system is closer to the width obtained for identicaldiffusion coefficients. As for the front speed, the departure from the dilution limit reducesthe effects induced by the difference between the diffusion coefficients. We have performed kinetic Monte Carlo simulations of the master equation associatedwith a chemical system involving two species A and B. The two species have two differentdiffusion coefficients, D A and D B , and are engaged in the autocatalytic reaction A + B → D B of the consumed species. The speed shift obtained for two differentdiffusion coefficients is the same as in the case D A = D B . The main result deducedfrom the master equation is that the front speed sensitively depends on the diffusion co-efficient D B . For D B larger than D A , the front speed decreases as D B increases and issignificantly smaller than the prediction of the linear cutoff theory. The speed decreaseobtained for large values of D B /D A is related to the number N B ε of B particles at theposition of the most advanced A particle in the leading edge of the front. When speciesB diffuses faster that species A, N B ε is significantly smaller than the steady-state value N .We carefully derived the nontrivial expression of the master equation in a concentratedsystem with cross-diffusion. The transition rates are deduced from the diffusion fluxes inthe linear domain of irreversible thermodynamics. The transition rates associated withdiffusion depend on the number of particles not only in the departure cell but also inthe arrival cell. Qualitatively, the conclusions drawn for a dilute solution and D A = D B remain valid, but the front properties deduced from the master equation with cross-diffusion depart less from those obtained for D A = D B . The dependence of the frontproperties on D B /D A in a concentrated system are softened with respect to the dilute case.Cross-diffusion mitigates the impact of the difference between the diffusion coefficients. This publication is part of a project that has received funding from the European Union’sHorizon 2020 (H2020-EU.1.3.4.) research and innovation program under the Marie Sklodowska-Curie Actions (MSCA-COFUND ID 711859) and from the Polish Ministry of Science andHigher Education for the implementation of an international cofinanced project.18 eferences [1] R. A. Fisher, Annals of Eugenics , 355 (1937).[2] A. N. Kolmogorov, I.G. Petrovsky, and N.S. Piskunov, Bulletin of Moscow StateUniversity Series A: Mathematics and Mechanics , 1-25 (1937).[3] W. van Saarloos, Phys. Rep. Mathematical Biology (Springer, Berlin, 1989).[5] V. Mendez, D. Campos, and F. Bartumeus,
Stochastic Foundations in MovementEcology: Anomalous diffusion, invasion fronts and random searches (Springer, Berlin,2014).[6] E. Bouin, V. Calvez, N. Meunier, S. Mirrahimi, B. Perthame, G. Raoul, R. Voituriez,C. R. Acad. Sci. Paris, Ser. I , 761 (2012).[7] J. Fort, N. Isern, A. Jerardino, and B. Rondelli, p. 189-197, in Simulating Prehistoricand Ancient Worlds, Eds. J. A. Baracelo and F. Del Castillo, Springer, Cham (2016).[8] R. Mancinelli, D. Vergni, A. Vulpiani, Physica D , 175 (2003).[9] D. Froemberg, H. Schmidt-Martens, I. M. Sokolov, and F. Sagues, Phys. Rev. E ,011128 (2008).[10] X. Cabré and J.-M. Roquejoffre, C. R. Acad. Sci. Paris, Ser. I , 1361 (2009).[11] F. El Adnani and H. Talibi Alaoui, Topol. Methods Nonlinear Anal. , 43 (2010).[12] G. Morgado, B. Nowakowski, and A. Lemarchand, Phys. Rev. E , 022205 (2019).[13] V. K. Vanag and I. R. Epstein, Phys. Chem. Chem. Phys. , 897-912 (2009).[14] D. G. Leaist, Phys. Chem. Chem. Phys. , 4732–4739 (2002).[15] V. K. Vanag, F. Rossi, A. Cherkashin, and I. R. Epstein, J. Phys. Chem. B ,9058–9070 (2008). 1916] F. Rossi, V. K. Vanag, and I. R. Epstein, Chem. A Eur. J. , 2138–2145 (2011).[17] M. A. Budroni, L. Lemaigre, A. De Wit and F. Rossi, Phys. Chem. Chem. Phys. ,1593 (2015).[18] M. A. Budroni, J. Carballido-Landeira, A. Intiso, A. De Wit, and F. Rossi, Chaos , 064502 (2015).[19] L. Desvillettes, Th. Lepoutre, and A. Moussa, SIAM J. Math. Anal. , 820 (2014).[20] L. Desvillettes and A. Trescases, J. Math. Anal. Appl. , 32 (2015).[21] L. Desvillettes, T. Lepoutre, A. Moussa, and A. Trescases, Commun. Part. Diff. Eq. , 1705 (2015).[22] Advances in Reaction-Cross-Diffusion Systems , A. Jüngel, L. Chen, and L. Desvil-lettes, Eds, Nonlinear anal. , 1-492 (2017).[23] E. Daus, L. Desvillettes, H. Dietert, J. Differ. Equ. , 3861 (2019).[24] A. Moussa, B. Perthame, and D. Salort, J. Nonlinear Sci. , 139 (2019).[25] H. P. Breuer, W. Huber, and F. Petruccione, Physica D , 259 (1994).[26] A. Lemarchand, A. Lesne, and M. Mareschal, Phys. Rev. E , 4457 (1995).[27] E. Brunet and B. Derrida, Phys. Rev. E , 2597 (1997).[28] J. S. Hansen, B. Nowakowski and A. Lemarchand, J. Chem. Phys. , 034503(2006).[29] D. Panja, Phys. Rep. , 87-174 (2004).[30] J. G. Conlon and C. R. Doering, J. Stat. Phys. , 421 (2005).[31] J. Mai, I. M. Sokolov, V. N. Kuzovkov, and A. Blumen, Phys. Rev. E , 4130(1997). 2032] G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems (Wiley,New York, 1977).[33] C.W. Gardiner,
Handbook of Stochastic Methods for Physics, Chemistry and the Nat-ural Sciences (Springer, Berlin, 1985).[34] D. T. Gillespie, J. Chem. Phys. , 2340 (1977).[35] Wave front for a reaction-diffusion system and relativistic Hamilton-Jacobi dynamicsS. Fedotov, Phys. Rev. E , 5040 (1999).[36] S. Mirrahimi, G. Barles, B. Perthame, and P.E. Souganidis, SIAM J. Math. Anal. , 4297 (2012).[37] C. Bianca and A. Lemarchand, Physica A , 1 (2015).[38] S. R. de Groot and P. Mazur, Non-Equilibrium Thermodynamics (North-Holland,Amsterdam, 1962).[39] L. Signon, B. Nowakowski, and A. Lemarchand, Phys. Rev. E93