# Single-channel speech enhancement using learnable loss mixup

###### Abstract

Generalization remains a major problem in supervised learning of single-
channel speech enhancement. In this work, we propose learnable loss mixup
(LLM), a simple and effortless training diagram, to improve the generalization
of deep learning-based speech enhancement models. Loss mixup, of which
learnable loss mixup is a special variant, optimizes a mixture of the loss
functions of random sample pairs to train a model on virtual training data
constructed from these pairs of samples. In learnable loss mixup, by
conditioning on the mixed data, the loss functions are mixed using a non-
linear mixing function automatically learned via neural parameterization. Our
experimental results on the VCTK benchmark show that learnable loss mixup
achieves 3.26 PESQ, outperforming the state-of-the-art.

Index Terms: speech enhancement, noise reduction, deep neural network,
convolutional neural network, supervised learning, regression, mixup, loss
mixup, learnable loss mixup

##  1 Introduction

Speech enhancement, the problem of extracting a clean speech signal from its
noisy version, plays a crucial role in speech applications and has been
extensively investigated in the literature. Classical approaches employ signal
processing algorithms [1, 2, 3, 4] and produce satisfactory results in
stationary noise environments. However, they fail in more challenging
scenarios, such as nonstationary noise, and have been surpassed by deep
supervised learning methods [5, 6, 7, 8, 9, 10, 11]. Deep supervised learning
for speech enhancement uses pairs of clean and noisy speech signals to train
an expressive deep neural network (DNN) able to robustly predict clean speech
in highly non-stationary noise environments. They have become the workhorse
for modern speech enhancement systems.

State-of-the-art supervised learning-based speech enhancement in the
literature typically trains the model by solving an Empirical Risk
Minimization (ERM) problem [12]. This approach, however, has a major drawback
of learning the network behavior only on the training samples. This limitation
leads to undesirable behavior of the model beyond the training set. A
straightforward and common remedy is to extend the training set which,
unfortunately, is expensive and laborious to acquire and label in speech
applications. Another common approach is data augmentation which, however, is
domain-dependent and thus requires expert knowledge. Improving generalization,
the network’s ability to produce robust predictions in unseen noise
environments, thus remains a major challenge in speech enhancement.

Our contribution. To tackle this problem, we introduce loss mixup (LM) with a
special variant called learnable loss mixup (LLM). During training, loss mixup
constructs virtual samples by mixing random noisy training samples. It, then,
trains the network on the augmented data by minimizing a mixture of the loss
functions of the individual mixing samples. In essence, loss mixup is not only
a data augmentation method but also an instance of Vicinal Risk Minimization
(VRM) [13], a training regime optimizing the loss function within the vicinity
of the training set, suitable for regression problems such as speech
enhancement. The advantages of loss mixup are threefold. First, it augments
the training data with richer noise and interference conditions via virtual
noisy samples. This data augmentation procedure is domain-independent and can
be implemented in a few lines of code without expert knowledge. Moreover, by
monitoring the network in the vicinity of training samples, loss mixup
encourages smooth regression curves and thus stabilizes the network behavior
outside the training set. Finally, optimizing the loss mixture instead of the
target mixture retains the individual clean targets of the mixing noisy
samples, avoiding the degradation of the network in regression settings [14].

Manually tuning the mixing parameter, sampled from a mixing distribution in
our framework, for the loss mixture is nontrivial. Therefore, we propose
learnable loss mixup to automatically learn the mixing distribution via
gradient descent by reparameterizing the mixing distribution parameters into
loss function variables. Furthermore, we introduce conditions allowing
reparameterization by neural networks which can be trained efficiently via
backpropagation. On the VCTK speech enhancement dataset [15], we show that
learnable loss mixup significantly outperforms standard ERM training. It
achieves 3.26 PESQ, surpassing the previous state-of-the-art by 0.06 points.

Related work. Our method resembles the sample mixing mechanism of mixup [16],
a new training method aiming to improve the generalization of supervised
learning for classification problems. Similar to our loss mixup, mixup, which
we call label mixup in this paper, adopts the Vicinal Risk Minimization
approach by generating virtual samples from random training samples. However,
unlike loss mixup, label mixup trains the network on virtual targets obtained
from the labels of the respective mixing samples. As label mixing is a form of
label smoothing [17], label mixup is beneficial to classification tasks and
has achieved both better generalization and increased model calibration in
classification applications across various data domains such as images, audio,
text, and tabular data [16, 18]. Nevertheless, applying label mixup to
regression problems, such as speech enhancement, is inappropriate since mixing
clean targets produces noisy labels which degrade the network performance in
regression settings [14]. As mentioned above, loss mixup overcomes this issue
by retaining the targets of the mixing samples and optimizing the mixture of
the respective individual loss functions. This insight is supported by our
ablation study, in the experiment section, confirming the superiority of loss
mixup over label mixup in speech enhancement. To our knowledge, we are the
first to propose the sample mixup principle beyond non-classification tasks.

Outline. The rest of the paper is organized as follows. We introduce loss
mixup in Section 2 and draw connection between loss mixup and label mixup. In
Section 3, we present learnable loss mixup by detailing the reparameterization
trick and conditions allowing an efficient neural parameterization of the
mixing function. We discuss experimental results and their implications in
Section 4, and, finally, conclude our work in Section 5.

##  2 Loss mixup

We seek an enhanced version of an observed noisy speech signal $\bm{x}$ via a
learning-based approach. As commonly the case in practice, the noisy speech is
assumed to follow an additive noise model: [tex/equations.tex]eq:noisymodel
where $\bm{s}$ and $\bm{n}$ are the unknown clean speech and the unknown
additive noise, respectively. We train a deep neural network (DNN) $f$ to
extract the clean speech from the noisy speech without noise information.

We adopt a supervised learning regime that trains the DNN on a dataset
$\mathcal{D}=\\{(\bm{x}_{j},\bm{s}_{j})\\}_{j}^{M}$ of $M$ pairs of noisy and
clean speech. Let $\ell{}$ be a loss function enforcing the consistency
between the enhanced speech and the corresponding clean speech, current state-
of-the-art methods typically train the network by solving the empirical risk
minimization (ERM) problem: [tex/equations.tex]eq:erm where
$f\left(\bm{x}_{j}\right)$ is the prediction for noisy sample $\bm{x}_{j}$.
Despite being the workhorse of modern speech enhancement systems, ERM
possesses a major drawback: It causes the network to memorize the training
data by limiting the training on a finite set of samples. Data memorization,
in turn, leads to poor network generalization.

###  2.1 Loss mixup training

To alleviate the generalization issue, we propose loss mixup (LM) training, a
generalized version of mixup [16], for speech enhancement. Given two random
training samples $(\bm{x}_{j},\bm{s}_{j})$ and $(\bm{x}_{k},\bm{s}_{k})$, and
a mixing parameter $\lambda\in[0,1]$, loss mixup constructs virtual noisy
samples [tex/equations.tex]eq:mixup-sample and trains the network by
minimizing [tex/equations.tex]eq:loss-mixup-loss Importantly, the virtual
noisy samples are constructed during the training stage only. During
inference, for a noisy speech signal $\bm{x}$, the enhanced speech is obtained
from the network’s output $f(\bm{x})$.

Similar to label mixup, the success of loss mixup is highly sensitive to the
mixing parameter $\lambda$ [19, 18, 20]. For instance, Zhang et. al. [16]
proposed for $\lambda$ to be drawn from a symmetric Beta distribution
$\beta(\alpha,\alpha)$ which varies dramatically for different values of
$\alpha$ as demonstrated in Fig. LABEL:fig:beta_visualization. Moreover,
several studies [16, 19] indicate that the optimal mixing distribution is
dataset dependent, and manually tuning the mixing distribution to fit a given
dataset is expensive. This challenge motivates our learnable loss mixup
variant in Section 3.

[tex/figures.tex]beta-visualization

###  2.2 Connection to label mixup

Before introducing learnable loss mixup in the next section, here, we draw a
connection between loss mixup and label mixup and highlight that applying them
to speech enhancement yields different model performance.
[tex/theorems.tex]thm:label-v-loss [tex/theorems.tex]proof:label-v-loss
[tex/theorems.tex]corr:mixup-classification

Corollary LABEL:eqn:mixup_classification implies that loss mixup and label
mixup are equivalent for classification and regression tasks utilizing the
cross-entropy and mean squared error loss. However, for application-specific
loss functions such as the log-spectral distance (LSD), a common loss function
in speech applications, loss mixup and label mixup lead to different model
performance for speech enhancement. Our ablation study in Section 4 confirms
this insight.

##  3 Learnable loss mixup

###  3.1 The reparameterization trick

We follow the approach in [16] that samples $\lambda$ from a mixing
distribution $p(\alpha)$. As discussed in the previous section, manually
tuning the mixing distribution is impractical. Furthermore, directly learning
the hyperparameters $\alpha$ of the mixing distribution is unfeasible since
they are absent from the loss function. Therefore, we reparameterize the mixed
loss function (LABEL:eq:loss-mixup-loss) by sampling $\lambda$ from the
uniform distribution $\mathcal{U}(0,1)$ and shaping $\lambda$ with a
parameterized mixing function $\phi$: [tex/equations.tex]eq:learnable-loss-
mixup-loss

Learnable loss mixup (LLM), which trains the network by minimizing
(LABEL:eq:learnable-loss-mixup-loss) instead of (LABEL:eq:loss-mixup-loss),
allows using gradient descent to learn the hyperparameters of the mixing
distribution. In particular, assume $\lambda$ follows an arbitrary
distribution $p(\alpha)$ and let $\phi_{\alpha}$ be the inverse cumulative
distribution function (CDF) of $p(\alpha)$, the reparameterization trick
enables us to calculate Monte Carlo estimates of the loss function
differentiable with respect to (w.r.t.) $\alpha$:

$\begin{split}\nabla_{\alpha}\mathbb{E}_{\lambda\sim
p(\alpha)}\mathcal{L}(\lambda)&=\nabla_{\alpha}\mathbb{E}_{\lambda\sim\mathcal{U}(0,1)}\mathcal{L}(\phi_{\alpha}(\lambda))\\\
&=\mathbb{E}_{\lambda\sim\mathcal{U}(0,1)}\nabla_{\alpha}\mathcal{L}(\phi_{\alpha}(\lambda)).\\\
\end{split}$ (1)

The optimal mixing distribution $p(\alpha)$, then, is obtained by minimizing
the last term of (1) w.r.t. $\alpha$ via gradient descent.

Eq. (1) implies that any inverse CDF is a candidate for the reparametrization
of a mixing distribution. This insight motivates us to define below desired
properties of the mixing function $\phi$. They allow $\phi$ to be
parameterized by neural networks which can be trained efficiently via
backpropagation.

###  3.2 Representation theorem for the mixing function

######  Definition 3.1.

$\phi:[0,1]\rightarrow[0,1]$ is a _mixing function_ if

  1. 1.

$\phi(1)=1$,

  2. 2.

$\forall\lambda\in[0,1]:1-\phi(\lambda)=\phi(1-\lambda),$

  3. 3.

and $\phi$ is monotonically increasing.

The first part of the definition implies loss mixup reverts to ERM when
$\lambda=1$. The second part enforces symmetry: loss mixup produces the same
loss when $(\bm{x}_{j},\bm{s}_{j})$ and $(\bm{x}_{i},\bm{s}_{i})$ are switched
and $\lambda$ becomes $1-\lambda$. The third part ensures that $\phi$ is an
inverse CDF. Our definition implies a general representation for mixing
functions that reduces the parametrization of $\phi$ to a simpler function
$\rho$.

[tex/theorems.tex]thm:representation

[tex/theorems.tex]proof:representation

