# Modeling in food across the scales: towards a universal mass transfer simulator of small molecules in food

O. Vitrac, Murielle Hayert

► To cite this version:

O. Vitrac, Murielle Hayert. Modeling in food across the scales: towards a universal mass transfer simulator of small molecules in food. SN Applied Sciences, 2020, 2 (9), pp.1509. 10.1007/s42452-020-03272-2. hal-02915283

HAL Id: hal-02915283

https://hal.inrae.fr/hal-02915283

Submitted on 10 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Modeling in food across the scales: towards a universal mass transfer simulator of small molecules in food

O. Vitrac¹ Ⓒ · M. Hayert¹

Received: 20 December 2019 / Accepted: 28 July 2020 / Published online: 13 August 2020

© The Author(s) 2020 OPEN

## Abstract

Multiscale modeling in food is the cutting-edge strategy to revisit food structure and food composition to meet specific targets such as bioavailability, oral perception, or to evaluate the contamination of food by chemicals. A special implementation of Langevin dynamics is proposed to describe mass transfer in structured food. The concepts of random walks over discrete times and physicochemical interactions are connected via an exact solution of the Fokker-Planck equation across interfaces. The methodology is illustrated on the calculation of effective diffusivities of small solutes in emulsions in relationship with their polydispersity, the volume fraction of dispersed phase $d = [0.1, 0.4]$, the ratio of diffusion coefficients between the two phases, $r_D = [10^{-2}, 10^2]$, and the partition coefficients between the continuous and disperse phases, $K = [10^{-2}, +\infty[$. Simulated diffusion paths are detailed in 2D emulsions and the effective diffusivities compared with the core-shell model of Kalnin and Kotomin (J Phys A Math Gen 31(35):7227-7234, 1998). The same effects are finally tabulated for 3D emulsions covering the full range of food applications. The methodology is comprehensive enough to enable various extensions such as chemisorption, adsorption in the surfactant layer, local flows, flocculation/creaming.

**Keywords** Multiscale modeling · Simulation · Langevin dynamics · Mathematical homogenization · Diffusion coefficient · Emulsion · Statistical physics

## List of symbols

$\{\alpha\}_{i=1..n}$ Reciprocal correlation times in Eq. (1) (SI unit: s⁻¹)

**A** $n \times n$ matrix associated with low frequency deterministic motions (SI unit: s⁻¹)

**B** $n \times n$ unitriangular matrix associated with the classical decomposition of **A** (–)

$\delta t$ Time step (SI unit: s)

$d$ Globule density or globule volume fraction (–)

$D_0$ Diffusion coefficient in a homogeneous medium (SI unit: m² s⁻¹)

$D^c$ Scalar diffusivity in the continuous phase (SI unit: m² s⁻¹)

**D**$^c$ $n \times n$ diffusion tensor in the continuous phase (SI unit: m² s⁻¹)

$D^d$ Scalar diffusivity in the dispersed/discontinuous phase (SI unit: m² s⁻¹)

**D**$^d$ $n \times n$ diffusion tensor in the dispersed/discontinuous phase (SI unit: m² s⁻¹)

$\emptyset$ Dimensionless diameter of the globule (–)

$F^\ddagger$ Free energy at the transition state

$F_c$ Free sorption energy in the continuous phase at constant volume (SI unit: J)

$F_d$ Free sorption energy in the dispersed/discontinuous phase at constant volume (SI unit: J)

$F_o$ Fourier number (–)

**I** $n \times n$ identity matrix (–)

$k_B$ Boltzmann constant (SI unit: J K⁻¹)

$k_{c \to d}$ Frequency of the transition $c \to d$ (SI unit: s⁻¹)

$k_{d \to c}$ Frequency of the transition $d \to c$ (SI unit: s⁻¹)

$K = \frac{p_{\text{eq}}^d}{p_{\text{eq}}^c}$ Partition coefficient between the dispersed and the continuous phase (–)

$l$ Hopping length (SI unit: m)

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<footnote>O. Vitrac, olivier.vitrac@agroparistech.fr | ¹UMR 0782 SayFood Paris-Saclay Food and Bioproduct Engineering Research Unit, INRAe, AgroParisTech, Université Paris-Saclay, 91300 Massy, France.</footnote>
$l_0$ Reference length scale (SI unit: m)

$\Lambda$ $n \times n$ diagonal matrix coding for the diffusion tensor (SI unit: $m^2 s^{-1}$)

$msd$ Mean-square displacement (SI unit: $m^2$)

$M$ Grid step to place globules (-)

$\mu_{\emptyset}$ Mean dimensionless diameter (-)

$\mu_{ex}^c$ Excess chemical potential of one single solute (SI unit: J)

$n$ Number of dimensions (-)

$\xi$ $n \times 1$ random vector (-)

$p_{eq}^c$ Probability density at equilibrium in the continuous phase (-)

$p_{eq}^d$ Probability density at equilibrium in the dispersed/discontinuous phase (-)

$r_D = \frac{D^c}{D^d}$ Ratio of diffusivities between the continuous and the dispersed/discontinuous phases (-)

$\mathbf{r}$ Position (SI unit: m)

$s$ Sampling factor in Eq. (12) (-)

$\sigma_{\emptyset}$ Dimensionless standard deviation of globule diameters (-)

$t$ Time (SI unit: s)

$T$ Absolute temperature (SI unit: K)

$\Phi$ Volume ratio of concentric regions in Eq. (14) (-)

# 1 Introduction

Mass transfer of small molecules plays a significant role in evaluating the shelf-life, the flavor, the mouthfeel, and the intestinal absorption of food constituents. The presented approach focuses on organic solutes with sensorial, nutritional, or sanitary interest, such as flavors, antioxidants or chemical contaminants originating from food contact materials. Similar descriptions could also apply to water vapor in dry foods. The state of the art offers a broad range of calculation and simulation methods to estimate diffusive fluxes in food [1] with a particular food geometry or to evaluate diffusivities [2–6] or thermodynamic properties [7–12] from molecular structures. There is, however, a lack of methodologies to describe the effect of food structure and the spatial organization of phases at intermediate scales above the molecular scale and below the scale where the effective medium approximation can be assumed. Although Fickian diffusion is strictly limited to ideal mixtures [13], continuous approximations averaged over a representative volume are physically acceptable at the cost of proper identification of driving forces and transport coefficients [14]. Such descriptions are commonly used to simulate mass transfer of water and organic solutes in complex food products or materials [15–19], but they fail to capture the effects of temperature (thermodynamic effects), the evolution of the microstructure during

processing, ingestion or digestion, and more generally for any fractal arrangement of phases. As examples, they are poorly adapted when the dispersed phase is colloidal and is moving to the continuous phase by Brownian motion, sedimentation or aggregation, or when the mean free path of molecules is comparable to the length scale of the details, or when the slowdown of random walks is controlled by rare traps [20]. Particle methods and other meshless methods offer an alternative by avoiding space discretization and Lagrangian formulation. For fluids, smoothed particle hydrodynamics [21–24] and moving particle semi-implicit scheme [25] are the most popular. Alternative techniques of integration of kinematic equations on intricate spatial scales comprise high order spectral approximations [26], lattice gas models [27] and kinetic Monte-Carlo [28, 29]. Describing mass transfer of solutes instead of advection requires, however, defining particles commensurable to molecules. They carry properties of real molecules, including mass, charge, position, and velocity. Representing all molecules of a food portion will remain inaccessible for several decades, but efficient and accurate descriptions can be recovered by following the trajectories of a smaller collection of solutes governed by Langevin dynamics (LD). The approach is very generic, and this study emphasizes on the correct implementation of thermodynamical concepts at interfaces such as local thermodynamical equilibrium, detailed balance, time and microscopic reversibility. A distinction between all concepts has been reviewed in Gorban [30]. In LD simulations, interactions with surrounding molecules are captured via an implicit solvent model and a proper stochastic equation, enabling the incorporation of temperature effects. In an infinite medium, the representation is equivalent to the Fokker-Planck equation for the evolution of the probability density of a population of molecules and hence similar to the convection-diffusion equation. No bias is therefore expected. This work illustrates how such a general methodology can be devised and implemented. With this respect and without a loss of generality, a simple case-study is presented. It aims at calculating (in few minutes or hours on a multicore processor) the effective solute diffusivities in 2D and 3D emulsions near their percolation threshold.

The chosen problem captures enough complexity—the coupling between transport and thermodynamics around complex interfaces while enabling visualization of 2D trajectories. The formulation shows how properties broadly available in the literature, such as intrinsic diffusion coefficients, activity or partition coefficients, association/disassociation coefficients with macromolecules, can be mapped onto LD simulations. There is additionally significant literature on the subject to be used for validation and discussion. Earlier results, including the Maxwell-Garnett mixing
rule, have been recapitulated by Crank [31]. Subsequent works [32, 33] considered obstacle effects (e.g., impenetrable globules), but neglected the partitioning of the diffusing species between the continuous phase and the globules. Ronis [34] studied the effect of partitioning for a continuous phase separating two globules and evidenced the role of trapping in globules on effective diffusivities. The repeated ring approximation failed, however, to provide a practical solution for solutes with a low chemical affinity for the continuous phase. The core-shell (or coated sphere) model of Kalnin and Kotomin [35], discussed further in Kalnin et al. [36], provided the first extension of the Maxwell-Garnett equation when partitioning occurs. The free geometry parameter, Φ, introduced in Eq. (29) in Kalnin et al. [36] has however no physical meaning and needs to be adjusted according to the arrangement and morphology of the dispersion.

LD simulations solve all previous issues by enabling an averaging over various arrangements and packaging of globules. The principle of simulations and the theoretical treatment of mass transfer across interfaces are developed in section two. In particular, this work generalizes the kinetic approximation of Lejay and Maire [37] at interfaces, by introducing a proper treatment of sorption isotherms. The methodology to determine the effective diffusivities from long-term dynamics are detailed in section three. The following section recapitulates the scaled results obtained for a wide range of conditions, including the effect of diffusion coefficients ratio, partition coefficients, volume fraction of the dispersed phase, and emulsion morphology. The reliability of our results is finally compared to the analytical model of Kalnin and Kotomin [35].

## 2 Langevin dynamics (LD) formulation in structured foods

Duderstadt and Martin [38], transport theory in condensed matter differs from the approaches encountered in classical physics because it involves particles, not a continuum theory of matter, such as fluid dynamics. In particle transport, the random nature of particle interaction events requires replacing densities and mass fluxes by a field of probability densities or distribution functions. The exact number of particles in a specific region of space and time cannot be predicted. Concentrations and particle densities used at macroscopic scale represent, indeed, only the average number of particles per elemental volume unit within a certain period of time. Though it is broadly used at macroscopic scale, the first moments (averages of particle densities or fluxes) are not sufficient to describe diluted systems at microscopic scale, and two classes of problem, so-called direct or indirect problems, address

the issue. The direct problem consists in calculating the distribution of particles in the medium given the composition and geometry of the host medium and the location of eventual sources of particles in the medium. The inverse problem is more difficult and aims at evaluating the effective characteristics of the medium through which the particles have propagated. The underlying random walk process has been described either as a succession of random steps on discrete positions or times or as a continuous process in time or space. Continuous representations are preferred to describe anomalous diffusion. In this work, solute transport is described as a special case of skewed Brownian motion random walks. Following the same ideas as Walsh [39] and Yor [40], the diffusion process is generalized by replacing random variates of cartesian coordinates by polar ones, where each particle moves along rays emanating from its previous position. This scheme is well adapted to describe random walks with discontinuous diffusion coefficients [41]. In short, particle positions are continuous, but they evolve over discrete times.

### 2.1 Equations of motion

Solutes dispersed at low concentration are described as self-moving particles due to thermal fluctuations and without interacting with themselves. This transport mechanism homogenizing solute density in any homogeneous and condensed phase is known as tracer diffusion or molecular diffusion. At macroscopic scale, the distribution of the diffusing species obeys to the general Fokker-Planck equation [42]. At microscopic scale, the motions of small solutes, e.g. not massive particles, are essentially controlled by viscous forces, so that their positions result from the random collisions with surrounding molecules. The energy of the motion is randomly pumped in or pumped out the system via a stochastic formulation of equations of motions. By assuming the conservation of momentum and angular momentum, the n-dimensional stochastic trajectory is captured by a special case of Langevin equation:

$$ \frac{d\mathbf{r}(t)}{dt} = -\mathbf{A}\mathbf{r}(t) + \mathbf{B}\mathbf{\Lambda}^{1/2}\xi(t) \quad (1) $$

where $\mathbf{r}(t)$ is the position vector in an Eulerian reference frame. $\mathbf{A}$ is a $n \times n$ matrix associated with low-frequency deterministic motions, whose eigenvalues $(\alpha_i)_{i=1,n}$ represent the corresponding reciprocal correlation times. $\xi(t)$ is a random vector, whose components are normal variates, with zero means and unitary standard deviations. The last term on the right-hand side is related to the symmetric positive-definite diffusion tensor, $\mathbf{D}$, via its unique Cholesky decomposition:
$$2\mathbf{D} = \mathbf{B}\mathbf{\Lambda}^{1/2}(\mathbf{B}\mathbf{\Lambda}^{1/2})^{T} \quad (2)$$

where $\mathbf{A}$ is a diagonal matrix of eigenvalues and where the columns of $\mathbf{B}$ are the corresponding unit eigenvectors.

In practice, Eq. (1) is straightforwardly integrated with a fixed time step, $\delta t$, via an Itô–Taylor expansion [43]:

$$\mathbf{r}(t + \delta t) = (\mathbf{I} - \mathbf{A}\delta t)\mathbf{r}(t) + \mathbf{B}\mathbf{\Lambda}^{1/2}\sqrt{\delta t}\xi(t) \quad (3)$$

where $\mathbf{I}$ is the identity matrix.

When $\delta t$ is chosen much smaller than $1/\alpha_1$, where $\alpha_1$ is the largest eigenvalue of $\mathbf{A}$, the visited distance by molecular diffusion is much larger than the distance visited by any other deterministic cause. $\mathbf{B}\mathbf{\Lambda}^{1/2}\sqrt{\delta t}\xi(t)$ is the random increment of the Wiener process associated to trace diffusion. Since all increments increase mainly as $\sqrt{\delta t}$ the integration scheme (3) has a low global order of $1/2$. The convergence is slow and long-term iterations are required to sample the trajectory in space and frequency. Since only the second moments of the stochastic trajectory are desired, the choice of the initial starting positions do not modify the final result. Hence a parallel integration of Eq. (1) over $m$ independent trajectories including each $n$ steps offer considerable acceleration by several magnitude orders.

In this work, Eq. (1) applies to both the continuous phase, denoted “c”, and the dispersed phase in it, denoted “d”. The corresponding diffusion tensors are noted $\mathbf{D}^c$ and $\mathbf{D}^d$, respectively. For simplicity, diffusivities were considered isotropic and diagonal matrices replaced diffusion tensors with uniform elements, $\mathbf{D}^c\mathbf{I}$ and $\mathbf{D}^d\mathbf{I}$ with $\mathbf{D}^c$ and $\mathbf{D}^d$ representing isotropic scalar diffusivities with SI units in $m^2 s^{-1}$.

## 2.2 An intuitive description of the free energy landscape on diffusivities

Diffusion on a free-energy landscape is not intuitive but plays a significant role in dense and multiphasic media such as food. It introduces two important consequences on trajectories: it can make some positions more probable and it modifies waiting times before a translation can be observed (i.e., a change in macroscopic state). For instance, heterogeneous dynamics are responsible for the broad mass dependence of solute diffusivities observed in solid polymers [3]. A simplified interpretation of the role of free energy landscape on diffusivities is presented in this section. The random walk along a periodic two-well potential involving two macrostates (or positions) denoted “c” (continuous phase) and “d” (dispersed phase) separated by a distance $l/2$ is presented in Fig. 1, and corresponds to the stationary sequence (trajectory) $\cdots \to c \to d \to c \to d \cdots$. In this 1D case and the transport across interfaces, it is essential to note that the waiting times are random not the positions. The thermodynamical equilibrium between both states is determined by the partition coefficient between the two macrostates, $K$, defined as:

$$K = \frac{p_{eq}^d}{p_{eq}^c} = \exp\left(\frac{\mu_{ex}^c - \mu_{ex}^d}{k_B T}\right) \quad (4)$$

where $\{p_{eq}^k\}_{k=c,d}$ and $\{\mu_{ex}^c\}_{k=c,d}$ are respectively the solute probability density at equilibrium and its excess chemical potentials in phase $k$. A more general definition of the local thermodynamical equilibrium at the interface between the macrostates/positions $c$ and $d$ can be envisioned within the framework of transition state theory (TST) modified by Kramers [44] and Mel'Nikov [45] and experimentally evidenced for transport of isolated molecules by Vestergaard et al. [46]. According to Kramers and

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Random walk along a periodic double-well potential in one dimension corresponding to Eqs. (6) and (7)</figcaption></figure>
in the limits of weak friction (i.e. no inertia), the escape rate of a Brownian particle out of a deep potential well is subjected to many recrossings. Additional discussions on symmetrizable (reversible) jump processes can be found in Chen [47]. At the exact position of the interface $\partial\Omega$, the probabilities to move forward or backward are the same. In other words, the flux of probabilities of transitions from c to d and in the opposite direction should be equal:

$$p_{eq}^{c} k_{c \to d} = p_{eq}^{d} k_{d \to c} \qquad (5)$$

The trajectories of solutes at the interface are not ballistic and are induced by the thermal fluctuations of the total free energy of the entire system at constant volume, F. A translation (whatever the direction: from A to B or from B to A) occurs when the system acquires an energy larger than the critical value $F^\ddagger$ separating the two wells. The height of this energy barrier is greater than the free energies associated with the sorption states c and d, denoted $\{F^k\}_{k=A,B}$. By considering that the thermal fluctuations are Poissonian events (independent and stationary), the escape rates from A to B are hence exponentially distributed with an expectation, in the canonical ensemble, given by:

$$k_{A \to B} = \frac{k_B T}{h} \exp \left( \frac{F^A - F^\ddagger}{k_B T} \right) = k_0 \exp \left( \frac{F^A - F^\ddagger}{k_B T} \right) \qquad (6)$$

where h is the Plank constant and $k_B$ the Boltzmann constant.

By noticing that the Wiener increments are identical for the sequences $d \to c \to d$ and $c \to d \to c$, the effective diffusivity on the double-well surface energy reads:

$$2D = l^2 k_{d \to c \to d} = \frac{l^2}{(k_{d \to c})^{-1} + (k_{c \to d})^{-1}} = \frac{k_0 l^2}{\frac{K}{\exp(-\frac{F_\Delta}{k_B T})} + \frac{1}{\exp(-\frac{F_\Delta}{k_B T})}} = \frac{k_0 l^2}{1 + K} \exp(-\frac{F_\Delta}{k_B T}) = \frac{D_0}{1 + K} \qquad (7)$$

$$\begin{gathered} w i t h \quad F_{\Delta} = F^{\ddagger} - \max(F^{c}, F^{d}) \geq 0 \quad a n d \\ K = \max(p_{eq}^{c}, p_{eq}^{d}) / \min(p_{eq}^{c}, p_{eq}^{d}). \end{gathered}$$

Equation (7) is illustrative of a pure thermodynamic control on diffusivities, which should not be neglected in simulation and require proper implementation. Effective diffusivities, D, deviate from the homogeneous medium as soon as large K values are met. More discussions on the effect of surface energies can be found in Vitrac and Hayert [29]. In food systems, partitioning between phases is very likely due to specific interactions with the host medium at molecular scale [10, 11] or due to the difference in size between the host molecules in phases c and d [8, 12].

## 2.3 Interface conditions for discrete times: partition function for a single hopping event

Equation (7) demonstrates that that changes in chemical affinity near interfaces increase correlation times (see the interpretation of the thermodynamic factor in Eq. 10 of Fang and Vitrac [4]). Integrating Eq. (1) over discrete times while enforcing local thermal equilibrium via Eq. (5) and variable residence times is particularly challenging. Classical Riemann–Stieltjes integral of microscopic displacements along time cannot be used across interfaces, because almost all paths of Brownian motion fail to have bounded variation. Itô calculus already applied to derive Eq. (3), provides an alternative by defining the trajectory as a limit in probability of Riemann sums, but the probability density function must be twice continuously differentiable in space. This requirement is again not verified close to any interface $\partial\Omega$, where abrupt changes in diffusion coefficients or chemical affinity occur. To reach a uniform probability at equilibrium, particle tracking methods used in porous or fractured media (i.e. $D_c/D_d > 1$ and $K=1$) rely on semi-reflecting barriers as reviewed by Delay et al. [48] and Salamon et al. [49]. Reflection is considered, however, only on the side of the high dispersive medium [50] with a reflection rate proportional to $(\sqrt{D^c} - \sqrt{D^d}) / (\sqrt{D^c} + \sqrt{D^d})$ [48]. One corollarily is that a different number of particles cross the interface in a given amount of time [51] from c to d, and from d to c in violation of the principles of local thermodynamical equilibrium. The verification at any time and in a probabilistic sense of mass conservation and of a local partitioning of particles, as envisioned by Kramers, is challenging and not subjected to abundant literature.

First-passage properties offer an alternative to deriving a thermodynamically consistent numerical scheme at the microscopic scale (i.e. at both microscopic space and time scales). First-passage and first-escape times are defined as the time to hit the interface and to leave the current phase, respectively. By assuming that the local direction of translation is not affected by the proximity with interfaces (i.e. no long-range interactions such as Coulombic ones), the probability to hit the interface within a given time is given by the distribution of the particle along the line connecting the current position (coordinate along the line $x = x_0$, at $t = t_0$) and the interface $\partial\Omega$ ($x = 0$). By indexing 1 the original half-plane $x \ge 0$ and 2, its complementary part on the other side of the interface, the probability density
function to find a particle at position x along the line after the dimensionless time $Fo = \frac{D_1(t-t_0)}{x_0^2}$ is given for the two half-segments on both sides of $\partial\Omega$ as:

$$ 
\begin{aligned}
p_{x \ge 0}^{(\text{Fo},x)} &= A \cdot G(\text{Fo}, x_0) + B \cdot G(\text{Fo}, -x_0) \quad \text{with } x \ge 0 \\
p_{x \le 0}^{(\text{Fo},x)} &= C \cdot G\left(\frac{D_2}{D_1}\text{Fo}, \sqrt{\frac{D_1}{D_2}}x_0\right) \quad \text{with } x < 0
\end{aligned} 
\quad (8) $$

with $G(\text{Fo}, x) \frac{1}{2x_0\sqrt{\pi Fo}} \exp\left[-\frac{(1-x/x_0)^2}{4Fo}\right]$ the Gaussian kernel.

In domain 1, the particle may either diffuse or being reflected. Domain 2 contains the fraction of transmitted particle (i.e. crossing $\partial\Omega$). Constants A, B and C are inferred from the condition of local equilibrium $\frac{p_{x \le 0}^{(\text{Fo},x=0)}}{\partial x} = K$, the conservation of flux $D_1 \frac{\partial p_{x \ge 0}^{(\text{Fo},x=0)}}{\partial x} = D_2 \frac{\partial p_{x \le 0}^{(\text{Fo},x=0)}}{\partial x}$ and the mass balance conservation law $\int_0^\infty p_{x \le 0}^{(\text{Fo},x_0)} dx + \int_0^\infty p_{x \ge 0}^{(\text{Fo},x_0)} dx = \text{pr}(1 \to 1 | \partial\Omega, \text{Fo}, x = x_0) + \text{pr}(1 \to 2 | \partial\Omega, \text{Fo}, x = x_0) = 1$. By denoting $u = \left(\frac{x}{x_0} + 1\right)^2 / 4\text{Fo}$, $v = \left(\frac{x}{x_0} - 1\right)^2 / 4\text{Fo}$ and $w = \frac{1}{2} \left(\frac{D_1}{D_2}\right)^2 \frac{1}{\text{Fo}}$, the distribution of probability reads piece-wisely across the interface:

On finite time scales, remaining solutes (i.e. in a probabilistic sense as we discuss only one event) include solutes which did not have “enough” time to cross and reflected ones. Both populations are indiscernible. On very long time scales, all remaining solutes have been already transmitted at least once and were back. The previous result of TST (Eq. 5) is now extended to semi-infinite geometries (i.e. the system is aperiodic). The new partition function between the remaining (i.e. reflected) and crossing probabilities depends not only on K but also on the square root ratio of diffusion coefficients (i.e. ratio of transmission rates) between the two half-spaces separated by the boundary $\partial\Omega$. The thermodynamical equilibrium is exactly and only verified at the interface $\partial\Omega$ (i.e. the ratio of probabilities is set by K), but it cannot be reached macroscopically since both media are semi-infinite. When particles crossed media with different properties (diffusivity or chemical affinity), they are not only dispersed, but they are also refracted and even diffracted when they cross multiple interfaces.

## 2.4 Choice of time steps and generalization to 2D and 3D structures

As discussed by Balescu [52, see Chapter 9] the sampling of random walks in inhomogenous media is significantly

$$ 
\begin{aligned}
p_{x \ge 0}^{(\text{Fo},x_0)} &= \frac{e^{w/2-v} \left(1 + K \sqrt{\frac{D_1}{D_2}}\right) + e^{w/2-u} \left(1 - K \sqrt{\frac{D_1}{D_2}}\right)}{2Z} && \text{for } x \ge 0 \\
p_{x \le 0}^{(\text{Fo},x_0)} &= \frac{K e^{\frac{2D_1^2-D_2^2}{4D_2^2 Fo}} e^{-\frac{D_1 \left(x-x_0 \sqrt{\frac{D_1}{D_2}}\right)^2}{4D_2 Fo x_0^2}}}{Z} && \text{for } x \le 0
\end{aligned} 
\quad (9) $$

with $Z = \sqrt{\pi} \left( x_0 e^{w/2} \sqrt{Fo} + K \sqrt{\frac{D_2}{D_1}} Fo x_0 \left( e^{\frac{2D_1^2-D_2^2}{4D_2^2 Fo}} \left[ 1 - \text{erf}\left(\frac{1}{2} \frac{D_1}{D_2} \frac{1}{\sqrt{Fo}}\right) \right] + e^{w/2} \text{erf}\left(\frac{1}{2\sqrt{Fo}}\right) \right) \right)$

Equation (9) enables to calculate the partition function associated with the considered hopping event. By interpreting it as the ratio between the fraction of remaining and transmitted solutes for a particle initially located at $x_0$ from $\partial\Omega$ after a dimensionless duration $Fo$, one gets:

$$ 
\begin{aligned}
& \frac{\text{pr}(1 \to 1 | \partial\Omega, \text{Fo}, x = x_0)}{\text{pr}(1 \to 2 | \partial\Omega, \text{Fo}, x = x_0)} \\
& = \frac{1}{K} \sqrt{\frac{D_1}{D_2}} \frac{1 + K \sqrt{\frac{D_1}{D_2}} \text{erf}\left(\frac{1}{2\sqrt{Fo}}\right)}{1 - \text{erf}\left(\frac{D_1}{D_2} \frac{1}{2\sqrt{Fo}}\right)} \exp\left[-\left(\frac{D_1^2}{D_2^2} - 1\right) \frac{1}{Fo}\right] \\
& \to Fo \to \infty \frac{1}{K} \sqrt{\frac{D_1}{D_2}}
\end{aligned} 
\quad (10) $$

improved (i.e., better ergodicity) when ordinary dynamics such as Eq. (1) are transformed into statistical dynamics. In this perspective, the concept of “sampled time” should be envisioned as a “time-ordered path” to explore the configuration space accessible to particles rather than the true physical time (see Chapter 4 in Landau and Binder [53]). A not necessary but sufficient condition to ensure a Gibbsian equilibrium state is to verify at any time the detailed balance condition, also known as microreversibility, across $\partial\Omega$ [54]. The term of Gibbsian equilibrium is used because statistical mechanics is associated since its foundation with two partly incompatible definitions of equilibrium. Chemical and food engineers are familiar with the Gibbsian equilibrium. The equilibrium is defined from the stationary distribution of
ensembles maximizing entropy. The Boltzmannian equi-
librium imposes only the macro-variables to be almost
constant on the equilibrium region, such as considered
previously for semi-infinite media (see for example the
discussion in Werndl and Frigg [55]). In shorts, Eq. (10)
applied to each "jumping event" along a line in n = 1, 2, 3
dimensions enforces Gibbsian equilibrium irrespectively
the implementation of the interface condition. Since the
proposed algorithm of first passage near interfaces is
analogous to the Fermat's principle in optics, any step of
duration δt intercepting the interface was split into two
times δt₁ + δt₂ = δt. δt₁ was the ballistic time required
to hit the interface and δt₂ the complimentary time to
complete partial specular reflection or partial refraction.
Because only one possibility was considered during δt₂
(particles were not duplicated), reflection was chosen
when a random number picked uniformly between 0
and 1 was lower than pr(1 → 1|∂Ω, Fo, x = x₀) and refrac-
tion otherwise. Because refractive indices were associ-
ated with the phase velocity of the diffusion front, they
were chosen proportionally to √D₁/D₂. Implementating
reflection in a billard ball fashion (i.e. with equal indicent
and reflected angles across the local normal) and reflec-
tion according to the Snell-Descartes law guaranteed
that the normal probability flux was conserved across
the interface for any initial incidence of the particle in
2D and 3D.

On time scales commensurable to δt, any trajectory
looks ballistic (jump process) with a random direction.
Since an analytical equation was used, increasing δt did
not provoke any loss of information, while interfaces could
be assumed plane and separating two semi-infinite half-
spaces. The distribution of probability in a convex region
is obtained by sampling several recrossings from different
directions. For convex objects such as globules, the inte-
gration time step must be chosen small enough so that
the probability of leaving the smallest object (e.g., globule)
from its center and without making a turn is almost zero.
The escape probability for small objects increases indeed
dramatically with highly curved interfaces. Figure 2 com-
pares exact analytical solutions and the corresponding
escape rates for a single macroscopic jump between two
infinite and semi-infinite spaces [Eq. (9) in 1D] and for
spherical globules (numerical simulation). The comparable
condition corresponds to D_c/D_d = D_g/D_A = 4 and K = 2 for a sol-
ute initially located at a distance x₀ of the nearest interface.
Size effects are shown for radii varying from x₀ to 20x₀.

Both descriptions succeed in reproducing the diffrac-
tion ring (two local maxima): one corresponding to the
initial position and an other emanating from the outside of
the interface where the chemical affinity and diffusion are
higher. The escape rates of the half-space (one exit point

corresponding to one single ballistic event in the simula-
tion) and of the entire object (many exit points) do not
coincide and should not be confused together (Fig. 2b).
Only the domain associated to low Fo values is sought in
simulations. Additionally, the calculations illustrate that
diffusion has broken apparent symmetry around fractal
and convex objects with poor penetrability. Mass diffusion
evolves rapidly to surface diffusion around impermeable
crowded or percolating objects. The complications associ-
ated with the arrangement of spheres was discussed early
in Chapter 12 of Crank's reference textbook [31].

**3 Simulated conditions**

**3.1 Food description**

The presented case-study addresses generically the diffu-
sion, sorption/desorption, and permeation in and across
a food including a dispersed phase, as shown in Fig. 3. The
continuous phase is assumed to be dense, viscous and
associated with homogeneous diffusion properties. The
discontinuous phase consists of globules with fixed posi-
tions. The eventual surfactant layer is neglected. Solutes
(aroma, contaminants, etc.) can diffuse freely within both
phases. The effective diffusivities are averaged on several
2D and 3D periodic systems sampling the real size distri-
bution, arrangement, and packing of globules. A dimen-
sionless formulation applicable to many emulsionned food
products (dairy products, seasonings, etc.) is proposed. In
its present form, application to foams is discarded, and all
phases need to be dense. The fundamental reason is that
the particle velocity should be comparable between all
phases: few µm/min or less in liquids and solids against
several km/min in gases.

**3.2 Scaled emulsions**

2D and 3D emulsions were represented as dispersed glob-
ules mapped onto a scaled periodic cell over [0, 1] × [0, 1[
and [0, 1] × [0, 1] × [0, 1[ respectively. Periodic condi-
tions enabled to reproduce an infinite medium with a
periodicity l₀. The size of each globule was sampled from
a log-normal distribution with arithmetic average, μØ,
and standard deviation, σØ. Globules were inserted ran-
domly according to a generic algorithm independent of
the shape of the object to be inserted. The positions of
all phases remained fixed during the simulation. To facili-
tate the sampling of the free space, a M × M or M × M × M
fine grid (with M > 5 · 10³) was used to mark sub-cell as
occupied or free space. For each new globule of diam-
eter Ø to be inserted, a non-centered coarse grid with a
step-size equal to 2 · Ø was superposed to the fine one.
<figure><img src="image_6.png" /><figcaption>**Fig. 2** Comparison of escape rates from **a** interfaces, **b**, **c** from globules for a solute a solute rate initially located at $x_0$ for $D_B/D_A = 4$ and $K = 2$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3** Principle of the idealization of the **a** real food structure (e.g., dairy product) as **b** equivalent emulsions, including one continuous phase (liquid or gel, not shown) and a globular dispersed phase (corresponding volume fractions: $d = 0.4, 0.3, 0.2$)</figcaption></figure>
All sub-cells of the coarse grid that contained enough free space were randomly numbered and iteratively tested. Within each sub-cell, all possible centers laying within the relative positions [-Ø/2, Ø/2] x [-Ø/2, Ø/2] or [-Ø/2, Ø/2] x [-Ø/2, Ø/2] x [-Ø/2, Ø/2] (the sub-cell center was assumed to be located at 0,0 or 0,0,0) were checked for insertion so that the domain enclosed by the globule did not overlap any occupied site on the fine grid. When no feasible position was encountered on the coarse grid, the coarse grid was randomly translated up to 200 times. When no insertion was possible, a new globule was randomly picked. The whole process stopped when the set density in globules, defined by the fraction of total occupied sites, was achieved within a tolerance about 1%. All the structures that did not meet the target were discarded.

The tested conditions are summarized in Table 1. Two extreme morphologies were studied: a monodisperse emulsion, denoted E0, and a polydisperse one, denoted E1. The effect of tortuosity was tested by varying the volume fraction of dispersed phase, d, between 0.1 and 0.4. To reduce possible bias due to the arrangement of globules, all trajectories were averaged over 10–20 different structures. Examples of emulsions and their corresponding size distributions are given in Fig. 4. According to the set density in globules, each morphology included between 70 and 800 globules in 2D. Calculations in 3D comprised more than 10⁴ globules.

### 3.3 Tested range of physico-chemical properties

The diffusion coefficient in the dispersed phase, $D^d$, was chosen as reference as it enables to emphasize the effect of entropic trapping. Its value was set lower than $I_0^2(4M^2\delta t)^{-1}$ so that more than 95% of elemental trajectories during $\delta t$ had a length lower than the step of the fine grid. This condition ensured that even small globules were correctly sampled. The final granularity was far beyond the possibilities order spectral methods, such as finite elements. Simulations were performed for ratios $r_D = D^c/D^d$ ranging between $10^{-2}$ and $10^2$. The range covered most of the food applications (oil in water, water in oil, yogurt, sauce, etc.). Partition coefficients K, defined as the continuous-to-dispersed phase ratio of concentration

at thermodynamical equilibrium were sampled over the whole scale: from 0 (no affinity for the continuous phase) to +∞ (no affinity for the dispersed phase). The interface conditions previously described was implemented using a search procedure, which did not presuppose any shape for the dispersed phase (e.g., disc or sphere). Interfaces and their surface normal were detected with the following tier approach: tier (1) coarse search on the fine grid, tier (2) continuous search by finding the intersection of trajectory with the parametric equation of the object contour, tier (3) recursive search of reflections using steps 1 and 2. An example of trajectory (one particle) is detailed in Fig. 5. The whole scheme can manage up to ten successive reflections and diffusions during $\delta t$ if needed. Simulations are performed at constant temperature.

### 3.4 Identification of diffusion coefficients

The effective diffusion tensor D was estimated from the asymptotic rate of increase of the covariance matrix of displacements with diffusion time:

$$D = \frac{1}{2} \left( \lim_{\tau \to \infty} \frac{\partial}{\partial \tau} \left\langle \left( \mathbf{r}_{(t+\tau)} - \mathbf{r}_{(t)} \right) \left( \mathbf{r}_{(t+\tau)} - \mathbf{r}_{(t)} \right)^{\mathrm{T}} \right\rangle_{t,i=1..m} \right) \quad (11)$$

where $\langle \rangle_{t,i=1..m}$ is the ensemble average operator over all simulated times, trajectories and sampled emulsions. It is worth noticing that $t$ and $\tau$ are not equivalent. The first is indeed the simulation time (i.e., the physical time counted from the origin), whereas the latter is so-called diffusion time and represents a lag. Since only stationary trajectories vs $\tau$ are studied (i.e., disregarding the origin of times), it is a common practice to denote lags $t$ in figures when it represents $\tau$.

The isotropic diffusivity is finally given by the expectation of the quadratic form:

$$ \begin{aligned} D &= \frac{\mathrm{tr}(\mathbf{D})}{n} = \frac{1}{2n} \left( \lim_{\tau \to \infty} \frac{\partial}{\partial \tau} \left\langle \left( \mathbf{r}_{(t+\tau)} - \mathbf{r}_{(t)} \right)^{\mathrm{T}} \left( \mathbf{r}_{(t+\tau)} - \mathbf{r}_{(t)} \right) \right\rangle_{t,i=1..m} \right) \\ &= \frac{1}{2n} \lim_{\tau \to \infty} \frac{\partial}{\partial \tau} \mathrm{msd}(\tau) \end{aligned} \quad (12) $$

<table><caption>**Table 1** Tested morphologies of emulsions</caption><thead><tr><th rowspan="2">Morphology</th><th colspan="2">Globule diameter (×l<sub>0</sub>)</th><th colspan="4">Number of repetitions</th></tr><tr><th>μ<sub>Ø</sub></th><th>σ<sub>Ø</sub></th><th>d=0.1</th><th>d=0.2</th><th>d=0.3</th><th>d=0.4</th></tr></thead><tbody><tr><td>E<sub>0</sub> (monodisperse)</td><td>4 × 10<sup>-2</sup></td><td>3 × 10<sup>-8</sup></td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><td>E<sub>1</sub> (polydisperse)</td><td>2 × 10<sup>-2</sup> ± 4 × 10<sup>-3</sup></td><td>2.8 × 10<sup>-2</sup> ± 6 × 10<sup>-3</sup></td><td>20</td><td>20</td><td>20</td><td>20</td></tr></tbody></table>

Mean and standard deviations are based on sampled values in presented simulations
<figure><img src="image_8.png" /><figcaption>Fig. 4 Example of 2D emulsions on a periodic[0, 1] × [0, 1] domain with different volume fractions in dispersed phase, a, e d = 0.1, b, f d = 0.2, c, g d = 0.3, d, h d = 0.4; a–d monodisperse emulsions (E₀)</figcaption></figure>

and polydisperse emulsions (E₁). Histograms plot the distributions of the surface area of globules fitted with a continuous log-normal distribution

where tr() is the trace operator and n is the number of dimensions. msd(τ) represents the average mean-square displacement (variance assuming zero mean) on time scale τ.

The asymptotic regime is reached when the displacements are time-independent and with a white noise spectrum. In the absence of correlations, msd(τ) increases as a power law τ$^{α}$, with an exponent α close to unity. Confinement inside or between packed globules increases considerably the time required to reach independent displacements. For the particular value K = 0, confinement occurs over all time scales, and the asymptotic diffusion coefficient is zero since $\lim_{\tau \to \infty} msd(\tau) = dl_0^2$.

Simulation lengths were adapted to get msd values s times larger than $dl_0^2$, with s ≥ 2. According to Eqs. (11) and (13), the necessary condition reads:

$$m \cdot t > t_{\min} = s \frac{dl_0^2}{2nD^c} = sdM^2\delta t \quad (13)$$

s values up to 20 and m values up to 500 were applied in this study. All simulations were performed for products mt longer than $10^7 \delta t$. When m is large, starting from a distribution of particles reproducing an equilibrium distribution accelerates convergence.
<figure><img src="image_9.png" /><figcaption>Fig. 5 Details of a trajectory close to an interface: **a** reflection ($d \to d|_{\partial\Omega}$), **b** transmission ($d \to c|_{\partial\Omega}$). The depicted condition corresponds to $K = 10$, $r_D = 2$, $d = 0.3$. The positions of both inserts are identified by double arrows and coordinates</figcaption></figure>

# 4 Results and discussion

## 4.1 Typical 2D trajectories and mean-square displacements

### 4.1.1 Trajectories

The effect of $r_D$ and $K$ on the trajectory of $m = 5$ diffusants is illustrated in Fig. 6 for $r_D = 2$ and $K = \{0.1, 1, 10\}$ for a polydisperse emulsion with $d = 0.3$. The corresponding cumulated diffusion time was set to $m \cdot t = 5 \cdot 10^6 \cdot \delta t$ to get a partly ergodic system but without reaching equilibrium (i.e., a single particle could not visit all globules; as a result, the trajectories can be distinguished). Corresponding residence times averaged over all diffusants were calculated on a $50 \times 50$ grid and plotted as isocontours.

Raising the affinity for the dispersed phase ($K \to 0$) confined the trajectories inside the globules (Fig. 6a) and increased the residence time in the globules (Fig. 6b). For $K = 1$, the trajectories spread over all the simulated without concentrating in any particular phase (Fig. 6c). For all conditions and irrespectively of the $K$ value, the residence times were longer inside the globules than in the continuous phase (Fig. 6d), in agreement with the twice lower diffusion coefficient in the disperse phase. In the very long term, residence times were not uniform between globules. As exemplified in Fig. 2b, c, dwelling times in large objects (slow escape rates) than in small ones (high escape rates).

The condition $K = 10$ filled more uniformly the space between globules (Fig. 6e) while visiting less frequently the dispersed phase (Fig. 6f). In regions where globules were concentrated, residence time increased due to more tortuous paths.

From previous examples, trajectories exhibited strong correlations either in globules or between globules. Locally low diffusivities, tightly packed globules, and high chemical affinity for the dispersed phase caused confinement alone or in combination. Persistent entropic trapping associated with negative correlations is known to lead to significant subdiffusive regimes on intermediate time scales (i.e., before displacements become decorrelated).

### 4.1.2 Mean-square displacements

Dimensionless mean-square displacements (msd) versus dimensionless diffusion time ($tD^dI_0^{-2}$) are plotted for a wide range of combinations $r_D \times K \times d \times E$ in Fig. 7. Msd values of each curve were averaged over the trajectories of 100 diffusants. Since longterm diffusion can propagate only through the continuous phase, the asymptotic diffusion of the continuous phase, $r_D$, is shown as reference. On a log-log scale, the asymptotic regime is identified by a unitary slope ($\alpha \approx 1$) and the effective diffusivity by its intercept. Negative correlations (subdiffusive regime associated with $\alpha(t) \le 1$) slowed down the increase of msd and translated the curve to the right.

In decreasing order of significance, $r_D$ and $K$ had the main effects. Increasing $r_D$ translated msd curves in the same fashion (Fig. 7d-f), whereas the effect of $K$ (for $K \ne 1$) was not monotonic and modified the local slope $\alpha(t)$ on intermediate time scales (Fig. 7a-c). Large $K$ values increased diffusion rates on all time scales. By contrast, the effect of small $K$ values was not symmetric and affected more significantly long-term diffusivities with drops down to two magnitude orders in tested conditions.

In details, $K$ and $r_D$ could have either cooperative or antagonist effects according to the considered condition. This interpretation was already suggested by Eq. (10). When the diffusion is slow in the continuous phase, the diffusion through the dispersed phase compensates the slowdown partly. This effect is very likely when the concentration in globules is high and when the values of $K$ are low but not too small. By contrast, large $K$ values tend to reroute all particles between and around globules as they are mainly seen as obstacles. Close to the percolation threshold, the dispersed phase becomes almost a co-continuous phase, and the roles of each phase can be reversed. When diffusion flux occurs mainly across the dispersed phase, the cross-section of contact and passage between globules limits the effective diffusion coefficient.
<figure><img src="image_10.png" /><figcaption>**Fig. 6** **a, c, e** Typical trajectories during $10^6 \cdot \delta t$, **b, d, f** corresponding residence times in the emulsion plotted in Fig. 4g ($d = 0.3$) for $r_D = 2$ and **a, b** $K = 0.1$, **c, d** $K = 1$, **e, f** $K = 10$. Five independent trajectories are plotted, residence times are averaged over the five trajectories trajectories and calculated on a $50 \times 50$ grid</figcaption></figure>

## 4.2 Individual effects of main parameters

### 4.2.1 Effect of $K$

The complex effect of $K$ on the asymptotic diffusion coefficient is analyzed in Fig. 8. A search algorithm was applied to identify the hydrodynamic regime while discarding the very terminal region with poor sampling. Confidence intervals were calculated from the repetitions for different emulsions (Table 1). Uncertainty was ranged between a few percents and a factor 2 according to tested conditions. The poorest estimates were obtained when low $r_D$ and $K$ values were combined together. Their combination increased the duration of the sub-diffusive regime drastically and increased the total simulation cost dramatically:

both simulation length and the number of particles needed to be augmented. In practice, the effect of very low $r_D$ and $K$ values could not be explored with a reasonable cost beyond a practical lower limit of $10^{-2}$. However our methodology and code remain valid beyond this limit. It was verified that the asymptotic diffusion coefficient was 0 when $K=0$ whatever the applied $r_D$ value. By contrast, there was no computational limitation for applying large $K$ values. Infinite $K$ values were simulated by assuming a value of $10^{99}$. This condition made the globules impenetrable and asymptotic diffusion coefficients $D/D^d$ converged to $r_D$.

For the same $r_D$ value, a variation of $K$ between $10^{-2}$ and $+\infty$ generated a variation of $D/D^d$ up to five folds. A greater effect was even observed for values lower than
$10^{-2}$. The effect of $K$ varied significantly according to the value of $r_D$. When the diffusivity was higher in the dis-continuous phase, the asymptotic diffusion coefficient did not vary monotonously with $K$ (Fig. 8d) and reached a maximum close to $K = 1$. The maximum was exactly located at $K = 1$ when $r_D = 1$ (Fig. 8e). The maximum was obtained for $K$ values slightly lower than unity for $r_D < 1$ and above unity for $r_D > 1$ (Fig. 8f-j). $K \to \infty$ did not lead to any maximum.

It was noticeable that the effect of $d$ was more pro-nounced for values of $K$ and $r_D$ above unity. For $r_D$ values higher than 1, the maximum was obtained for the low-est concentration in globules and was always lower than the $r_D$ value. For $r_D$ values lower than 1, the trend was the opposite, the maximum was obtained for the highest $d$ values and reached values higher than $r_D$.

### 4.2.2 Effect of $r_D$

Mean square displacement presented in Fig. 8 suggested that $r_D$ was the main scaling coefficient of $D/D^d$ values. Figure 9 plots on a log-log scale the values of $D/D^d$ versus $r_D$ to identify the value of the critical exponent $\alpha$ which would correspond to the asymptotic model $D/D^d \propto r_D^\alpha$. Besides, the plotted curve $y = x$ makes it possible to identify the particular conditions, which verified $D/D^d > r_D$.

The critical exponent $\alpha$ was lower than 1 for almost all conditions except when $r_D$ and $K$ values were together lower than 2. $\alpha$ values slightly higher than 1 yielded effective diffusion coefficients higher than in the continuous phase only for intermediate $K$ slightly lower than 1 (Fig. 9b-d). As a rule of thumb, $\alpha$ decreased with $r_D$. The decay varied from values close to 1 down to 0.25 in this study. The decaying rate was minimal for large $K$ values. It was almost 0 when $K \to \infty$.

The results simulated with a complete reflection against globules ($K \to +\infty$) are compared with the well-known tortuous model $D = (1-d)D^c$ [56] in Fig. 10 by plotting $r_D^{-1} \cdot D/D^d = D/D^c$ with $r_D$. The plots should theoretically depend on $d$ and not on $r_D$. The tortuous model was globally well verified while exhibiting significant dispersion between repetitions. A lack of convergence was noticed for low $d$ and $r_D$ values, that, when of reflections, was not sufficiently sampled.

### 4.3 Comparison of simulated results with the model proposed by Kalnin and Kotomin [35]

A more general attempt of validation was thought by comparing our simulations with the predictions of Eq. (29) in Kalnin and Kotomin [35]. In this model, diffusion occurs in concentric spheres/discs. The polar/spherical topology does not match the descriptions in Cartesian coordinates, but it

offers an acceptable description when all globules are non-contact. In our notations, the analytical model in 2D reads:

$$ \frac{D}{D^d} = \frac{r_D}{1 - \Phi - \frac{\Phi}{K}} \cdot \left[ 1 + 2 \cdot \left( \frac{1}{K} - r_D \right) \cdot \frac{\Phi}{r_D + \frac{1}{K} - \left( \frac{1}{K} - r_D \right) \cdot \Phi} \right] \quad (14) $$

Equation (14) accepts a free parameter associated with the volume ratio of concentric regions, $0 \le \Phi \le 1$. It was fitted to our simulated results obtained for the same $d$ value and morphology. Simulated and analytical estimates are compared in Fig. 11 for all tested conditions where a good convergence was achieved. The agreement was globally correct, but negative deviations were observed for inhomogeneous emulsions ($E_1$) and positive ones for the smallest $D/D^d$ values. Figure 11b provides an interpretation of the deviations by showing that $\Phi$ cannot be captured with the sole effect of $d$, but that it depends also non-linearly of $K$ or $r_D$. Confidence intervals on $\Phi$ were obtained indeed by adjusting independently $\Phi$ to the results simulated with similar $K$ and $r_D$ values. The comparison of the width of confidence intervals with the effect of $d$ shows that contribution of $d$ is as high as the other variations due to $K$ or $r_D$. As a result, $\Phi$ should be envisioned as a nonlinear function of $d, K, r_D$ and of the morphology. The presented results extend therefore the model and the conclusions of Kalnin and Kotomin [35] on more realistic morphologies.

### 4.4 Diffusivities in 3D emulsions

Simulations on 3D emulsions generalize the results obtained on 2D ones. They are more expensive to calculate and require computers with multiple cores (up to 256 cores were used) to reach convergence in a couple of hours. Asymptotic diffusivities are plotted in Fig. 12 as isocontours. Contour plots exhibited a half square shape confirming the predominance effect of $r_D$ (rate of diffusion in the continuous phase) when $K$ values were greater than unity $10^{-1}$. The effect of $K$ was significant mainly for values lower than 1. Oscillations were observed for some isocontours and were related to the non-linear combination of the effects of $r_D$ and $K$ around unity, as already shown in 2D emulsions (Fig. 8). Concentration in globules and polydispersity modified the location of the maximum effect of $K$ scale at constant $r_D$ values.

### 5 Conclusions

Mass transfer in food is known to interact with the food structure. Effective properties are not only scarce, but they do not explicit the complex relationship with local composition and microstructure. Multiscale modeling of
<figure><img src="image_11.png" /></figure>
►Fig. 7 Typical mean square displacements (msd) versus diffusion time. **a–c** Effect of *K* for morphology *E*₁, *d* = 0.3 and **a** *r*<sub>D</sub> = 0.1, **b** *r*<sub>D</sub> = 1, **c** *r*<sub>D</sub> = 10; **d–f** effect of *r*<sub>D</sub> for morphology *E*₁, *d* = 0.3 and **d** *K* = 0.1, **e** *K* = 1, **f** *K* = 10; **g–i** effect of *d* for morphology *E*₁, *r*<sub>D</sub> = {0.1, 1, 10} and **g** *K* = 0.1, **h** *K* = 1, **i** *K* = 10, **j–l** effect of morphology for *d* = 0.3, *r*<sub>D</sub> = {0.1, 1, 10}, **j** *K* = 0.1, **k** *K* = 1, **l** *K* = 10. Theoretical msd curves corresponding to an asymptotic diffusion coefficient *D*/*D*<sup>d</sup> = *r*<sub>D</sub> are plotted in dashed lines in a–c and g–l

mass transfer is with this respect much more preferable
as it opens the possibility of optimizing properties and
exploring virtually new trends regardless of the existence
of any underlying theory and of the food feasibility. This
work proposes a comprehensive framework to simulate
mass transfer over discrete times using Langevin dynam-
ics (LD) capable of mimicking the random walk of diluted
solutes in polyphasic and multiconstituent food systems.
The parameterization of the simulation is simplified by the
use of intrinsic transport properties respectively to each
phase, such as diffusion coefficients, activity coefficients
or sorption isotherms. These properties, as well as detailed
microstructures, are either broadly available or can be
measured independently. At each time step, the continu-
ous and probabilistic interpretation of the first "interface
hit" enforce both mass conservation and local thermody-
namical equilibrium across the picked jumping direction.
The integration scheme can handle potentially any kind of
interfaces or frontiers, including interfaces between two
fluids, one fluid and a solid, two solids or interactions with
macromolecules. The only condition is that the normal
direction at the "hit" point is known. In addition, the anal-
ogy with the Kramers interpretation of chemical reaction
paths opens the extension to diffusion-reaction problems
and chemisorption.

A detailed case-study is shown covering hundreds of
configurations met in emulsionned food systems when the
globules are diluted or concentrated near their percolation
threshold: K = [10⁻², ∞ [ × r_D = [10⁻², 10²[
×d = [0.1, 0.4] × [E₀, E₁]. The use of a dimensionless formu-
lation enabled to cover almost all practical cases far
beyond the domain of applicability of the Maxwell–Gar-
nett equation and its extensions. Beyond their analogy

with microscopic observations of emulsions, trajectories
in 2D systems offer of visual interpretation of the geomet-
ric construction of skewed trajectories in crowded environ-
ments. Transitions between free and confined diffusion
under thermodynamical control are directly observable.
To accelerate convergence, it is preferable to implement
the methodology on multicore processors with small sys-
tems rather than simulating a single very large system. In
its 2D and 3D versions, the simulated systems might be
seen as oversimplified emulsions without surfactants and
relative displacements of the globules respectively to the
continuous phase. It is, however, straightforward to add a
third phase to represent globules as core–shell objects and
to implement proper interactions with the considered sol-
ute. Mass transport in micelles and liposomes can be rep-
resented in the same fashion. Adding the contribution of
sedimentation/creaming, coalescence and Oswald ripen-
ing requires integrating the Stokesian dynamics of the
globules. By replacing a series of static and randomly gen-
erated emulsion by a series snapshots representing the
evolution of the emulsion over time, fundamental results
proven on bond-percolation models [57–59] can be used
to derive dynamic diffusivities in a straightforward man-
ner. The time–frequency frequency effective diffusion
coefficient, *D*<sub>dynamic</sub> (−i*ω*, *λ*), associated with the renewal
emulsion
rate *λ* is indeed the analytical continuation of the “static”
diffusivity as *D*<sub>dynamic</sub> (−i*ω*, *λ*) = *D*<sub>static</sub> (−i*ω* + *λ*), with
*D*<sub>static</sub> (−i*ω*) = *D*<sub>dynamic</sub> (−i*ω*, *λ* = 0). By covering the
microstructure scale, LD simulations fill the gap between
modeling at the atomistic/molecular scale, where compo-
sitional effects hold, and conventional deterministic or
Bayesian modeling where decisions are taken based on
shelf-life, bioavailability and food safety. Tackling the cou-
pling of these models across scales is the next frontier for
systems (reactions, structures) evolving on different time
scales. We will present the coupling between transport
and oxidation reactions in a companion paper.
<figure><img src="image_12.png" /><figcaption>**Fig. 8** Effects of *K* on the asymptotic diffusion coefficient for typical *r*<sub>D</sub> values: **a** 10<sup>-2</sup>, **b** 10<sup>-1</sup>, **c** 0.25, **d** 0.5, **e** 1, **f** 2, **g** 4, **h** 10, **i** 20, **j** 100. The effects of *d* and morphology are also depicted. 95% confidence intervals based on repetitions are depicted as error bars</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Effect of $r_D$ on asymptotic diffusivities for typical K values: **a** $10^{-2}$, **b** $10^{-1}$, **c** 0.25, **d** 0.5, **e** 1, **f** 2, **g** 4, **h** 10, **j** 100. The effects of *d* and *morpology* are also depicted</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 10** Scaling relationship of the asymptotic diffusivities in 2D emulsions for $K \to +\infty$. The corresponding theoretical “obstacle” model $D/D^c = 1 - d$ is plotted as dotted lines</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 11** a Comparison between simulated results and the fitted algebraic model proposed by Kalnin and Kotomin [35]. The curve $y=x$ is plotted in dotted line. b Values of the fitted parameter $\Phi$ in</figcaption></figure>

the Kalnin and Kotomin model by assuming it varies only with $d$. 80% confidence intervals plotted as error bars were calculated by considering possible effects of $K$ and $r_D$
<figure><img src="image_16.png" /><figcaption>**Fig. 12** Asymptotic diffusivities in 3D emulsions with *K* and *r*<sub>D</sub> for *a*, *e* *d* = 0.1, *b*, *f* *d* = 0.2, *c*, *g* *d* = 0.3, *d*, *h* *d* = 0.4, *a*-*d* monodisperse emulsions *E*<sub>0</sub>, *e*-*h* polydisperse emulsions *E*<sub>1</sub></figcaption></figure>

## Compliance with ethical standards

**Conflict of interest** On behalf of all authors, the corresponding author states that there is no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Datta AK (2016) Toward computer-aided food engineering: mechanistic frameworks for evolution of product, quality and safety during processing. J Food Eng 176:9–27

2. Durand M, Meyer H, Benzerara O, Baschnagel J, Vitrac O (2010) Molecular dynamics simulations of the chain dynamics in mono-disperse oligomer melts and of the oligomer tracer diffusion in an entangled polymer matrix. J Chem Phys 132(19):194902

3. Fang X, Domenek S, Ducruet V, Refregiers M, Vitrac O (2013) Diffusion of aromatic solutes in aliphatic polymers above glass transition temperature. Macromolecules 46(3):874–888

4. Fang X, Vitrac O (2017) Predicting diffusion coefficients of chemicals in and through packaging materials. Crit Rev Food Sci Nutr 57(2):275–312

5. Vitrac O, Lezervant J, Feigenbaum A (2006) Decision trees as applied to the robust estimation of diffusion coefficients in polyolefins. J Appl Polym Sci 101(4):2167–2186

6. Zhu Y, Welle F, Vitrac O (2019) A blob model to parameterize polymer hole free volumes and solute diffusion. Soft Matter 15(43):8912–8932

7. Gillet G, Vitrac O, Desobry S (2009) Prediction of solute partition coefficients between polyolefins and alcohols using a generalized Flory–Huggins approach. Ind Eng Chem Res 48(11):5285–5301

8. Gillet G, Vitrac O, Desobry S (2010) Prediction of partition coefficients of plastic additives between packaging materials and food simulants. Ind Eng Chem Res 49(16):7263–7280

9. Kadam A, Karbowiak T, Voilley A, Bellat J-P, Vitrac O, Debeaufort F (2014) Sorption of *n*-hexane in amorphous polystyrene. J Polym Sci Part B Polym Phys 52(19):1252–1258

10. Nguyen P-M, Guiga W, Dkhissi A, Vitrac O (2017) Off-lattice Flory–Huggins approximations for the tailored calculation of activity coefficients of organic solutes in random and block copolymers. Ind Eng Chem Res 56(3):774–787

11. Nguyen P-M, Guiga W, Vitrac O (2017) Molecular thermodynamics for food science and engineering. Food Res Int 88(Part A):91–104

12. Vitrac O, Gillet G (2010) An off-lattice Flory–Huggins approach of the partitioning of bulky solutes between polymers and interacting liquids. Int J Chem Reactor Eng. https://doi.org/10.2202/1542-6580.2094

13. Nauman EB, He DQ (2001) Nonlinear diffusion and phase separation. Chem Eng Sci 56(6):1999–2018

14. Busch K, Soukoulis CM (1995) Transport properties of random media: a new effective medium theory. Phys Rev Lett 75(19):3442–3445

15. Chen XD (2007) Moisture diffusivity in food and biological materials. Dry Technol 25(7–8):1203–1213

16. Doulia D, Tzia K, Gekas V (2000) A knowledge base for the apparent mass diffusion coefficient (DEFF) of foods. Int J Food Prop 3(1):1–14

17. Mittal GS (1999) Mass diffusivity of food products. Food Rev Int 15(1):19–66

18. Panagiotou NM, Krokida MK, Maroulis ZB, Saravacos GD (2004) Moisture diffusivity: literature data compilation for foodstuffs. Int J Food Prop 7(2):273–299

19. Zogzas NP, Maroulis ZB, Marinos-Kouris D (1996) Moisture diffusivity data compilation in foodstuffs. Dry Technol 14(10):2225–2253

20. Sznitman A (2001) On a class of transient random walks in random environment. Ann Probab 29(2):724–765

21. Gingold RA, Monaghan JJ (1977) Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc 181(3):375–389
22. Lucy LB (1977) A numerical approach to the testing of the fission hypothesis. Astron J 82:1013-1024

23. Monaghan JJ (1988) An introduction to SPH. Comput Phys Commun 48(1):89-96

24. Monaghan JJ (1994) Simulating free surface flows with SPH. J Comput Phys 110(2):399-406

25. Koshizuka S, Oka Y (1996) Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl Sci Eng 123(3):421-434

26. Zhu J, Chen L-Q, Shen J, Tikare V (2001) Computing the effective diffusivity using a spectral method. Mater Sci Eng A 311(1-2):135-141

27. Küntz M, van Mier J, Lavallée P (2001) A lattice gas automaton simulation of the nonlinear diffusion equation: a model for moisture flow in unsaturated porous media. Transp Porous Med 43(2):289-307

28. Vauvre J-M, Patsioura A, Kesteloot R, Vitrac O (2015) Multiscale modeling of oil uptake in fried products. AIChE J 61(7):2329-2353

29. Vitrac O, Hayert M (2007) Effect of the distribution of sorption sites on transport diffusivities: a contribution to the transport of medium-weight-molecules in polymeric materials. Chem Eng Sci 62(9):2503-2521

30. Gorban AN (2014) Detailed balance in micro- and macrokinetics and micro-distinguishability of macro-processes. Results Phys 4:142-147

31. Crank J (1975) The mathematics of diffusion, 2nd edn. Clarendon Press, Oxford

32. Belova IV, Murch GE (2003) The effective diffusivity in two-phase material. Defect Diffus Forum 218-220:79-87

33. Riley MR, Muzzio FJ, Buettner HM, Reyes SC (1994) Monte Carlo calculation of effective diffusivities in two- and three-dimensional heterogeneous materials of variable structure. Phys Rev E 49(4):3500-3503

34. Ronis D (1987) Diffusion in random two-phase media. Phys Rev A 36(4):1908-1928

35. Kalnin JR, Kotomin E (1998) Modified Maxwell-Garnett equation for the effective transport coefficients in inhomogeneous media. J Phys A Math Gen 31(35):7227-7234

36. Kalnin JR, Kotomin EA, Maier J (2002) Calculations of the effective diffusion coefficient for inhomogeneous media. J Phys Chem Solids 63(3):449-456

37. Lejay A, Maire S (2010) Simulating diffusions with piecewise constant coefficients using a kinetic approximation. Comput Methods Appl Mech Eng 199(29-32):2014-2023

38. Duderstadt JJ, Martin WR (1979) Transport theory. Wiley, New York

39. Walsh JB (1978) The local time of the Brownian sheet, in Temps locaux. Astérisque 52-53:47-61

40. Yor M (1997) The Brownian filtration, Tsirel'son's examples, and Walsh's Brownian motions. In: Yor M (ed) Some aspects of Brownian motion: part II: some recent martingale problems. Birkhäuser, Basel, pp 103-123

41. Lejay A (2006) On the constructions of the skew Brownian motion. Probab Surv 3:413-466

42. Risken H (1989) The Fokker-Planck equation: methods of solution and applications, 2nd edn. Springer, Berlin

43. Kloeden PE, Platen E (1992) Numerical solution of stochastic differential equations. Springer, Berlin

44. Kramers HA (1940) Brownian motion in a field of force and the diffusion model of chemical reactions. Physica 7(4):284-304

45. Mel'Nikov VI (1991) The Kramers problem: fifty years of development. Phys Rep 209(1-2):1-71

46. Vestergaard CL, Mikkelsen MB, Reisner W, Kristensen A, Flyvbjerg H (2016) Transition state theory demonstrated at the micron scale with out-of-equilibrium transport in a confined environment. Nat Commun 7(1):10227

47. Chen M-F (2004) From Markov chains to non-equilibrium particle systems. World Scientific Publishing, Singapore

48. Delay F, Ackerer P, Danquigny C (2005) Simulating solute transport in porous or fractured formations using random walk particle tracking. Vadose Zone J 4(2):360

49. Salamon P, Fernández-García D, Gómez-Hernández JJ (2006) A review and numerical assessment of the random walk particle tracking method. J Contam Hydrol 87(3-4):277-305

50. Hoteit H, Mose R, Younes A, Lehmann F, Ackerer P (2002) Three-dimensional modeling of mass transfer in porous media using the mixed hybrid finite elements and the random-walk methods. Math Geol 34(4):435-456

51. Labolle EM, Quastel J, Fogg GE (1998) Diffusion theory for transport in porous media: transition-probability densities of diffusion processes corresponding to advection-dispersion equations. 34(7):1685-1693

52. Balescu R (1997) Statistical dynamics: matter out of equilibrium. Imperial College Press, London

53. Landau DP, Binder K (2000) A guide to Monte Carlo Simulations in statistical physics. Cambridge University Press, Cambridge

54. Marsh CA, Coveney PV (1998) Detailed balance and H-theorems for dissipative particle dynamics. J Phys A Math Gen 31(31):6561-6568

55. Werndl C, Frigg R (2017) Mind the gap: Boltzmannian versus Gibbsian equilibrium. Philos Sci 84(5):1289-1302

56. Cussler EL (1997) Diffusion: mass transfer in fluid systems, 2nd edn. Cambridge University Press, Cambridge

57. Druger SD, Nitzan A, Ratner MA (1983) Dynamic bond percolation theory: a microscopic model for diffusion in dynamically disordered systems. I. Definition and one-dimensional case. J Chem Phys 79(6):3133-3142

58. Dürr O, Volz T, Dieterich W, Nitzan A (2002) Dynamic percolation theory for particle diffusion in a polymer network. J Chem Phys 117(1):441-447

59. Harrison AK, Zwanzig R (1985) Transport on a dynamically disordered lattice. Phys Rev A 32(2):1072-1075

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Terms and Conditions

Springer Nature journal content, brought to you courtesy of Springer Nature Customer Service Center GmbH (“Springer Nature”).

Springer Nature supports a reasonable amount of sharing of research papers by authors, subscribers and authorised users (“Users”), for small-scale personal, non-commercial use provided that all copyright, trade and service marks and other proprietary notices are maintained. By accessing, sharing, receiving or otherwise using the Springer Nature journal content you agree to these terms of use (“Terms”). For these purposes, Springer Nature considers academic use (by researchers and students) to be non-commercial.

These Terms are supplementary and will apply in addition to any applicable website terms and conditions, a relevant site licence or a personal subscription. These Terms will prevail over any conflict or ambiguity with regards to the relevant terms, a site licence or a personal subscription (to the extent of the conflict or ambiguity only). For Creative Commons-licensed articles, the terms of the Creative Commons license used will apply.

We collect and use personal data to provide access to the Springer Nature journal content. We may also use these personal data internally within ResearchGate and Springer Nature and as agreed share it, in an anonymised way, for purposes of tracking, analysis and reporting. We will not otherwise disclose your personal data outside the ResearchGate or the Springer Nature group of companies unless we have your permission as detailed in the Privacy Policy.

While Users may use the Springer Nature journal content for small scale, personal non-commercial use, it is important to note that Users may not:

1. use such content for the purpose of providing other users with access on a regular or large scale basis or as a means to circumvent access control;

2. use such content where to do so would be considered a criminal or statutory offence in any jurisdiction, or gives rise to civil liability, or is otherwise unlawful;

3. falsely or misleadingly imply or suggest endorsement, approval, sponsorship, or association unless explicitly agreed to by Springer Nature in writing;

4. use bots or other automated methods to access the content or redirect messages

5. override any security feature or exclusionary protocol; or

6. share the content in order to create substitute for Springer Nature products or services or a systematic database of Springer Nature journal content.

In line with the restriction against commercial use, Springer Nature does not permit the creation of a product or service that creates revenue, royalties, rent or income from our content or its inclusion as part of a paid for service or for other commercial gain. Springer Nature journal content cannot be used for inter-library loans and librarians may not upload Springer Nature journal content on a large scale into their, or any other, institutional repository.

These terms of use are reviewed regularly and may be amended at any time. Springer Nature is not obligated to publish any information or content on this website and may remove it or features or functionality at our sole discretion, at any time with or without notice. Springer Nature may revoke this licence to you at any time and remove access to any copies of the Springer Nature journal content which have been saved.

To the fullest extent permitted by law, Springer Nature makes no warranties, representations or guarantees to Users, either express or implied with respect to the Springer Nature journal content and all parties disclaim and waive any implied warranties or warranties imposed by law, including merchantability or fitness for any particular purpose.

Please note that these rights do not automatically extend to content, data or other material published by Springer Nature that may be licensed from third parties.

If you would like to use or distribute our Springer Nature journal content to a wider audience or on a regular basis or in any other manner not expressly permitted by these Terms, please contact Springer Nature at

onlineservice@springernature.com