Asymptotic Expansion Method with Respect to Small Parameter for Ternary Diffusion Models

Ternary diffusion models lead to strongly coupled systems of PDEs. We choose the smallest diffusion coefficient as a small parameter in a power series expansion whose components fulfill relatively simple equations. Although this series is divergent, one can use its finite sums to derive feasible numerical approximations, e.g. finite difference methods (FDMs).


Introduction
The subject of divergent series dates back to the distant past and is linked with the names of outstanding mathematicians such as Euler, Poincaré, Borel, Padé, and Birkhoff. each of the PDEs depends on the gradient of all other variables, so these models exhibit the phenomenon of crossdiffusion. This makes the PDE system strongly coupled. In addition, the coefficients depend on the unknown functions and the parabolicity of this system is conditional, namely it can be proved only if solutions remain in an admissible set. The invariance of this set is connected to mass conservation laws. Any violation of this property leads to serious analytical and computational problems. We overcome these problems by constructing conservative numerical schemes, i.e. preserving a discrete version of mass. In spite of mathematical and numerical difficulties, three dimensional interdiffusion problems are still of interest. Interdiffusion becomes an important topic in electrochemistry and biochemistry, e.g. molecular channels (nano-channels) [15].

Preliminaries
Let Ω be a bounded domain in ℝ n with a smooth boundary Ω. The closure of Ω is denoted by Ω . In this manuscript, we investigate the system of diffusion equations on [0, T] × Ω, where D 1 , D 2 , D 3 > 0 and D is the drift velocity given by The system (1) is considered with the initial conditions and the Neumann boundary conditions where denotes the outward normal vector to the boundary Ω.
Our analysis of the strongly coupled system (1) starts from the simplest case w ≡ 0. Consider two chemical substances and their normalized densities u, v ≥ 0, u + v ≡ 1 described by the system of differential equations It is clear that the assumptions on initial functions u 0 , v 0 ≥ 0 and u 0 + v 0 ≡ 1 imply u, v ≥ 0, u + v ≡ 1. Indeed, consider new dependent variables p, q defined as linear combinations of the functions u, v: Then we obtain the following system Since u 0 + v 0 ≡ 1 it is seen that p ≡ 1 solves the first equation of system (6). Based on the relation v ≡ 1 − u we obtain from (4) partial differential equation Therefore, the implication 0 ≤ u 0 ≤ 1 ⇒ 0 ≤ u ≤ 1 follows easily from the maximum principle.
There are various numerical techniques for approximations of solutions of the ternary diffusion models (1)-(3). In particular, we refer to the papers [16,17]. The authors of [16] analyze the method of lines for diffusion equations. The paper [17] deals with iterative methods for this problem. The authors of both papers study stability and convergence of the constructed method in L 2 and in the Sobolev space W 1,∞ . In this paper, as a tool for approximation of solutions to the initial boundary value problems for systems of equations obtained by an asymptotic expansion, we use a finite-difference methodology. In the construction of difference schemes we apply general ideas presented in [18,19].

Small Parameter Method
Suppose that initial functions u 0 , v 0 , w 0 ∈ C 2 (Ω, ℝ) in (2) satisfy the relation u 0 + v 0 + w 0 = 1 for x ∈Ω. In the sequel we assume that the diffusion coefficients D 1 , D 2 , D 3 in the system (1) satisfy the relation D 1 > D 2 > D 3 > 0. If u, v, w satisfy (1)-(3), then u + v + w ≡ 1 and the system (1) reduces to the system of two equations It is more convenient to use the same change of variables u, v to p, q as in (5): For these functions we obtain the system of differential equations where : = D 3 . Observe that if = 0 then (7) is equivalent to (6). We expand the functions p and q in the series with respect to the small parameter : .
Once we substitute these expansions to the system (7) and compare the coefficients we get the recurrence relation for k ≥ 0. We represent the initial conditions for the unknown functions p [k] and q [k] as follows: It is easy to see that If we compare the right-hand sides of (10) and (8), we obtain the initial conditions for the unknown coefficients p [k] , q [k] for k ≥ 1: Using a small parameter expansion of functions p and q we replace the original system of equations (7) by a separate systems of equations (9) for each pair of functions (p [k] , q [k] ), k ≥ 0. The complexity of these systems is rising with increasing values of k. Therefore, the expansions (8) are practically limited to several terms. According to this we are mainly interested in the asymptotic nature of this truncated expansion for → 0, ≪ 1, then the possible divergence of the whole series is in general not important.
The study of asymptotic expansions was first made by Poincaré. According to Poincaré's definition, a series ∑ ∞ n=0 p [n] (t, x; ) n is said to be asymptotic to a function p(t, x; ), i.e.
iff for every N ≥ 0 and → 0, ≪ 1, we have This definition is equivalent to the following conditions Taking the first limit under consideration, we have become solutions of the same initial boundary value problem which leads to the conclusion that the limit is equal to 1 for every N. The same result we obtain for the expansion of q(t, x; ). Therefore, the expansions (8) are asymptotic if only ≪ 1. Even if the series do not converge, the asymptotic expansions give a good approximation of solutions if we take under consideration only few terms of these expansions with a very small value of [11,20]. This conclusion is confirmed by our numerical analysis in one dimensional case presented in Example 1. It turns out that our method is very effective for multidimensional couples, even in the case of complicated geometries.