Note that the reparameterization form (LABEL:eq:representation) of $\phi$
produces a regularization effect given certain choices of $\rho$.
Specifically, Fig. 1 shows the shapes of the mixing function for various
$\rho$ satisfying Theorem LABEL:thm:representation. Convex $\rho$ flattens
$\phi$ near the endpoints, reflecting the principle behind common
regularization methods, such as dropout, that small changes in the input
preserve the output. Conversely, concave $\rho$ flattens $\phi$ near the
midpoint, which reduces the model’s dependence on noisy data by forcing both
small and large mixing to produce similar losses. As oppose to manually tuning
a global mixing distribution that is optimal on average for the whole dataset,
learnable loss mixup automatically learns the optimal $\rho$, and thus $\phi$,
for each mixed input, combining the advantages of both regularization regimes.

Figure 1: Shape of mixing function $\phi$ for different $\rho$.

Importantly, Theorem LABEL:thm:representation permits flexible neural
parameterization of the mixing function $\phi$. Below, we introduce an
instance of such parameterization allowing efficient learning of the mixing
function.

###  3.3 Neural parameterization of the mixing function

We propose a computational efficent parameterization of $\phi$ by conditioning
on the virtual noisy samples. In particular, given a virtual samples
$\tilde{\bm{x}}$ as in (LABEL:eq:mixup-sample), assume our speech enhancement
network produces an embedding $\text{embed}\left(\tilde{\bm{x}}\right)$, we
use [tex/equations.tex]eq:embed-mlp to parameterize $\phi$, where $C>1$,
$\sigma$ is the sigmoid function, and MLP is a Multilayer Perception. An
example of speech enhancement networks producing such embedding
$\text{embed}\left(\tilde{\bm{x}}\right)$ is the UNet described in Table
LABEL:table:unet.

The function in (LABEL:eqn:embed_mlp) satisfies Theorem
LABEL:thm:representation and enables computational efficiency via
backpropagation as $\phi$ in this case is conditioned on an embedding of the
mixed samples. Furthermore, when $C>1$, $\rho$ is Eq. (LABEL:eqn:embed_mlp) is
either convex or concave and thus, following the insight in the previous
subsection, regularizes the network to prevent overfitting.

##  4 Experimental results

In this section, we empirically evaluate and benchmark learnable loss mixup
against state-of-the-art speech enhancement methods on a public speech
enhancement dataset.

###  4.1 Dataset

The VCTK dataset [15] is a popular benchmark for single-channel speech
enhancement [21, 22, 23, 24, 25, 26]. It contains 30 speakers (28 training, 2
test), 15 noise types (10 training, 5 test), and 8 signal-to-noise ratios (4
training, 4 test). We randomly split 1% of the training set for validation,
and evaluate the models on the test set. The noise types and speakers in the
test set are different from those in the training data.

###  4.2 Objective metrics

We use four popular objective measures for speech enhancement to evaluate the
enhanced speech: CSIG, CBAK, COVL, and PESQ. They are computed by comparing
the enhanced speech with the corresponding clean reference of each test
sample. The first three metrics predict the Mean Opinion Score (MOS) that
would result from human perceptual trials: CSIG predicts the signal distortion
MOS; CBAK produces the MOS of background-noise intrusiveness; and COVL
computes the MOS of the overall signal quality. They produce MOS values from
$1$ to $5$. The last metric, PESQ standing for Perceptual Evaluation of Speech
Quality, is a broadly used objective measure for speech quality. Its score is
in the range between -0.5 to 4.5. For all of the metrics, the higher value
corresponds to the better quality of the enhanced speech. In this work, we
focus on PESQ as it is most perceptually consistent with our listening test.

###  4.3 Experimental setup

We predict the log power spectrum of clean speech from that of the
corresponding noisy signal and combine the prediction with the noisy phase to
produce the enhanced speech. To obtain the spectrograms, we downsample the
data to 16kHz, and apply a 512-point STFT with a Hanning window of size 512
and hop length 256. The last frequency bin is removed to form spectrogram
inputs of size 64 x 256 x 1. We train a UNet on the log-spectral distance
(LSD) loss [26] using the Adam optimizer with learning rate 10-4,
$\beta_{1}=0.5$, $\beta_{2}=0.9$, for 450 epochs. In addition, we use batch
size 64 and apply 0.1 $\ell_{2}$-regularization. The UNet, summarized in Table
LABEL:table:unet, has $56$ layers with $20,322,818$ parameters in total. Each
decoder layer operates on a concatenation of the preceding layer and a skip
connection from its corresponding layer in the encoder. Every convolution
layer in the UNet, except the last, is followed by batch normalization,
dropout with 0.1 drop rate, and Leaky ReLU with 0.2 negative slope. For $\rho$
in (LABEL:eqn:embed_mlp), we set $C=5$, the UNet’s bottleneck layer as the
embedding, and an MLP with one hidden layer of width 512.

[tex/tables.tex]mixup-unet

For mixup methods, the virtual noisy samples are generated during the training
process only. To compute the metrics on the test set, for each test sample
($\bm{x}$, $\bm{s}$), we compare the prediction $f(\bm{x})$ against the target
speech $\bm{s}$. The final metric on the test set is the average of the
metrics of the test samples.

###  4.4 Results

We benchmarked our proposed framework against classical and state-of-the-art
speech enhancement methods including MMSE-GAN [23], D+M [24], UNet [25],
Affinity Minimization [26], T-GSA [21], DEMUCS [27], and RHR-Net [22]. Table
LABEL:tab:benchmark-results shows the numerical results of the aforementioned
objective metrics on the test dataset for all benchmarked frameworks. For
reference, we also report the objective measures computed for the noisy test
signals. The results indicated that learnable loss mixup outperformed the
benchmarked methods in all perceptual metrics excepts CBAK. In particular, it
surpassed the state-of-the-art by at least $0.06$ in PESQ.
[tex/tables.tex]benchmark-results [tex/tables.tex]ablation-results

We further conducted an ablation study on VCTK over 5 random seeds, and list
our findings in Table LABEL:table:ablation-results. Here, learnable label
mixup refers to label mixup in which the mixup parameter $\lambda$ in
(LABEL:eq:label-mixup-loss) is substituted by the mixup function $\phi$
parameterized by (LABEL:eqn:embed_mlp). Using ERM as a baseline, we found that
learnable loss mixup , favoring VRM, boosted the performance of our UNet model
significantly from 3.18 to 3.26 PESQ. Furthermore, consistent with our earlier
discussions, label mixup, despite adopting VRM, is inappropriate for speech
enhancement: Learnable label mixup fared significantly worse than ERM at 3.10
PESQ. Finally, loss mixup surpassed ERM slightly by 0.02 points at 3.20 PESQ,
demonstrating the benefits of loss mixup even when the mixing is non-optimal.

##  5 Conclusions

We presented loss mixup, a training framework resembling data augmentation and
Vicinal Risk Minimization, to improve the generalization of deep supervised
learning for speech enhancement. Furthermore, we derived neural
parameterization of the mixing function, resulting in learnable loss mixup, to
efficiently learn the optimal mixing function via gradient descent.
Experimental results showed that our framework outperforms state-of-the-art
speech enhancement methods in the literature. Our results suggested that loss
mixup is a promising approach to boost the generalization of supervised
learning for regression problems.

## References

  * [1] S. F. Boll, “Suppression of acoustic noise in speech using spectral subtraction,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 27, no. 2, pp. 113–120, 1929.
  * [2] J. S. Lim and A. V. Oppenheim, “All-pole modeling of degraded speech,” _IEEE Trans. on Acoustics, Speech, and Signal Processing_ , vol. 26, no. 3, pp. 197–210, 1978.
  * [3] ——, “Enhancement and bandwidth compression of noisy speech,” _Proceedings of the IEEE_ , vol. 67, no. 12, pp. 1586–1604, 1979.
  * [4] Y. Ephraimand and D. Malah, “Speech enhancement using a minimum mean-square error log-spectral amplitude estimator,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 33, no. 2, pp. 443–445, 1985.
  * [5] S. Pascual, A. Bonafonte, and J. Serra, “Segan: Speech enhancement generative adversarial network,” _Interspeech_ , 2017.
  * [6] D. S. Williamson and D. Wang, “Time-frequency masking in the complex domain for speech dereverberation and denoising,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 25, no. 7, 2017.
  * [7] D. Rethage, J. Pons, and X. Serra, “A wavenet for speech denoising,” _ICASSP_ , 2018.
  * [8] F. G. Germain, Q. Chen, and V. Koltun, “Speech denoising with deep feature losses,” _Arxiv_ , 2018.
  * [9] M. H. Soni, N. Shah, and H. A. Patil, “Time-frequency masking- based speech enhancement using generative adversarial network,” _ICASSP_ , 2018.
  * [10] K. Tan and D. Wang, “A convolutional recurrent neural network for real-time speech enhancement,” _Interspeech_ , 2018.
  * [11] A. E. Bulut and K. Koishida, “Low-latency single channel speech enhancement using u-net convolutional neural networks,” _ICASSP_ , 2020.
  * [12] M. Mohri, A. Rostamizadeh, and A. Talwalkar, _Foundations of Machine Learning_. The MIT Press: The MIT Press, 2012.
  * [13] O. Chapelle, J. Weston, L. Bottou, and V. Vapnik, “Vicinal risk minimization,” in _NIPS_ , 2020.
  * [14] D. Tanaka, D. Ikami, T. Yamasaki, and K. Aizawa, “Joint optimization framework for learning with noisy labels,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 5552–5560.
  * [15] C. Valentini-Botinhao, X. Wang, S. Takaki, and J. Yamagishi, “Speech enhancement for a noise-robust text-to-speech synthesis system using deep recurrent neural networks.” 2016.
  * [16] H. Zhang, M. Cisse, Y. N. Dauphin, and D. Lopez-Paz, “mixup: Beyond empirical risk minimization,” _arXiv preprint arXiv:1710.09412_ , 2017.
  * [17] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the inception architecture for computer vision,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 2818–2826.
  * [18] S. Thulasidasan, G. Chennupati, J. A. Bilmes, T. Bhattacharya, and S. Michalak, “On mixup training: Improved calibration and predictive uncertainty for deep neural networks,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 13 888–13 899.
  * [19] H. Guo, Y. Mao, and R. Zhang, “Mixup as locally linear out-of-manifold regularization,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, 2019, pp. 3714–3722.
  * [20] V. Verma, A. Lamb, C. Beckham, A. Najafi, I. Mitliagkas, D. Lopez-Paz, and Y. Bengio, “Manifold mixup: Better representations by interpolating hidden states,” in _International Conference on Machine Learning_ , 2019, pp. 6438–6447.
  * [21] J. Kim, M. El-Khamy, and J. Lee, “T-gsa: Transformer with gaussian-weighted self-attention for speech enhancement,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6649–6653.
  * [22] J. Abdulbaqi, Y. Gu, S. Chen, and I. Marsic, “Residual recurrent neural network for speech enhancement,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6659–6663.
  * [23] M. H. Soni, N. Shah, and H. A. Patil, “Time-frequency masking-based speech enhancement using generative adversarial network,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 5039–5043.
  * [24] J. Yao and A. Al-Dahle, “Coarse-to-fine optimization for speech enhancement,” _arXiv preprint arXiv:1908.08044_ , 2019.
  * [25] A. E. Bulut and K. Koishida, “Low-latency single channel speech enhancement using u-net convolutional neural networks,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6214–6218.
  * [26] D. Tran and K. Koishida, “Single-channel speech enhancement by subspace affinity minimization,” in _INTERSPEECH 2020_. IEEE, 2020.
  * [27] A. Defossez, G. Synnaeve, and Y. Adi, “Real time speech enhancement in the waveform domain,” _arXiv preprint arXiv:2006.12847_ , 2020.


# Spatially heterogeneous dynamics and locally arrested density fluctuations
from first-principles.

