Non-Markovian out-of-equilibrium dynamics: A general numerical procedure to construct time-dependent memory kernels for coarse-grained observables

TL;DR

This paper introduces a numerical method to reconstruct non-stationary memory kernels from experimental or simulation data, enabling better understanding of non-equilibrium dynamics in coarse-grained systems.

Contribution

The authors develop a series expansion approach to compute time-dependent memory kernels from two-time correlation functions, applicable to non-stationary systems.

Findings

Memory kernel duration exceeds particle attachment timescales in crystallization.

Method successfully applied to Lennard-Jones melt crystallization.

Provides insights into non-Markovian effects in out-of-equilibrium processes.

Abstract

We present a numerical method to compute non-equilibrium memory kernels based on experimental data or molecular dynamics simulations. The procedure uses a recasting of the non-stationary generalized Langevin equation, in which we expand the memory kernel in a series that can be reconstructed iteratively. Each term in the series can be computed based solely on knowledge of the two-time auto-correlation function of the observable of interest. As a proof of principle, we apply the method to crystallization from a super-cooled Lennard Jones melt. We analyze the nucleation and growth dynamics of crystallites and observe that the memory kernel has a time extent that is about one order of magnitude larger than the typical timescale needed for a particle to be attached to the crystallite in the growth regime.

Equations

$$\frac{d{A}_t}{dt}$$

$$\frac{\partial C(t^{\prime },t)}{\partial t}$$

$$C(t^{\prime },t)=C(t^{\prime },t^{\prime })+{\int }_{t^{\prime }}^{t}d\tau C(t^{\prime },\tau )\mathcal{J}(\tau ,t)$$

$$\mathcal{J}(t^{\prime },t)=\omega (t^{\prime })+{\int }_{t^{\prime }}^{t}d\tau K(t^{\prime },\tau )$$

