Double diffusivity model under stochastic forcing

TL;DR

This paper enhances the double diffusivity model by incorporating stochastic effects into the internal length gradient framework, improving the alignment of theoretical predictions with experimental observations in nanopolycrystals.

Contribution

It introduces stochasticity into the ILG-based diffusion model, bridging the gap between deterministic theory and experimental relaxation times.

Findings

Stochastic ILG model aligns better with experimental relaxation times.

Inclusion of boundary layer fluctuations improves model accuracy.

Provides a benchmark for gradient-based continuum models in fluctuating environments.

Abstract

The "double diffusivity" model was proposed in the late 1970s, and reworked in the early 1980s, as a continuum counterpart to existing discrete models of diffusion corresponding to high diffusivity paths, such as grain boundaries and dislocation lines. Technically, the model pans out as a system of coupled {\it Fick type} diffusion equations to represent "regular" and "high" diffusivity paths with "source terms" accounting for the mass exchange between the two paths. The model remit was extended by analogy to describe flow in porous media with double porosity, as well as to model heat conduction in media with two non-equilibrium local temperature baths e.g. ion and electron baths. Uncoupling of the two partial differential equations leads to a higher-ordered diffusion equation, solutions of which could be obtained in terms of classical diffusion equation solutions. Similar equations could also be derived within an "internal length" gradient (ILG) mechanics formulation applied to diffusion problems, {\it i.e.}, by introducing nonlocal effects, together with inertia and viscosity, in a mechanics based formulation of diffusion theory. This issue becomes particularly important in the case of diffusion in nanopolycrystals whose deterministic ILG based theoretical calculations predict a relaxation time that is only about one-tenth of the actual experimentally verified timescale. This article provides the "missing link" in this estimation by adding a vital element in the ILG structure, that of stochasticity, that takes into account all boundary layer fluctuations. Our stochastic-ILG diffusion calculation confirms rapprochement between theory and experiment, thereby benchmarking a new generation of gradient-based continuum models that conform closer to real life fluctuating environments.

Equations

$$\frac{\partial {\rho }_1}{\partial t}=D_1{\nabla }^2{\rho }_1-{\kappa }_1{\rho }_1+{\kappa }_2{\rho }_2+{\eta }_1(\mathbf{x},t)$$

$$\frac{\partial {\rho }_2}{\partial t}=D_2{\nabla }^2{\rho }_2+{\kappa }_1{\rho }_1-{\kappa }_2{\rho }_2+{\eta }_2(\mathbf{x},t).$$

$$⟨{\eta }_i(\mathbf{x},t)⟩$$

$$⟨{\eta }_i(\mathbf{x},t){\eta }_j({\mathbf{x}}^{\prime },t^{\prime })⟩$$

$${\partial }_t\rho +\tau {{\partial }_t}^2\rho =D{\nabla }^2\rho +c^{\ast }D{\partial }_t{\nabla }^2\rho +cD{\nabla }^4\rho .$$

$${\partial }_t\rho +\tau {{\partial }_t}^2\rho =D{\nabla }^2\rho +c^{\ast }D{\partial }_t{\nabla }^2\rho +cD{\nabla }^4\rho +\eta (\mathbf{x},t),$$

$${\partial }_t\rho =D{\nabla }^2\rho ,$$

$$\mathcal{M}\left(\begin{array}[]{c}\hat{\rho}_{1}\\ \hat{\rho}_{2}\\ \end{array}\right)=\left(\begin{array}[]{c}\hat{\eta}_{1}\\ \hat{\eta}_{2}\\ \end{array}\right),$$

$$\mathcal{M}=\left(\begin{array}[]{cc}-i\omega+D_{1}k^{2}+\kappa_{1}&-\kappa_{2}\\ -\kappa_{1}&-i\omega+D_{2}k^{2}+\kappa_{2}\\ \end{array}\right).$$

$$\psi (\mathbf{x},t)=\int d^d\mathbf{k}\int d\omega \hat{\psi }(\mathbf{k},\omega )e^{(ik.x-\omega t)},$$

$$\left(\begin{array}[]{c}\hat{\rho}_{1}\\ \hat{\rho}_{2}\\ \end{array}\right)=\mathcal{M}^{-1}\left(\begin{array}[]{c}\hat{\eta}_{1}\\ \hat{\eta}_{2}\\ \end{array}\right).$$

$$\mathcal{M}^{-1}=\left(\begin{array}[]{cc}M_{1}&M_{2}\\ M_{3}&M_{4}\\ \end{array}\right),$$

$$M_1$$

$$M_2$$

$$M_3$$

$$M_4$$

$${{\rho }_1}^{\text{rms}}$$

$${{\rho }_2}^{\text{rms}}$$

$${{\rho }_{12}}^{\text{rms}}$$

$${({{\rho }_1}^{\text{rms}})}^2$$