J. Lira-Escobedo  Instituto de Física “Manuel Sandoval Vallarta”, Universidad
Autónoma de San Luis Potosí, Álvaro Obregón 64, 78000, San Luis Potosí,
México.  J.R. Vélez-Cordero  Investigadores CONACYT-Instituto de Física
“Manuel Sandoval Vallarta”, Universidad Autónoma de San Luis Potosí, Álvaro
Obregón 64, 78000, San Luis Potosí, México.  Pedro E. Ramírez-González
<EMAIL_ADDRESS>Investigadores CONACYT-Instituto de Física “Manuel
Sandoval Vallarta”, Universidad Autónoma de San Luis Potosí, Álvaro Obregón
64, 78000, San Luis Potosí, México.

###### Abstract

We present a first-principles formalism for studying dynamical heterogeneities
in glass forming liquids. Based on the Non-Equilibrium Self-Consistent
Generalized Langevin Equation theory, we were able to describe the time-
dependent local density profile during the particle interchange among small
regions of the fluid. The final form of the diffusion equation contains both,
the contribution of the chemical potential gradient written in terms of a
coarse-grained density and a collective diffusion coefficient as well as the
effect of a history-dependent mobility factor. With this diffusion equation we
captured interesting phenomena in glass forming liquids such as the cases when
a strong density gradient is accompanied with a very low mobility factor
attributable to the denser part: in such circumstances the density profile
falls into an arrested state even in the presence of a density gradient. On
the other hand, we also show that above a certain critical temperature,which
depends on the volume fraction, any density heterogeneity relaxes to a uniform
state in a finite time, known as equilibration time. We further show that such
equilibration time varies little with the temperature in diluted systems but
can change drastically with temperature in concentrated systems.

###### pacs:

64.70.P, 82.70.Dd

##  I Introduction

The study of structural heterogeneities constitute a common area out of many
different processes occurring near or far from equilibrium. From the
technological viewpoint they are important, for example, in the stability of
bubbly flows Hoefsloot and Krishna (1993); Tyagi and Buwa (2017), food
products de Kruif (2012); Zhu  _et al._ (2019), water purification processes
involving coagulation and flocculation Wilts, Herzberger, and Long (2018) and
in the transport of crude oil along the extraction pipes avoiding asphaltene
deposition Alhosani and Daraboina (2020), to mention a few. On the other hand,
spatial and dynamical heterogeneities have been observed in well controlled
laboratory experiments or numerical simulations involving the formation of
dendritic structures Rátkai, Pusztai, and Gránásy (2019), spontaneous
emulsification Solans, Morales, and Homs (2016) or in glass-forming liquids
near its glass transition temperature or $T_{g}$ Xinyi  _et al._ (2019). Among
the several attempts that have been made in embracing some of these phenomena
in a common foundational theoretical description, the Dynamic Density
Functional Theory (DDFT) arises as one of the most useful theoretical
frameworks Wittkowski and Löwen (2011). Classical Density Functional Theory
(DFT) has been widely applied to study equilibrium properties of uniform and
homogeneous fluids Wu and Li (2007) and it has helped to understand the
complex behavior of non-uniform liquids as wellEvans (1979, 1992) . The
dynamic extension of DFT, initially based on phenomenological considerations
Dieterich, Frisch, and Majhofer (1990) and later derived from first principles
Marconi and Tarazona (2000, 1999); Archer and Evans (2004), has allowed to
include time-dependent problems for the evolution of the density profile. This
dynamic formalism of DFT allows the study of slow dynamics in liquids but it
is unable to describe the rapid increase of the relaxation time in supercooled
liquids Wu and Li (2007). On the other hand, taking into account DFT arguments
together with some other microscopic considerations, Random First Order
Transition (RFOT) theory has emerged as an alternative framework that has
already proved to be very useful in the context of glass-transitionsLubchenko
and Wolynes (2004, 2007).

Another existent theoretical framework for time-dependent problems that allows
the study of glass transitions is the Non-Equilibrium Self-Consistent
Generalized Langevin Equations (NE-SCGLE) theory Ramírez-González and Medina-
Noyola (2010a, b). Such a formalism is based on functional derivatives (or
correlation functions) and its most general version can describe the dynamics
and kinetics of non-uniform systems by solving a closed set of time-dependent
equations for the density profile, correlations functions and particle
mobility. The NE-SCGLE theory, like RFOT, starts with a proper Helmholtz free
energy defining the interparticle interaction potential and then proceeds with
the computation of the structural evolution of liquids under relaxation
processes. It has been shown, for example, that within the NE-SCGLE treatment,
dynamic heterogeneity and aging come naturally from solving such a closed set
of equations Ramírez-González and Medina-Noyola (2010a); Lira-Escobedo  _et
al._ (2019, 2021). In contrast, for other formalism such as RFOT, these
features of glasses are only observed if one introduces strong correlations
between particles Lubchenko and Wolynes (2007). In Lira-Escobedo et al. Lira-
Escobedo  _et al._ (2019, 2021) we have addressed some qualitative difference
and similarities between the NE-SCGLE theory and RFOT. Unfortunately, solving
the equations belonging to the most general version of the NE-SCGLE theory is
a difficult task and so far, it has been solved for simplified but
illustrative problems Olais-Govea, López-Flores, and Medina-Noyola (2015);
Rivas-Barbosa  _et al._ (2018); Lira-Escobedo  _et al._ (2019); Mendoza-Méndez
_et al._ (2017); Ramírez-González and Medina-Noyola (2010b); Sánchez-Díaz,
Ramírez-González, and Medina-Noyola (2013); Olais-Govea  _et al._ (2019). In
particular, the NE-SCGLE theory has been applied in uniform systems under
isochoric conditions and applying temperature jumps (instantaneous change) to
study the structural relaxation. The equilibrium counterpart of this
simplified version of the NE-SCGLE theory, known as SCGLE theory, shares
several similarities with the Mode Coupling Theory (MCT) Elizondo-Aguilera and
Voigtmann (2019). However, the non-equilibrium extension of this theory has
demonstrated to go beyond the limits of MCTLira-Escobedo  _et al._ (2021). A
thoroughly comparison between MCT and SCGLE and an analysis addressing the
differences between MCT and NE-SCGLE theories can be found elsewhere Elizondo-
Aguilera and Voigtmann (2019); Lira-Escobedo  _et al._ (2021). The simplified
version of the NE-SCGLE theory has also been extended to include cooling rate
effects with an exponential approach Lira-Escobedo  _et al._ (2019) and also
it has been extended to study the relaxation of the fictive temperature and
therefore the Kovacs’ kinetics glass transition signaturesLira-Escobedo  _et
al._ (2021).

The aim of this work is to propose a numerical method to solve the most
general version of the NE-SCGLE theory capable to compute the relaxation of
the density profile, including the possibility of dynamic arrest, in a system
that is non-uniform or heterogeneous from the beginning and which is subject
to an instantaneous quench. These initial heterogeneities can appear, in turn,
by small wavelength fluctuations of the density profile induced by the thermal
bath or by long wavelength disturbances given a certain artificial
preparation. In this way, the long-time heterogeneity is seen as a consequence
of the slow dynamics of the density correlation functions. Our formalism is
based on the functional Taylor’s expansion of the chemical potential
introduced by Evans et al. Evans (1979). Using this procedure, we were able to
reproduce well-known results of the non-linear diffusion equation such as the
linear dependence of the chemical potential on the density gradient and the
dependence of the density time derivative on the collective diffusion
coefficient Dieterich, Frisch, and Majhofer (1990); Dieterich, Fulde, and
Peschel (1980); Sadighi and Ganji (2007). In order to apply the theoretical
framework in a concrete example, we chose to study the diffusion of confined
soft-spheres under isothermal conditions. We demonstrate that the density-
relaxation has a strong dependence with the initial density-profile and
temperature. In this sense, we find that a system, initially containing
arrested zones, is unable to relax towards a state with uniform density even
in the presence of density gradients. Our method to calculate the density-
relaxation, given its simplicity, could be expanded to study more complicated
and realistic situations such as including cooling rates or considering fluxes
of heat or matter.

The paper is organized as follows: In Sec. II we introduce a simplified model
to study non-uniform systems. A local chemical state equation, needed in our
picture of non-uniform systems, is derived in Sec. III. In Sec. IV we derive a
method to calculate the local time-dependent mobility for a non-uniform
system. In Sec. V we finally apply our formalism to study the relaxation of
the particle density under isothermal conditions. Lastly, the conclusions and
the numerical methods of this work are presented in Sec. VI and in App. A.

##  II A simplified picture of a non-homogeneous fluid and the NE-SCGLE
equations.

A rather convenient form of visualizing a non-homogeneous fluid is to consider
a two dimensional patchwork model of a supercooled liquid as shown in Fig. 1.
Each cell is described by the local density $n_{i}(t)$ (the figure actually
shows the volume fraction $\varphi_{i}$ instead of the particle density),
where $i=1,2,...,N$, being $N$ the total number of cells and $t$ refers to the
waiting (laboratory) time. Let us assume that the $N$ is large enough so that
we can consider the local density as a continuum variable. Thus, the local
density $n_{i}(t)$ becomes $n(\mathbf{r},t)$, where $\mathbf{r}$ is the
position vector. In what follows we present the NE-SCGLE formalism in order to
describe such a model for heterogeneous fluids.

Figure 1: Color map of the volume fraction (particle density) of a
heterogeneous system. Each cell has a volume fraction of $\varphi_{i}$. The
amplification shows two adjacent cells each of them with a volume fraction of
$\varphi^{0}_{1}$ and $\varphi^{0}_{2}$. $x^{*}$ and $y^{*}$ are defined by
$x/L$ and $y/L$, where $L$ is the length of the system. This figure was built
using R software R Core Team (2020).

The most general version of the NE-SCGLE theory Ramírez-González and Medina-
Noyola (2010a) allows the description of the relaxation of the mean value of
the density profile $\bar{n}(\mathbf{r},t)$ as well as its fluctuations around
this mean using the non-uniform single time matrix covariance
$\sigma(k;\mathbf{r},t)$ and the non-uniform correlation functions
$C(k,\tau;\mathbf{r},t)$ and $C_{S}(k,\tau;\mathbf{r},t)$, where
$\bar{n}(\mathbf{r},t)$, $\sigma(k;\mathbf{r},t)$, $C(k,\tau;\mathbf{r},t)$
and $C_{S}(k,\tau;\mathbf{r},t)$ are coupled by the non-uniform time-dependent
mobility $b(\mathbf{r},t)$ and the non-uniform time-dependent irreducible
memory function $\Delta\zeta(\tau;\mathbf{r},t)$, where $k$ refers to the
Fourier space. NE-SCGLE theory further assumes that the system is isotropic,
thus, the two-point correlation functions depend only on the distance between
such points, $|\mathbf{x}-\mathbf{x}^{\prime}|$. In this way, the inverse
Fourier transform of the two-point correlation functions become
$\sigma(|\mathbf{x}-\mathbf{x}^{\prime}|;\mathbf{r},t)=\overline{\delta
n(\mathbf{r}+\mathbf{x},t)\delta n(\mathbf{r}+\mathbf{x}^{\prime},t)}$,
$C(|\mathbf{x}-\mathbf{x}^{\prime}|,\tau;\mathbf{r},t)=\overline{\delta
n(\mathbf{r}+\mathbf{x},t+\tau)\delta n(\mathbf{r}+\mathbf{x}^{\prime},t)}$
and its self component
$C_{S}(|\mathbf{x}-\mathbf{x}^{\prime}|,\tau;\mathbf{r},t)$ calculated as

$f(|\mathbf{x}-\mathbf{x^{\prime}}|)=\frac{1}{(2\pi)^{3}}\int
d^{3}ke^{-i\mathbf{k}\cdot|\mathbf{x}-\mathbf{x^{\prime}}|}f(k)\;,$ (1)

written for any function $f(x)$. The time-evolution equation for the average
density profile, $\bar{n}(\mathbf{r},t)$, is given by the temporally and
spatially local but non-linear diffusion equation Ramírez-González and Medina-
Noyola (2010a)