Numerical Analysis
This section is devoted to a numerical analysis of solutions to the problem (9), (10), obtained by the small parameter expansion method. We construct a conservative difference scheme and present the numerical behavior of ternary solid solutions based on this method.

Finite Difference Method
We discretize the system of partial differential equations (9). Consider a uniform partition of We convey that h will represent an approximation to the exact value of the func- where 1 is the jth coordinate. We employ standard and linear difference Denote by ∇ h the discrete equivalent of the gradient opera- (2) jj then we define With these conventions, the finite-difference scheme will be given by the following discrete system with the discrete initial conditions The discrete Neumann boundary conditions have the form for i ∈ {0, … , K}.

Numerical Simulations
In this section we explore the capability of the constructed method (11)-(13) to provide good approximations to the exact solution of the PDE problem (1)-(3). We present some examples of tracer and up-hill diffusion in ℝ 1 and tracer diffusion in ℝ 2 . Examples in one-dimensional case are similar to those presented in [16]. Therefore, one can compare the results of the method of lines obtained in [16] with the results obtained by our numerical scheme. We also provide more complicated two-dimensional cases which are absent from the literature.
Example 1 (Tracer at interface) Let n = 1. In the numerical analysis we take L = 1 and the diffusion parameters D 1 = 0.18, D 2 = 0.08, D 3 = . The initial distribution u 0 is a small perturbation of the Heaviside function and v 0 is a symmetry of u 0 about the ordinate axis. The solution w is initiated by a smoothed small peak where A denotes the characteristic function of a set A, and for x ∈ [ − 1, 1]. We analyze the results for the PDE system (1) with the initial distributions given by (14) with respect to the small parameter . Illustrative simulations are provided with: h 0 = 2 ⋅ 10 −5 and h = 2 ⋅ 10 −2 . Based on formula (5) we can derive We consider expansions (8) for k = 0, 1, 2: Numerical results for the respective u [0] , v [0] , w [0] calculated by (15) are presented in Fig. 1. Note that we obtain an expected behavior of densities u, v, w regardless of the value selection. • k = 1, two terms in each expansions: p = p [0] + p [1] , q = q [0] + q [1] . In Fig. 2 we observe the distributions of u [1] , v [1] , w [1] computed by (15) with different values of . Figure 3 presents the behavior of partial sums u [0] + u [1] , v [0] + v [1] , w [0] + w [1] with different values of . Although oscillations of u [1] , v [1] , w [1] increase with decreasing values of , we obtain reasonable approximations of u, v, w by means of partial sums u [0] + u [1] ,v [0] + v [1] , w [0] + w [1] , respectively. Notice that these partial sums take values in the interval [0, 1] for all . • k = 2, three terms in each expansions: [2] . Figure 4 presents the behavior of u [2] , v [2] , w [2] computed by (15). Numerical results of partial sums [2] are presented in Fig. 5. Similarly as in the case k = 1, it can be observed that the oscillations of u [2] , v [2] , w [2] increase. Unlike the case k = 1, the partial sums for u, v, w take values in [0, 1] for sufficiently small . If = 0.01, the partial sums for u, v are not monotone. In Fig. 6 one can observe additional results for partial sums [2] at times t = 0.1 and t = 0.2 with fixed = 0.001. We present these results to illustrate the dynamics of approximate distributions.
We conclude that the number of terms in the expansions (8) is closely related to the size of the parameter . Once we use more terms in the expansions, we need to decrease the value of to get reasonable approximation of u, v, w by means of respective partial sums, i.e. taking values in the interval [0, 1] and preserving monotonicity of partial sums for u, v.
Example 2 (Up-hill diffusion) We consider the onedimensional case and we set L = 1 similarly as in Example 1. We take D 1 = 0.2, D 2 = 0.08, D 3 = 0.001 = . We put some monotone functions as initial distributions of substances u, v and we assume a constant initial distribution w: where Results of numerical simulations for partial sums [2] computed with h 0 = 2 ⋅ 10 −5 , h = 2 ⋅ 10 −2 are given in Fig. 7(b) at time t = 0.1. Figure 7a presents initial functions u 0 , v 0 , w 0 .
No numerical analysis for tracer diffusion in twodimensional case have been done so far in the literature because of the computational complexity, not to mention any efficient comparison of numerical results with reallife experiments. Due to our approach we constructed a numerical method which allows to perform computations  [1] , v [1] , w [1] for different values of at time t = 0.8 Fig. 3 (Tracer at interface) the case k = 1: partial sums u [0] + u [1] , v [0] + v [1] , w [0] + w [1] for different values of at time t = 0.8 Fig. 4 (Tracer at interface) the case k = 2: functions u [2] , v [2] , w [2] for different values of at time t = 0.8  [2] , v [0] + v [1] [2] , w [0] + w [1] + 2 w [2] for different values of at time t = 0.8 even for multi-dimensional cases. We present two examples of numerical simulations for tracer diffusion in ℝ 2 .
In Fig. 9a we observe approximate distributions of partial sums for u, v, w with k = 2 and the step sizes h 0 = 2 −15 , h = 2 −3 at time t = 0.2. Figure 9b presents enlargement of w.
In the above example initial distribution of u and v is symmetrical with respect to the plane x = 0. Now we present two-dimensional example where the symmetry of initial distribution is broken.  [2] , v [0] + v [1] [2] , w [0] + w [1] + 2 w [2] at time t = 0.1 where and are the same as in Example 3. We perform simulations with parameters: the step sizes h 0 = 2 −15 , h = 2 −3 and the length of series expansions k = 2. In Fig. 11(a) the approximate distributions of partial sums for u, v, w at time t = 0.2 are presented. Figure 11b shows enlargement of w.
We conduct a brief analysis of marker positions in time, see [21]. By 'marker' we understand a small amount of substance placed in the diffusion couple, this substance moves according to locally not-balanced diffusive fluxes. Recall that diffusive fluxes are proportional to gradients of u, v, w. In the above tracer diffusion examples the marker is represented by the substance w. Referring to Example 1 the trajectory of the marker position is described by the ODE initial value problem where D is the drift velocity. Since the diffusion coefficient D 3 is small, it can be omitted, then the solution of the Cauchy problem (16) can be approximated by the solution of the following ODE problem To solve the problems (16) and (17) numerically, we first use a polynomial fitting for discrete functions u, v, w with respect to spatial variable. Then approximate solutions to the above ODE problems are computed by a second-order Runge-Kutta method. The results of computations with step sizes: h 0 = 2 × 10 −5 , h = 2 × 10 −2 are illustrated in Fig. 12. A mesh refinement leads to better accuracy, i.e. the solution to problem (17) gives a very good approximation