$${({{\rho }_2}^{\text{rms}})}^2$$

$$⟨{\rho }_1(\mathbf{k},\omega ){{\rho }_1}^{\ast }(-\mathbf{k},-\omega )⟩$$

$$⟨{\rho }_2(\mathbf{k},\omega ){{\rho }_2}^{\ast }(-\mathbf{k},-\omega )⟩$$

$$<\rho_{1}({\bf k},\omega)\>{\rho_{2}}^{*}(-{\bf k},-\omega)>=\bigg{[}\frac{2\gamma_{1}\kappa_{1}(D_{2}k^{2}+\kappa_{2})+2\gamma_{2}\kappa_{2}(D_{1}k^{2}+\kappa_{1})}{\zeta_{n}}\bigg{]}+i\bigg{[}\frac{2\gamma_{1}\omega^{2}\kappa_{1}+2\gamma_{2}\omega^{2}\kappa_{2}}{{\zeta_{n}}}\bigg{]}.$$

$$\Omega_{0}\approx\frac{1}{2^{1/4}}{\big{[}2D^{2}k^{4}+4Dk^{2}(\kappa_{1}+\kappa_{2})+2{(\kappa_{1}+\kappa_{2})}^{2}\big{]}}^{1/4}.$$

$${{\rho }_1}^{\text{rms}}$$

$${{\rho }_2}^{\text{rms}}$$

$${[{\rho }_i^{\text{rms}}(\mathbf{r})]}^2$$

$${\rho }_1^{\text{rms}}(\mathbf{r})$$

$${\rho }_2^{\text{rms}}(\mathbf{r})$$

$${{\rho }_1^{\text{rms}}(\mathbf{r})}_{k\to \text{large}}$$

$${{\rho }_2^{\text{rms}}(\mathbf{r})}_{k\to \text{large}}$$

$${\rho }_1^{\text{rms}}(\mathbf{r})$$

$${\rho }_2^{\text{rms}}(\mathbf{r})$$

$${[{{\rho }_i}^{\text{rms}}(T)]}^2$$

$$\displaystyle{\big{[}\rho_{1}^{\text{rms}}(T)\big{]}}^{2}$$

Full text

I Introduction

The subject of material science has traditionally dealt with hard and soft matter based objects, typically at observable macroscopic, and mesoscopic, and microscopic scales. Until very recently, most branches of traditional science used to be expressible within such length scales for which appropriate phenomenological theories have been well established over time. Atomistic or lattice based approaches have also been developed and related to molecular dynamics and quantum mechanical simulation codes have been advanced. However, the regime between the microscopic and atomic scale, i.e., the nanoscale regime, especially the length scale between 5 and 100 nm, still remains a major modeling challenge.

In particular, the advent of nanocrystalline and ultrafine grain nanocomposite materials, along with their tremendous possibilities in material engineering implementations have challenged all existing conventional and well established theoretical and experimental realizations. The cause of this can be easily attributed to the fact that at manometric scales, most materials, both soft and hard, often have widely different physicomechanical and thermochemical properties to their macroscopic counterparts that prove elusive to existing scientific analyses.

A simple case in hand is that of diffusion in nanocrystals, and more generally, in nanocomposites, for which it has been shown that the diffusivity could be many orders of magnitude larger than the diffusion constant for bulk lattice diffusion of mesoscopic and microscopic scale materials Aifantis2011 . This is in line with other mechanical, electrical, magnetic and chemical properties of nanoscopic objects where nanoscale moduli differ significantly from their macroscale and microscale equivalents. A compromise between the atomistic and conventional continuum mechanical engineering modeling approaches for nanopolycrystals is reviewed in Aifantis2011 , within a Laplacian based internal length gradient (ILG) generalization of classical deformation (elasticity, plasticity) and diffusion theories. A more thorough analysis, including size effects in chemomechanics and electromechanics, as well as the role of stochasticity, due to internal stress fluctuations, can be availed in detail from Aifantis2016 . This article refrains from explicit discussion of the role of stochasticity in diffusion at nanoscales, a task earmarked for the present article.

