Long-time asymptotics of non-degenerate non-linear diffusion equations
LLong-time asymptotics of non-degenerate non-linear diffusion equations
Ivan C. Christov, a) Akif Ibraguimov, b) and Rahnuma Islam c)1) School of Mechanical Engineering, Purdue University, West Lafayette,Indiana 47907, USA Department of Mathematics & Statistics, Texas Tech University, Lubbock,Texas 79409, USA (Dated: 22 May 2020)
We study the long-time asymptotics of prototypical non-linear diffusion equations.Specifically, we consider the case of a non-degenerate diffusivity function that is a(non-negative) polynomial of the dependent variable of the problem. We motivatethese types of equations using Einstein’s random walk paradigm, leading to a partialdifferential equation in non-divergence form. On the other hand, using conservationprinciples leads to a partial differential equation in divergence form. A transformationis derived to handle both cases. Then, a maximum principle (on both an unboundedand a bounded domain) is proved, in order to obtain bounds above and below for thetime-evolution of the solutions to the non-linear diffusion problem. Specifically, thesebounds are based on the fundamental solution of the linear problem (the so-calledAronson’s Green function). Having thus sandwiched the long-time asymptotics ofsolutions to the non-linear problems between two fundamental solutions of the linearproblem, we prove that, unlike the case of degenerate diffusion, a non-degeneratediffusion equation’s solution converges onto the linear diffusion solution at long times.Select numerical examples support the mathematical theorems and illustrate theconvergence process. Our results have implications on how to interpret asymptoticscalings of potentially anomalous diffusion processes (such as in the flow of particulatematerials) that have been discussed in the applied physics literature. a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected] a r X i v : . [ m a t h . A P ] M a y . INTRODUCTION Almost two centuries ago, Robert Brown observed the apparently random motion ofpollen particles on the surface of a liquid layer (Brown, 1828). Since then, what has becomeknown as “Brownian motion” (Frey and Kroy, 2005) continues to offer scientific insights intomicroscopic phenomena, including in frontier areas such as microrheological measurementsof complex fluids (Zia, 2018). However, Brown’s work did not yield a working theory. Three-quarters of a century later, Einstein proposed the first complete mathematical descriptionof this phenomenon (Einstein, 1905). Specifically, Einstein showed that the spread of thepollen particles obeys a diffusion process, when viewed macroscopically in the sense of aprobability distribution of where the particles might be found. The diffusion process arisesfrom the random motion (walk) of the pollen particles caused by their endless collisions withthe thermally agitated molecules of the fluid in which they are suspended. The final piece insolving the puzzle of Brown’s experiments was Einstein’s determination (contemporaneouslywith Sutherland (1905) and von Smoluchowski (1906)) of the diffusivity in terms of the fluidproperties (viscosity, temperature) via the previous result of Stokes on viscous fluid drag(Stokes, 1851).In this classical example of diffusion, the process is governed by a linear equation, specif-ically a linear evolutionary parabolic partial differential equation (PDE) for the probability,concentration or another related quantity that describes the collection of particles at themacroscale. The diffusivity is constant and set by the Stokes–Einstein formula. Since thisclassical work, however, diffusion equations have been derived (and analyzed) as the gov-erning equation of various other phenomena as well, ranging from flows through porousmedia (Barenblatt, 1952; Philip, 1970) to the motion of free interfaces bounding thin liq-uid films (Oron, Davis, and Bankoff, 1997) to high-temperature shock wave phenomena(Zel’dovich and Raizer, 1967). In contrast to Brownian motion and Einstein’s theory, thelatter examples lead to non-linear diffusion problems in which the diffusivity is a functionof the dependent variable (say, concentration or probability). Interestingly, almost all ofthese examples feature degenerate non-linear diffusion, i.e., the diffusivity vanishes when thedependent variable (or its gradient (Celik et al. , 2017)) vanishes (V´azquez, 2007).A more recent example of a diffusion process concerns the flows of granular materials,such as sand. These materials are macroscopic and the thermal fluctuations of “molecules”2re irrelevant. The granular material must be driven by external forces to flow, which givesrise to collective diffusion. Specifically, during flow, particles (grains) collide with each other.Although these collisions are deterministic, they occurs so often and in such variety that onemay consider the end result to be random particle velocity fluctuations, much like to thoseimparted on pollen by the fluid’s molecules in Brown’s 1828 experiments. For a collectionof identical (in shape, size and density) particles, one might expect that the diffusivity isindependent of concentration (Lacey, 1954; Savage, 1993), and indeed experiments supportthis claim (Cahn et al. , 1966; Zik and Stavans, 1991). However, granular materials are oftenhighly heterogeneous, i.e., they are mixtures of particles of different shapes, size and density(Ottino and Khakhar, 2000; Umbanhowar, Lueptow, and Ottino, 2019). In the case ofparticles of different sizes, depending on the mixture proportions, the positive time intervalbetween small-small or small-large binary collisions depends on the number of particles of agiven type locally. The small-small and small-large collisions are not identical as the largerparticles impact a bigger force onto the smaller ones, therefore produce a different free jump.Nonetheless, the free jump frequency distribution itself could be assumed independent ofthe particle concentration because the jumps are entirely set by the collision physics (mass,velocity, coefficient of restitution, etc.).Then, there are three main parameters that are involved in Einstein’s random walkparadigm (Einstein, 1905): the frequency distribution of free jumps ϕ ( δ ), the length ofa free jump δ , and time interval τ within which particles perform a free jump. A thoughtexperiment along the lines of Einstein can incorporate non-linearity in the random walkmodel. The influence of particle concentration can be taken into account through the inter-val of free jumps. In § III, we will discuss how the diffusivity will depend, through τ , on howmany large (or small particles) are in the vicinity of a spatial location. This observationleads to concentration-dependent diffusion of poly-disperse granular materials (Fischer et al. ,2009). Interestingly, however, the non-constant diffusivity in these cases is not degenerate(Ristow and Nakagawa, 1999; Dury and Ristow, 1999), as even in the absence of largeparticles nearby, the small particles still collide and, thus, “diffuse” (Christov and Stone,2012).Thus, we have provided ample motivation that in many applications involving the flow ofliquids, gases, and even particulates, the dependent variable in the problem, such as density,pressure or concentration, can be governed by a non-linear evolutionary parabolic PDE, i.e., a3iffusion equation. In the case of degenerate diffusivity, a large mathematical literature existsdiscussing the qualitative properties of solutions, asymptotics, and so on (DiBenedetto, 1993;V´azquez, 2007). On the other hand, the case of non-degenerate diffusion has not received asmuch attention. Nevertheless, in applications involving such PDE, one observes “regular”diffusive scalings at long times (in the sense of intermediate asymptotics (Barenblatt, 1996)),which has caused no small amount of controversy in interpreting experiment and simulationdata (Christov and Stone, 2012). Motivated by the need for a clear mathematical answer regarding the long-time behaviorof solutions (distinct from the short-time behavior studied by Christov and Stone (2012);Sekimoto and Fujita (2019)) of non-degenerate non-linear parabolic equations, in the presentwork, we establish new mathematical results showing that the long-time asymptotics ofsolutions are given by a Gaussian profile and its corresponding “normal” (as apposed to“anomalous”) scalings obtainable from the linear diffusion equation. Specifically, using themaximum principle ( §§ V and VI), we obtain estimates, from above and below, on the solutionof the non-linear problem (applicable to equations in both divergence and non-divergenceform, as shown in § IV), proving that the solution to the non-linear problem converges tothe Green’s function of the linear problem (suitably shifted). These mathematical resultsare illustrated via numerical simulations in § VII, and conclusions are stated in § VIII.
II. POSITION OF THE PROBLEM
As a generic model for all the previously mentioned phenomena, consider:˜ Lu = ∂u∂t − ∂∂x (cid:18) P ( u ) ∂u∂x (cid:19) = 0 , ( x, t ) ∈ ( −∞ , + ∞ ) × (0 , ∞ ) , (1)where P ( u ) is a polynomial with P ( u ) > ∀ u ≥
0, and u is a scalar quantity that ischaracteristic of the physical systems (to be made precise below). Equation (1) is subjectto a non-negative localized initial condition: u ( x,
0) = u ( x ) . (2) Here, by “scalings” we mean that the solution u ( x, t ) can, at any t sufficiently large, be trans-formed/collapsed as u ( x, t ) (cid:55)→ κ t n U (cid:0) κ xt − n (cid:1) , to a good approximation, into a universal profile U ( · ),for some suitable dimensional constants κ , and scaling exponents n , .
4o fully pose the problem, appropriate growth conditions should be satisfied at | x | = ∞ bythe solution and the initial condition:lim | x |→∞ u ( x, t ) | x | γ = 0 ∀ t ≥ , (3)for some γ >
0. The solution, which is initially non-negative, should remain non-negativefor all times: u ( x, t ) ≥ , ∀ ( x, t ) ∈ ( −∞ , + ∞ ) × [0 , ∞ ) . (4)Observe that solutions to Eqs. (1)–(3) obey a conservation principle (see also § IV B): (cid:82) + ∞−∞ u ( x, t ) d x = (cid:82) + ∞−∞ u ( x ) d x > ∀ t ≥
0, as easily shown by direct integration andapplication of the decay condition as | x | → ∞ .Since we are studying a generic mathematical problem, we do not concern ourselves withthe units. Physically, this just means that we have made x dimensionless by some charac-teristic domain length x c , we have made P dimensionless by some characteristic diffusivity P c (such that, for P (cid:55)→ P/P c , we may take P (0) = 1 now), we have made u dimensionlessby some characteristic scale u c (say, such that for u (cid:55)→ u/u c , (cid:82) + ∞−∞ u ( x, t ) d x = 1 ∀ t ≥ t (cid:55)→ t/t c dimensionless by the characteristic diffusion time t c = x c /P c .Here, the subscript ‘ c ’ stands for ‘characteristic.’In Eq. (1), we consider the case in which P ( u ) is a polynomial such that P ( u ) > ∀ u ≥ P ( u ) has no real roots. Hence, Eq. (1) is non-degenerate, whichis specifically the case of interest here. A particular example we are interested in, based onprevious studies (see, e.g., Ristow and Nakagawa, 1999; Christov and Stone, 2012), is P ( u ) = a + a u with a >
0, which is of the form assumed above. Equation (1) is usuallyderived from conservation principles (see § IV B) in conjunction with the physical law thatthe flux is proportional to the gradient of the scalar function u . Examples of the latter laware Darcy’s, Fourier’s, Fick’s, etc. (Bird, Stewart, and Lightfoot, 2002).At the same time, this process can be also be modeled using the Einstein paradigm( § III), which he used to derive the linear diffusion equation describing Brownian motionmacroscopically. Therefore, we next discuss this derivation, in detail, for the non-linearcase. We will show that the governing PDE has non-divergence form and can be derivedfrom a probabilistic model via a thought experiment. To connect the non-divergence-formequation (obtained from Einstein’s paradigm) to the divergence-form equation (1) (obtained5rom a conservation principle), in § IV, we show there exists a closed-form mapping betweenthe solutions of the two.
III. EINSTEIN PARADIGM: “FROM RANDOM MOTION OFPARTICLES TO DIFFUSION”
In this section, we explain how non-linearity can be incorporated into Einstein’s randomwalk model of Brownian motion. In the celebrated work of Einstein (1905, 1956), a mathe-matical model is derived on the basis of a thought experiment (or,
Gedankenexperiment inEinstein’s own terminology (Perkowitz, 2010)). On the basis of its generality, this model canthen be applied to a number of physical processes, in which a random walk occurs (see, e.g.,Gardiner, 2009, § Assumption 1.
There exists a time interval τ , which is very small compared to the ob-servable time intervals but large enough that the motions, performed by particles during twoconsecutive time intervals τ , can be considered as mutually independent events. Assumption 2.
The distance traveled during the time interval τ , without undergoing acollision, is called the “free jump” and has a finite size, δ . Assumption 3.
The particles are not allowed to interact chemically (i.e., they cannot ag-glomerate, or breakup, or react with a solvent).
Let the total number of particles present in the system be N . The number d N of particlesexperiencing a displacement that lies between δ and δ + d δ in the time interval τ is given byd N = N ϕ ( δ ) d δ, (5)where ϕ is the probability density function of particle jumps such that (cid:90) δ max δ min ϕ ( δ ) d δ = 1 . (6)6instein assumed that ϕ is localized, i.e., it differs from zero only in a range of δ valuesabout δ = 0. This assumption seems natural, but it does not have always be true for allphysical diffusion processes. Remark 1.
Note that in the Einstein paradigm, τ , δ and ϕ are characteristics of the physicalprocess. In general, these three parameters can be functions of both the spatial variable x and the time variable t , as well as other physical quantities (depending on the problem).Furthermore, at this point, nothing prevents τ , δ and ϕ from also being functions of thedependent variable (and its derivatives), such as the number of particles N .In the present work, however, we will assume that the length of the free jumps δ and theirfrequency distribution ϕ ( δ ) are fixed (by the underlying physics) w.r.t. N . Then, the onlyparameter involved that can depend on N is the time interval τ , which leads to the non-linearnature of the diffusion process below. Assumption 4.
Let f ( x, t ) be the number of particles per unit volume. Then, the numberof particles found at time t + τ between two planes perpendicular to the x -axis, with abscissas x and x + d x , is given by f ( x, t + τ ) · d x = (cid:18)(cid:90) δ max δ min f ( x + δ, t ) ϕ ( δ )d δ (cid:19) · d x. (7)Next, by Caratheodory’s theorem, there exists a function ψ ( x, t ) such that f ( x, t + τ ) = f ( x, t ) + τ ψ ( x, t + τ ) , (8)where lim τ → ψ ( x, t + τ ) = ∂f ( x, t ) ∂t . (9)However, we shall not formally take the limit. Instead, since τ (cid:28) t (see Remark 3 below),we make the approximation ψ ( x, t + τ ) ≈ ∂f ( x, t ) ∂t . (10)Nest, using the Taylor expansion of f ( x + δ, t ) in powers of δ , we obtain f ( x + δ, t ) = f ( x, t ) + δ ∂f∂x + δ ∂ f∂x + · · · . (11)Thus, we obtain from Eqs. (7) and (11): f ( x, t )+ τ ∂f∂t = f ( x, t ) (cid:90) δ max δ min ϕ ( δ ) d δ (cid:124) (cid:123)(cid:122) (cid:125) =1 by Eq. (6) + ∂f∂x (cid:90) δ max δ min δϕ ( δ ) d δ + ∂ f∂x (cid:90) δ max δ min δ ϕ ( δ ) d δ + · · · , (12)7r τ ∂f∂t = (cid:90) δ max δ min δϕ ( δ ) d δ · ∂f∂x + (cid:90) δ max δ min δ ϕ ( δ ) d δ · ∂ f∂x + · · · . (13)Now, we impose restrictions on the jump size distribution ϕ , to make precise Einstein’sassumption on the localized nature of ϕ : Assumption 5. | δ min | (cid:29) | δ min | and | δ max | (cid:29) | δ max | . Assumption 6. ϕ ( δ ) = ϕ ( − δ ) . Due to Assumption 5, we keep only second-order (in δ ) terms in the right-hand side ofEq. (13). Due to Assumption 6, the first and further odd moments vanish. Now, we definethe diffusivity (diffusion coefficient) D via the second moment of ϕ :1 τ (cid:90) δ max δ min δ ϕ ( δ ) d δ = D, (14)whence Eq. (13) leads to the well-known linear diffusion equation for the function f countingthe number of particles per unit volume: ∂f∂t = D ∂ f∂x . (15) Remark 2.
According to our interpretation of the Einstein paradigm (Remark 1), the diffu-sion coefficient D in Eq. (14) can be a function of x and t , or even f , via the time interval τ . In other words, unlike previous works that apply the Einstein paradigm to non-linear dif-fusion (Boon and Lutsko, 2007; Lenzi et al. , 2019), the non-linearity in the present contextcomes into play via τ and its possible direct dependence on f specifically. Remark 3. “Einstein’s derivation is really based on a discrete time assumption, that impactshappen only at times 0, τ , 2 τ , 3 τ , ...” (Gardiner, 2009, p. 5). In other words, in the Einsteinparadigm (Remark 1), τ is considered to be the finite time between particle collisions. Thismicroscopic time scale is assumed to be small compared to the observational (macroscopic)time scale, τ (cid:28) t , but it is not taken to zero. Therefore, here, Eq. (14) is interpreted asstated, not as a limiting process. V. NON-LINEAR PARABOLIC EQUATIONS IN DIVERGENCE ANDNON-DIVERGENCE FORMA. Non-linear model arising from the Einstein paradigm (equation innon-divergence form)
First, we establish that the Einstein paradigm can be used to obtain a non-linear parabolicdiffusion equation for a Brownian-like process.
Assumption 7.
Let the number of particles per volume, f ( x, t ) , be linearly proportional toa scalar function, such as the concentration v ( x, t ) in the medium, which is assumed to behomogeneous and isotropic. We postulate that time-interval of free jumps τ (cid:37) a − > as v (cid:38) . In other words, we consider the concentration-dependent diffusivity function D ( v ) = a + F ( v ) , (16) where F ( v ) is a non-decreasing homogeneous function. Under Assumption 7 and for general for x ∈ R d , Eq. (15) takes the form Lv = ∂v∂t − D ( v )∆ v = 0 , (17)where D ( v ) is the diffusion coefficient for concentration v at the point ( x, t ), and ∆ v = (cid:80) di =1 ∂ v∂x i is the Laplacian operator applied to v . Taking the scalar function v ≥ non-degenerate . B. Non-linear model arising from the conservation law principle (equation indivergence form)
Another way to take into account of non-linearity is to employ the more traditionalconservation law for the density ρ ( x, t ) of a substance of interest with attendant flux vector (cid:126)J ( x, t ) (i.e., the conservation of mass or “continuity” equation (Bird, Stewart, and Lightfoot,2002; Dafermos, 2016)). This equation, for a finite control volume, is written in divergenceform as ∂∂t (cid:90) V ρ d x + (cid:73) ∂ V (cid:126)J · (cid:126)ν d s = 0 . (18)Here, V ⊂ R d is the control volume, ∂ V ⊂ R d − is its boundary (with unit normal vector (cid:126)ν ), and t ∈ R is the time variable. 9 ssumption 8. To allow for the consideration of different physical phenomena governed bythe same equations, let the density ρ = Au + B be a linear function of some scalar u ( x, t ) ,where A > and B ≥ are suitable dimensional constants. For example, u can be theconcentration. Assumption 9.
Suppose that Fick’s law (Fick, 1855a,b) holds a.e. That is, the flux vector (cid:126)J is proportional to the gradient of concentration: (cid:126)J = − P ∇ ρ = − P A ∇ u, (19) where the proportionality factor P is precisely the diffusivity (in general, a tensor of secondrank) (Bird, Stewart, and Lightfoot, 2002). Assumption 10.
In an isotropic medium, the diffusivity P in Eq. (19) is a scalar. Specifi-cally, let P ( · ) be a non-decreasing function of u : P ( u ) = a + G ( u ) , where G ( u ) ≥ . Remark 4.
Assumption 10 is simply a restatement of Assumption 7 but for the case of adivergence-form equation derived from the conservation law principle.
Under the above assumptions, and applying Green’s theorem for the arbitrary controlvolume V , Eq. (18) can be transformed to an equation in divergence form at space-timepoint ( x, t ): ˜ Lu = ∂u∂t − ∇ · ( P ( u ) ∇ u ) = 0 . (20)Here, due to Assumption 10, the diffusivity is a non-decreasing function w.r.t. u , and wetake the constant a > u ≥ P ( u ) > ∀ u , and it follows that Eq. (20) is non-degenerate . C. Mapping between solutions of equation in divergence and non-divergenceform
The two governing diffusion equations, i.e., Eqs. (20) and (17), introduced above areobviously related. Now, we prove that, for any positive polynomial function P ( u ) withnon-negative coefficients satisfying Assumption 10, there exists a function D ( v ) satisfyingAssumption 7 s.t. P ( u ) = D ( v ). 10 heorem 1. Let P ( u ) be a polynomial P ( u ) = a + a u + · · · + a n u n (21) with all non-negative coefficients and at least one strictly positive coefficient. Consider thetransition formula v = (cid:90) u P ( ξ ) d ξ. (22) Then, from monotonicity of integral it follows that there exists a function F ( v ) s.t. P ( u ) = D ( v ) , (23) where D ( v ) = a + F ( v ) as given in Eq. (16) .Proof. For first-order polynomials, the construction of the function F ( v ) is explicit. Indeedfor n = 1, P ( u ) = a + a u . Let F ( v ) = − a + (cid:113) a + 2 a v. (24)Then, a direct substitution shows that: D ( v ) = P ( u ) . For the general case, this construction is not explicit.Next, we show that if u and v are related by the transition formula (22), and D ( v ) satisfiesEq. (23), then the left-hand side of Eq. (17) takes the form: Lv : ∂v∂t − P ( u )∆ v. (25)Observe that both functions u and v are involved in Eq. (25).Denote by C , the class of function that have continuous second derivatives w.r.t. x andcontinuous first derivatives w.r.t. t . Then, we prove: Lemma 1. If u, v ∈ C , , then the transition formula (22) implies that Lv = ∂v∂t − D ( v )∆ v = ∂v∂t − P ( u )∆ v = P ( u ) ˜ Lu, (26) where ˜ Lu is defined in (20) . roof. Let v = ψ ( u ), where ψ ( u ) = (cid:82) u P ( ξ ) d ξ and so, P ( u ) = ψ (cid:48) ( u ). Then, we compute: Lv = ∂v∂t − P ( u )∆ v = ψ (cid:48) ( u ) ∂u∂t − P ( u ) ψ (cid:48)(cid:48) ( u ) |∇ u | − P ( u ) ψ (cid:48) ( u )∆ u = ψ (cid:48) ( u ) (cid:20) ∂u∂t − P ( u ) ψ (cid:48)(cid:48) ( u ) ψ (cid:48) ( u ) |∇ u | − P ( u )∆ u (cid:21) = ψ (cid:48) ( u ) ˜ Lu.
Corollary 1.
Since ψ (cid:48) ( u ) > , if u is a solution of the equation ˜ Lu = 0 , then v is solutionof the equation Lv = 0 . V. MAXIMUM PRINCIPLE ON A BOUNDED DOMAIN
In this section, we prove a maximum principle for the solution, following Ilyin, Kalash-nikov, and Oleynik (2002) (see also Landis, 1998). Let R d be the d -dimensional real Eu-clidean space. We also consider the ( d + 1)-dimensional space R d +1 , in which the spatialcoordinates are augmented by time: ( x, t ) = ( x , x , . . . , x d , t ). Now, suppose that U ⊂ R d is bounded, and t >
0. We define the cylindrical region Ω = U × (0 , T ] ⊂ R d +1 , andits parabolic boundary Γ = ( U × { t = 0 } ) ∪ ( ∂U × (0 , T ]). We also consider the layer H ⊂ R d +1 ∩ (0 , T ].Again, let C , (Ω) be the class of continuous functions in ¯Ω that have two continuousderivatives in x and one continuous derivative in t inside the domain Ω. Henceforth, bothfunctions u ( x, t ) and v ( x, t ) are assumed to be in C , (Ω). Lemma 2.
For a given function u ∈ C , , if P ( u ) ≥ a > and Lv > in Ω , then v ( x, t ) ≥ min Γ(Ω) v in Ω .Proof. Suppose there is a point ( x , t ) ∈ Ω such that min ¯Ω v = v ( x , t ) < min Γ(Ω) v , thenthe minimum of v ( x, t ) attained at ( x , t ) ∈ Ω, for x ∈ U and t ≤ T . Therefore, at thepoint ( x , t ), ∇ v = 0, ∂v∂t ≤ v ≥ Lv ≤ x , t ), which is a contradiction. Lemma 3.
For a given function u ∈ C , , if P ( u ) ≥ a > ∀ u and Lv ≥ in Ω , then v ≥ min Γ(Ω) v in Ω . roof. Consider the function w ( x, t ) = Kt + v ( x, t ) for ( x, t ) ∈ Ω and t < T , K >
0. Then, Lw = ∂w∂t − P ( u )∆ w = K + ∂v∂t − P ( u )∆ v > Lw >
0, we obtain w ( x, t ) ≥ min Γ(Ω) w in Ω for K >
0. Thus, Kt + v ( x, t ) ≥ min Γ(Ω) v in Ω, which implies KT + v ( x, t ) ≥ min Γ(Ω) v in Ω for K . Taking K → v ≥ min Γ(Ω) v .From standard maximum principle follows the comparison lemma: Lemma 4. If P ( u ) ≥ a > ∀ u (Assumption 10) and Lv = ∂v ∂t − P ( u )∆ v ≥ Lv = ∂v ∂t − P ( u )∆ v in Ω , and v ≥ v on Γ , then v ≥ v in Ω .Proof. It is sufficient to apply the maximum principle to the function w ( x, t ) = v ( x, t ) − v ( x, t ) using the properties of the functions v and v . VI. MAXIMUM PRINCIPLE ON AN UNBOUNDED DOMAIN
Let r = (cid:0) (cid:80) di =1 x i (cid:1) / . As in the previous section, we follow Ilyin, Kalashnikov, andOleynik (2002) and take into account that the coefficients of the operator L are given by P ( u ). Lemma 5.
Suppose the function u ( x, t ) is continuous in Ω = R d × [0 , T ] such that − m , m > . Under theseassumptions, if Lv ≥ and v | t =0 > , then v ( x, t ) ≥ for all ( x, t ) ∈ Ω .Proof. Consider, the auxiliary function w ( x, t ), s.t. w ( x, t ) = mr (cid:0) r + Kt (cid:1) e αt + v ( x, t ) , on the auxiliary domain consisting of the cylinder Q r = { ( r, t ) | r ≤ r , ≤ t ≤ T } , wherethe constants K > α > ∀ r >
0. Indeed, Lw = mr e αt (cid:0) K + αr + Kαt − dP ( u ) (cid:1) + Lv.
Since Lv ≥ P ( u ) = (cid:15) ( r ) r + C , then for r ≥
1, we Lw ≥ mr e αt (cid:0) K + Kαt + αr − d · o( r + 1) (cid:1) > .
13n the other hand, for r < P ( u ) ≤ M for some M . It follows that, if α > dMr , then Lw ≥ mr e αt ( K + Kαt + αr − dM (cid:1) > . Consider now w ( x, t ) in Q r . For t = 0, w ( x, ≥ v ( x, r = r , w ( x, t ) ≥ m + v ≥
0. Then, according to Lemma 2, the inequality w ( x, t ) ≥ R d +1 ∩ (0 < t ≤ T ) is contained in Q r for asufficiently large r . Now, we have w ( x, t ) = mr ( r + Kt ) e αt + v ( x, t ) ≥
0. Taking the limitas r → ∞ , v ( x, t ) ≥ ∀ ( x, t ) ∈ Ω as desired.Similarly, one can prove:
Lemma 6.
Suppose the function u ( x, t ) is such that < a ≤ P ( u ) ≤ (cid:15) ( r ) r + C with lim r →∞ (cid:15) ( r ) = 0 . Let v ∈ C , and Lv ≤ , then v ( x, t ) ≤ max v ( x, for all ( x, t ) ∈ Ω . Now, from Theorem 1 and Lemma 6, it follows:
Corollary 2.
Assume P ( u ) is bounded and ˜ Lu = 0 , then v = (cid:82) u P ( ξ ) d ξ ≤ (cid:82) max u ( x, P ( ξ ) d ξ . However, we would like to make the bound on v more precise. First we prove: Lemma 7.
Let w ( x, t ) = CF s,β ( x, t ) for t ≥ , where the constant C > and F s,β ( x, t ) = ( t + 1) − s e − | x | β ( t +1) (27) is the barrier function. If < a ≤ P ( u ) ≤ β and s ≤ a d β , then Lw ≥ .Proof. First, consider the function F s,β ( x, t ) with constant s >
0, for t (cid:29)
1. Second,compute Lw = C (cid:32) − st + 1 + | x | β ( t + 1) − P ( u ) | x | β ( t + 1) + P ( u ) d β ( t + 1) (cid:33) ( t + 1) − s e − | x | β ( t +1) = C (cid:32) | x | β ( t + 1) (cid:18) β − P ( u ) (cid:19) + 12 β ( t + 1) (cid:18) dP ( u ) − sβ (cid:19)(cid:33) F s,β ( x, t ) . Now, the conjecture of the lemma follows directly from this expression along with the as-sumptions.
Lemma 8.
Let w ( x, t ) = F s,β ( x, t ) as in Eq. (27) . Assume that s and β satisfy the sameconditions as in the previous lemma. If Lv ≤ and w ≥ v on Γ , then v ( x, t ) ≤ w ( x, t ) . roof. Consider, ζ ( x, t ) = w ( x, t ) − v ( x, t ) . Then, Lζ = Lw − Lv.
Since Lv ≤ Lζ ≥ Lw.
Then, due to the conditions on s and β , Lζ ≥
0. Also, w ( x, ≥ v ( x,
0) implies ζ ( x, ≥
0. Then by Lemma 3, ζ ( x, t ) = F s,β ( x, t ) − v ( x, t ) ≥ v ( x, t ) ≤ F s,β ( x, t )in Ω for s > β .From the above results and the positivity of the function u , it follows that 0 < a ≤ P ( u ) ≤ const . Therefore, one can apply Aronson’s estimate (Aronson, 1967) for the Greenfunction of the divergence-form problem ˜ Lu = 0 with bounded coefficients. Namely: Theorem 2.
There exists G ( x, t ) , Green function of the Cauchy problem for the equation ˜ Lu = 0 , s.t. C t − d e − | x | β − t ≤ G ( x, t ) ≤ C t − d e − | x | β + t for some constants C and C , at every ( x, t ) ∈ Ω ( t > ). Remark 5.
In Aronson’s estimate (Aronson, 1967) discussed above, β ∓ depend on the ellip-ticity constant and the spatial bounds of the divergence-form equation’s coefficients. Specif-ically, if we write ˜ L ( · ) = ∂∂t ( · ) − ∇ · ( a ( x, t ) ∇ ( · )) , then the constants β ∓ can be such that β − ≤ a ( x, t ) ≤ β + . Of course, in this representation needed to apply Aronson’s estimate, thecoefficient a ( x, t ) depends on u ( x, t ) directly. Then, since u ( x, t ) → as t → ∞ , it followsthat β − = β + = a (recall a = P (0) ) in the limit as t → ∞ . Thus, finally, we have achieved our main conclusion: Let u ( x, t ) be a solution of theCauchy problem for ˜ Lu = 0, with ˜ L as defined in Eq. (20), and 0 < a ≤ P ( u ) ≤ (cid:15) ( r ) r + C with lim r →∞ (cid:15) ( r ) = 0. Then, starting from a compact initial condition (such as that givenby Eq. (30) below, or other choices suitable for executing the proof with the barrier functionfrom Eq. (27)), there exists constants c and c such that c e − | x | a t ≤ t d u ( x, t ) ≤ c e − | x | a t as t → ∞ . (28)15his mathematical result holds for any non-degenerate polynomial diffusivity function P ( u )as in Eq. (21), and thus significantly constrains the long-time asymptotic scalings thatsolutions to non-degenerate non-linear diffusion equations can exhibit. Determining, orat least constraining, the possible scaling behaviors of such solutions was our motivatingscientific question, which we have now answered. VII. NUMERICAL EXPERIMENTS TO ILLUSTRATE THE MAINMATHEMATICAL RESULTS
In this section, we illustrate our mathematical results (in d = 1 dimensions) with selectednumerical simulations. Specifically, we show that the Green function c t − e − | x | a t , as inEq. (28), does indeed describe, quantitatively, the long-time asymptotic behavior of thesolutions of non-degenerate non-linear parabolic equations.There are many numerical methods that one can use to solve the scalar parabolic equa-tion (1) subject to the initial condition (2) (Strikwerda, 2004). For simplicity, we use the pdepe subroutine of Matlab x ∈ [ − x max , + x max ] (0 < x max < ∞ ), the asymptotic de-cay condition (3) must be replaced with an appropriate BC at x = ± x max . In our numericalexamples, we choose the interval to be large enough ( x max = 200 or larger, depending on thefinal simulation time T ), so that this boundary condition does not influence the diffusionprocess of a localized initial condition. Then, we impose the “natural” (Neumann) boundaryconditions ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) x = ± x max = 0 ∀ t ∈ [0 , T ] . (29)At least 10 000 x -grid points are used for the discretization, and time integration is performedby the adaptive, variable-order multistep stiff solver ode15s in Matlab (Shampine andReichelt, 1997). 16 . Compact initial condition
First, we take the initial condition, u ( x,
0) = u ( x ), to be a box of unit area: u ( x ) = x , | x | ≤ x , , | x | > x . (30)We take x = 1 without loss of generality, but we do note that various constants (in bounds,etc.) will depend on x . Then, we solve numerically the initial-boundary-value problemconsisting of Eqs. (1), (29), (30) on the finite space-time domain [ − x max , + x max ] × (0 , T ].Figure 1(a) shows an example time-evolution of this non-linear diffusion process. Clearly,the long-time numerical solution to the non-linear problem with P ( u ) = 1 + u (note P (0) = a = 1) converges, visually at least, as t → ∞ to the fundamental (Gaussian) solution u G ( x, t ) = 1 (cid:112) πP (0) t e − x P (0) t (31)of the linear problem with P ( u ) = 1 and an initial condition of unit area, (cid:82) + ∞−∞ u ( x ) d x = 1.Note that, since we are interested in the long-time asymptotics (specifically, the scaling ofthe solution), we neglect details arising from the fact the initial condition (30) (see Kleinsteinand Ting, 1971; Witelski and Bernoff, 1998; Christov and Stone, 2012) is not a point source(Dirac δ ) for which, specifically, the exact solution to the linear problem is the fundamentalsolution in Eq. (31).On rescaling the solution u ( x, t ) of the non-linear diffusion problem using the “nor-mal” (linear) diffusion scalings as u ( x, t ) (cid:55)→ (4 πP (0) t ) / U ( ξ ) with ξ = x/ (4 P (0) t ) / , inFig. 1(b,i), we observe the curves begin to approach the Gaussian profile U ( ξ ) = e − ξ . Thebound, Eq. (28) proved in § VI, is difficult to evaluate point-wise numerically because β inthe barrier function (27) (and β ∓ in Theorem 2) can depend, in particular, on t . However,in the rescaled coordinates, it is possible to evaluate the bound in an L sense by computingthe norm of the difference between the rescaled solution of the non-linear problem and theGaussian profile U ( ξ ) = e − ξ . This “error” decays algebraically in time, as Fig. 1(b,ii) shows.Note that the algebraic decay in t is expected from previous estimates of the convergencerate of solutions to nonlinear parabolic equations, starting from arbitrary initial data, to-wards their self-similar intermediate asymptotics (Kleinstein and Ting, 1971; Witelski andBernoff, 1998; Bernoff and Witelski, 2010). The numerical observation that the norm of the17
10 -5 0 5 1000.10.20.30.4 (a) -2 0 200.20.40.60.81 10 -1 -3 -2 -1 (b) FIG. 1. (a) Time-evolution of u ( x, t ) (visualized by 100 solution curves with color changing frompurple/dark to green/light, from t = 0 to t = T = 10) starting from a unit box initial condition (30)(dash-dotted curve) at t = 0, with P ( u ) = 1+ u . (b,i) At long times, when rescaled in the “normal”diffusive way, the solution appears to converge the Gaussian/fundamental solution (31) of the linearproblem with P ( u ) = 1 (dashed curve). (b,ii) The L norm of the difference between the numericalsolution of the non-linear problem and the Gaussian profile, having both been expressed in rescaledvariables as in (b,i), decays algebraically in time (dashed line is a reference slope of t − ). In allplots, the abscissas have been truncated for clarity. “error” decays in time indicates that the bound proved in § VI is accurate, and the solutionof the non-linear problem converges (as t → ∞ ) to the Gaussian profile.Figure 2 shows the equivalent of Fig. 1 but with P ( u ) = 1+ u +10 u . Clearly, as predictedby the mathematical theory, the long-time asymptotics are similar but the constants in thebound (28), and prior theorems and lemmas, change. Furthermore, due to 10 u dominating18
10 -5 0 5 1000.10.20.30.4 (a) -2 0 200.20.40.60.81 10 -1 -2 -1 (b) FIG. 2. Same as Fig. 1 but with the diffusivity P ( u ) = 1 + u + 10 u being a quadratic polynomial. u at early times, the approach to the ultimate long-time asymptotics takes longer, thusthe convergence in Fig. 2(b) looks worse than in Fig. 1(b) for the same integration time-interval: t ∈ (0 , T = 10]. B. Non-compact initial condition with algebraic decay
Second, we take u ( x ) to have slow algebraic decay at infinity as per Eq. (3): u ( x ) = (cid:18) γ − x ( γ −
1) + 1] (cid:19) , | x | ≤ x , ( | x | − x + 1) − γ , | x | > x . (32)Note that the function was chosen so that we still have unit “mass,” i.e., (cid:82) + ∞−∞ u ( x ) d x =1 < ∞ , which restricts the decay rate to γ >
1. Again, x = 1 without loss of generality.19
30 -20 -10 0 10 20 3000.10.20.3 (a) -2 0 200.20.40.60.81 10 -4 -3 -2 -1 (b) FIG. 3. Same as Fig. 1 but starting from the non-compact initial condition in Eq. (32).
Let us take γ = 3 for the remainder of this numerical experiment. Figure 3 shows theequivalent of Fig. 1 but starting from the initial condition in Eq. (32). Clearly, the long-timeasymptotics are similar as in the previous two examples in § VII A starting from compactinitial conditions. Of course, the constants in the various theorems and lemmas change,and there are quantitative differences in the convergence process. Importantly, to reach“sufficiently large” t (cid:29) T . Here, we had touse x max = 2 000, 100 000 x -grid points and integrate up to T = 100 to observe the long-timeasymptotics. Values of γ closer to 1 require even more computational care and resources.This example is interesting as our construction, leading up to the main result in § VI, wasbased on a barrier function, Eq. (27), with exponential decay as | x | → ∞ . However, theinitial condition (and, hence, the solution) in this example does not satisfy this condition,having only algebraic decay. Yet, the numerical results suggest that the final bound obtained20n Eq. (28), holds just as well for this example, a case that violates the assumptions of thetheorems and lemmas proved above. VIII. CONCLUSION
In this paper, we discussed how Einstein’s random walk paradigm can be employed toderive non-linear parabolic equations in non-divergence form, as a continuum description ofa random-walk diffusion process. Then, by proving that a mapping between divergence-formand non-divergence-form parabolic equations exists, we connected the derivation from theEinstein paradigm to the traditional derivation of the diffusion equation from a conservationlaw, i.e., the continuity equation along with Fick’s (or Fourier’s, etc.) law.The mapping theorem enabled us to obtain accurate estimates for the Green function ofthe Cauchy problem for the divergence-form non-linear diffusion equation. These estimateswere used to prove bounds on the solution of the non-linear diffusion equation, from bothabove and below, and establish the long-time asymptotics of the non-linear equation’s so-lutions. Specifically, we proved that the solution to the non-degenerate non-linear problemconverges to the fundamental solution (Gaussian distribution) of the linear diffusion prob-lem, leading to the bound in Eq. (28), which is valid from above and from below. Numericalsimulations for some example non-linear equations quantitatively support the mathematicalresults proved via Aronson’s estimate.Importantly, the present work sheds light on the issue of anomalous diffusion scalingsin certain areas of applied physics. Specifically, by proving that the fundamental solutionof the linear diffusion problem is the long-time asymptotic behavior (from below and fromabove) of solutions to any non-degenerate parabolic equation with strictly positive polyno-mial diffusivity, from arbitrary initial data, we are led to suggest that any “anomaly” inthe scalings might be an artefact of the short time-duration of an experiment. Ultimately,the physics must justify whether degeneracy of the governing equation is expected (or not)and why. Indeed, for degenerate non-linear parabolic equations, a wealth of different scalingfunctions and transformations (different from the one defined in Eq. (27)), leading to boundsdifferent from Eq. (28), are allowed (see, e.g., Barenblatt, 1996; Witelski and Bernoff, 1998).Non-degeneracy, on the other hand as we proved above, imposes strict restrictions on theasymptotic scaling behavior of a non-linear diffusion equation’s solutions.21nterestingly, the numerical example in § VII B suggests that our mathematical resultsalso hold for initial conditions and solutions that lack the exponential decay required by thebarrier function in Eq. (27). Therefore, in future work, it would be of interest to attempt,or to determine whether it is even possible, to generalize the proof in § VI using a barrierfunction with algebraic decay as | x | → ∞ .Finally, it should be noted that, in an interesting paper, Bricmont, Kupiainen, and Lin(1994) used the re-normalization group (RG) method (Wilson, 1983) to prove certain resultsabout the asymptotic behaviour of non-linear diffusion equations. Specifically, taking thedifference between the fundamental solution of the linear problem (Gaussian distribution)and the solution of the non-linear PDE in divergence form, they are able to prove convergenceat a point ( x, t ) = ( √ t, t ). The RG method is generic, and it does not use the maximumprinciple. On the other hand, our approach is based on maximum principle leading to anAronson estimate. Importantly, our approach is applicable to both non-linear equations inboth divergence and in non-divergence form; in the latter case, the coefficients may evendepend on the spatial variable. As a result, our approach leads to an accurate estimate,from both above and below, of the solution of the non-linear equation in terms of Aronson’sGreen function. That is to say, the estimates proved in this work cannot be obtained viathe RG method of Bricmont, Kupiainen, and Lin (1994). ACKNOWLEDGEMENTS
Acknowledgment is made to the donors of the American Chemical Society PetroleumResearch Fund for partial support of I.C.C., under ACS PRF award
DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the correspondingauthor upon reasonable request. 22
EFERENCES
Aronson, D. G., “Bounds for the fundamental solution of a parabolic equation,” Bull. Amer.Math. Soc. , 890–896 (1967).Barenblatt, G. I., “On some unsteady fluid and gas motions in a porous medium,” Prik.Mat. Mekh. (PMM) , 67–78 (1952), in Russian.Barenblatt, G. I., Scaling, Self-similarity, and Intermediate Asymptotics , Cambridge Textsin Applied Mathematics, Vol. 14 (Cambridge University Press, New York, 1996).Bernoff, A. J.and Witelski, T. P., “Stability and dynamics of self-similarity in evolutionequations,” J. Eng. Math. , 11–31 (2010).Bird, R. B., Stewart, W. E., and Lightfoot, E. N., Transport Phenomena , 2nd ed. (JohnWiley & Sons, New York, NY, 2002).Boon, J. P.and Lutsko, J. F., “Nonlinear diffusion from Einstein’s master equation,” EPL , 60006 (2007).Bricmont, J., Kupiainen, A., and Lin, G., “Renormalization group and asymptotics of solu-tions of nonlinear parabolic equations,” Commun. Pure Appl. Math. , 893–922 (1994).Brown, R., “A brief account of microscopical observations made in the months of june, julyand august 1827, on the particles contained in the pollen of plants; and on the generalexistence of active molecules in organic and inorganic bodies,” Phil. Mag. Ser. 2 , 161–173(1828).Cahn, D. S., Fuerstenau, D. W., Healy, T. W., Hogg, R., and Rose, H. E., “Diffusionalmechanism of solid–solid mixing,” Nature , 494–496 (1966).Celik, E., Hoang, L., Ibragimov, A., and Kieu, T., “Fluid flows of mixed regimes in porousmedia,” J. Math. Phys. , 023102 (2017).Christov, I. C.and Stone, H. A., “Resolving a paradox of anomalous scalings in the diffusionof granular materials,” Proc. Natl Acad. Sci. USA , 16012–16017 (2012).Dafermos, C. M., Hyperbolic Conservation Laws in Continuum Physics , 4th ed., Grundlehrender mathematischen Wissenschaften, Vol. 325 (Springer-Verlag, Berlin/Heidelberg, 2016).DiBenedetto, E.,
Degenerate Parabolic Equations , Universitext (Springer-Verlag, New York,NY, 1993).Dury, C. M.and Ristow, G. H., “Axial particle diffusion in rotating cylinders,” Gran. Matt. , 151–161 (1999). 23instein, A., “ ¨Uber die von der molekularkinetischen theorie der w¨arme geforderte bewegungvon in ruhenden fl¨ussigkeiten suspendierten teilchen,” Ann. Phys. (Leipzig) , 549–560(1905).Einstein, A., Investigations on the Theory of the Brownian Movement (Dover Publications,Mineola, NY, 1956) edited by R. F¨urth, Translated by A. D. Cowper.Fick, A., “On liquid diffusion,” Phil. Mag. Ser. 4 , 30–39 (1855a).Fick, A., “Ueber diffusion,” Ann. Phys. Chem. (Leipzig) , 59–86 (1855b).Fischer, D., Finger, T., Angenstein, F., and Stannarius, R., “Diffusive and subdiffusive axialtransport of granular material in rotating mixers,” Phys. Rev. E , 061302 (2009).Frey, E.and Kroy, K., “Brownian motion: a paradigm of soft matter and biological physics,”Ann. Phys. (Leipzig) , 20–50 (2005).Gardiner, C. W., Stochastic Methods: A Handbook for the Natural and Social Sciences , 4thed., Springer Series in Synergetics, Vol. 13 (Springer-Verlag, Berlin, 2009).Ilyin, A. M., Kalashnikov, A. S., and Oleynik, O. A., “Linear second-order partial differentialequations of the parabolic type,” J. Math. Sci. , 435–542 (2002).Kleinstein, G.and Ting, L., “Optimum one-term solutions for heat conduction problems,”Z. Angew. Math. Mech. (ZAMM) , 1–16 (1971).Lacey, P. M. C., “Developments in the theory of particle mixing,” J. Appl. Chem. , 257–268(1954).Landis, E. M., Second Order Equations of Elliptic and Parabolic Type , Translations of Math-ematical Monographs, Vol. 171 (American Mathematical Society, Providence, RI, 1998).Lenzi, E. K., Lenzi, M. K., Ribeiro, H. V., and Evangelista, L. R., “Extensions and solutionsfor nonlinear diffusion equations and random walks,” Proc. R. Soc. A , 20190432(2019).Oron, A., Davis, S. H., and Bankoff, S. G., “Long-scale evolution of thin liquid films,” Rev.Mod. Phys. , 931–980 (1997).Ottino, J. M.and Khakhar, D. V., “Mixing and segregation of granular materials,” Annu.Rev. Fluid Mech. , 55–91 (2000).Perkowitz, S., “Gedankenexperiment,” in Encyclopædia Britannica Online (EncyclopædiaBritannica, Inc., 2010).Philip, J. R., “Flow in porous media,” Annu. Rev. Fluid Mech. , 177–204 (1970).24istow, G. H.and Nakagawa, M., “Shape dynamics of interfacial front in rotating cylinders,”Phys. Rev. E , 2044–2048 (1999).Savage, S. B., “Disorder, diffusion, and structure formation in granular flow,” in Disorderand Granular Media , edited by A. Hansen and D. Bideau (Elsevier, Amsterdam, 1993) pp.255–285.Sekimoto, K.and Fujita, T., “Symmetry in self-similarity in space and timeshort time tran-sients and power-law spatial asymptotes,” Symmetry , 1489 (2019).Shampine, L. F.and Reichelt, M. W., “The MATLAB ODE suite,” SIAM J. Sci. Comput. , 1–22 (1997).Skeel, R. D.and Berzins, M., “A method for the spatial discretization of parabolic equationsin one space variable,” SIAM J. Sci. Stat. Comput. , 1–32 (1990).von Smoluchowski, M., “Zur kinetischen theorie der brownschen molekularbewegung undder suspensionen,” Ann. Phys. (Leipzig) , 756–780 (1906).Stokes, G. G., “On the effect of the internal friction of fluids on the motion of pendulums,”Trans. Cambridge Philos. Soc. , 8–106 (1851).Strikwerda, J. C., Finite Difference Schemes and Partial Differential Equations , 2nd ed.(Society for Industrial and Applied Mathematics, Philadelphia, PA, 2004).Sutherland, W., “A dynamical theory of diffusion for non-electrolytes and the molecularmass of albumin,” Phil. Mag. Ser. 6 , 781–785 (1905).Umbanhowar, P. B., Lueptow, R. M., and Ottino, J. M., “Modeling segregation in granularflows,” Annu. Rev. Chem. Biomol. Eng. , 129–153 (2019).V´azquez, J. L., The Porous Medium Equation: Mathematical Theory (Oxford UniversityPress, Oxford, UK, 2007).Wilson, K. G., “The renormalization group and critical phenomena,” Rev. Mod. Phys. ,583–600 (1983).Witelski, T. P.and Bernoff, A. J., “Self-similar asymptotics for linear and nonlinear diffusionequations,” Stud. Appl. Math. , 153–193 (1998).Zel’dovich, Y. B.and Raizer, Y. P., Physics of Shock Waves and High-Temperature Hydro-dynamic Phenomena , Vol. II (Academic Press, New York, 1967) chapter X § , 371–405 (2018).Zik, O.and Stavans, J., “Self-diffusion in granular flows,” Europhys. Lett.16