$\frac{\partial\bar{n}(\mathbf{r},t)}{\partial
t}=D_{0}\nabla\cdot\bigg{\\{}\bar{n}(\mathbf{r},t)b(\mathbf{r},t)\nabla\bigg{[}\beta\mu[\mathbf{r};\bar{n}(\mathbf{r},t)]\bigg{]}\bigg{\\}},$
(2)

where $D_{0}$ is the short-time diffusion coefficient and
$\beta\mu[\mathbf{r};\bar{n}(\mathbf{r},t)]$ the chemical potential with the
brakets indicating functional dependence on $\bar{n}(\mathbf{r},t)$. The time-
equation for the non-uniform single time matrix covariance
$\sigma(k;\mathbf{r},t)$ reads as follows:

$\displaystyle\frac{\partial\sigma(k;\mathbf{r},t)}{\partial t}$
$\displaystyle=-2k^{2}D_{0}\bar{n}(\mathbf{r},t)b(\mathbf{r},t)\mathcal{E}(k;\bar{n}(\mathbf{r},t))\sigma(k;\mathbf{r},t)$
$\displaystyle+2k^{2}D_{0}\bar{n}(\mathbf{r},t)b(\mathbf{r},t),$ (3)

where
$\mathcal{E}(k;\bar{n}(\mathbf{r},t))\equiv\mathcal{E}(k;\bar{n}(\mathbf{r},t),T(t))$
(assuming that the temperature profile $T(\mathbf{r},t)$ remains uniform) is
the Fourier Transform of the second functional derivative matrix
$\mathcal{E}(|\mathbf{x}-\mathbf{x}^{\prime}|;\bar{n}(\mathbf{r},t),T(t))$ of
the Helmholtz free energy density-functional $\mathcal{F}[n;T]$,

$\mathcal{E}(|\mathbf{x}-\mathbf{x}^{\prime}|;\bar{n}(\mathbf{r},t),T(t))\equiv\bigg{(}\frac{\delta^{2}\mathcal{F}[n;T/k_{B}T]}{\delta
n(\mathbf{\mathbf{r}+\mathbf{x}})\delta
n(\mathbf{\mathbf{r}+\mathbf{x}^{\prime}})}\bigg{)}\bigg{|}_{T,n},$ (4)

evaluated at $n=\bar{n}(\mathbf{r},t)$ and $T=T(t)$. The relaxation of the
non-uniform correlation functions $C(k,\tau;\mathbf{r},t)$ and
$C_{S}(k,\tau;\mathbf{r},t)$ are given by:

$\displaystyle\frac{\partial C(k,\tau;\mathbf{r},t)}{\partial\tau}$
$\displaystyle=\lambda(k;\mathbf{r},t)\Delta\zeta(\tau;\mathbf{r},t)\sigma(k,\mathbf{r},t)$
$\displaystyle-\frac{k^{2}D_{0}C(k,\tau;\mathbf{r},t)\bar{n}(\mathbf{r},t)}{\sigma(k;\mathbf{r},t)}$
$\displaystyle-\lambda(k;\mathbf{r},t)\frac{\partial}{\partial\tau}$
$\displaystyle\int_{0}^{\tau}d\tau\Delta\zeta(\tau-\tau^{\prime};\mathbf{r},t)C(k,\tau^{\prime};\mathbf{r},t),$
(5)

$\displaystyle\frac{\partial C_{S}(k,\tau;\mathbf{r},t)}{\partial\tau}$
$\displaystyle=\lambda(k;\mathbf{r},t)\Delta\zeta(\tau;\mathbf{r},t)$
$\displaystyle-k^{2}D_{0}C_{S}(k,\tau;\mathbf{r},t)$
$\displaystyle-\lambda(k;\mathbf{r},t)\frac{\partial}{\partial\tau}$
$\displaystyle\int_{0}^{\tau}d\tau\Delta\zeta(\tau-\tau^{\prime};\mathbf{r},t)C_{S}(k,\tau^{\prime};\mathbf{r},t),$
(6)

in which the function $\lambda(k;\mathbf{r},t)$ is an interpolating function
defined by

$\lambda(k;\mathbf{r},t)=\frac{1}{1+\bigg{(}\frac{k}{k_{c}(\mathbf{r},t)}\bigg{)}^{2}},$
(7)

where, $k_{c}(\mathbf{r},t)=ck_{max}(\mathbf{r},t)$, with $c\in\mathbb{R}$ and
$k_{max}(\mathbf{r},t)$ being the main peak of the local time-dependent
structure factor
$S(k;\mathbf{r},t)=\sigma(k;\mathbf{r},t)/\bar{n}(\mathbf{r};t)$ at
$\mathbf{r}$ and $t$. Parameter $c$ has been set to 1.109 in all the studies
using NE-SCGLE in order to reproduce the arrest scenario of hard spheres
Juárez-Maldonado  _et al._ (2007).

The irreducible non-uniform memory function $\Delta\zeta(k,\tau;\mathbf{r},t)$
is then computed by: Medina-Noyola (1987)

$\displaystyle\Delta\zeta(\tau;\mathbf{r},t)=\frac{D^{0}}{3(2\pi)^{3}\bar{n}(\mathbf{r},t)}$
$\displaystyle\times\int
d\mathbf{k}k^{2}\bigg{[}\frac{\sigma(k;\mathbf{r},t)/\bar{n}(\mathbf{r},t)-1}{\sigma(k;\mathbf{r},t)}\bigg{]}^{2}C(k;\tau,\mathbf{r},t)C_{S}(k;\tau,\mathbf{r},t).$
(8)

Finally, the time-evolution of the non-uniform mobility $b(\mathbf{r},t)$ is
given by:

$b(\mathbf{r},t)=\frac{1}{1+\int_{0}^{\infty}d\tau\Delta\zeta(\tau;\mathbf{r},t)}.$
(9)

Eqs. (2) - (9) constitute a closed set of equations for the dynamics and
kinetics of a non-uniform system and they are the most general version of the
NE-SCGLE theory for colloidal dynamics. The only external input is the
chemical state equation $\beta\mu[\mathbf{r};\bar{n}(\mathbf{r},t)]$. A
complete description of the theory discussed here can be found elsewhere
Ramírez-González and Medina-Noyola (2010a).

Notice that Eq. (2) reduces to the mean-evolution equation of the dynamic DFT
when the mobility $b(\mathbf{r},t)$ is taken as constant Dieterich, Frisch,
and Majhofer (1990); Marconi and Tarazona (2000, 1999); Archer and Evans
(2004). Let us emphasize that in both, the dynamic DFT and in the NE-SCGLE
theory, the local and non-linear equation for the density profile is obtained
starting from the Langevin equation Ramírez-González and Medina-Noyola
(2010a). Thus, in principle, the NE-SCGLE theory extends the picture of DFT
for colloidal dynamics by including a non-uniform time-dependent mobility.

##  III Local chemical state equation and local diffusion equation

The most general state chemical equation of a heterogeneous fluid can be
written as:Ramírez-González and Medina-Noyola (2010a)

$\beta\mu[\mathbf{r},n(\mathbf{r},t)]=\beta\mu^{*}(\beta)+\mbox{ln}n(\mathbf{r},t)-c[\mathbf{r};n(\mathbf{r},t)]+\beta\psi(\mathbf{r}),$
(10)

where $\psi(\mathbf{r})$ is the external potential acting on the particles at
position $\mathbf{r}$. The terms
$\beta\mu^{*}(\beta)+\mbox{ln}n(\mathbf{r},t)$ are the ideal-gas contribution
to the chemical potential and $-c[\mathbf{r};n(\mathbf{r},t)]$ represents the
deviations of the ideal gas due to the interparticle interactions. Thus, in
the absence of external fields, we can expand
$\beta\mu[\mathbf{r},n(\mathbf{r},t)]$ as a functional Taylor’s series around
a uniform density $n_{u}$, namely

$\beta\mu[\mathbf{r},n(\mathbf{r},t)]=\beta\mu[n_{u}]+\int
d\mathbf{r^{\prime}}\frac{\delta\beta\mu[n(\mathbf{r},t)]}{\delta
n(\mathbf{r^{\prime}},t)}\bigg{|}_{n_{u}}\tilde{n}(\mathbf{r^{\prime}},t)\\\
+\cdots,$ (11)

where $\tilde{n}(\mathbf{r},t)=n(\mathbf{r},t)-n_{u}$ is the deviation from
the uniform density. Keeping up to linear terms in the expansion and deriving
Eq. (10) with respect to $n(\mathbf{r},t)$ to compute the derivative inside
the integral of Eq. (11), we obtain that:

$\beta\mu[\mathbf{r};n(\mathbf{r},t)]]=\beta\mu[n_{u}]+\\\ \int
d\mathbf{r}^{\prime}\bigg{[}\frac{\delta(|\mathbf{r}-\mathbf{r^{\prime}}|)}{n_{u}}-c[\mathbf{r},\mathbf{r^{\prime}};n_{u}]\bigg{]}\tilde{n}(\mathbf{r^{\prime}},t),$
(12)

where $c[\mathbf{r},\mathbf{r^{\prime}};n_{u}]$ is the Ornstein–Zernike direct
correlation function McQuarrie (2000). Defining
$\mathbf{r}=\mathbf{r}+\mathbf{x}$ and
$\mathbf{r}^{\prime}=\mathbf{r}+\mathbf{x}^{\prime}$ and assuming that the
system is isotropic we have that:

$\beta\mu[\mathbf{r},\mathbf{x};n(\mathbf{r},t)]=\beta\mu[n_{u}]+\\\ \int
d\mathbf{x}^{\prime}\bigg{[}\frac{\delta(|\mathbf{x}-\mathbf{x^{\prime}}|)}{n_{u}}-c[|\mathbf{x}-\mathbf{x^{\prime}|;\mathbf{r},n_{u}}]\bigg{]}\tilde{n}(\mathbf{x^{\prime}},\mathbf{r},t),$
(13)

where
$\tilde{n}(\mathbf{x^{\prime}},\mathbf{r},t)=n(\mathbf{r}+\mathbf{x}^{\prime},t)-n_{u}$.
Thus, in order to calculate the chemical potential at $\mathbf{r}$-position we
need to integrate over the $\mathbf{x}$ coordinate:

$\beta\mu[\mathbf{r};n(\mathbf{r},t)]=\frac{1}{V}\int_{V}d\mathbf{x}\beta\mu[\mathbf{r},\mathbf{x};n(\mathbf{r},t)],$
(14)

being $V$ the volume of the system.

Introducing Eq. (13) into Eq. (14) we therefore obtain that:

$\beta\mu[n(\mathbf{r},t)]=\frac{\beta\mu[n_{u}]}{V}\int_{V}d\mathbf{x}+\\\
\frac{\tilde{n}(\mathbf{r},t)}{V}\int_{V}\int_{V}d\mathbf{x}d\mathbf{x^{\prime}}\bigg{[}\frac{\delta(|\mathbf{x}-\mathbf{x^{\prime}}|)}{n_{u}}-c[|\mathbf{x}-\mathbf{x}^{\prime}|;n_{u}]\bigg{]},$
(15)

where we have assumed that
$\tilde{n}(\mathbf{x^{\prime}},\mathbf{r},t)\approx\tilde{n}(\mathbf{r},t)$
and neglected the dependence of $c[|\mathbf{x}-\mathbf{x}^{\prime}|;n_{u}]$ on
$\mathbf{r}$ since it is evaluated at $n_{u}$. It is important to remark that
the first assumption considers that the fluctuation differences in the
vicinity of $\mathbf{r}$ are so small that they can be approximated up to the
zero order term or just considering the position centered at $\mathbf{r}$.
Now, the second term of the last equation can be written as