Nanoscale diffusion or diffusion in nanopolycrystals in particular, is viewed as a two-phase process, one through the bulk and the other across the grain boundary (GB) space. In the case representing an admixture of two families of grain boundaries (low and high GBs), depending on the proportion of the two types of GBs, it may be necessary to account for a third type of diffusion, a variant of the paradigmatic stochastic nanodiffusion akc-eca-2016 dynamics. For decreasing grain size, the density of the triple grain boundary junction (TJ) increases significantly and then, TJ has also to be identified as a separate family of high-diffusivity paths. In any case, independently of how one will identify the types of the families of paths available for diffusion, the mathematical model that we discuss here distinguishes between only two types of paths, fast and slow. These two paths with varying rate kinetics are differentiated by assigning two different diffusivities, thereby allowing mass transfer of diffusion species between them. For example, in the case of very small grain sizes ($\sim$10 nm), where a large number of triple junctions is present, the two types of diffusion may be identified with GB and TJ, while bulk diffusion may be neglected as being comparably at a much slower rate. For severe plastic deformation (SPD) fabricated polycrystals with larger grain sizes ($\sim$100 nm), the two types of high diffusivity paths may be identified with low-angle (equilibrium) and high-angle (nonequilibrium) grain boundaries which are densely populated with dislocation/disclination defects. For microscale polycrystals with even larger grain sizes ($\sim$1 $\mu$m), the two families of diffusion paths are represented by the bulk/grain interior and the surrounding/grain boundary space. With the above interpretation, we focus on modeling stochastically enhanced counterpart of the deterministic continuum model proposed and analyzed in Aifantis1979JApp ; Aifantis1979Acta ; Aifantis1980Quarter ; Aifantis1980Acta ; Aifantis1980Engg ; Beskos1980 ; Aifantis1980Mech ; Konstantinidis1998 ; Konstantinidis1999 ; Konstantinidis2001 ; Forest2010 ; Kalampakas2012 .

More specifically, we elaborate on the coupled system of the two partial differential equations proposed in Aifantis1979JApp , solved in Aifantis1980Quarter , and thereafter extended in Aifantis1980Acta ; Aifantis1980Engg to explain double porosity in media and through Beskos1980 ; Aifantis1980Mech to describe heat transport in materials with two temperatures. All of these theoretical enclaves were eventually implemented in Konstantinidis1998 ; Konstantinidis1999 ; Konstantinidis2001 to interpret experimental measurements of diffusion in polycrystals and nanocrystalline aggregates. Some thermodynamic aspects of the deterministic double diffusivity model have been relatively recently discussed in Forest2010 , where the approach is probability based. Phenomenological coefficients have recently been provided in Kalampakas2012 .

From the above discussion, it follows that a vital contributing factor that has remained unattended in the legion of double diffusivity studies is that of thermally induced fluctuations arising out of the structural difference between lattice and grain boundary spacings, especially close to the boundary layers as also arising out of material imperfections and structural randomness. In line with the well established legion of stochastically forced flow models Forster1977 ; DeDominicis1979 ; Yakhot1986 ; Chattopadhyay2000 representing archaic dynamical randomness generated close to the sheared boundary layers, as also due to structural imperfections, like crack propagation, this article will explore this realistic limit of double-diffusion, thereby accounting for all modes of randomness. This will be structured within the well-knit Langevin formulation of stochastic dynamics Risken ; Barabasi . Phenomenologically, this can be seen as an external stochastic force that is randomly redistributing the relevant spatial structure, for example, that of the high-diffusivity paths, in which a stochastic increase (decrease) in temperature extends (contracts) the interlayer grain boundary distance between two nanosized grains and thus alters the configuration of the structural defects (dislocations, disclinations). Such a multi ensembled stochastic reorganization of the lattice space distribution is sometimes known to create new universality classes Barabasi ; Cross1993 as well.

II The Models

Crack-fracture propagation through wavefronts and diffusion in complex media have long been known to be complementary physical realizations that can be addressed using a combination of continuum mechanics and numerical simulation Barabasi . ILG, on the other hand, has been established as a powerful theoretical tool to address stress-strain deformation aspects in materials, as well as in related thermomechanical, electromechanical and chemomechanical processes as also in complex heterogeneous media Aifantis2011 ; Aifantis2016 ; Aifantis2014 . The role of randomness and some combined deterministic gradient-stochastic models have been considered in Barabasi ; Aifantis2016 ; Aifantis2014 for “higher-ordered” deformation models. In contrast, the role of stochasticity has not been considered for higher-ordered diffusion models such that the double-diffusivity equations or the ${4}^{\text{th}}$ ordered diffusion equations resulting from their uncoupling, following on from the ILG mechanics.

Our theoretical structure will combine two models separately in which Model 1 will define the spatiotemporal dynamics of each individual phase with respect to the concentration of the diffusing species in the phase concerned. With reference to Aifantis2011 , where the deterministic version of this model referred to bivariate dynamics, we will call this the “non-conserved double-diffusive model” (NDD). Model 2, on the other hand, will focus on the higher-ordered diffusion equations resulting from uncoupling the NDD equations which also hold for the total concentration, i.e., the sum of the two concentrations in the “slow” and “high” diffusivity paths. Once again, keeping in mind the origin of its deterministic analog, we will call this the “conserved double-diffusive model”.

A key outcome of this analysis is the fundamentally different dynamical structure functions (two-point spatiotemporal functions) estimated for these two classes of models. This is highly non-trivial, and a consequence of symmetry violation due to the stochastic forcing term in these models, since the deterministic versions of these two classes of models converged to the same “universality class”.

