{'fraction_non_alphanumeric': 0.03168836411343223, 'fraction_numerical': 0.007720254314259764, 'mean_word_length': 4.685312679288583, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tensorial mobilities for accurate solution of transport problems in models with diffuse interfaces

The general problem of two-phase transport in phase-field models is analyzed: the flux of a conserved quantity is driven by the gradient of a potential through a medium that consists of domains of two distinct phases which are separated by diffuse interfaces. It is shown that the finite thickness of the interfaces induces two effects that are not present in the analogous sharp-interface problem: a surface excess current and a potential jump at the interfaces. It is shown that both effects can be eliminated simultaneously only if the coefficient of proportionality between flux and potential gradient (mobility) is allowed to become a tensor in the interfaces. This opens the possibility for precise and efficient simulations of transport problems with finite interface thickness.


Introduction
Phase-field models have recently enjoyed a rapidly growing popularity as a compact and elegant simulation tool for moving boundary problems in such diverse fields as solidification [1,2], fluid dynamics [3] or solid-state transformations [4,5]. Their technical advantage resides in the implicit representation of interfaces by one or several phase fields, i.e. fields that are defined in the entire space, take constant values within the bulk of each domain, and exhibit smooth but steep variations in well-localized interfacial regions. The tedious procedure of front tracking is avoided by introducing equations of motion for the phase fields that are coupled to the relevant transport variables. The price to pay for this advantage is the introduction of a new length scale into the model: the thickness W of the phase-field front. For a given macroscopic problem, simulations with a phase-field model yield in general results that depend on the value of W .
In the field of crystal growth, great progress towards efficient and precise simulations has been made by reducing this dependence on the interface thick-ness [6,7,8,9,10]. This was made possible by a detailed analysis of the model equations using the method of matched asymptotic expansions, which is a systematic procedure to calculate the effective boundary conditions "seen" by the macroscopic transport field. Since this analysis is carried out within a perturbation approach, these boundary conditions are naturally expressed as a power series in W . Within the phase-field community the limit W → 0 is referred to as the sharp-interface limit. When the corrections due to the finite interface thickness are taken into account for choosing the model parameters, the accuracy of the phase-field method can be drastically improved. This procedure, which has been called thin-interface limit [6], has so far been worked out only for a few specific physical systems.
To be more precise, let us consider the problem of solidification, in which the relevant transport process that limits the growth of the crystal is the diffusion of heat and/or solute. Two cases are completely solved: the symmetric model, where the diffusion coefficients in the two phases are identical [6], and the onesided model, in which no diffusion takes place within the solid phase [8,9]. However, so far no method has been found to eliminate all thin-interface effects in the case of arbitrary diffusion coefficients in the two phases, despite some recent progress [7,11,2] (for a more detailed discussion, see [12]).
As will be pointed out here, part of this problem arises from the fact that for a truly two-sided model (with different diffusivity in each phase) even the stationary transport problem without interface motion exhibits thin-interface effects. This prevents a solution of the problem by the antitrapping approach, which has been successful for the one-sided model [8,9] and for the two-sided model with vanishing diffusion current in one phase [11].
In fact, such thin-interface effects are fairly general and arise in a whole class of problems, namely, two-phase transport through a complex structure. Examples of relevant physical situations are the conduction of electrical current or heat through a two-phase material with different conductivities, the magnetic flux through a two-phase material with different susceptibilities, or fluid flow though a porous medium with variations in permeabilities. At first sight, the advantages of using the phase-field method in these cases are less obvious since the interfaces do not move and therefore the problem of front tracking does not arise. However, the geometrical representation of complex-shaped surfaces can be difficult even without interfacial motion, especially in three dimensions. Additionally, the use of a stationary phase-field function makes it possible to prescribe a given boundary condition at the interface in a straightforward manner [13,14], and to impose arbitrary boundary conditions at the border of a physical domain of complex shape, see [15,16] and references therein.
The problem of representing a complex interface through a diffuse boundary gains additional relevance because the use of tomographic methods for structure determination becomes more and more widespread. In such methods, the structure representation takes the form of a matrix consisting of discrete pixels (or voxels in three dimensions) that contain binary or intensity data indicating whether a point in space is "filled" or "empty". From these data, the "true" structure (represented for example by discrete sharp surface elements) has to be reconstructed by image analysis techniques. The phase-field method is an easy and robust method to obtain a smoothened representation of such data [17,18]. It could be interesting to use directly this smoothened structure for accurate calculation of transport processes instead of going through the additional steps of determining the "sharp" surface geometry. However, it will be shown below that thin-interface effects are also present in the "simple" problem of two-phase transport, even without interface motion. Therefore, these effects must be quantified and if possible eliminated. As we will demonstrate below, two effects that depend on the interface thickness are present: transport along the surface, and an interface resistance. In the standard phase-field formulation, where the transport is described by a scalar coefficient whose value depends on the phase field, these two effects cannot be eliminated simultaneously. In contrast, if the transport coefficient is allowed to become a tensor inside the diffuse interfaces, there are enough degrees of freedom in the model to eliminate both effects.

Problem formulation
All the problems listed above have a common structure, namely, the flux of a conserved quantity is driven by a potential gradient, where V is a potential. The structure of Eq. (1) is standard in out-of-equilibrium thermodynamics: a linear relationship between a flux and a thermodynamic driving force (the potential gradient). The mobility coefficient M (φ) depends on the phase field φ. If the transport process is electric conduction, V is the electrostatic potential and M the conductivity; for diffusive mass transport V is the chemical potential and M the atomic mobility.
The transported quantity satisfies a conservation law, which is valid both in the bulk and at the surfaces, and for a time-independent solution (steady-state flow) reads ∇ · j = 0.
The problem specification is completed by a boundary condition for the potential at the interfaces. We assume continuity of the potential, where V + and V − are the values of the potential when the interface is approached from the two sides. This corresponds to a rapid exchange of the transported quantity between the two sides of the interface. Since we consider a fixed and immobile two-phase structure, the phase field is independent of time. We will assume that the two constant values that designate the two phases are φ = 0 and φ = 1, and that the field φ varies between these two limits continuously through a front region of width W . For the sake of concreteness, the reader may have in mind a sigmoid function such as [1 + tanh(x/W )] /2, but the explicit form of this function is not important. The only hypothesis we make is that for a straight interface, the profile of φ is odd with respect to the point φ = 1/2, that is, φ(x) = 1 − φ(−x) for an interface centered at x = 0. This is the case in all standard phase-field models.

Surface current
For simplicity of exposition, it is useful to focus on a concrete example. Consider the conduction of electric current through a two-phase material. Then, V is the electrostatic potential, and M (φ) is the phase-dependent electric conductivity. Furthermore, consider a straight interface normal to the x direction, centered at x = 0. A potential gradient along the y direction (along the interface) is imposed by sandwiching the material between two parallel plates located at ±L/2 that are held at constant potentials ±U . Since the phase field φ and hence the conductivity M are constant along any line of constant x (although interpolation type Table 1: Coefficients of the second order regression J(ǫ) = J 0 + c 1 ǫ + c 2 ǫ 2 for the three different interpolation methods. These coefficients are obtained from a quadratic regression to the data indicated by the symbols of Fig. 1.
their values differ for different values of x), Eq. (2) yields a constant potential gradient U/L directed along the y direction. Therefore, the total current J that flows between the two plates is given by Since we have considered a sample that extends to infinity, this current is clearly infinite. However, we will be concerned only with the excess of this current with respect to the sharp-interface value. The latter is obtained as the current that would flow if φ(x) was a step function, that is, the space between the two plates is filled with material 1 of conductivity M 1 for x < 0, and with material 2 of conductivity M 2 for x > 0. This yields The difference between these two expressions is the excess of current δJ due to the variation of conductivity over a zone of finite thickness. This excess is localized in the interface, and can therefore be interpreted as a additional current along the surface. It can be written as the product of the potential gradient and a surface conductivity This surface transport coefficient has two obvious properties: (i) for an interface profile of fixed functional form, φ(x) = f (x/W ), M s is proportional to the interface thickness (as can be shown by a simple change of variables), and (ii) it is strictly zero for any value of W if In this case, the excesses of current on the two sides of the interface exactly compensate. For a phase-field profile that satisfies φ(−x) = 1 − φ(x), this can be simply achieved by choosing This will be called direct interpolation in the following.

Surface resistance
Let us now analyze again a planar interface normal to the x direction, but this time crossed by a steady current J ⊥ along x. In this case, the continuity equation immediately yields that the current is constant (independent of x). Then, the potential V satisfies the simple equation Integration along x yields where V is the potential at x = 0 (an integration constant). In contrast, if the interface was sharp, the potential would simply be given by for x > 0 and x < 0, respectively. Of course, outside the diffuse interfaces, the slopes of V (x) are identical in Eqs. (11) and (12). Therefore, the asymptote of the diffuse-interface expression is of the form for x → ±∞. The constants V + and V − (the interface potentials "seen" from the region outside the diffuse interface) are readily obtained from the matching of this expression to Eq. (11), Of particular interest is the fact that these surface potentials can be different, in contradiction to the assumption of Eq. (3). The difference δV = V + − V − can be written as the product of the current J ⊥ and an interface resistance This resistance is often called Kapitza resistance and has been frequently observed in experiments and simulations [19,20,21,22,23]. Again, it is obvious that it is proportional to the interface thickness, and that it vanishes if the integral is exactly zero. This can be achieved for any value of W by the inter- This will be called inverse interpolation in the following.

Tensorial mobility
In summary, the two interface effects (surface current and surface resistance) can each be eliminated by a specific choice for the interpolation function of the mobility. Since these interpolations are mutually exclusive, it seems as if necessarily one of the two effects must remain nonzero. However, the current is a vector quantity, and the two effects are linked to distinct components of the current vector: the excess surface conductivity is relevant only for the components parallel to an interface, whereas the surface resistance modifies the boundary conditions for the normal component. Inside the two phases, where each medium is isotropic, the Curie principle requires to choose a scalar mobility to relate the current and the potential gradient. However, in the presence of an interface, isotropy of space is broken and a tensorial transport coefficient is permitted. In- deed, the gradient of the phase field can be readily used to define the interface normal n = ∇φ/|∇φ|, which provides a second direction that is independent of the potential gradient. Then, we can define the transport coefficient by where 1 is the unit tensor, with two independent interpolation functions M ⊥ (φ) and M || (φ). If we interpolate M ⊥ according to Eq. (16) and M || according to Eq. (9), both thin-interface effects are eliminated. Hence, for this interpolation the transport problem defined by Eq. (1) becomes with components (in two dimensions) Here, we have designated by M ij the elements of the symmetric tensor M (φ). The simple calculations developed above are valid only for planar interfaces. However, for a sufficiently smooth interface (that is, with a local radius of curvature R satisfying R ≫ W ), a local curvilinear coordinate system can be defined in which the above relations remain valid at least up to second order in ǫ = W/R. It should be mentioned that a similar strategy has been used recently to develop efficient phase-field models for surface diffusion [24].
whereas the tensorial interpolation leads to the following equation Details about the discretization and the method used to solve these equations can be found in the appendix. The thin-interface effects arising in the spherical inclusion problem are quantified by measuring the total flux at x = 1, and by plotting J as a function of the ratio ǫ = W/R between the interface width W and the radius of the disk R, see Fig. 1. As shown in Fig. 1, the three interpolations converge to the same value of J 0 when ǫ → 0. The estimation of J 0 is given by a quadratic regression in ǫ of the simulation data. The coefficients of the regressions are listed in Table 1. Within the truncation error [that is O h 2 ∼ 10 −6 ] the three interpolations give the same value of J 0 = 1.752750 ± 10 −6 , as shown in the second column of Table 1.
In addition, we have estimated the rate of convergence of J(ǫ) to the sharpinterface limit J 0 for the three interpolations of the mobility. As shown in Fig.  2, the direct and inverse interpolations converge only linearly with ǫ to this limiting value, whereas the tensorial interpolation suppress linear thin-interface effects, which leads to a convergence that is almost quadratic in ǫ.
For further illustration, we consider a second geometrical configuration formed by two half-disks of radius R placed in a rectangular domain with the same boundary condition of the spherical inclusion problem, see Fig. 3. We choose a small value of the ratio between M inside and outside the half-disks, e.g. M 2 /M 1 = 0.1, and in Fig. 3 we plot the streamlines for this configuration. In this configuration, the flux is constricted in the narrow space between the two half disks. The difference between the three interpolations is largest for streamlines which are locally almost tangent to the half-disks; note the divergence of the different lines close to the tips of the half-disks. As shown by Eq. (8), choosing M according to the direct interpolation cancels exactly the surface conductivity M s for a flat interface. Therefore, thin-interface errors affecting streamlines which are parallel to the boundary of variation of the transport coefficient are more pronounced in case of inverse interpolation (red dashed lines) than in case of direct interpolation (black dash-dotted lines), see Fig. 3. In Fig. 4, we show again the values of the total current as a function of the ratio ǫ, and again the tensorial interpolation performs much better than the other two, as expected.

Conclusion
We have investigated a phase-field model with a mobility tensor, in which normal and parallel components of the flux are interpolated with distinct functions of the phase field. Contrary to phase-field models with scalar mobilities, this method makes it possible to eliminate at the same time the additional surface diffusion and the surface resistance, which are both linked to the finite thickness of the interfaces. This opens the possibility to perform accurate simulations of two-phase transport problems with enlarged interface thickness, which can lead to dramatic savings in computation time.
The complete elimination of thin-interface effects to first order in the interface thickness is also a prerequisite for the development of a quantitative crystal growth model for arbitrary ratio of the diffusion coefficients in the two phases. Indeed, a tensorial diffusion coefficient can remove both surface diffusion and the Kapitza resistance [12] from such models. Even if a second order asymptotic analysis of a time dependent version of the tensorial problem (17) has been recently performed [26], this methodology is still not capable of removing all thin-interface effects from phase-field models with a non-stationary φ. The obstacle that needs to be overcome for the successful development of such a model is to find a coupling of the transport equation to an evolution equation for the phase field that yields the correct boundary conditions for the transport field at moving interfaces. We hope to be able to report on this problem in the near future.
where h −1 △ + k is the standard forward difference operator acting on the k direction. Finally, Eq. (21) is discretized by employing the backward difference Evidently, the discretization of the elliptic problem arising from the tensorial interpolation M is more involved. In fact, two parts of Eq. (22) mix x and y derivatives, i.e. ∂ x (M xy ∂ y V ) and ∂ y (M xy ∂ x V ). In order to guarantee second order accuracy and maximum compactness of the discrete stencil it is convenient to place the M xy component of the tensor M on a third grid, whose nodes are shifted by h/2 in the two directions [28]. The tensorial interpolation produces an additional contribution to j x and to j y As before, by applying backward differentiation we obtain the discrete version of Eq. (22) 1 The MAC arrangement ensures a second order accuracy of the truncation error of the two elliptic problems [27,28]. Eq. (30) and (33) are two linear systems of equations where the unknowns are the values of the potential V i,j . These linear systems can be easily solved through any iterative method, for example by the Successive Over Relaxation (SOR) [29]. To find an accurate solution of these problems we iterate the SOR algorithm until the maximum residue at the nodes of the primary grid is comparable with the truncation error, i.e. ρ ≤ h 2 .