$\int_{V}\int_{V}d\mathbf{x}d\mathbf{x}^{\prime}\bigg{[}\frac{\delta(|\mathbf{x}-\mathbf{x^{\prime}}|)}{n_{u}}-c[|\mathbf{x}-\mathbf{x}^{\prime}|;n_{u}]\bigg{]}=\\\
V\int_{V}d\mathbf{y}\bigg{[}\frac{\delta(|\mathbf{y}|)}{n_{u}}-c[|\mathbf{y}|;n_{u}]\bigg{]},$
(16)

making the substitution $\mathbf{y}=\mathbf{x}-\mathbf{x^{\prime}}$. The
integral in $\mathbf{y}$ can be readily found using the properties of the
Fourier transform:

$\int_{V}d\mathbf{y}\bigg{[}\frac{\delta(|\mathbf{y}|)}{n_{u}}-c[|\mathbf{y}|;n_{u}]\bigg{]}=\frac{1}{n_{u}}-C^{eq}(k=0;n_{u}),$
(17)

where, for any value of k, the general expression reads as:

$\frac{1}{n_{u}}-C^{eq}(k;n_{u})=\int_{V}d\mathbf{y}e^{-i\mathbf{k}\cdot|\mathbf{y}|}\bigg{[}\frac{\delta(\mathbf{y})}{n_{u}}-c[n_{u},y]\bigg{]}.$
(18)

The function $C^{eq}(k;n_{u})$ is related to the equilibrium structure factor
McQuarrie (2000) by $S^{eq}(k;n_{u})=(1-n_{u}C^{eq}(k;n_{u}))^{-1}$, thus

$\frac{1}{n_{u}}-C^{eq}(k=0;n_{u})=\frac{1}{n_{u}S^{eq}(k=0;n_{u})}.$ (19)

Combining Eqs. (15) to (19) we have that

$\beta\mu[n(\mathbf{r},t)]=\beta\mu[n_{u}]+\frac{1}{n_{u}S^{eq}(k=0;n_{u})}\tilde{n}(\mathbf{r},t).$
(20)

With the result of the Eq. (20) we have that the diffusion equation reads as
follows:

$\frac{\partial\bar{n}(\mathbf{r},t)}{\partial
t}=D^{*}_{0}\nabla\cdot[b(\mathbf{r},t)\bar{n}(\mathbf{r},t)\nabla\bar{n}(\mathbf{r},t)],$
(21)

where $D_{0}^{*}=D_{0}/n_{u}S^{eq}(k=0;n_{u})=D_{coll}(\mathbf{r})/n_{u}$,
being $D_{coll}(\mathbf{r})=D_{0}/S^{eq}(k=0;n_{u})$ the collective diffusion
coefficient that describe long-wavelength density fluctuations behavior
Dieterich, Fulde, and Peschel (1980).

It is instructive to compare the present diffusion expression, Eq. (21), with
other theoretical views that goes from the simple Fickian expression and
passing through the Onsager phenomenological coefficients up to the well-known
Maxwell-Stefan diffusion for multicomponent systems Erminia and Celestino
(2010). As discussed by Wang and LeVan Y and M. Douglas (2008), all these
three theoretical descriptions are somehow equivalent and their selection is
more a matter of applicability (evidently, Fick’s diffusion coefficient is,
for example, reducible in the sense that it can be further decomposed in terms
including Maxwell-Stefan or Onsager coefficients). Comparing with the present
formulation, it is evident that it also shares elements with these other
expressions and equivalences can be defined among them. Notice that Eq. (21)
has exactly the same form of the the standard non-linear diffusion equation
for the relaxation of density with diffusion coefficient
$D(\bar{n}(\mathbf{r},t))=D^{*}_{0}(\mathbf{r})b(\mathbf{r},t)\bar{n}(\mathbf{r},t)$Dieterich,
Frisch, and Majhofer (1990); Sadighi and Ganji (2007). Also, the present
formulation uses the gradient of the chemical potential, as it is done in the
Onsager and Maxwell-Stefan expressions, and also contains a non-local mobility
which is analogous to the friction terms capture in Maxwell and Stefan’s
analysis. In any case, although the present approach was done for
monocomponent systems, it has the advantage that is free of phenomenological
coefficients and its components came directly from the microestructure and
interparticle interactions information.

Eq. (21) together with the set of Eqs. (3)-(9) represent a simplified version
of the NE-SCGLE theory for the relaxation of the average of the density
profile , $\bar{n}(\mathbf{r},t)$, in a heterogeneous systems. However, the
solution of those equations is a difficult task and we require some other
approximations.

##  IV Calculation of the local time-dependent mobility.

Let us mention that $\bar{n}(\mathbf{r},t)$ along with its fluctuations
$\sigma(k,\mathbf{r},t)$ describe the overall macroscopic relaxation in the
sense that if we know $b(\mathbf{r},t)$ for all $\mathbf{r}$ and $t$, we do
not need to solve for the microscopic information provided by Eqs. (5) - (9).
As a first step, lets assume that we are only interested in the relaxation of
$\bar{n}(\mathbf{r},t)$. Furthermore, lets assume that any historical
relaxation path of the fluctuations $\sigma(k,\mathbf{r},t)$ visit, in their
sequence of transitions, the exactly same values as its equilibrium
counterparts obtained in an ordered set of equal $\bar{n}$’s (i.e., we assumed
instantaneous local relaxation for the fluctuations). Formally, we can
decouple $\sigma(k,\mathbf{r},t)$ from the others NE-SCGLE equations by
writing
$\sigma(k,\mathbf{r},t)=\sigma^{eq}(k;\bar{n}(\mathbf{r},t))+\Delta\sigma(k,\mathbf{r},t)$,
where
$\Delta\sigma(k,\mathbf{r},t)\equiv\sigma(k,\mathbf{r},t)-\sigma^{eq}(k;\bar{n}(\mathbf{r},t))$,
and then neglect the non-equilibrium effects provided by
$\Delta\sigma(k,\mathbf{r},t)$. With this assumption we no longer need the
time equation for $\sigma(k,\mathbf{r},t)$. In particular, we can write Eq.
(21) as:

$\frac{\partial\varphi(\mathbf{r},t)}{\partial
t}=D^{*}_{0}\nabla\cdot[b(\varphi(\mathbf{r},t))\varphi(\mathbf{r},t)\nabla\varphi(\mathbf{r},t)],$
(22)

where $\varphi(\mathbf{r},t)$ is the time-dependent volume fraction,
$\varphi(\mathbf{r},t)=\pi a^{3}n(\mathbf{r},t)/6$ for spheres of radius $a$.
Henceforth, we focus on the relaxation of the local time-dependent mobility
$b(\varphi(\mathbf{r},t))$ which can be calculated by solving the equations
for the equilibrium intermediate scattering functions,
$F^{eq}(k,\tau;[\varphi(\mathbf{r},t)]),F^{eq}_{S}(k,\tau;[\varphi(\mathbf{r},t)])$
and the equilibrium memory kernel
$\Delta\zeta^{eq}(\tau;[\varphi(\mathbf{r},t)])$ given by:

$\frac{\partial
F^{eq}(k,\tau;[\varphi(\mathbf{r},t)])}{\partial\tau}=\lambda(k)\Delta\zeta^{eq}(\tau;[\varphi(\mathbf{r},t)])S^{eq}(k,[\varphi(\mathbf{r},t)])\\\
-\frac{k^{2}D_{0}F^{eq}(k,\tau;[\varphi(\mathbf{r},t)])}{S^{eq}(k,[\varphi(\mathbf{r},t)])}\\\
-\lambda(k)\frac{\partial}{\partial\tau}\int_{0}^{\tau}d\tau\Delta\zeta^{eq}(\tau-\tau^{\prime};[\varphi(\mathbf{r},t)])F^{eq}(k,\tau^{\prime};[\varphi(\mathbf{r},t)]),$
(23)

$\frac{\partial
F^{eq}_{S}(k,\tau;[\varphi(\mathbf{r},t)])}{\partial\tau}=\lambda(k)\Delta\zeta^{eq}(\tau;[\varphi(\mathbf{r},t)])\\\
-k^{2}D_{0}F^{eq}_{S}(k,\tau;[\varphi(\mathbf{r},t)])\\\
-\lambda(k)\frac{\partial}{\partial\tau}\int_{0}^{\tau}d\tau\Delta\zeta^{eq}(\tau-\tau^{\prime};[\varphi(\mathbf{r},t)])F^{eq}_{S}(k,\tau^{\prime};[\varphi(\mathbf{r},t)]),$
(24)

$\Delta\zeta^{eq}(\tau;[\varphi(\mathbf{r},t)])=\frac{D_{0}}{6\pi^{2}n(\mathbf{r},t)}\int_{0}^{\infty}dkk^{4}\times\\\
\bigg{[}\frac{S^{eq}(k;[\varphi(\mathbf{r},t)])-1}{S^{eq}(k;[\varphi(\mathbf{r},t)])}\bigg{]}^{2}\\\
\times
F^{eq}(k,\tau;[\varphi(\mathbf{r},t)])F^{eq}_{S}(k,\tau;[\varphi(\mathbf{r},t)]),$
(25)

where $S^{eq}(k,[\varphi(\mathbf{r},t)])$ is the equilibrium structure factor
and the brakets $[\dots]$ inside $F^{eq}$, $F^{eq}_{S}$, $\Delta\zeta^{eq}$
and $S^{eq}$ mean functional dependence on the local volume fraction
$\varphi(\mathbf{r},t)$. $\lambda(k;\mathbf{r},t)$ is again given by:

$\lambda(k;\mathbf{r},t)=\frac{1}{1+\bigg{(}\frac{k}{k_{c}(\mathbf{r},t)}\bigg{)}^{2}},$
(26)

where now $k_{c}(\mathbf{r},t)=ck_{max}(\mathbf{r},t)$ is defined in terms of
the $S^{eq}(k;\varphi(\mathbf{r},t))$ profile (here again $c=1.109$ to
reproduce the dynamics of hard-spheres Juárez-Maldonado  _et al._ (2007)). The
mobility is obtained by:

$b(\varphi(\mathbf{r},t))=\bigg{[}1+\int_{0}^{\infty}d\tau\Delta\zeta(\tau;[\varphi(\mathbf{r},t)])\bigg{]}^{-1}.$
(27)

Let us emphasize that the zero order approximation of $\sigma(k,\mathbf{r},t)$
does not mean that the fluctuations are not important in the relaxation of the
density $\bar{n}(\mathbf{r},t)$. In fact, the fluctuations
$\sigma^{eq}(k,\bar{n}(\mathbf{r},t))=S^{eq}(k;\bar{n}(\mathbf{r},t))/\bar{n}(\mathbf{r},t)$
enter in our picture through the equilibrium mobility
$b(\bar{n}(\mathbf{r},t))$ which contains the information of the interparticle
forces. Thus, our approach introduces a self-consistent method to describe the
relaxation of the density and the equilibrium mobility. The self-consistent
scheme provided by NE-SCGLE theory is the main difference with other
approaches, for instance within dynamic DFT Marconi and Tarazona (1999, 2000)
where a constant mobility is assumed. On the other hand, analytic functions
for the local diffusion equations are required to apply nonlinear diffusion
equations Dieterich, Frisch, and Majhofer (1990); Sadighi and Ganji (2007).

Given the structure of the Eqs. (22) - (27) we can calculate separately the
equilibrium mobility $b(\varphi)$ for any value of the density $\varphi$ and
then use interpolation formulas to compute $b(\varphi)$ inside Eq. (22) and
solve for the time evolution of $\varphi(\mathbf{r},t)$. Given a initial
condition, $\varphi(\mathbf{r},t=0)$, we can built an initial mobility profile
$b(\varphi(\mathbf{r},t=0))$ and then use Eq. (22) to calculate the change in
the density for the marching time $t>0$. In App. A we provide a numerical
method to solve Eq. (22).

##  V A fundamental model of dynamical heterogeneties: isothermal density
relaxation of soft spheres.