$$\mathcal{J}(t^{\prime },t)=j_0(t^{\prime },t)+{\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_0(t^{\prime },\tau )\mathcal{J}(\tau ,t)$$

$${\mathcal{S}}_0(t^{\prime },t)$$

$$j_0(t^{\prime },t)$$

$$\mathcal{J}(t^{\prime },t)=j_0(t^{\prime },t)$$

$$+{\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_1(t^{\prime },\tau )\mathcal{J}(\tau ,t)$$

$$\mathcal{J}(t^{\prime },t)=j_0(t^{\prime },t)+{\int }_{t^{\prime }}^{t}d\tau \mathcal{S}(t^{\prime },\tau )j_0(\tau ,t)$$

$${\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_n(t^{\prime },\tau ){\mathcal{S}}_m(\tau ,t)={\mathcal{S}}_{n+m+1}(t^{\prime },t),$$

$$\tilde{A}(\mathbf{\Gamma },t)\equiv \left[A(\mathbf{\Gamma })-{\mu }_A(t)\right]/{\sigma }_A(t)$$

$$\frac{d{\tilde{A}}_t}{dt}={\int }_0^{t}d\tau \tilde{K}(\tau ,t){\tilde{A}}_{\tau }+{\tilde{\eta }}_t$$

$$\frac{d{A}_t}{dt}=$$

$$ϵ(s):=\frac{1}{s(T-s)}{\int }_0^{s}d{t}^{\prime }{\int }_s^{T}dt\left|1-\frac{C(t^{\prime }s)C(s,t)}{C(t^{\prime },t)}\right|,$$

$$\mathcal{J}(t^{\prime },t)=j_0(t^{\prime },t)+{\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_0(t^{\prime },\tau )\mathcal{J}(\tau ,t)$$

$${\mathcal{S}}_0(t^{\prime },t)$$

$$j_0(t^{\prime },t)$$

$$\mathcal{J}(t^{\prime },t)=$$

$$+{\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_0(t^{\prime },\tau ){\int }_{\tau }^{t}d{\tau }^{\prime }{\mathcal{S}}_0(\tau ,{\tau }^{\prime })\mathcal{J}({\tau }^{\prime },t)$$

$$=$$

$$+{\int }_{t^{\prime }}^{t}d{\tau }^{\prime }\left[{\int }_{t^{\prime }}^{{\tau }^{\prime }}d\tau {\mathcal{S}}_0(t^{\prime },\tau ){\mathcal{S}}_0(\tau ,{\tau }^{\prime })\right]\mathcal{J}({\tau }^{\prime },t)$$

$$=$$

$$+{\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_1(t^{\prime },\tau )\mathcal{J}(\tau ,t)$$

$${\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_1(t^{\prime },\tau )\mathcal{J}(\tau ,t)=$$

$$+{\int }_{t^{\prime }}^{t}d\tau$$

$$={\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_1$$

$$+{\int }_{t^{\prime }}^{t}d{\tau }^{\prime }$$

$$={\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_1(t^{\prime },\tau )j_0(\tau ,$$

$${\mathcal{S}}_n(t^{\prime },t)={\int }_{t^{\prime }}^{t}d{\tau }_1{\int }_{t^{\prime }}^{{\tau }_1}d{\tau }_2\cdots {\int }_{t^{\prime }}^{{\tau }_{n-1}}d{\tau }_n{\mathcal{S}}_0(t^{\prime },{\tau }_n){\mathcal{S}}_0({\tau }_n,{\tau }_{n-1})\cdots {\mathcal{S}}_0({\tau }_2,{\tau }_1){\mathcal{S}}_0({\tau }_1,t)$$

$$\mathcal{J}(t^{\prime },t)=j_0(t^{\prime },t)+\sum _{n=0}^{\infty }{\int }_{t^{\prime }}^{t}d\tau {\mathcal{S}}_n(t^{\prime },\tau )j_0(\tau ,t)$$

$$C(t^{\prime },t)=C(t^{\prime },t^{\prime })\exp \left({\int }_{t^{\prime }}^{t}d\tau \gamma (\tau )\right)$$

$${\mathcal{S}}_0(t^{\prime },t)$$

$${\mathcal{S}}_1(t^{\prime },t)$$

$$=\gamma (t^{\prime }){\int }_{t^{\prime }}^{t}d\tau \gamma (\tau )e^{{\int }_{t^{\prime }}^{t}d\tau \gamma (\tau )}$$

$${\mathcal{S}}_2(t^{\prime },t)$$

$$=-\gamma (t^{\prime }){\int }_{t^{\prime }}^{t}d\tau {\int }_{t^{\prime }}^{\tau }d{\tau }^{\prime }\gamma ({\tau }^{\prime })\gamma (\tau )e^{{\int }_{t^{\prime }}^{t}d\tau \gamma (\tau )}$$

$${\mathcal{S}}_n(t^{\prime },t)=(-1{)}^{n+1}$$

$$\times {\int }_{t^{\prime }}^t$$

Full text

Introduction - The dynamics of complex many-body systems are often modelled in terms of the effective dynamics of a small set of relevant observables. Depending on the context, these observables are called ’reaction coordinates’ or ’order parameters’Peters (2016a); Rohrdanz et al. (2013); Sittel and Stock (2018). E.g. in biophysics one might be interested in the evolution of the shape of a protein during a folding experiment, but not in the motion of every single water molecule. Or in materials science, one might model the dynamics of a phase transition in terms of a suitable order parameter without resolving the details of the microscopic motion of the atoms. In the 1960’s Mori and Zwanzig developed a projection operator formalism to derive the equation of motion of such averaged observables, the Generalized Langevin Equation (GLE) Zwanzig (1961); Mori (1965). The GLE is valid only if the density of microstates is stationary. Based on the same formalism, but for non-stationary densities of microstates Grabert (1982), we have recently derived a general structure for the equation of motion of reaction coordinates Meyer et al. (2019, 2017). The structure applies to any dynamical process for which the microscopic equations of motion are deterministic and for any phase-space observable, even if it contains an explicit dependence on time. The resulting non-stationary Generalized Langevin Equation (nsGLE) is thus the equation that needs to be solved, if one wishes to derive reaction coordinate dynamics outside of thermal equilibrium.

The GLE and the nsGLE contain an effective friction term that is non-local in time and integrates over the history of the process. The function that controls these history effects is called “the memory kernel”. The non-locality in time makes the analysis of the GLE and the nsGLE mathematically inconvenient. In applications the memory kernel is therefore often approximated by a Dirac delta-distribution such that a time-local Langevin Equation is recovered Huopaniemi et al. (2006); Hörger et al. (2008); Espanol and Löwen (2009); Shea and Kreuzer (2011); Knoch and Speck (2015). This assumption, that we refer to as the “Markovian approximation”, is however in practice often not verified before it is used, which potentially leads to inaccurate or wrong results. To go beyond this approximation, there are two possible routes. The first one consists in searching (or constructing) a set of new reaction coordinates for which a Markovian description is correct Peters and Trout (2006); Rohrdanz et al. (2011); Pérez-Hernández et al. (2013); Peters (2016b); Chodera and Noé (2014). This method is often useful, but the reaction coordinates constructed may either be abstract quantities, which are not guaranteed to be accessible in experiment, or a large set of them might be required, which implies that the Markovian Langevin equation needs to be solved in a high dimensional space Hegger and Stock (2009); Schaudinnus et al. (2015, 2016). The second route consists in keeping the original observable despite its non-Markovian dynamics, and finding a way to evaluate the corresponding memory kernel. Various methods in this spirit do already exist, but they are restricted to processes with stationary microstate distributions, i.e. they assume that the memory kernel is invariant under translations in time. We present here a general method to compute memory kernels for arbitrary processes and for arbitrary observables from MD simulations or experimental data, i.e. a procedure to analyze the out-of-equilibrium dynamics of reaction coordinates and to construct their equation of motion.

The main idea of the procedure is to recast the nsGLE in a form in which the memory kernel is expressed as a sum of convolution products. Each term is determined recursively from the previous ones, and it probes the behaviour of memory at successively longer times. The discrepancies from the Markovian limit can be assessed by analyzing the functional shape of the computed terms and their relative timescales. The method neither requires modeling nor approximation and uses as a single input the two-time auto-correlation function of the observable under study, which is easily accessible in simulations or experiments. We also note that the computational time of the procedure scales linearly with the number of terms computed in the expansion.

As a proof of principle, we apply this new technique to the process of crystallization of a super-cooled Lennard-Jones melt. Crystallization processes are usually described in the framework of Transition State Theory, which is based on a Markovian approximation. We test this assumption here and reconstruct the memory kernel of the nsGLE which governs the evolution of nucleation and growth of crystalline clusters. We observe significant memory effects.

Numerical method - Consider a system of $N\gg 1$ degrees of freedom $\left\{\Gamma_{i}\right\}_{i\leq N}$ that evolve according to deterministic microscopic equations of motion (e.g. Hamilton’s equations of motion), and a phase-space observable $A(\mathbf{\Gamma})$ that is fully determined by the location $\mathbf{\Gamma}$ in phase-space. Next, take an ensemble, i.e. a large number of copies of the system, and allow these copies to be initialized according to any phase-space distribution, in particular a non-stationary one. The microscopic equations of motion as well as the observable $A$ can also be explicitly time-dependent. These formal definitions are general enough to encompass a very broad spectrum of processes. We showed in ref. Meyer et al. (2019) that for any such observable and for any dynamical process, regardless of how far from equilibrium it evolves, one can always define functions $\omega(t)$, $K(t^{\prime},t)$ and $\eta_{t}$ such that the equations of motion for $A_{t}$ and for its auto-correlation function $C(t^{\prime},t)=\left\langle A^{*}(t^{\prime})A(t)\right\rangle$ are

[TABLE]

where $K(\tau,t)$ is the memory kernel and the average is taken over the ensemble of non-equilibrium trajectories. The time-dependence as a subscript denotes the dependence on a single trajectory, whereas the time-dependence between parentheses indicates a fixed function of time independent of the trajectory. In particular we showed $\omega(t)=\text{d}\left(\ln\sqrt{C(t,t)}\right)/\text{d}t$. If the timescale on which $A_{t}$ evolves is much longer than the typical microscopic timescale, $\eta_{t}$ can be interpreted as a noise.

Several methods are already available to infer memory kernels from simulation data, e.g. Fourier-Laplace analysis Schnurr et al. (1997); Mokshin et al. (2005); Shin et al. (2010); Córdoba et al. (2012), projection operator analysis Darve et al. (2009); Carof et al. (2014); Li et al. (2017), parametrization techniques Lei et al. (2016); Ma et al. (2016), and iterative numerical inference schemes Jung et al. (2017). However, these methods are applicable only to stationary processes where the memory kernel $K(t^{\prime},t)$ effectively depends only on the difference $t^{\prime}-t$. We propose here a method that applies in both the stationary and the non-stationary case. The basic idea is to use a measured (or simulated) auto-correlation function $C(t^{\prime},t)$ to construct the memory kernel by means of an iterative procedure:

We formally integrate eqn. (2) into

[TABLE]

where we have defined

[TABLE]

such that $K(t^{\prime},t)=\partial_{t}\mathcal{J}(t^{\prime},t)$. Taking the derivative of eq. (3) with respect to $t^{\prime}$ and rearranging terms allows to write

[TABLE]

where we have defined

[TABLE]

(Note that the first time-derivative of the auto-correlation function $\partial_{t^{\prime}}C(t^{\prime},t)$ is all that is required for the following steps. We do not need to take any further time-derivatives of the input data.)

Next we substitute $\mathcal{J}(t,\tau)$ on the right-hand side of eq. (15) by eq. (15) itself.

[TABLE]

where we have defined $\mathcal{S}_{1}(t^{\prime},t)=\int_{t^{\prime}}^{t}d\tau\mathcal{S}_{0}(t^{\prime},\tau)\mathcal{S}_{0}(\tau,t)$. This construction is similar to a typical diagrammatic expansion in liquid state theory, e.g. in the context of the Ornstein-Zernike equation Hansen and McDonald (1990). $\mathcal{S}_{1}(t^{\prime},t)$ takes into account correlations at intermediate times $\tau$ similar to the way correlations between two particles in a liquid would be mediated by other particles.

We iterate the recursive substitution of $\mathcal{J}(\tau,t)$ below the integral and finally obtain

[TABLE]

where $\mathcal{S}(t^{\prime},t):=\sum_{n=0}^{\infty}\mathcal{S}_{n}(t^{\prime},t)$, and the functions $\mathcal{S}_{n}$ are defined recursively via the identity

[TABLE]

which is valid for any $(n,m)\in\mathbb{N}^{2}$. Note that in the limit $\omega=0$ we have $\mathcal{J}(t^{\prime},t)\stackrel{{\scriptstyle\omega=0}}{{\to}}-\mathcal{S}(t^{\prime},t)$. In general, higher orders in the expansion have impact on longer times, and the number of terms needed for the sum $\mathcal{S}$ to converge at a certain time yields information about the strength and the time-extend of the memory effects. Details about this statement as well as about the Markovian limit are given in appendix.

As all terms on the right-hand side of eqn. (21) are expressed in terms of the auto-correlation function $C(t^{\prime},t)$, we propose the following numerical scheme to compute $K(t^{\prime},t)$:

Carry out a set of simulations (or experiments) and measure the observable $A_{t}$ for each trajectory. 111The number of trajectories must be chosen such that the initial phase-space distribution $\rho_{0}(\mathbf{\Gamma})$ is probed with sufficient precision. 2. 2.

Compute the two-time auto-correlation function $C(t^{\prime},t)=N_{\text{traj}}^{-1}\sum_{i\in\text{traj}}A_{t^{\prime}}^{*(i)}A_{t}^{(i)}$. 3. 3.

Compute $\mathcal{S}_{0}$ and $j_{0}$ using eqns. (16) and (17). 4. 4.

Compute $\mathcal{S}_{n+1}$ recursively using eqn. (10) with $m=0$. Stop whenever for any times $t^{\prime}$ and $t$, $\mathcal{S}_{n}(t^{\prime},t)\ll\sum_{j=0}^{n}\mathcal{S}_{j}(t^{\prime},t)$. 5. 5.

Compute $\mathcal{J}$ using eqn. (21) and finally $K(t^{\prime},t)=\partial_{t}\mathcal{J}(t^{\prime},t)$.

This method is general enough to be applied to any dynamical process and any observable, even far from equilibrium, because it relies only on the structure of equs. (1) and (2). Note that once the memory kernel has been constructed, the corresponding nsGLE can be solved to predict the dynamics of the process. One has thus obtained an effective coarse-grained description in terms of one coordinate rather than $N$.

Side Remark - If one intends to study the fluctuations of $A_{t}$ independently from the evolution of the average trajectory, eqn. (1) may not be convenient. Instead one can use a modified, explictly time-dependent phase-space observable $\tilde{A}(\mathbf{\Gamma},t)$

[TABLE]

where $\mu_{A}(t)=\left\langle A(t)\right\rangle$ is the time-dependent average of $A$ and $\sigma_{A}(t)^{2}=\left\langle\left|A(t)\right|^{2}\right\rangle-\left\langle A(t)\right\rangle\left\langle A^{*}(t)\right\rangle$ is its time-dependent variance. The obervable $\tilde{A}$ measure deviations from the average trajectory, normalized by the variance of the process, it is thus a unitless number indicating whether a particular trajectory is delayed or advanced compared to the average one. We can define the corresponding auto-correlation function as $\tilde{C}(t^{\prime},t)=\left\langle\tilde{A}^{*}(t^{\prime})\tilde{A}(t)\right\rangle$ (the notation $\tilde{\cdot}$ from this point on refers to any function related to $\tilde{A}$), which characterizes how fluctuations decorrelate irrespective of the amplitude of the variable $A$. The structure of both equ. (1) and equ. (2) remains valid for $\tilde{A}$ and $\tilde{C}$. We thus introduce a new memory kernel $\tilde{K}$, a fluctuating force $\tilde{\eta}_{t}$ with $\left\langle\tilde{\eta}(t)\right\rangle=0$, and a drift function $\tilde{\omega}$ which turns out to vanish because $\tilde{C}(t,t)=1$ is constant, yielding

[TABLE]

This is then used to rewrite the equation of motion for $A_{t}$. Using $\sigma_{A}(t)d\tilde{A}_{t}/dt=d\Delta A_{t}/dt-\Delta A_{t}\dot{\sigma_{A}}(t)/\sigma_{A}(t)$, where $\Delta A_{t}=A_{t}-\mu_{A}(t)$, we can rewrite eqn. (12)

[TABLE]

with $\Omega(t)=d\ln\left(\sigma_{A}(t)\right)/dt$, $\mathcal{K}(t^{\prime},t)=\frac{\sigma_{A}(t)}{\sigma_{A}(t^{\prime})}\tilde{K}(t^{\prime},t)$ and $f_{t^{\prime}}(t^{\prime},t)=\sigma_{A}(t)\tilde{\eta}_{t^{\prime}}(t^{\prime},t)$. Hence there is a direct correspondence between the equation of motion for $\tilde{A}$ and the one for $A$. Note that $\omega\neq\Omega$ and $K\neq\mathcal{K}$ such that eqns. (1) and (13) are not identical, but their mathematical structure is the same. When studying memory effects and deriving an equation of motion for an observable, one can chose between these equivalent descriptions. In both cases, the method to numerically infer memory kernels is valid.

Lennard-Jones crystallization - As a test of the numerical procedure we analyzed the dynamics of crystallization. We carried out molecular dynamics simulations of $N=32,000$ particles interacting via a 6-12 Lennard-Jones potential. We used a cubic box with periodic boundary conditions and we ran the dynamics in the NVT ensemble, using a Nosé-Hoover thermostat to control the temperature. We first equilibrated the liquid phase at density $\rho=1$ and temperature $T=2$ (in Lennard-Jones reduced units), and then instantaneously quenched the temperature to $T=0.75$, for which the equilibrium phase is known to be a crystal. We then let the system evolve freely (i.e. we did not use any biasing scheme to speed up sampling). In order to monitor the formation and growth of crystallites, we used the method introduced by ten Wolde et. al Rein ten Wolde et al. (1996) based on orientational bond order parameters Steinhardt et al. (1983). A crystalline cluster is defined as a set of neighbouring crystal-like particles. The observable that we used as a reaction coordinate is the number $N_{c}$ of particles in the largest crystalline cluster. We observed the evolution of $N_{c}(t)$, see fig. 1, upper panel, where we also show a typical snapshot of a crystallite. A total of 4019 trajectories were used for the analysis.

We computed for each trajectory the modified variable $\tilde{N}_{c}$ as defined in eqn. (11) (see fig. 1, lower panel), and its auto-correlation function $\tilde{C}(t^{\prime},t)$. We then applied the numerical method presented above. As we used $\tilde{C}(t^{\prime},t)$, the drift term $\tilde{\omega}$ vanishes, which implies $j_{0}=-\mathcal{S}_{0}$ and $\mathcal{J}=-\mathcal{S}$. We computed $S_{0}(t^{\prime},t)=\partial_{t^{\prime}}\tilde{C}(t^{\prime},t)$, and we then iteratively applied the recursion relation eq. (10) with $m=0$, i.e. $\int_{t^{\prime}}^{t}d\tau\mathcal{S}_{0}(t^{\prime},\tau)\mathcal{S}_{n}(\tau,t)=\mathcal{S}_{n+1}(t^{\prime},t)$. We show in fig. (2) the functions $\mathcal{S}_{n}$, as well as the function $\mathcal{J}=-\sum_{n}S_{n}$. If the nucleation process were Markovian, $\mathcal{J}$ would be a step function of $t$, constant for $t<t^{\prime}$ and $t>t^{\prime}$ and discontinuous at $t=t^{\prime}$. Although we observe this discontinuity, the behaviour around $t\neq t^{\prime}$ is far from being constant. In addition, we note that $\tilde{N}_{c}(t)$ vanishes for each trajectory in the limit $t\to 0$, which implies that $\mathcal{S}_{0}(t^{\prime},t)\sim 0$ for $t\to 0$, and thus $J(t^{\prime},t)=-\mathcal{S}(t^{\prime},t)\sim 0$ for $t\to 0$. This is observed numerically, and suggests a non-trivial behaviour for the memory kernel $K$ at short times. The algorithm needed 15 iterations in $\mathcal{S}$ to properly converge. The functions $S_{n}(t^{\prime},t)$ probe the behaviour of memory effects for different times. Each new order $n$ presents extrema at increasing values of $|t^{\prime}-t|$. The summation of all the orders results in the smooth function $\mathcal{S}=-\mathcal{J}$ (see fig. 2).

We checked the validity of the function $\mathcal{J}$ obtained in this way by using it as an input to compute the right-hand side of eqn. (3), and we compared the result to the left-hand side, i.e. $\tilde{C}(t^{\prime},t)$ itself. As shown in figure 3, the overlap is very good 222The discrepancies at very short times are due to numerical errors in the computation of the derivative of $\tilde{C}(t^{\prime},t)$. Due to the initial vanishing width of the crystllite size distribution, the number of trajectories needed to eliminate statistical fluctuations at short times is very large. Here our aim is to demonstrate the validity of our method and not capture the details of the nucleation process, therefore we did not run additional simulations to resolve the early stage dynamics.. This test confirms that the method presented is able to reconstruct the dynamics of a reaction coordinate, and that it can be used to develop numerical coarse-graining schemes for dynamics out of thermal equilibrium.

Once this check on $\mathcal{J}(t^{\prime},t)$ was performed, we computed its derivative with respect to $t$, yielding $K(t^{\prime},t)$ which is also shown in fig. (2). In addition we show vertical lines indicating the mean first passage times for various values of $N_{c}$. A committor analysis of the trajectories showed that the critical crystallite contains ca. 80 particles. Hence we observe memory both, in the nucleation and in the growth regime. The average time needed for a particle to be attached to a cluster surface area of $\sigma_{LJ}^{2}$ was $\tau_{att}\simeq 5$ (in LJ units), which is about one order of magnitude smaller that the time extent of the memory kernel. We thus conclude that memory effects are not negligible in the Lennard-Jones crystallization process.

Conclusion - We have introduced a numerical method to construct memory kernels for any process for which the non stationary Generalized Langevin Equation is relevant, regardless of how far from equilibrium the system is. This procedure requires little computational effort and no modeling assumption since its only input is the two-time auto-correlation function of the observable under study. The method can also be applied to a modified version of the nsGLE, where the memory kernel and the fluctuating force contribute only to the fluctuations about the average of the observable of interest. We investigated the crystallization process as a proof of principle, and we have shown that the procedure allows to reconstruct the dynamics of the problem. In this particular example, we demonsrate that memory effects play a significant role for nucleation and growth dynamics.

Acknowledgement - We thank G. Stock for valuable feedback on the manuscript.

Supplemental material

Quantifying Non-Markovianity

To reconstruct the full memory kernel $K(t^{\prime},t)$ is an interesting way of assessing memory effects, however we can come up with a simpler measure of non-Markovianity. The Markovian limit is reached if $K(t^{\prime},t)=-\gamma(t)\delta(t^{\prime}-t)$, yielding $C(t^{\prime},t)=C(t^{\prime},t^{\prime})e^{-\int_{t^{\prime}}^{t}d\tau\gamma(\tau)}$ for any times $t^{\prime}<t$. Therefore, for any time $s\in[t^{\prime},t]$ a Markovian process fulfills $C(t^{\prime},t)=C(t^{\prime},s)C(s,t)$. We test the validity of this decomposition by defining a unitless quantity $\epsilon(s)$:

[TABLE]

where $T$ is the total time of the process. For a given value of $s$, $\epsilon(s)$ measures the deviation from the exponential decomposition at any point in the 2D-domain $\left[0,s\right]\times\left[s,T\right]$ and averages over it. By definition, $\epsilon(s)$ vanishes in the Markovian case. The dependence on $s$ allows to quantify how non-Markovian behaviour evolves throughout the process. In fig. 4 we show $\epsilon(s)$ for the Lennard-Jones crystallization process, using $\tilde{C}(t^{\prime},t)$ as input.

Details of the derivation for $\mathcal{J}(t^{\prime},t)$

We recall the starting point of our derivation for $\mathcal{J}$, that is

[TABLE]

where we have defined

[TABLE]

By substituting $\mathcal{J}(t,\tau)$ on the right-hand side of eq. (15) by eq. (15) itself, we obtain

[TABLE]

where we have defined $\mathcal{S}_{1}(t^{\prime},t)=\int_{t^{\prime}}^{t}d\tau\mathcal{S}_{0}(t^{\prime},\tau)\mathcal{S}_{0}(\tau,t)$. We can continue and insert eqn. (15) into the right-hand side of eqn. (Details of the derivation for $\mathcal{J}(t^{\prime},t)$), and find

[TABLE]

where we have defined $\mathcal{S}_{2}(t^{\prime},t)=\int_{t^{\prime}}^{t}d\tau\mathcal{S}_{1}(t^{\prime},\tau)\mathcal{S}_{0}(\tau,t)$. By defining more generally

[TABLE]

we can iterate the recursive substitution of $\mathcal{J}(\tau,t)$ below the integral and finally obtain

[TABLE]

Markovian limit

In the Markovian limit, we have $K(t^{\prime},t)=\gamma(t^{\prime})\delta(t^{\prime}-t)$, with $\gamma(t)<0$, yielding

[TABLE]

for $t^{\prime}<t$. Let us compute the first orders of $\mathcal{S}_{n}(t^{\prime},t)$. We first have

[TABLE]

We infer

[TABLE]

The nested integral on the right-hand side, together with the $(-1)^{n}$ factor, consist in the $n$-th order term in the Taylor expansion of the the exponential function $\exp\left(-\int_{t^{\prime}}^{t}d\tau\gamma(\tau)\right)$ :

[TABLE]

Again, this is valid for $t^{\prime}<t$. In the case $t^{\prime}>t$, the same calculation would yield $\mathcal{S}(t^{\prime},t)=\gamma(t^{\prime})$. In other words, for any $t$ and $t^{\prime}$, we would have $\mathcal{S}(t^{\prime},t)=-\mathcal{J}(t^{\prime},t)=\gamma(t^{\prime})\left[2\Theta(t^{\prime}-t)-1\right]$, where $\Theta$ is the Heaviside function. The derivative with respect to $t$ would then yield $K(t^{\prime},t)=\gamma(t^{\prime})\delta(t^{\prime}-t)$, which is consistent with what we started from.

An interesting property of the functions $\mathcal{S}_{n}(t^{\prime},t)$ is the fact that their derivatives with respect to $t$ can be written as

[TABLE]

In other words, we have $\partial_{t}\mathcal{S}_{n}(t^{\prime},t)=0\Leftrightarrow\mathcal{S}_{n}(t^{\prime},t)=\mathcal{S}_{n-1}(t^{\prime},t)$: the extremum of the function $\mathcal{S}_{n}(t^{\prime},t)$ is located where it crosses the previous order. This is a graphical check that can easily be performed to test the Markovian assumption.

Stationary case

In the stationary limit, $\gamma(t)$ becomes a constant. We then obtain, for $t^{\prime}<t$

[TABLE]

At constant $t^{\prime}$, $\mathcal{S}_{n}(t^{\prime},t)$ presents an extremum at $t^{*}_{n}=t^{\prime}-n/\gamma$, and $\mathcal{S}_{n}(t^{\prime},t_{n}^{*})=-\gamma n^{n}e^{n}/n!\stackrel{{\scriptstyle n\to\infty}}{{\propto}}n^{-1/2}$. In a stationary process, one can hence test the Markovian assumption by checking if the extrema of the functions $\mathcal{S}_{n}(t^{\prime},t)$ are located with an equal spacing and if their amplitude decay as $1/\sqrt{n}$.