II.1 Forced non-conserved double-diffusion (NDD) model

In line with the narrative presented earlier, we will arrive at the first stochastic gradient nanomechanics (SGNM) model of nanodiffusion starting from a deterministic Fick’s diffusion like model (Model 1, as referred to in Aifantis2011 ). Defining $\rho_{1}$ and $\rho_{2}$ as the concentrations of the diffusing species, for example, in the intercrystalline (IC) and TJ space, respectively, the stochastically forced Fick’s diffusion model can be written as

[TABLE]

In Eqs. (1a) and (1b) above, $D_{1}$ and $D_{2}$ refer to the diffusion constants in the IC and TJ phases respectively, with $\kappa_{1}$ and $\kappa_{2}$ denoting the respective formation or depletion rate of diffusive substances at these phases. A negative sign before $\kappa_{i}\>\rho_{i}$ ($i$=1, 2) would indicate depletion, i.e., loss of diffusion species “jumping” to the other “phase” while a positive sign would indicate replenishment of that particular phase. In the given form above, the deterministic concentration abides a certain conservation law (detailed below) with $\eta_{i}({\bf x},t)$ (i=1,2) being the white thermal noise perturbing this dynamics and is defined as follows

[TABLE]

where $d$ represents the spatial dimension and $\gamma_{i}$ ($i$=1,2) are the noise strengths corresponding to $\eta_{1}$ and $\eta_{2}$ respectively.

II.2 Forced conserved couble-diffusion (CDD) model

Using a series of transformations, $\tau=\dfrac{1}{\kappa_{1}+\kappa_{2}}$, $D=\tau(\kappa_{1}D_{2}+\kappa_{2}D_{1})$, $c^{*}=\tau(D_{1}+D_{2})/D$ and $c=-\tau D_{1}D_{2}/D$, as detailed in Aifantis2011 , the deterministic model corresponding to Eqs. (1a) and (1b) (for the noiseless case $\eta_{i}=0$) can be easily shown to abide by a conservation law in a scaled variable set defined as follows:

[TABLE]

This fourth-order equation containing, in addition to the classical Fick’s law, a second time derivative inertial or “telegrapher” term ($\tau\nabla^{2}\rho$); a third order mixed spatiotemporal or pseudoparabolic term ($c^{*}D\partial_{t}\nabla^{2}\rho$), and a fourth order spatial or biharmonic term ($cD\nabla^{4}\rho$) is obtained by uncoupling Eqs. (1a) and (1b) with $\eta_{2}=0$. It holds for both individual concentrations $\rho_{1}$ and $\rho_{2}$, as well as for its sum $\rho=\rho_{1}+\rho_{2}$, i.e., the total concentration. The boundary layer fluctuation and random structural imperfection perturbed stochastic model can be derived from the deterministic model defined in Eq. (3) as a stochastically forced model with an additive (uncorrelated white) noise. The corresponding root-mean-square spatiotemporal “width” will define the spatiotemporal evolution of the interface separating the two phases. The resultant model (Model 2, as referred to in Aifantis2011 ) is defined below

[TABLE]

where $\rho({\bf x},t)$ is the mass density of the separation width of the two phases, in which $\tau=\frac{1}{\kappa_{1}+\kappa_{2}},\>D=\tau(\kappa_{1}D_{2}+\kappa_{2}D_{1}),\>c^{*}=\tau\frac{(D_{1}+D_{2})}{D},\>c=-\tau\frac{D_{1}D_{2}}{D}$, and $\eta({\bf x},t)$ is the stochastic fluctuation.

For the special case for which $D_{1}=D_{2}=D$ in Eqs. (1a) and (1b) give

[TABLE]

that is a simple diffusion equation in the variable $\rho=\rho_{1}+\rho_{2}$. One must note that both spatiotemporal and reflection conservation as inherent to Eq. (4) will be lost in the presence of a non-zero noise which is our starting model. Instances of such general applications of double-diffusion model in analyzing heterogeneous growth process have also been explored Showalter2004 ; Klein2012 .

In the following sections, we will separately analyze the single phase dynamics of the respective IC and TJ concentrations (Model 1) by estimating autocorrelation functions of the variables $\rho_{1}$ and $\rho_{2}$ respectively, from Eqs. (1a) and (1b). These individual autocorrelations will then be compared against the autocorrelation evaluated from Eq. (4) that defines the dynamics of the total concentration. The following sections will then estimate the spatiotemporal dynamics of the interface in presence of noise from estimation of the spatiotemporal correlation functions. This later part of the analysis will be compared against separate “thin” and “thick”-film conditions as detailed in Aifantis2011 to establish the importance of the stochastic contribution.

III Phase Evolution Dynamics of the NDD Model