In order to develop a fundamental picture of the dynamical heterogeneties in
glass-forming liquids, let us consider again the two dimensional patchwork
model of a supercooled liquid as shown in Fig. 1. At this point we have all
the theoretical elements required for the description of the relaxation
process of the local volume fraction $\varphi(\mathbf{r},t)$. In order to
obtain quantitative results, let us consider two adjacent cells in the patch
work model of Fig 1, each of them characterized by the initial local volume
fraction $\varphi^{0}_{1}$ and $\varphi^{0}_{2}$ respectively, like the ones
showed in the amplification of Fig. 1. Hence, rather than describe all the
diffusion throughout the patchwork, we describe in detail the diffusion along
two adjacent cells; this process can be further spread over all the system. In
this manner we can analyze in a simpler way the main predictions of the NE-
SCGLE theory.

Let us apply the numerical scheme shown in App. A (Eq. (45)) to a one-
dimensional heterogeneous mono-component fluid made of soft spheres of
diameter $a$ immersed in a solvent and interacting by the 6-12 ($\nu=6$)
truncated Lennard-Jones potential, (6-12TLJP) i.e., $\beta
u(r)=\epsilon[(\frac{a}{r})^{2\nu}-2(\frac{a}{r})^{\nu}+1]$ for $r\leq a$ and
$\beta u(r)=0$ for $r\geq a$. In order to calculate the structure factor for
this potential we used the equivalence between soft spheres and hard spheres,
$S^{eq}_{SS}(ka;\varphi,T,\nu)=S^{eq}_{HS}(\lambda_{1}ka;\lambda_{1}^{3}\varphi)$,
where $\lambda_{1}$ is computed with the blip functionHansen and McDonald
(1990),
$\lambda_{1}(T^{*},\nu)=\bigg{\\{}1-3\int_{0}^{1}dxx^{2}\mbox{exp}\bigg{[}-\frac{1}{T^{*}}(\frac{1}{x^{2\nu}}-\frac{2}{x^{\nu}}+1)\bigg{]}\bigg{\\}}$
with $T^{*}$ being the reduced temperature. The structure factor
$S^{eq}_{HS}(\lambda_{1}ka;\lambda_{1}^{3}\varphi)$ was obtained using the
Percus-Yevick closure for the Ornstein-Zernike equation along with the Verlet-
Weis approximation for dense liquidsPercus and Yevick (1958); Verlet and Weis
(1972). The fluid zone and the arrest zone are divided by a critical
temperature $T_{c}^{*}=T_{c}^{*}(\varphi)$. This critical temperature
separates the states that fully equilibrate in any temperature history of
those which are unable to equilibrate regardless the thermal history Lira-
Escobedo  _et al._ (2019, 2021). The transition temperature
$T^{*}_{c}(\varphi)$ was obtained with the equilibrium dynamic arrest
criterion derived from an asymptotic analysis of the memory function
$\Delta\zeta(\tau,t)$ when $\tau\rightarrow\infty$ (Juárez-Maldonado and
Medina-Noyola, 2008a; Sánchez-Díaz, Vizcarra-Rendón, and Juárez-Maldonado,
2009; Juárez-Maldonado and Medina-Noyola, 2008b). The dynamic arrest criterion
reads as follows:

Figure 2: Critical temperature $T_{c}^{*}$ as function of the volume fraction
$\varphi$. The inset shows the equilibrium mobility $b(\varphi)$ for different
temperatures $T^{*}$. The point of arrest occurs at $\varphi_{c}$ where
$b(\varphi_{c})=0$.

$\frac{1}{\gamma}=\frac{1}{6\pi^{2}n}\int_{0}^{\infty}dkk^{4}\\\
\times\frac{[S^{eq}(k;n,T)-1]^{2}\lambda^{2}(k)}{[\lambda(k)S^{eq}(k;n,T)+k^{2}\gamma][\lambda(k)]+k^{2}\gamma]},$
(28)

where $\gamma$ is a localization length. If $\gamma$ is infinite, the system
remains in a fluid-like state, whereas if it becomes finite, the system is
considered arrested. Fig. 2 shows the critical temperature $T_{c}^{*}$ as
function of the volume fraction $\varphi$. Inversely, we can also define the
critical volume fraction $\varphi_{c}=\varphi_{c}(T_{c}^{*})$ as function of
temperature $T^{*}$, having $\varphi_{c}$ the same meaning of $T_{c}^{*}$. The
hard-sphere system is obtained in the limit $T^{*}\rightarrow 0$, giving a
value of $\varphi_{c}=0.583$, in complete agreement with previous results
Pérez-Ángel  _et al._ (2011); Brambilla  _et al._ (2009).

The first input for solving Eq. (22) is the equilibrium mobility $b(\varphi)$
which is calculated separately by means of Eqs. (23) - (27). See App. A for
more details in the numerical solution of Eq. (22). The inset of Fig. 2 shows
the equilibrium mobility $b(\varphi)$ as function of the volume fraction
$\varphi$ for different temperatures of the system of soft-spheres described
above. For the dilute case, given that the interparticle forces are too small,
the mobility $b(\varphi)$ is almost temperature independent and always close
to unity, i.e., the dynamics are controlled by the short-time diffusion
coefficient which only senses the viscosity of the solvent. As the density
increases, the mobility $b(\varphi)$ decreases until it reach $b(\phi)=0$.
This condition occurs at $\varphi_{c}$ and is equivalent to saying that the
system is arrested. From the inset of Fig. 2 it can be seen that as the
temperature increases, the critical volume fraction $\varphi_{c}$ also
increases, as predicted by Eq. (28) in the main plot of Fig. 2.

Once we have calculated the equilibrium mobility $b(\varphi)$ by solving the
discrete form of Eqs. (23) - (27) Juárez-Maldonado and Medina-Noyola (2008b),
the numerical solution of Eq. (45) as described in App. A is straightforward.
Let us consider an initial density profile, $\varphi(x^{*},t^{*}=0)$, given by
the following piece-wise function