Conclusions and Remarks
Before we close this paper with concluding remarks we would like to mention other applications of asymptotic expansion methods, in addition to those given in the introduction, to point out the importance of these techniques in a different branches of knowledge. Some applications can be observed in the fluid and solid mechanics theory. In [22] the asymptotic expansion method was used as a tool to derive a new two-dimensional shallow water model from the time-averaged non-dimensional Navier-Stokes equations. Asymptotic simplification procedures for linear and nonlinear wave propagation problems that contain large parameters are examined in [20]. A divergent asymptotic series found applications also in the boundary layer theory [23] and in the nonlinear wave theory [24]. We use the asymptotic expansion techniques to the ternary diffusion problem. Modeling of interdiffusion phenomena is not yet unified, and various methods differ from each other in the arbitrary choice of the reference velocity, i.e. the drift or convection velocity. In computational solid mechanics and physics it is defined as being equal to the mass average velocity [13]. In gases and fluids the drift is defined based on the volume average velocity [25]. It is still common to neglect drift (convection) and assume Fickian diffusion [26]. Such simplified approach to the mass transport was used there to a three-dimensional model of interdiffused quantum dots. The above methods do not allow considering all interdiffusion effects [26]. The use of widely accepted Onsager method [27,28] is narrow by very limited data on transport coefficients. In material sciences and chemistry the drift definition bases on the Darken