The primary focus in this section will be the evaluation of the individual phase dynamics of the IC and TJ phase variables $\rho_{1}$ and $\rho_{2}$ of the NDD model. Our starting point here will be the Fourier transformation of the core model presented in Eqs. (1a) and (1b), defined by variables $({\bf k},\omega)$ of Eqs. (1a) and (1b) to arrive at the following matrix form

[TABLE]

where the matrix $\mathcal{M}$ is defined as

[TABLE]

Here $\hat{\rho}_{1}$ and $\hat{\rho}_{2}$ represent the Fourier transformed version of the $(\rho_{1},\rho_{2})$ variables in the $({\bf k},\omega)$ space (3+1 dimensional) while the Fourier transforms themselves abide by the following generic form

[TABLE]

where $\psi$ generically represents either $\rho_{i}$ or $\eta_{i}$. From Eq. (6), we get

[TABLE]

This inverse of the matrix $\mathcal{M}$ takes the form

[TABLE]

where

[TABLE]

As like in any stochastically driven model DeDominicis1979 ; Chattopadhyay2000 , in the stochastically ensemble averaged state, the measurables will be the (Brownian) root-mean-square averaged quantities of their deterministic equivalents. The corresponding mean energy dissipation rate in such an ensemble-averaged state is thus defined as the time averaged kinetic term $\frac{1}{2}{\langle\partial_{t}\rho_{i}({\bf x},t)*\partial_{t}\rho_{i}({\bf x},t)\rangle}$ ($i$=1, 2) that necessitates evaluation of the following two autocorrelation functions:

[TABLE]

and a complementary set of cross-correlation functions given by

[TABLE]

In the above, the superscript “rms” stands for the root-mean-squares of the respective quantities under consideration post the stochastic (Brownian) average that is indicated by the “$<.>$” sign, while $\rho_{i}^{*}$ is the complex conjugate of $\rho_{i}$. The cross-correlations of the quantities with their complex conjugates emphasize the importance of the attenuation term in the dynamics (complex quantity); the fluctuation-dissipation theorem is always implicitly assumed in such analyses. In this work, we will assume noise cross-correlation to be zero (that is uncorrelated) and hence $\rho_{12}^{\text{rms}}=0$. The rms quantities are the ones of our interest, as this is what an experimental measurement will see, an allusion to the classical Brownian dynamics Barabasi .

Eqs. (9a) and (9b) can be explicitly written as

[TABLE]

III.1 Phase autocorrelation and crosscorrelation

Equation (8) can be solved to obtain the autocorrelation functions in the $k-\omega$ space:

[TABLE]

where the quantity ${{\zeta_{n}}}=(\omega^{4}+D_{2}^{2}k^{4})[\omega^{4}+{(D_{1}k^{2}+\kappa_{1})}^{2}]+2[\kappa_{1}\omega^{4}+D_{2}k^{2}(\omega^{4}+D_{1}k^{2}(D_{1}k^{2}+\kappa_{1}))]\kappa_{2}+(\omega^{4}+D_{1}^{2}k^{4})\kappa_{2}^{2}$ defines the pole structure and hence possible discontinuities in the spectral dynamics.

The Fourier transformed cross-correlation is even more interesting in that it shows a complex form in which the complex part represents attenuation. The structure looks as follows

[TABLE]

For $\gamma_{1}=\gamma_{2}$, as shown in Figure 1, the $\rho_{1}$ and $\rho_{2}$ plots merge with each other.

The following section alludes to the derivation of the spectral dynamics in the Fourier transformed space that could be easily encapsulated with the descriptions of the autocorrelation and cross-correlation functions, the former representing the “self-energy” structure while the latter relates to the “spatiotemporal” interaction of sites.

In order to solve for both auto and cross-correlation functions, we will need to evaluate Eqs. (12b) and (13) around the poles which are defined through the eighth degree polynomial equation ${{\zeta_{n}}}(\omega)=0$, where the roots of the equation will be given by $\omega=\Omega_{i}$ ($i$ = 1, 2,…, 8). These poles are given in details in Appendix I. The nature of the complexity of the analytical structure can be gauged even from the simplified special case of $D_{1}=D_{2}$, represented by the equation ${{\zeta_{n}}}={(\omega+i\Omega_{0})}^{4}{(\omega-i\Omega_{0})}^{4}=0$, where

[TABLE]

Within the ambits of this assumption ($D_{1}=D_{2}=D$), we arrive at the root-mean-squared autocorrelation functions for spatial dimension $d=3$ as follows

[TABLE]