$\displaystyle\varphi(x^{*},t^{*}=0)=\left\\{\begin{array}[]{cc}\varphi_{1}^{0}&\hskip
14.22636pt0\leq x^{*}\leq 1/2\\\ &\\\ \varphi_{2}^{0}&\hskip
14.22636pt1/2<x^{*}\leq 1\\\ \end{array}\right.,$ (32)

where $\varphi_{1}^{0}<\varphi_{2}^{0}<1$, $t^{*}=t/\tau_{M}$, $x^{*}=x/L$,
$L$ is the length of the system and $\tau_{M}$ is the macroscopic
characteristic time defined by $\tau_{M}=L^{2}/D^{*}_{0}$. Physically,
$\tau_{M}$ is the time required by a tagged particle to diffuse the available
volume. We used this configuration as a model for two adjacent cell each with
a density $\varphi^{0}_{1}$ and $\varphi^{0}_{2}$ (see the amplification in
Fig. 1). For reference, we refer to the density in $0<x^{*}<1/2$ as
$\varphi_{1}(x^{*},t^{*})$ and for $1/2<x^{*}<1$ as
$\varphi_{2}(x^{*},t^{*})$. By studying the relaxation of a system with an
initial density profile given by Eq. (32), we can understand the relaxation of
a heterogeneous system like the one showed in Fig. 1 since all cell pairs
composing the system showed in Fig. 1 should relax in a similar manner.

Figure 3:  The sigmoid function, $\varphi_{s}(x^{*})$, is plotted for several
values of $r$. As $r$ increases, the sigmoid function $\varphi_{s}(x^{*})$
approaches to $\varphi(x^{*},t^{*}=0)$. Figure 4:  Convergence of the value of
$\varphi_{ef}$ as the nodal numbers, $N_{n}$, increases. The parameters used
were $\varphi^{0}_{1}=0.5$, $\varphi^{0}_{2}=0.62$ and $r=100$. From
$N_{n}=500$ on-wards we have that
$\varphi_{ef}\approx(\varphi_{1}+\varphi_{2})/2=0.56$.

To avoid the point at $x^{*}=1/2$ where the derivative of
$\phi(x^{*},t^{*}=0)$ is undefined, the function $\varphi(x^{*},t=0)$ was
smoothed out according to:

$\varphi_{s}(x^{*})=\frac{\varphi_{1}^{0}e^{r/2}+\varphi_{2}^{0}e^{rx^{*}}}{e^{r/2}+e^{rx^{*}}},$
(33)

where $r$ controls the slope of the function. In Fig. 3 the function
$\varphi_{s}(x^{*})$ is plotted according to Eq. (33) for several values of
$r$ with fixed $\varphi_{1}^{0}=0.5$, $\varphi_{2}^{0}=0.62$ as an example.
Fig. 4 shows, for the same $\varphi^{0}_{1}$ and $\varphi_{2}^{0}$ values but
choosing $r=100$, the convergence of the value of
$\varphi_{ef}=\int_{0}^{1}\varphi_{s}(x^{*})dx^{*}$ as the number of nodal
points used in the discretization of Eq. 22 , $N_{n}$, increases. Notice that
$\varphi_{ef}\approx(\varphi_{1}^{0}+\varphi_{2}^{0})/2$. From Fig. 4 we can
conclude that a suitable choice of $N_{n}$, for the above parameters, is
$N_{n}=500$.

To see how the system relaxes for $t>0$ for the given initial condition, we
monitored the decay of the difference of the average of the densities in each
side of the box defined by

$\varphi^{*}(t)=\frac{\braket{\varphi_{2}(x^{*},t^{*})}-\braket{\varphi_{1}(x^{*},t^{*})}}{\braket{\varphi_{2}(x^{*},t^{*}=0)}-\braket{\varphi_{1}(x^{*},t^{*}=0)}},$
(34)

where

$\braket{\varphi_{1}(x^{*},t^{*})}=2\int_{0}^{\frac{1}{2}}\varphi(x^{*},t^{*})dx^{*}$
(35)

$\braket{\varphi_{2}(x^{*},t^{*})}=2\int_{\frac{1}{2}}^{1}\varphi(x^{*},t^{*})dx^{*}.$
(36)

so that a value of $\varphi^{*}(t)=0$ indicates that a uniform value of the
density has been achieved among the cells.

Fig.5 presents the time evolution of $\varphi^{*}(t^{*})$ for different pairs
of the initial fractions $\varphi^{0}_{1}$ and $\varphi^{0}_{2}$ as well as
for several values of $T^{*}$. The mobility $b(\varphi(x^{*},t^{*}))$, used in
our calculations, was linearly interpolated between the equilibrium mobilities
$b(\varphi_{i-1})$ and $b(\varphi_{i})$, previously calculated with the NE-
SCGLE framework, Eqs. (23)-(27), for some $i=1,\cdots,M$, being $M$ some
integer. To find $i$, the condition
$\varphi_{i-1}<\varphi(x^{*},t^{*})<\varphi_{i}$ must be satisfied; thus

Figure 5:  Relaxation of $\varphi^{*}(t^{*})$ for different initial condition
as function of the temperature $T^{*}$. (a) $\varphi^{0}_{1}=0.2$,
$\varphi^{0}_{2}=0.4$, (b) $\varphi^{0}_{1}=0.2$, $\varphi^{0}_{2}=0.62$ and
(c) $\varphi^{0}_{1}=0.5$, $\varphi^{0}_{2}=0.62$.

$b(\varphi(x^{*},t^{*}))=b(\varphi_{i-1})+\frac{b(\varphi_{i-1})-b(\varphi_{i})}{\varphi_{i-1}-\varphi_{i}}(\varphi(x^{*},t^{*})-\varphi_{i-1}).$
(37)

In general, we can see in Fig. 5 that $\varphi^{*}(t^{*})$ drop to zero (the
cells become uniform) after a certain time for all the chosen parameters
except when $T^{*}=0$. Please notice that since the values of $T^{*}$ are
defined at the beginning of the simulation, physically this means that
whatever the values of $\varphi^{0}_{1}$ and $\varphi^{0}_{2}$ are, the cells
are driven instantaneously to the selected temperature at $t^{*}=0$. When time
starts to march, the zone with higher density and low mobility
($\varphi^{0}_{2}$) starts to diffuse towards the zone with lower densities
and higher mobility until a uniform density is reached. During this
relaxation, the decrease in the density at $1/2<x^{*}<1$ induces an increase
in the average mobility, whereas in the region $0<x^{*}<1/2$ occurs the
inverse process, i.e., the mobility decreases as $\varphi_{1}(x^{*},t^{*})$
goes up. Also notice that an increase of the initial density gradient does not
necessarily translates into a faster density relaxation, as one will expect
due to an increase of the driving force (compare the curves between Figs.5(a)
and (b) where, for the first case the initial volume fraction difference is
0.2 while for the latter is 0.42. This effect is due to the fact that the
mobility of the particles is slower in the denser zone for the case denoted in
Fig. 5(b). The effect of the temperature, on the other hand, is to accelerate
or delay the relaxation of $\varphi^{*}(t^{*})$, as one can see in Fig. 5:
reduction of $T^{*}$ yields a shift in the equilibration time, $t^{*}_{eq}$,
to a higher values, where $t^{*}_{eq}$ is defined as the time that takes
$\varphi^{*}(t^{*})$ to reach zero, see also Fig. 6. This figure illustrates
that for the diluted system composed by the pair $\varphi^{0}_{1}=0.2$,
$\varphi^{0}_{2}=0.4$ (blue line), the shift in $t^{*}_{eq}$ as a function of
temperature is almost negligible. A similar trend happens in the pair
$\varphi^{0}_{1}=0.2$, $\varphi^{0}_{2}=0.62$ in which the difference in
$t^{*}_{eq}$ is more evident but still small. In contrast, for the dense
system having $\varphi^{0}_{1}=0.5$, $\varphi^{0}_{2}=0.62$ (red line), the
magnitude of $t^{*}_{eq}$ changes significantly as the reduced temperature
increases up to a point ($T^{*}\approx 100$) where the curve converges with
the diluted cases. This difference in diffusion between diluted and dense
systems is due to the fact that the former, in overall, have greater mobility
than the latter.

Figure 6: Equilibration time $t^{*}_{eq}$ as function of temperature $T^{*}$
for the initial condition of Fig. 5.

Another characteristic that it is worth noticing in our systems is that the
curves shown in Fig. 5 have the same shape and are only shifted horizontally
for any value of $T^{*}$ except zero as well as for any value of the initial
density gradient. This highlights again the relevance of the mobility
$b(\varphi(x^{*},t^{*}))$ in the relaxation of the density
$\varphi(x^{*},t^{*})$. In previous works using the NE-SCGLE theory, the
mobility has shown to play an important role in finding stationary arrested
states out of equilibrium Ramírez-González and Medina-Noyola (2010b); Sánchez-
Díaz, Ramírez-González, and Medina-Noyola (2013); Olais-Govea, López-Flores,
and Medina-Noyola (2015); Lira-Escobedo  _et al._ (2019, 2021); the value
$b(\varphi(x^{*},t^{*}))=0$ for example means that the particles localized at
$x^{*}$ and being at time $t^{*}$ cannot longer diffuse, leading to a system
trapped in an heterogeneous configuration even if non-zero density gradients
are present. Hence, the system is found in a non-equilibrium arrested
heterogeneous state. Given that we are considering an isothermal process
attached to equilibrium (because we are using the equilibrium version of the
NE-SCGLE theory to calculate $b(\varphi)$), for temperatures $T^{*}$ lower
than $T^{*}_{c}(\varphi^{0}_{2})$ the system will remain static as time
marches with no particle diffusion. In Fig. 5 (b) and (c) the horizontal
dotted pink lines represent the evolution of $\varphi^{*}(t^{*})$ for
$T^{*}=0$. With Eq. (28) we have determined that
$T^{*}_{c}(\varphi^{0}_{2}=0.62)=0.02$, thus, for this particular temperature,
$T^{*}<T^{*}_{c}(\varphi^{0}_{2})$ and the density of the configuration
remains unchanged. Thus, in general, for a high heterogeneous system at
temperature $T^{*}$, no matter the size of the clusters characterizing the
initial heterogeneity, the clusters of density $\varphi(\mathbf{r},t)$ will
relax only if $T^{*}_{c}(\varphi(\mathbf{r},t))<T^{*}$.

##  VI Conclusions

In this work we have shown that the NE-SCGLE theory, a theory that studies
structural relaxation, can be used to solve the spatial and temporal evolution
of density heterogeneities including the possibility of dynamic arrest. In
this way, the present theory can manage the full range of volume fractions,
$\varphi(\mathbf{r},t)$, and reduced temperatures, $T^{*}$, conforming the
dynamic phase diagram: from the low density, free-diffusing regime controlled
by the short-time self-diffusion coefficient and in which all heterogeneities
fade away (uniformity), up to the high density, low mobility regime controlled
by the collective diffusion and in which dynamic arrest can take place.

In order to proceed with the computations we used a model of interacting
spherical soft colloids and a solution strategy that consists in computing the
so-called non-uniform mobility, $b(\mathbf{r},t)$, and the mean volume
fraction in two separate steps. In the first one we computed the mobility
using the NE-SCGLE apparatus and which, in turn, uses the microestructural
information of the colloidal model via correlation functions. The key
assumption proposed here was that the noise or fluctuations around the
transient values of the density are essentially the same as the values around
the corresponding stationary value of $\bar{n}(\mathbf{r})$. Therefore, with
the NE-SCGLE equations we can tabulate the values of $b(\mathbf{r},t)$ and use
them as input in the diffusion equation via interpolation formulas.

In the second step we only require to focus efforts in solving the diffusion
equation and handle the expression of the chemical potential which, again,
recapitulates the microscopic information via the entropic and enthalpic part
(particle interactions) of the total potential. After Taylor expanding the
chemical state function around a uniform density value and assuming that the
spatial gradients of the density fluctuations are small, we end up with an
expression in terms of its value at the uniform density but corrected by the
equilibrium structure factor which takes care of the cumulative many-body
interaction of the particles. These terms are further collected in a
generalized diffusion coefficient. The final form of the diffusion equation
was solved using the finite element formulation (weak form) which has the
benefit in reducing the degree of differentiation of the original equation and
which can be further used as a basis to handle problems in different
applications.

In order to show the efficacy of our method to solve the density relaxation,
we studied in particular the density-decay of a initially non-uniform system
(soft-spheres interacting with the truncated Lennard-Jones potential) under
isothermal conditions. To illustrate how the density relaxes, we divided the
system into two boxes each one having different initial densities and then
monitor the time evolution of the difference between the average density of
the two boxes. For a given set of initial density profiles, we observed that
the system relaxes towards a uniform density in a time period named the
equilibration time. The overall method showed consistency in the sense that
this equilibration time increases as the temperature of the system decreases,
as expected. We further observed that for the same value of $T^{*}$, a diluted
system can diffuse faster than a denser system even if the latter has a higher
initial density gradient: such an effect is a direct consequence of a lower
mobility factor in dense systems. Finally, the present methodology can also
capture instances ($T^{*}<T^{*}_{c}$) where the initial density gradient is
accompanied by a very small mobility factor: in this case the initial
heterogeneity is seen to be time-arrested.

Finally, it is worth saying that although our method to study the diffusion of
glass forming colloids is apparently simple, there is still plenty of room to
bring even closer these kind of theoretical developments to real industrial
applications. The formalism presented here can be expanded easily to other
more complicated situations than the isothermal case, for instance,
considering cooling/heating rates acting over the system. Future efforts
should be directed as well in extending the theory to consider non-uniform
temperature landscapes (heat flow) as well as fluid flow. About the latter, we
anticipate that flow not only changes the mean or stationary quantities such
as the structure factor (case of slow or laminar flow) but also changes the
way the density fluctuations relaxes (case of transient or turbulent flows).

## Acknowledgements

The authors acknowledge the financial support provided by CONACyT México
through grants: Cátedras CONACyT-1631 and CB-2015-01-257636. Besides we
appreciate the technical assistance of Jose Limon and the computational
infrastructure provided by LANIMFE .

##  Appendix A Numerical solution of Eq. (22)

In this appendix we will show a numerical scheme to solve Eq. (22) using the
Finite Element formulation. First, lets write Eq. (22) in dimensionless form
as follows:

$\frac{\partial\varphi(\mathbf{r}^{*},t^{*})}{\partial
t^{*}}=\nabla^{*}\cdot[b(\varphi(\mathbf{r}^{*},t^{*}))\varphi(\mathbf{r}^{*},t^{*})\nabla^{*}\varphi(\mathbf{r}^{*},t^{*})],$
(38)

where $t^{*}=t/\tau_{M}$, $\mathbf{r}^{*}=\mathbf{r}/L$, $L=V^{1/3}$ is
characteristic length and $\tau_{M}$ is the macroscopic characteristic time
defined by $\tau_{M}=L^{2}/D^{*}_{0}$. Physically, $\tau_{M}$ is the time
required by a tagged particle to diffuse the available volume.

Eq. (38) is a non-linear second order partial differential equation. The most
easy, fast and general way for solving Eq. (38) is using numerical methods
such as finite differences or finite element. In the finite element approach,
the variational formulation (weak form) of Eq. (38) in one dimension is given
byLarson and Bengzon (2013)

$\int_{0}^{1}\frac{\partial\varphi(x^{*},t^{*})}{\partial
t^{*}}\varphi_{i}(x^{*})dx^{*}=\bigg{[}b(\varphi(x^{*},t^{*}))\varphi(x^{*},t^{*})\varphi_{i}(x^{*})\bigg{]}_{0}^{1}\\\
-\int_{0}^{1}b(\varphi(x^{*},t^{*}))\varphi(x^{*},t^{*})\frac{\partial\varphi(x^{*},t^{*})}{\partial
x^{*}}\frac{d\varphi_{i}(x^{*})}{dx^{*}}dx^{*}$ (39)

where the region of integration $\mathcal{I}=[0,1]$ is a domain where the
function $\phi(\mathbf{x^{*}},t)$ is piecewise defined and $\varphi_{i}(x)$
are the test functions in one dimension defined in the $i$-interval of a
partition of $\mathcal{I}$Larson and Bengzon (2013). Although the test
functions, $\varphi_{i}(x)$, are in general non orthogonal, we can write any
continuous piece-wise linear function by constructing an expansion based on
$\varphi_{i}(x)$. Thus,

$\varphi(x^{*},t^{*})=\sum_{j=0}^{N_{n}-1}\eta_{j}(t^{*})\varphi_{j}(x^{*}),$
(40)

where $N_{n}$ are the total number of nodes in the partition of $\mathcal{I}$.
Substituting the expansion of $\varphi(x^{*},t^{*})$ and considering that
there is no flux of matter at the outer boundaries, i.e., the first term of
the r.h.s. of Eq. (39) is zero, hence,

$\sum_{j=0}^{N_{n}-1}\frac{d\eta_{j}(t^{*})}{dt^{*}}\int_{0}^{1}\varphi_{j}(x^{*})\varphi_{i}(x^{*})dx^{*}=\\\
-\sum_{j=0}^{N_{n}-1}\eta_{j}(t^{*})\int_{0}^{1}b(n(x^{*},t^{*}))\varphi(x^{*},t^{*})\frac{d\varphi_{j}(x^{*})}{dx^{*}}\frac{d\varphi_{i}(x^{*})}{dx^{*}}dx^{*},$
(41)

where we have only used the expansion of $\varphi(x^{*},t^{*})$ in the
derivative of the r. h. s. of Eq. (39) since we have a non-linear differential
equationLarson and Bengzon (2013). In matrix form the last equation can be
written as

$A\dot{\eta}(t)=-B(t)\eta(t)$ (42)

where

$A_{ij}=\int_{0}^{1}\varphi_{i}(x^{*})\varphi_{j}(x^{*})dx^{*},$ (43)

and

$B_{ij}(t)=\int_{0}^{1}b(\varphi(x^{*},t^{*}))\phi(x^{*},t^{*})\frac{d\varphi_{i}(x^{*})}{dx^{*}}\frac{d\varphi_{j}(x^{*})}{dx^{*}}dx^{*},$
(44)

with $i,j=1,\cdots,N_{n}-1$. In order to obtain the discrete time form of Eq.
(42) we can make $\dot{\eta}(t)\approx(\eta_{l+1}-\eta_{l})/\Delta t$, where
$\eta_{l}=\eta(t_{0}+l\Delta t)$, with $l=0,1...,$ being $\Delta t$ the time
step and $t_{0}=0$. Applying this approximation for $\dot{\eta}(t)$ and doing
$B_{l}=B_{l-1}$ (it is understood that $B_{l}=B(t_{0}+l\Delta t)$), the
discrete form of Eq. (42) ends in the following form:

$\bigg{[}A+\Delta tB_{l}\bigg{]}\eta_{l+1}=A\eta_{l}.$ (45)

## References

  * Hoefsloot and Krishna (1993) H. Hoefsloot and R. Krishna, “Influence of gas density on the stability of homogeneous flow in bubble columns,” Ind. Eng. Chem. Res. 32, 747–750 (1993).
  * Tyagi and Buwa (2017) P. Tyagi and V. Buwa, “Dense gas-liquid-solid flow in a slurry bubble column: measurement of dynamic characteristics, gas volume fraction and bubble size distribution,” Chem. Eng. Sci. 173, 326–362 (2017).
  * de Kruif (2012) C. K. de Kruif, “Concluding remarks: the future of soft matter and food structure,” Faraday Discuss 158, 523–527 (2012).
  * Zhu  _et al._ (2019) Q. Zhu, Y. Pan, X. Jia, J. Li, M. Zhang, and L. Yin, “Review on the stability mechanism and application of water-in-oil emulsions encapsulating various additives,” Comp. Rev. Food Sci. F. 18, 1660–1675 (2019).
  * Wilts, Herzberger, and Long (2018) E. Wilts, J. Herzberger, and T. Long, “Addressing water scarcity: cationic polyelectrolytes in water treatment and purification,” Polym. Int. 67, 799–814 (2018).
  * Alhosani and Daraboina (2020) A. Alhosani and N. Daraboina, “Unified model to predict asphaltene deposition in production pipelines,” Energy Fuels 34, 1720–1727 (2020).
  * Rátkai, Pusztai, and Gránásy (2019) L. Rátkai, T. Pusztai, and L. Gránásy, “Phase-field lattice boltzmann model for dendrites growing and moving in melt flow,” npj Computational Materials 5, 113–10 (2019).
  * Solans, Morales, and Homs (2016) C. Solans, D. Morales, and M. Homs, “Spontaneous emulsification,” Curr. Opin. Colloid In. 22, 88–93 (2016).
  * Xinyi  _et al._ (2019) W. Xinyi, X. Wen-Sheng, Z. Hao, and D. Jack F., “Universal nature of dynamic heterogeneity in glass-forming liquids: A comparative study of metallic and polymeric glass-forming liquids,” J. Chem. Phys. 151, 184503 (2019).
  * Wittkowski and Löwen (2011) R. Wittkowski and H. Löwen, “Dynamical density functional theory for colloidal particles with arbitrary shape,” Molecular Physics 109, 2935–2943 (2011).
  * Wu and Li (2007) J. Wu and Z. Li, “Density-functional theory for complex fluids,” Annu. Rev. Phys. Chem. 58, 85–112 (2007).
  * Evans (1979) R. Evans, “The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids,” Advances in physics 28, 143–200 (1979).
  * Evans (1992) R. Evans, “Density functionals in the theory of nonuniform fluids,” Fundamentals of inhomogeneous fluids 1, 85–176 (1992).
  * Dieterich, Frisch, and Majhofer (1990) W. Dieterich, H. Frisch, and A. Majhofer, “Nonlinear diffusion and density functional theory,” Zeitschrift für Physik B Condensed Matter 78, 317–323 (1990).
  * Marconi and Tarazona (2000) U. M. B. Marconi and P. Tarazona, “Dynamic density functional theory of fluids,” Journal of Physics: Condensed Matter 12, A413 (2000).
  * Marconi and Tarazona (1999) U. M. B. Marconi and P. Tarazona, “Dynamic density functional theory of fluids,” The Journal of chemical physics 110, 8032–8044 (1999).
  * Archer and Evans (2004) A. J. Archer and R. Evans, “Dynamical density functional theory and its application to spinodal decomposition,” The Journal of chemical physics 121, 4246–4254 (2004).
  * Lubchenko and Wolynes (2004) V. Lubchenko and P. G. Wolynes, “Theory of aging in structural glasses,” The Journal of chemical physics 121, 2852–2865 (2004).
  * Lubchenko and Wolynes (2007) V. Lubchenko and P. G. Wolynes, “Theory of structural glasses and supercooled liquids,” Annu. Rev. Phys. Chem. 58, 235–266 (2007).
  * Ramírez-González and Medina-Noyola (2010a) P. Ramírez-González and M. Medina-Noyola, “General nonequilibrium theory of colloid dynamics,” Phys. Rev. E 82, 061503 (2010a).
  * Ramírez-González and Medina-Noyola (2010b) P. Ramírez-González and M. Medina-Noyola, “Aging of a homogeneously quenched colloidal glass-forming liquid,” Phys. Rev. E 82, 061504 (2010b).
  * Lira-Escobedo  _et al._ (2019) J. Lira-Escobedo, D. Varela-Varela, P. Mendoza-Méndez, and P. Ramírez-González, “First-principles prediction of multiple stationary states in glass-forming liquids,” The Journal of chemical physics 151, 234501 (2019).
  * Lira-Escobedo  _et al._ (2021) J. Lira-Escobedo, P. Mendoza-Méndez, M. Medina-Noyola, G. McKenna, and P. Ramírez-González, “On a fundamental description of the kovacs’ kinetic signatures in glass-forming systems,” The Journal of Chemical Physics 155, 014503 (2021).
  * Olais-Govea, López-Flores, and Medina-Noyola (2015) J. M. Olais-Govea, L. López-Flores, and M. Medina-Noyola, “Non-equilibrium theory of arrested spinodal decomposition,” The Journal of Chemical Physics 143, 174505 (2015), https://doi.org/10.1063/1.4935000 .
  * Rivas-Barbosa  _et al._ (2018) R. Rivas-Barbosa, E. Lázaro-Lázaro, P. Mendoza-Méndez, T. Still, V. Piazza, P. E. Ramírez-González, M. Medina-Noyola, and M. Laurati, “Different routes into the glass state for soft thermo-sensitive colloids,” Soft matter 14, 5008–5018 (2018).
  * Mendoza-Méndez  _et al._ (2017) P. Mendoza-Méndez, E. Lázaro-Lázaro, L. E. Sánchez-Díaz, P. E. Ramírez-González, G. Pérez-Ángel, and M. Medina-Noyola, “Crossover from equilibration to aging: Nonequilibrium theory versus simulations,” Phys. Rev. E 96, 022608 (2017).
  * Sánchez-Díaz, Ramírez-González, and Medina-Noyola (2013) L. E. Sánchez-Díaz, P. Ramírez-González, and M. Medina-Noyola, “Equilibration and aging of dense soft-sphere glass-forming liquids,” Phys. Rev. E 87, 052306 (2013).
  * Olais-Govea  _et al._ (2019) J. M. Olais-Govea, L. López-Flores, J. B. Zepeda-López, and M. Medina-Noyola, “Interference between the glass, gel, and gas-liquid transitions,” Scientific reports 9, 1–7 (2019).
  * Elizondo-Aguilera and Voigtmann (2019) L. F. Elizondo-Aguilera and T. Voigtmann, “Glass-transition asymptotics in two theories of glassy dynamics: Self-consistent generalized langevin equation and mode-coupling theory,” Physical Review E 100, 042601 (2019).
  * Dieterich, Fulde, and Peschel (1980) W. Dieterich, P. Fulde, and I. Peschel, “Theoretical models for superionic conductors,” Advances in Physics 29, 527–605 (1980).
  * Sadighi and Ganji (2007) A. Sadighi and D. Ganji, “Exact solutions of nonlinear diffusion equations by variational iteration method,” Computers & Mathematics with Applications 54, 1112–1121 (2007).
  * R Core Team (2020) R Core Team,  _R: A Language and Environment for Statistical Computing_, R Foundation for Statistical Computing, Vienna, Austria (2020).
  * Juárez-Maldonado  _et al._ (2007) R. Juárez-Maldonado, M. Chávez-Rojo, P. Ramírez-González, L. Yeomans-Reyna, and M. Medina-Noyola, “Simplified self-consistent theory of colloid dynamics,” Physical Review E 76, 062502 (2007).
  * Medina-Noyola (1987) M. Medina-Noyola, “The generalized langevin equation as a contraction of the description. an approach to tracer diffusion,” Faraday Discuss. Chem. Soc. 83, 21–31 (1987).
  * McQuarrie (2000) D. McQuarrie,  _Statistical Mechanics_ (University Science Books, 2000).
  * Erminia and Celestino (2010) L. Erminia and A. Celestino, “On the maxwell-stefan approach to diffusion: A general resolution in the transient regime for one-dimensional systems,” J. Phys. Chem. B 114, 151––164 (2010).
  * Y and M. Douglas (2008) W. Y and L. M. Douglas, “Mixture diffusion in nanoporous adsorbents: Equivalence of fickian and maxwell-stefan approaches,” J. Phys. Chem. B 112, 8600––8604 (2008).
  * Hansen and McDonald (1990) J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids_ (Elsevier, 1990).
  * Percus and Yevick (1958) J. K. Percus and G. J. Yevick, “Analysis of classical statistical mechanics by means of collective coordinates,” Physical Review 110, 1 (1958).
  * Verlet and Weis (1972) L. Verlet and J.-J. Weis, “Equilibrium theory of simple liquids,” Physical Review A 5, 939 (1972).
  * Juárez-Maldonado and Medina-Noyola (2008a) R. Juárez-Maldonado and M. Medina-Noyola, “Alternative view of dynamic arrest in colloid-polymer mixtures,” Phys. Rev. Lett. 101, 267801 (2008a).
  * Sánchez-Díaz, Vizcarra-Rendón, and Juárez-Maldonado (2009) L. E. Sánchez-Díaz, A. Vizcarra-Rendón, and R. Juárez-Maldonado, “Erratum: Ionic and wigner glasses, superionic conductors, and spinodal electrostatic gels: Dynamically arrested phases of the primitive model [phys. rev. lett. 103, 035701 (2009)],” Phys. Rev. Lett. 103, 179902 (2009).
  * Juárez-Maldonado and Medina-Noyola (2008b) R. Juárez-Maldonado and M. Medina-Noyola, “Theory of dynamic arrest in colloidal mixtures,” Phys. Rev. E 77, 051503 (2008b).
  * Pérez-Ángel  _et al._ (2011) G. Pérez-Ángel, L. E. Sánchez-Díaz, P. E. Ramírez-González, R. Juárez-Maldonado, A. Vizcarra-Rendón, and M. Medina-Noyola, “Equilibration of concentrated hard-sphere fluids,” Physical Review E 83, 060501 (2011).
  * Brambilla  _et al._ (2009) G. Brambilla, D. El Masri, M. Pierno, L. Berthier, L. Cipelletti, G. Petekidis, and A. B. Schofield, “Probing the equilibrium dynamics of colloidal hard spheres above the mode-coupling glass transition,” Physical review letters 102, 085703 (2009).
  * Larson and Bengzon (2013) M. G. Larson and F. Bengzon,  _The finite element method: theory, implementation, and applications_ , Vol. 10 (Springer Science & Business Media, 2013).