It is interesting to note that for $D_{1}=D_{2}$, the deterministic versions of Eqs. (1a) and (1b) involving $\rho_{1},\>\>\rho_{2}$ obey Eq. (4) for a special form of the constants while $\rho_{1}+\rho_{2}$ still abides by Fick’s law; however, if $D_{1}\neq D_{2}$, then both $\rho_{1},\>\rho_{2}$, and also $\rho_{1}+\rho_{2}$ obey Eq. (4). The answer lies in the average energy dissipation rates of each of the two phases (IC and TJ) in the thermally driven system; these rates are equal to each other with exact values calibrated against the system parameters ($D,\kappa_{1},\kappa_{2}$) involved. We should emphasize here that the simplifying assumption ($D_{1}=D_{2}=D$) used in arriving at the above result in no way sullies the implication of this analysis. For all realistic experimental observations Konstantinidis1998 ; Konstantinidis1999 concerning double diffusivity, including its application in explaining the oxygen diffusivity in barium superconductors Konstantinidis2001 , the two diffusive constants typically differ by about three orders of magnitude whose exact correlation forms can be analyzed using the representations in the Appendix. Figure 1 uses identical noise strengths $\gamma_{1}=\gamma_{2}=1$.

For anisotropic noise, the equivalent representation is provided in Fig. 2.

III.2 Spatial Correlation of Phases

In this section, our attention will be focused on evaluating how the concentration of each phase changes with spatial distance in the dynamical equilibrium limit for spatial dimension $d=3$. Mathematically, this implies evaluation of the respective spatial correlation functions of each phase for all times and then taking ensemble averages over all noise realizations.

By definition, we have

[TABLE]

where $J_{0}(x)$ represents the zeroth order Bessel function for the scalar variable $x$ and $i$=1, 2. The above Eq. (16) leads to the following density expression for both phases as given below

[TABLE]

where $k_{0}$ and $k_{m}$ refer to the minimum and maximum of the wave vector $k$ (measured as the inverse of the characteristic system length). In our numerical solution, as shown in Fig. 4, we have chosen $k_{0}=0$ and $k_{m}=10000$. The plots shown in Fig. 3 remain mostly unaffected by the specific choice of $k_{0}$ and $k_{m}$. This can be easily seen from a cross-check of the integration kernels given in Eqs. (17a) and (17b) in the limit $k\to$ large, as given below:

[TABLE]

Within a very short interval, contributions to the correlation functions from such large-spatial separations can be seen to decay to zero as shown in Fig. 3. This confirms the convergence of the integrals in Eqs. (17a) and (17b). Figure 3 validates the convergence of Eq. (17a); a similar analysis could be repeated for Eq. (17b) in the limit of $k\rightarrow$ large to establish a similar convergence. It must be remembered, though, that $\rho_{1}^{\text{rms}}=\rho_{2}^{\text{rms}}$ is a result of our assumption $D_{1}=D_{2}$; the correlations will have different values for the different phases for $D_{1}\neq D_{2}$, the precise nature of which can be estimated from the expressions presented in the Appendix. Once again, for the isotropic case ($\gamma_{1}=\gamma_{2}$), as shown in Fig. 4, the $\rho_{1}$ and $\rho_{2}$ plots merge with each other.

In order to get a feel for the functional dependence of the spatial correlation function for the “nondegenerate” case $D_{1}\neq D_{2}$, the relevant correlation functions can be approximately shown to be as follows (details in the Appendix)

[TABLE]

where $D_{\text{av}}=\frac{D_{1}+D_{2}}{2}$ and \Omega_{\text{av}}=-\dfrac{\bigg{(-T_{1}-k^{2}D_{1}(2\kappa_{1}+T_{2})-(\kappa_{1}+\kappa_{2})(\kappa_{1}+\kappa_{2}+T_{2})-k^{2}D_{2}(2\kappa_{2}+T_{2})\bigg{)}}^{1/4}}{2^{1/4}}, with $T_{1}=k^{4}(D_{1}^{2}+D_{2}^{2})$ and T_{2}=\sqrt{{\bigg{(}k^{2}(D_{1}-D_{2})+\kappa_{1}\bigg{)}}^{2}+2\bigg{(}k^{2}(D_{2}-D_{1})+\kappa_{1}\bigg{)}\kappa_{2}+\kappa_{2}^{2}}.

For the special case of $\gamma_{1}=\gamma_{2}=1$, comparing with the “degenerate” case $D_{1}=D_{2}={10}^{-9}$, we arrive at a very similar functional behavior for the non-degenerate case as well.

In Figure 4, we show how the spatial correlation saturates with increasing separation distance, for the case $\gamma_{1}=\gamma_{2}$. In order to portray the situation for the case of anisotropic noise ($\gamma_{1}\neq\gamma_{2}$), below we plot Fig. 5. While the qualitative features remain unchanged, due to a large $\gamma_{2}=2\gamma_{1}$, the saturation level of $\rho_{2}^{\text{rms}}$ can be seen to be way above $\rho_{1}^{\text{rms}}$, although the crossover point remains roughly unchanged. This feature clearly suggests that noise anisotropy is not a qualitatively devolving feature of this dynamics.

Unlike the deterministic case (as in Aifantis2011 ), the dynamical equilibrium of the thermal noise driven two-phased system stabilizes to the same spatial concentration spread for both phases. This clearly suggests a difference at the qualitative level, as well as obvious quantitative differences (compared to Aifantis2011 ).

III.3 Temporal correlation of phases

It is well known that a fundamental consideration in multi-phase systems is the time evolution of the interface separating two different phases. Often such systems are known to be stochastically perturbed and hence nonequilibrium in nature, potentially rendering the relevant dynamics as oscillatory, or with oscillatory-rotatory instability leading to chaos Hu2003 . Phase control through synchronization driving such systems away from the chaotic bifurcation point has in fact benchmarked the hare-lynx model in ecology Blasius2000 .

By analogy with two-phase systems, we calculate below the theoretical quantities which, in principle, can be compared with the experimental set ups describing the temporal correlation dynamics. In the present case, the relevant variables necessary to model such a stochastically driven two-phase system are ${{\rho_{i}}^{\text{rms}}(T)}$, which are encapsulated in the following equations

[TABLE]

for $i=1,2$. Our next target is to estimate the two-point temporal correlation functions, as the first-order euphemism of the probability density function (also connected to the fluctuation-dissipation theorem Barabasi ).

Calculating as before, we obtain these two-point correlation functions as follows

[TABLE]

A remarkable feature of the temporal correlation function plot as shown in Fig. 7 is the relatively low effect of the noise (an)isotropy compared to the spatial cases. The other aspect of these temporal correlation functions shown in Eqs. (21a) and (21b) is the tremendous stability with regard to noise fluctuations. In our simulations, we sampled across a wide range of noise strengths ${10}^{-2}<D<{10}^{-9}$ to find that the results shown in Fig. 7 remain unaffected by the value of $D$, as long as the stability condition $k_{\text{min}}>\sqrt{\dfrac{D}{\lambda_{2}}}$ is obeyed, where $k_{\text{min}}$ is the minimum allowed value of the wave vector $k$.

In the analysis above, we have deliberately refrained from going into the details of the cross-correlation description, both for the spatial as well as the temporal cases. While qualitatively the presence of an imaginary part in the correlation function indicates attenuation, the subject will be dealt with in more detail separately.

IV Nanodiffusion Spatiotemporal Dynamics

In order to compare nanodiffusion transport between the present thermally driven model and the paradigmatic Aifantis model Aifantis2011 , we will now calculate the spatiotemporal dynamics of the concentration fields $\rho_{1},\>\rho_{2}$ and $\rho=\rho_{1}+\rho_{2}$, based on Model 2, starting from the conserved SGNM model previously defined in Eq. (4). The corresponding solution of the phase concentrations as shown in Aifantis2011 will be compared with its thermally driven counterpart for periodic boundary conditions. As always, in the context of stochastic Brownian-type statistics, individual dynamical variables give way to their corresponding r.m.s. counterparts. Starting from Eq. (4) and using $\lambda_{1}=c^{*}D$ and $\lambda_{2}=cD$, we can rewrite the model as

[TABLE]

Fourier transformation of the above Eq. (22) in the ${\bf k}-\omega$ space gives

[TABLE]

The above Eq. (23) has 4 poles at $\omega=\omega_{i}$ ($i$ = 1, 2, 3, 4), such that $\omega_{i}=\pm\sqrt{T_{1}\pm\dfrac{\sqrt{T_{2}}}{2\tau^{2}}}$, in which $T_{1}=-\frac{1}{2\tau^{2}}+\frac{Dk^{2}}{\tau}-\frac{\lambda_{1}k^{2}}{\tau^{2}}-\frac{\lambda_{1}^{2}k^{4}}{2\tau^{2}}-\frac{\lambda_{2}k^{4}}{\tau}$, $T_{2}={(1-2Dk^{2}\tau+2\lambda_{1}k^{2}+\lambda_{1}^{2}k^{4}+2k^{4}\tau\lambda_{2})}^{2}-4\tau^{2}(D^{2}k^{4}-2D\lambda_{2}k^{6}+\lambda_{2}^{2}k^{8})$. To simplify calculations through a reduced model, we study the case for an overdamped system where $\tau\to 0$.

In this limit, the two poles turn out as $\pm i\Gamma$, where $\Gamma=\dfrac{k^{2}(D-\lambda_{2}k^{2})}{1+\lambda_{1}k^{2}}$. This gives autocorrelation function as

[TABLE]

where $\gamma_{0}$ is the strength of the Gaussian white noise.

A comparison between Eqs. (15a) and (15b), and Eq. (24) allows us to compare the quantitative difference between individual concentrations with individually added noise as against their total concentration with an overall added noise from the perspective of energy dissipation rate. This will be separately evaluated through the spatial autocorrelation function and the temporal autocorrelation function [definitions as in equations (16), and (21a), and (21b)] that are defined as per Amit1 ; Amit2 ; Amit3 as follows

[TABLE]

Figure 8 shows a periodic stabilizing pattern which is distinctly different from either of Fig. 4 or Fig. 5. This is most remarkable since the deterministic description does not show any qualitative difference between the two phase model and its equivalent single phase description Konstantinidis1998 . The temporal description (Fig. 9) too shows qualitatively distinctive features compared to Fig. 7. Fig. 9 indicates a sharp rise towards a saturation concentration density, a case of finite-sized saturation effect, as opposed to the decaying pattern represented in the anisotropic case in Fig. 7. In order to magnify the quantitative impact, we used a larger noise value $D=10^{-1}$ in arriving at these plots; however, the choice of the actual value does not impact the qualitative feature in any way, as the model remains remarkably stable to noise perturbations [detailed earlier after Eqs. (19a) and (19b)]. In fact, the quantitative corrections are less than 1% with every order change in noise [using formulas in Eqs. (25a) and (25b)], reconfirming the noise amplitude independence of the dynamics. Here we must indicate, though, that while the dynamics is largely unaffected by changes in the noise amplitude, the noise distribution function is expected to be vitally important to the dynamics, a feature that is presently being studied for future publications.

This reduced model has a cut-off at $k=\sqrt{\frac{D}{\lambda_{2}}}$ which defines its validity regime. We postpone analysis of the full model, includinga non-zero $\tau$, as defined in Eq. (22) for a later work.

V Conclusion

In this article, we have provided an initial analysis of the continuum double diffusivity model under stochastic forcing. Stochasticity is introduced separately as a white thermal noise, either in a Fick class of equations describing the concentrations in each type of diffusion path, or in a higher-order equation for the total concentration. Comparisons between the two cases have been made by comparing the respective classes of spatiotemporal cross-correlations and autocorrelations.

The autocorrelation plots in Figs. 1 and 2 conform qualitatively to predictions based on the deterministic double diffusivity models Konstantinidis1998 ; Konstantinidis1999 ; Konstantinidis2001 ; although the amplitudes are higher due to additional energy inputs through stochastic forcing. It is pertinent to remember that the root-mean-squared forms of the respective autocorrelation forms are the dimensional equivalents of the corresponding quantities in the deterministic models in Aifantis2016 ; Konstantinidis1998 ; Konstantinidis1999 ; Konstantinidis2001 and hence could be compared on a term-by-term basis. As to the crosscorrelation terms in the spatiotemporal dynamics (Figs. 3 to 7), the stochastic model ushers in a new regime of description where stochasticity mediates off-diagonal, often asymmetric forcing across multiple variables even at the first Gaussian approximation order. While this is very much an expected part of real life nanodynamic processes, the deterministic double-diffusivity model failed to capture this aspect that we have successfully made now. The results are verifiable using experimental data.

An important aspect of this analysis is the relative independence of both spatial and temporal correlation functions to the stochastic fluctuations. Over a wide range of noise strengths (${10}^{-9}<D<{10}^{-2}$), the correlation functions showed no qualitative change and very little quantitative change, thereby confirming the stability of this model to noise perturbations. This indirectly explains why some past theories Konstantinidis1998 ; Konstantinidis1999 ; Konstantinidis2001 have managed to arrive at experimental results reasonably accurately for some cases while faltering in others.

In summary, we point out some interesting features that stochasticity brings into the double diffusivity model, comparing existing deterministic terms against the relative stochastic forcing. This is only a first step toward integrating the double diffusivity properties with real life fluctuations that could modulate the process. As this model and the corresponding higher-order diffusion equation have been shown to effectively interpret transport in heterogeneous media possessing more than one family of conduction paths, as well as transport phenomena at the nanoscale, more detailed analysis focusing on the precise nature of randomness, e.g., forcing through non-Gaussian noise, as well as also delving deeper into the microscopic dynamics of the process, will be pursued in the next set of publications.

VI Acknowledgments

The combined support of ERC (Belgium)-13 and ARISTEIA II projects funded by GSRT (Greece) of the Green Ministry of Education, sponsoring A.K.C’s visits to Thessaloniki, is gratefully acknowledged.

VII Appendix I:

Poles of the autocorrelation function

The poles of the autocorrelation function as defined in Eq. (12b) can be obtained as follows:

[TABLE]

where $T_{1}=k^{4}(D_{1}^{2}+D_{2}^{2})$ and T_{2}=\sqrt{{\bigg{(}k^{2}(D_{1}-D_{2})+\kappa_{1}\bigg{)}}^{2}+2\bigg{(}k^{2}(D_{2}-D_{1})+\kappa_{1}\bigg{)}\kappa_{2}+\kappa_{2}^{2}}.