1 Introduction

Transport coefficients are archetypal examples of off-equilibrium properties, which describe in fact entropy production and the approach to equilibrium in extended systems, thus giving a quantitative meaning and a conceptual framework to the very notion of irreversibility and the arrow of time. While non-equilibrium statistical mechanics is still a very active and largely unsettled field of research, the relaxation of small off-equilibrium fluctuations and the response of systems to small perturbations have been given a rigorous theoretical foundation back in the fifties by the Green–Kubo (GK) theory of linear response [1,2,3,4]. Among other achievements, this theory provides a rigorous and elegant way to cast the computation of transport coefficients into the evaluation of equilibrium time correlation functions of suitably defined fluxes, thus making it accessible to equilibrium molecular dynamics (EMD) simulations. This feat notwithstanding, the several conceptual subtleties underlying the linear-response theory of transport coefficients are often dodged or disguised with a clumsy notation that makes it difficult to fully appreciate their scope and impact on the design and use of computer simulation methodologies. Also due to this predicament, a number of misconceptions have affected the otherwise mature and fecund field of computer simulation of transport in condensed matter, thus unduly limiting its scope. Foremost among these misconceptions is that the intrinsic indeterminacy of any local representations of an extensive quantity, such as the density or atomic break-up of the energy (or charge and mass, for that matter), would undermine the uniqueness of the transport coefficients that are derived from them. While similar reservations should apply to classical and quantum ab initio simulations alike, they have in fact impacted mainly the latter, to the extent that until recently it had been believed that the GK theory of heat transport could not be combined with quantum simulation methods based on electronic-structure theory. This quandary has been recently overcome, mainly thanks to the introduction of a so-called gauge invariance principle of transport coefficients [5, 6], which basically states that, under well-defined conditions, the value of a transport coefficient is largely independent of the detailed form of the conserved densities and fluxes from which they are derived through the GK formulae. In full generality, gauge invariance implies that the value of a transport coefficients is unchanged if the flux from which it is calculated is modified by adding to it a vector process whose power spectrum vanishes at \(\omega =0\) (such a process is conveniently dubbed non-diffusive). Further difficulties may arise in the case of multi-component systems, where the interaction among different fluxes make the transport of one conserved quantity (such as, e.g., energy), depend on the dynamics of the other hydrodynamical variables (such as, e.g., the number of different molecular species), thus muddling the definition of heat conductivity in these systems. This difficulty is solved by defining, e.g., the thermal conductivity as the ratio between the energy flux and the temperature gradient, when all the other conserved fluxes vanish. In this case, a further invariance principle, dubbed convective invariance [7], states that the thermal conductivity results to be unchanged if the definition of the energy flux is altered by adding to it an arbitrary linear combination of the mass fluxes of the molecular species constituting the system. While the value of the transport coefficients enjoys the invariance properties mentioned above, the statistical properties of the flux time series from which they are derived do depend on the microscopic representation of the conserved densities and current densities, thus substantially affecting the statistical error of the transport coefficient being computed. This dependence opens the way to optimizing this representation, so as to minimize the resulting statistical errors. This freedom can be exploited in conjunction with the recently introduced cepstral analysis of the current spectra [8] to substantially reduce the length of the EMD simulations needed to evaluate a transport coefficient to a given target accuracy.

In this paper, we review the concepts of gauge and convective invariance in the classical theory of transport in condensed matter, with emphasis on their application to the computation of transport coefficients from the GK theory of linear response and EMD, as well as some concepts and tools for the spectral analysis of the current time series, which can be used in conjunction with these invariance principles to substantially reduce the statistical errors affecting the EMD estimate of various conductivities. In Sect. 2, we briefly lay down the GK linear-response theory of transport, aiming at establishing some general terminology and notations. In Sect. 3, we introduce the concept of gauge invariance, whereas convective invariance is discussed in Sect. 4. In Sect. 5, we discuss a newly introduced spectral method (dubbed cepstral analysis) to evaluate and systematically reduce the statistical error affecting the estimate of transport coefficients from EMD. In Sect. 6, we specialize our discussion to ab initio heat transport and report a general expression for the microscopic heat flux suitable to density functional theory. In Sect. 7, we show how gauge invariance can be combined with concepts from topology to reveal some unexpected features of charge transport in ionic conductors. Section 8 finally contains our conclusions.

2 Theory

Transport theory is essentially a dynamical theory of hydrodynamical variables, i.e. of the long-wavelength components of the densities of conserved, extensive, variables, which in short we refer to as conserved charges and densities. In a simple fluid, for instance, the conserved charges are the energy, the three components of total momentum, and the total number of molecules of each chemical species. The corresponding transport coefficients are the thermal conductivity, the viscosity, and the various diffusivities. Let \(\hat{Q}\) be one such conserved charge, and \(\hat{q}({\varvec{r}})\) the corresponding density: \(\hat{Q} = \int \hat{q}({\varvec{r}}) \mathrm{d}{\varvec{r}}\). Here and in the following we adopt the convention that a hat, as in \(\hat{A}\), indicates an implicit dependence on the system’s phase-space variables, \(\Gamma \): \(\hat{A} = A(\Gamma )\), where \(\Gamma =\{x,p\}\) is the set of all the atomic coordinates and momenta. When neither a hat nor the phase-space argument are present, A indicates the expectation of \(\hat{A}\) over some suitably defined phase-space distribution, \(\rho (\Gamma )\) (most frequently the canonical or micro-canonical equilibrium distributions): \(A = \langle \hat{A} \rangle \doteq \int A(\Gamma ) \rho (\Gamma ) \mathrm{d}\Gamma \). When it will be necessary to distinguish between the expectation of a quantity, A, and its phase-space representation, \(A(\Gamma )\), the latter will be referred to as a phase-space sample of the former. Sometimes, we will need to indicate explicitly the implicit time dependence of a phase-space variable (no pun intended). When this is done, we will mean the phase-space variable \(\hat{A}(t)=A(\Gamma ,t)\doteq A(\Gamma _t)\), as a function of time and of the initial condition, \(\Gamma =\Gamma _0\), of a phase-space trajectory, \(\Gamma _t\), as determined by Hamilton’s equations of motion.

Locality implies that for any conserved density, \(q({\varvec{r}})\), a (conserved) current density, \({\varvec{\jmath }}({\varvec{r}})\), can be defined, such that the two of them satisfy the continuity equation, \( \dot{q}({\varvec{r}},t) +~\nabla ~\cdot ~{\varvec{\jmath }}({\varvec{r}},t) = 0\), where \(\dot{q}\) indicates the time derivative of q. A current density, \({\varvec{\jmath }}({\varvec{r}},t)\), satisfying the continuity equation with a conserved density, \(q({\varvec{r}},t)\), will be said to be conjugate to q. Strictly speaking, the continuity equation holds for the expectations of conserved densities and current densities, as well as for their phase-space samples, \(\hat{q}({\varvec{r}},t) = q({\varvec{r}}, \Gamma _t)\) \(\hat{{\varvec{\jmath }}} ({\varvec{r}},t) = {\varvec{\jmath }}({\varvec{r}},\Gamma _t)\). For the sake of unburdening the notation as much as possible, we will sometimes overlook the distinction between phase-space samples and their expectations.

The Fourier transform of the continuity equation reads:

$$\begin{aligned} \dot{\tilde{q}}({\varvec{k}},t)+i{\varvec{k}}\cdot \tilde{{\varvec{\jmath }}}({\varvec{k}},t)=0, \end{aligned}$$
(1)

where \(\tilde{q}({\varvec{k}})=\int q({\varvec{r}})\mathrm {e}^{i{\varvec{k}}\cdot {\varvec{r}}} \mathrm{d}{\varvec{r}}\) is the Fourier transform of q, and similarly for \({\varvec{\jmath }}({\varvec{r}})\) and any other function of \({\varvec{r}}\). Equation (1) shows that the smaller the wave-vector, \(|{\varvec{k}}|\), i.e. the longer the wavelength, the slower the dynamics of conserved densities and fluxes. This means that, at sufficiently long wavelength, conserved densities and current densities are adiabatically decoupled from the (zillions of) other atomically fast degrees of freedom. Also, translational invariance implies that conserved densities at different wavevectors do not interact with each other. As a consequence, the dynamics of hydrodynamic variables is determined by a handful of equations that couple them with each other at fixed wavevector. When the intensive variables conjugate to the conserved quantities depend on position sufficiently slowly, the system can be thought of as locally in thermal equilibrium and conserved currents can then be connected to the thermodynamical affinities (i.e. to the gradients of the intensive variables) through the so-called constitutive equations [9]. By combining the constitutive equations with the continuity equations of all the conserved densities and currents, the Navier–Stokes equation of classical hydrodynamics can be derived [9].

Let us consider a system described by a Hamiltonian \(\hat{H}^\circ \) and subject to a time-dependent external perturbation, \(\hat{V}(t)=\sum _i\int v_i({\varvec{r}},t) \hat{q}_i({\varvec{r}})\mathrm{d}{\varvec{r}}\), where \(\{\hat{q}_i\}\) is a set of conserved densities. The perturbed Hamiltonian reads:

$$\begin{aligned} \hat{H} = \hat{H}^\circ + \sum _i\int v_i({\varvec{r}},t) \hat{q}_i({\varvec{r}})\mathrm{d}{\varvec{r}}. \end{aligned}$$
(2)

Here the \(v_i\)s are to be treated as strengths of the perturbation in linear-response theory. As such, they are not phase-space functions, and they dependence on time only explicitly. To first order in the vs, the expected value of the conserved current densities conjugate to the qs, \(\{\hat{{\varvec{\jmath }}}_i\}\), can be obtained from the GK theory of linear response [3, 4] as [10]:

$$\begin{aligned} \jmath _{i\alpha }({\varvec{r}},t) = -\frac{1}{k_\mathrm{B}T} \\ \quad \times \sum _j\int \mathrm{d}{\varvec{r}}'\int _{-\infty }^t \mathrm{d}t' \Bigl \langle \,{\widehat{\!\!{j}}}_{i\alpha }({\varvec{r}},t) \dot{\hat{q}}_j({\varvec{r}}',t') \Bigr \rangle v_j({\varvec{r}}',t'), \end{aligned}$$
(3)

where \(k_\mathrm{B}\) is the Boltzmann’s constant, T the system’s temperature, and the correlation function, \(\langle \cdot \rangle \), is defined for a pair of general phase-space variables, \(\hat{X}\) and \(\hat{Y}\), as the equilibrium expectation over the initial conditions of a molecular trajectory, \(\Gamma _t\), of the time-lagged product of the values of the variables:

$$\begin{aligned} \begin{aligned} \langle \hat{X}(t) \hat{Y}(t') \rangle&= \langle \hat{X}(t-t') \hat{Y}(0) \rangle \\&= \int X(\Gamma _{t-t'}) Y(\Gamma _0) \rho (\Gamma _0) \mathrm{d}\Gamma _0. \end{aligned} \end{aligned}$$
(4)

The dependence of the correlation function in Eq. (4) on the time difference is due to time-translation invariance ensuing from the equilibrium condition. By the same token, space-translation invariance makes the correlation function in Eq. (3) only depend on \({\varvec{r}}-{\varvec{r}}'\), turning the integral in \(\mathrm{d}{\varvec{r}}'\) into a convolution. Using the continuity equation to replace the time derivative of the density with the divergence of the conjugate current density, Eq. (3) can be cast into a linear relation between the Fourier transforms of the longitudinal component of the current density and the forces acting on the system:

$$\begin{aligned} \tilde{\jmath }_{\parallel i}({\varvec{k}},t) =-\sum _j \int _{-\infty }^t \chi _{\parallel ij}({\varvec{k}},t-t') \tilde{f}_{\parallel j}({\varvec{k}},t') \mathrm{d}t', \end{aligned}$$
(5)

where \(\tilde{g}_\parallel ({\varvec{k}})=\frac{1}{k}{\varvec{k}}\cdot \tilde{{\varvec{g}}}({\varvec{k}})\) indicates the Fourier transform of the longitudinal component of a generic vector field \({\varvec{g}}({\varvec{r}})\), \(\chi _{\parallel ij}({\varvec{k}},t)=\frac{1}{k_BT}\bigl \langle \tilde{\hat{\jmath }}_{\parallel i}({\varvec{k}},t) \tilde{\hat{\jmath }}_{\parallel j}(-{\varvec{k}},0) \bigr \rangle \) is the longitudinal susceptibility of the current densities, and \({\varvec{f}}_i({\varvec{r}})=-\nabla v_i({\varvec{r}})\) is the force field associated with the \(v_i\) perturbing potential. For the sake of streamlining the notation and without much loss of generality, we will restrict ourselves to longitudinal perturbations and response currents, and drop the “\(\parallel \)” suffix from currents and forces. If the external perturbation is independent of time, in the long-wavelength limit Eq. (5) results in the Onsager relation between particle fluxes and applied forces [11, 12]:

$$\begin{aligned} J_n = \sum _{nm} \Lambda _{nm}F_m, \end{aligned}$$
(6)

where the index \(n=(i,\alpha )\) denotes for short the combination of the indices i, for the conserved charge, and \(\alpha \), for the Cartesian component; the fluxes \(J_n=\frac{1}{\Omega } \tilde{\jmath }_n(0)\) and forces \(F_n=\frac{1}{\Omega } \tilde{f}_n(0)\) are the macroscopic averages of the current densities and force fields, respectively; \(\Omega \) is the system’s volume, and

$$\begin{aligned} \Lambda _{nm}=\frac{\Omega }{k_BT}\int _0^\infty \Bigl \langle \hat{J}_n(t) \hat{J}_m(0) \Bigr \rangle dt \end{aligned}$$
(7)

is the matrix of Onsager’s transport coefficients [11, 12]. In the case of a charged fluid, for instance, the steady state charge flux, \({\varvec{J}}\), induced by a stationary electric field, \({\varvec{E}}\), is given by \({\varvec{J}}=\sigma {\varvec{E}}\), where the static electrical conductivity is \(\sigma =\frac{\Omega }{k_BT}\int _0^\infty \bigl \langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}} (0) \bigr \rangle \mathrm{d}t\).

We can reformulate the GK expression of Onsager’s coefficients in another equivalent representation, the so-called Helfand–Einstein (HE) formula, which will be expedient in the following and is also better behaved statistically, based on the identity:

$$\begin{aligned} \int _0^{\mathcal {T}}\mathrm{d}t \int _0^{\mathcal {T}}\mathrm{d}t' \, f(t' -t) = 2{\mathcal {T}}\int _0^{\mathcal {T}}\mathrm{d}t \, \left( 1 - \frac{t}{{\mathcal {T}}} \right) f(t){,} \end{aligned}$$
(8)

valid for any even function, \(f(-t)=f(t).\) Let \(\,{\widehat{\!\!{\varvec{J}}}}(t)\) be a stationary stochastic process representing a conserved flux, so that \(f(t,t') = \langle \,{\widehat{\!\!{\varvec{J}}}}(t)~\cdot ~ \,{\widehat{\!\!{\varvec{J}}}}(t')\rangle \) only depends upon \(|t-t'|\). By applying the identity above, we obtain

$$\begin{aligned} \int _0^\infty \left\langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}}(0) \right\rangle \mathrm{d}t = \lim _{{\mathcal {T}}\rightarrow \infty } \frac{1}{2{\mathcal {T}}} \left\langle \left| \int _0^{\mathcal {T}}\,{\widehat{\!\!{\varvec{J}}}}(t) \mathrm{d}t \right| ^2 \right\rangle . \end{aligned}$$
(9)

This is called the HE formula, which gives a transport coefficient of some conserved charge, as the ratio of the mean-square dipole, \({\varvec{D}}({\mathcal {T}})=\int _0^{\mathcal {T}}{\varvec{J}}(t)\mathrm{d}t\), displaced by the conserved flux, \({\varvec{J}}\), in a time t and time itself. This argument was first exploited by Einstein in his celebrated paper on Brownian motion [13] to establish the relation between diffusivity and velocity auto-correlation functions, and later extended by Helfand to general transport phenomena [14]. A comparison between the numerical performance of the GK and the HE formulas is displayed in Fig. 1 in the case of charge transport in a molten salt. The better stability of the HE integral is evident: not only does the HE integral converge faster than the GK one, but the variance on the first, even though growing linearly with time, is much smaller than the one on the second.

Fig. 1
figure 1

Comparison between the Green–Kubo and the Helfand–Einstein integrals, see Eq. (9), of the autocorrelation function of the charge flux for an ab initio MD simulation of molten KCl. The shaded area represents the confidence interval estimated via standard block analysis

Full size image

To understand the better statistical behaviour of the HE representation of transport coefficients with respect to the GK one, we leverage the relation between conductivities and the spectral properties of the conserved fluxes, which will also be instrumental in our subsequent considerations on data analysis in Sect. 5. Let us define

$$\begin{aligned} \begin{aligned} \lambda&= \lim _{{\mathcal {T}}\rightarrow \infty } \lambda _\mathrm{GK}({\mathcal {T}}) \\&= \lim _{{\mathcal {T}}\rightarrow \infty } \lambda _\mathrm{HE}({\mathcal {T}}) \\&= \frac{1}{2}S(0), \end{aligned} \end{aligned}$$
(10)

where

$$\begin{aligned} \begin{aligned} \lambda _\mathrm{GK}({\mathcal {T}})&=\int _0^{\mathcal {T}}\left\langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}}(0) \right\rangle \mathrm{d}t, \\ \lambda _\mathrm{HE}({\mathcal {T}})&= \int _0^{\mathcal {T}}\left\langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}}(0) \right\rangle \left( 1-\frac{t}{{\mathcal {T}}} \right) \mathrm{d}t, \\ S(\omega )&=\int _{-\infty }^\infty \left\langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}}(0) \right\rangle \mathrm {e}^{i\omega t}\mathrm{d}t. \end{aligned} \end{aligned}$$
(11)

One has

$$\begin{aligned} \begin{aligned} \lambda _\mathrm{GK}({\mathcal {T}})&= \frac{1}{2}\int _{-\infty }^\infty \Theta _{GK}^{\mathcal {T}}(t) \left\langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}}(0) \right\rangle \mathrm{d}t \\ \lambda _\mathrm{HE}({\mathcal {T}})&= \frac{1}{2}\int _{-\infty }^\infty \Theta _{HE}^{\mathcal {T}}(t) \left\langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}}(0) \right\rangle \mathrm{d}t, \end{aligned} \end{aligned}$$
(12)

where

$$\begin{aligned} \begin{aligned} \Theta _\mathrm{GK}^{\mathcal {T}}(t)&= {\left\{ \begin{array}{ll} 1 &{} \text {for} ~ |t|\le {\mathcal {T}}\\ 0 &{} \text {otherwise} \end{array}\right. }, \\ \Theta _\mathrm{HE}^{\mathcal {T}}(t)&= {\left\{ \begin{array}{ll} 1 - {\frac{\displaystyle |t|}{\displaystyle {\mathcal {T}}}}&{} \text {for} ~ |t|\le {\mathcal {T}}\\ 0 &{} \text {otherwise} \end{array}\right. } . \end{aligned} \end{aligned}$$
(13)

Using the Parseval–Plancherel identity [15], one gets

$$\begin{aligned} \lambda _X({\mathcal {T}}) = \frac{1}{4\pi } \int _{-\infty }^{\infty } \tilde{\Theta }_X^{\mathcal {T}}(\omega ) S(\omega ) \mathrm{d}\omega , \end{aligned}$$
(14)

where \(X=\mathrm{GK}\text { or }\mathrm{HE}\) and \(\tilde{\Theta }_X^{\mathcal {T}}(\omega ) = \int _{-\infty }^\infty \Theta _X^{\mathcal {T}}(t) \mathrm {e}^{i\omega t}\mathrm{d}t \). The two functions

$$\begin{aligned} \begin{aligned} \tilde{\Theta }^{\mathcal {T}}_\mathrm{GK}(\omega )&= 2 {\mathcal {T}}\mathrm {sinc}(\omega {\mathcal {T}}) \\ \tilde{\Theta }^{\mathcal {T}}_\mathrm{HE}(\omega )&= {\mathcal {T}}\mathrm {sinc}^2(\tfrac{\omega {\mathcal {T}}}{2}) , \end{aligned} \end{aligned}$$
(15)

where \(\mathrm {sinc}(x) \doteq \sin (x)/x\) is the cardinal sine function, are displayed in Fig. 2: as \({\mathcal {T}}\rightarrow \infty \), they both tend to \(2\pi \delta (\omega )\), where \(\delta \) indicates the Dirac delta function. One sees that the statistical accuracy of the HE estimate of the transport coefficients is higher than GK’s. Using the findings of Sect. 5, in Appendix A we demonstrate that, in the large-T limit, one has indeed the following relation between the statistical uncertainties: \(\Delta \lambda _\mathrm{HE} \approx \frac{1}{\sqrt{3}} \Delta \lambda _\mathrm{GK}\).

Fig. 2
figure 2

Fourier transforms of the cutoff functions entering the finite-time GK and HE expressions for the transport coefficients as integrals over the entire real axis, see Eq. (13)

Full size image

3 Gauge invariance

3.1 General theory

The gauge-invariance principle of transport coefficients is the condition by which transport coefficients are largely insensitive to the specific definition of the fluxes. In fact, from a microscopic standpoint, any two conserved densities, \(\hat{q}'({\varvec{r}},t)\) and \(\hat{q}({\varvec{r}},t)\), whose integrals over a volume \(\Omega \) differ by a quantity that scales as the volume boundary, should be considered equivalent in the large volume limit, \(\Omega \rightarrow \infty \). For instance, two equivalent densities may differ by the divergence of a (bounded) vector field \(\hat{{\varvec{b}}}({\varvec{r}},t)\):

$$\begin{aligned} \hat{q}'({\varvec{r}},t)=\hat{q}({\varvec{r}},t) - \nabla \cdot \hat{{\varvec{b}}}({\varvec{r}},t). \end{aligned}$$
(16)

In this sense, \(\hat{q}({\varvec{r}},t)\) and \(\hat{q}'({\varvec{r}},t)\) can be thought of as different gauges of the same scalar field. Since \(\hat{Q}=\int _\Omega \hat{q}({\varvec{r}},t) \mathrm{d}{\varvec{r}}\) is also conserved, for any given gauge of the conserved density, \(\hat{q}({\varvec{r}},t)\), a conserved current density can be defined, \(\hat{{\varvec{\jmath }}}({\varvec{r}},t)\), to satisfy the continuity equation, Eq. (1). By combining Eqs. (16) and (1) we see that conserved current densities, as well as macroscopic fluxes, transform under a gauge transformation as

$$\begin{aligned} \begin{aligned} \hat{{\varvec{\jmath }}}'({\varvec{r}},t)&= \hat{{\varvec{\jmath }}}({\varvec{r}},t) + \dot{\hat{{\varvec{b}}}}({\varvec{r}},t), \\ \,{\widehat{\!\!{\varvec{J}}}}'(t)&= \,{\widehat{\!\!{\varvec{J}}}}(t) + \dot{\hat{{\varvec{B}}}}(t), \end{aligned} \end{aligned}$$
(17)

where \(\hat{{\varvec{B}}}(t)=\frac{1}{\Omega } \int \hat{{\varvec{b}}}({\varvec{r}},t)\mathrm{d}{\varvec{r}}\), the time dependence being implicit via the phase-space point \(\Gamma _t\), e.g. \(\,{\widehat{\!\!{\varvec{J}}}}(t)={\varvec{J}}(\Gamma _t)\), and where the dot indicates the Poisson brackets with the unperturbed Hamiltonian. We conclude that the macroscopic energy fluxes in two different energy gauges differ by the total time derivative of a bounded phase-space vector function. Nonetheless, both definitions must lead to the same transport coefficient. To prove such a gauge-invariance principle, let us consider a generic transport process with conserved flux represented by the stochastic stationary process \(\,{\widehat{\!\!{\varvec{J}}}}(t)\), and define the (generalised, conserved) charge displacement per unit volume \( \hat{{\varvec{D}}}({\mathcal {T}})=\int _0^{{\mathcal {T}}} \,{\widehat{\!\!{\varvec{J}}}}(t) \, \mathrm{d}t\). According to the HE formulation we can express the transport coefficient of the process, \(\sigma \), as

$$\begin{aligned} \sigma&= c \lim _{{\mathcal {T}}\rightarrow \infty } \frac{\langle | \hat{{\varvec{D}}}({\mathcal {T}}) |^2 \rangle }{2{\mathcal {T}}}, \end{aligned}$$
(18)

where the factor c is specific to the transport process considered. The addition of a time-bounded term \(\hat{{\varvec{B}}}({\mathcal {T}})\) to \( \hat{{\varvec{D}}}({\mathcal {T}})\), resulting in the new displacement \( \hat{{\varvec{D}}}^\prime ({\mathcal {T}}) \equiv \hat{{\varvec{D}}}({\mathcal {T}}) + \hat{{\varvec{B}}}({\mathcal {T}})\) produces a transport coefficient \(\sigma '\) which coincides with \(\sigma \). In fact, by direct calculation, we have

$$\begin{aligned} \sigma '&\equiv c \lim _{{\mathcal {T}}\rightarrow \infty } \frac{\langle |\hat{{\varvec{D'}}}({\mathcal {T}}) |^2 \rangle }{2{\mathcal {T}}} \nonumber \\&=c \lim _{{\mathcal {T}}\rightarrow \infty } \frac{ \overbrace{\langle |\hat{{\varvec{D}}}({\mathcal {T}}) |^2 \rangle }^{\mathcal {O}({\mathcal {T}})} + \overbrace{2 \langle \hat{{\varvec{B}}}({\mathcal {T}}) \hat{{\varvec{D}}}({\mathcal {T}}) \rangle }^{\mathcal {O}({\mathcal {T}}^{1/2})} + \overbrace{\langle | \hat{{\varvec{B}}}({\mathcal {T}}) |^2 \rangle }^{\mathcal {O}({\mathcal {T}}^0)}}{2{\mathcal {T}}} \nonumber \\ {}&=c \lim _{{\mathcal {T}}\rightarrow \infty } \frac{\langle |\hat{{\varvec{D}}}({\mathcal {T}}) |^2 \rangle }{2{\mathcal {T}}} \equiv \sigma . \end{aligned}$$
(19)

3.2 Some considerations on boundary conditions

To conclude this section, we spend a few words to examine the role of boundary conditions (BC) in practical molecular dynamics simulations of a finite sample of the system. First of all, the simulation box must be larger than the relevant correlation/diffusion lengths of the system, for equilibrium properties to be independent of specific BC adopted in the simulation. In EMD simulations, periodic BC (PBC) are preferred since (1) they minimize size effects, and (2) the limit \(\lim _{{\mathcal {T}}\rightarrow \infty } \int _0^{\mathcal {T}}\langle \,{\widehat{\!\!{\varvec{J}}}}(t) \,{\widehat{\!\!{\varvec{J}}}}(0) \rangle \, \mathrm{d}t\) commutes with thermodynamic limit, where \(L,N\rightarrow \infty \) while the density \(L^3/N\) is kept fixed. This commutation no longer holds in open BC (OBC) at thermodynamic equilibrium, where the asymptotic time limit, \({\mathcal {T}}\rightarrow \infty \), must be taken only after the thermodynamic limit is performed.

Differently stated, PBC have to be preferred with respect to OBC, since they are the only ones which can sustain a steady-state flux [16]. Nonetheless, this poses some issues in the definitions of the fluxes, since the textbook definition relying on the first moment of the time-derivative of the (periodic) charge density \(\hat{q}({\varvec{r}},t)\),

$$\begin{aligned} \,{\widehat{\!\!{\varvec{J}}}}(t) = \frac{1}{\Omega } \int _{\Omega } \dot{\hat{q}}({\varvec{r}},t) \, {\varvec{r}} \mathrm{d}{\varvec{r}}, \end{aligned}$$
(20)

cannot be employed since \({\varvec{r}}\) is ill-defined in extended, PBC-closed systems. Strictly speaking, Eq. (20) is ill-defined in PBC for the very same reason why macroscopic polarisation is so in insulators [17]. The formal meaning of this equation is that it should be considered as the leading order of a Taylor series of the Fourier transform of the time derivative of the conserved density in powers of its argument: \(\dot{\hat{\tilde{q}}}({\varvec{k}},t) =-i {\varvec{k}} \cdot \,{\widehat{\!\!{\varvec{J}}}}(t)+\mathcal {O}(k^2)\). Computer simulations performed for systems of any finite size, L, give access to the Fourier components of conserved (current) densities at finite wave-vectors whose minimum magnitude is \(|{{\varvec{k}}}_\mathrm{min}|=\frac{2\pi }{L}\). Let \(c_L(t)= \bigl \langle \hat{\tilde{{\varvec{\jmath }}}}({{\varvec{k}}}_\mathrm{min},t) \cdot \hat{\tilde{{\varvec{\jmath }}}}(-{{\varvec{k}}}_\mathrm{min},0) \bigr \rangle _L \) be the (spatial Fourier transform of the) current-current correlation function evaluated at \( {\varvec{k}} = {{\varvec{k}}}_\mathrm{min}\), which can be easily evaluated from a MD simulation. Unfortunately, \(c_L(t)\) cannot be directly used to estimate, not even by extrapolation, the values of the transport coefficients, because the corresponding GK integral vanishes for any finite system size. In practice, the flux to be used to evaluate transport coefficients from the GK formula is obtained from Eq. (20) by formal manipulations that make the flux boundary insensitive. The time correlation function \(C_L(t) \doteq \langle \,{\widehat{\!\!{\varvec{J}}}}(t) \cdot \,{\widehat{\!\!{\varvec{J}}}}(0) \rangle _L\), computed for a system of size L, is well-defined as well, and has the property that, for any given t,

$$\begin{aligned} \lim _{L\rightarrow \infty } \bigr (C_L (t) - c_L(t) \bigr )=0. \end{aligned}$$
(21)

The tricky thing is the convergence of the limit in Eq. (21) is not uniform and, while \(\int _0^\infty c_L(t)\mathrm{d}t=0\), \(\int _0^\infty C_L(t)\mathrm{d}t\) is the non-vanishing GK integral yielding the transport coefficient we are after. Therefore, definitions for the fluxes suitable to PBC must be properly designed not only from the speculative standpoint, but also to run meaningful simulations. In what follows, we shall largely use these PBC-based definitions together with the gauge invariance principle to draw general conclusions on transport properties which do not depend on the system size, and hold in the thermodynamic limit.

4 Convective invariance

In a multicomponent system, the (relevant) conserved charges are the energy and the particles number (or, equivalently, the masses) of each atomic species. Since the total-mass flux (i.e. the total linear momentum) is itself a constant of motion, for a K-species system the number of independent conserved fluxes is equal to K (energy, plus \(K-1\) partial masses) [18]. Further constraints may reduce the number of relevant conserved fluxes. For instance, energy flux becomes the only relevant conserved flux in solids or in one-component molecular liquids, as long as the molecules do not dissociate. In true multicomponent systems (molten salts, solutions, etc.), neither the mass fluxes of the single atomic species are constant of motion nor their integral is a bound quantity. This fact must be taken into account in practical simulations of heat transport since the thermal conductivity relates, by definition, a gradient of temperature to the induced energy flux in the absence of convection, i.e. when the non-equilibrium average of the macroscopic mass flux vanishes.

To make things simpler but also more quantitative, let us consider a two-component system, like, e.g. a molten salt. The conserved fluxes are the energy flux \({\varvec{J}}_E\) and the mass flux of one of the species \({\varvec{J}}_M\). The following system of phenomenological equations holds:

$$\begin{aligned} {\left\{ \begin{array}{ll} {\varvec{J}}_E &{}= \Lambda ^{EE} \nabla \left( \frac{1}{T}\right) + \Lambda _{EM} \nabla \left( \frac{\Delta \mu }{T}\right) \\ {\varvec{J}}_M &{}= \Lambda _{ME} \nabla \left( \frac{1}{T}\right) + \Lambda _{MM} \nabla \left( \frac{\Delta \mu }{T}\right) , \end{array}\right. } \end{aligned}$$
(22)

where \(\Delta \mu \) is the difference between the chemical potentials of the two species [19], and where \({\varvec{J}}_{M}\equiv {\varvec{J}}_{M_1} = - {\varvec{J}}_{M_2}\), \({\varvec{J}}_{M_l}\) being the mass flux of species l, the last step following from the conservation of linear momentum. When we set the mass flux \({\varvec{J}}_M = 0\),

$$\begin{aligned} \nabla \left( \frac{\Delta \mu }{T}\right) = - \frac{\Lambda _{ME}}{\Lambda _{MM}} \nabla \left( \frac{1}{T}\right) \end{aligned}$$
(23)

which can be substituted in the first equation of the system to finally find

$$\begin{aligned} {\varvec{J}}_E&= -\kappa \nabla T \end{aligned}$$
(24)
$$\begin{aligned} \kappa&\equiv \frac{1}{T^2} \left[ \Lambda _{EE} - \frac{(\Lambda _{EM})^2}{\Lambda _{MM}} \right] , \end{aligned}$$
(25)

where we employed the symmetric property \(\Lambda _{EM}=\Lambda _{ME}\). In the light of GK theory, the thermal conductivity \(\kappa \) is obtained by removing from the GK integral of the energy flux a term which represents the contributions of convection/mass diffusion to heat flow [20].

It is straightforward to verify that a change in the definition of the microscopic energy flux by any multiple of the mass flux,

$$\begin{aligned} \,{\widehat{\!\!{\varvec{J}}}}_{E}^\prime = \,{\widehat{\!\!{\varvec{J}}}}_E + c \,{\widehat{\!\!{\varvec{J}}}}_M, \quad c \in \mathbb {R}, \end{aligned}$$
(26)

does not affect \(\kappa \), even if such change does change each of the Onsager coefficients in Eq. (25). We dub this peculiar property the “convective invariance” principle. This can be easily extended to more than two species, \(K>2\), thanks to standard linear algebra techniques. In such case, \(\kappa \) becomes

$$\begin{aligned} \kappa = \frac{1}{T^2} \left( \Lambda _{EE} - \sum _{l,m=1}^{K-1} \Lambda _{EM_l} (\mathrm {L}^{-1})^{l m} \Lambda _{EM_m} \right) , \end{aligned}$$
(27)

where \(\mathrm {L}=\{\Lambda ^{M_l M_m}\}\) is the square matrix of Onsager coefficients of the mass fluxes. The convective-invariance principle reads

$$\begin{aligned} \,{\widehat{\!\!{\varvec{J}}}}_{E}^\prime = \,{\widehat{\!\!{\varvec{J}}}}_E + \sum _{l=1}^{K-1} c_l \,{\widehat{\!\!{\varvec{J}}}}_{M_l}, \quad c_l \in \mathbb {R} \quad \Rightarrow \quad \kappa '=\kappa . \end{aligned}$$
(28)

Any linear combination of the mass fluxes can be added to the energy flux without affecting the thermal conductivity. This statement has a direct, crucial consequence concerning ab initio calculations: the heat conductivity cannot in fact depend on whether atomic cores contribute to the definition of the atomic energy, as they would in an all-electron calculation, or not, as they would when using pseudo-potentials. In the latter case, the energy of isolated atoms would depend on the specific form of pseudo-potential adopted, which is to a large extent arbitrary, while the heat conductivity in all cases should not. Thanks to convective invariance, shifting the zero of energy of each species by a quantity \(\delta E_l\) would result in a change \(\,{\widehat{\!\!{\varvec{J}}}}_{E}^\prime = \,{\widehat{\!\!{\varvec{J}}}}_E + \sum _{l=1}^{K-1} \frac{\delta E_l}{M_l} \,{\widehat{\!\!{\varvec{J}}}}_{M_l}\) which does not affect \(\kappa \), just like physical intuition would suggest. Finally, from a more practical way, convective invariance also avoids the calculation of partial enthalpies to dispose of the spurious self-energy effects, a rather tedious and cumbersome task [21,22,23].

5 Cepstral analysis

5.1 Wiener–Khinchin theorem

Cepstral analysis is a powerful spectral method introduced in the ’60s for the analysis of time series, mainly in the field of speech recognition and sound engineering [24]. To deploy its power to extract the transport coefficient from the time series of the relevant conserved fluxes, we must shift to a Fourier-space representation of stochastic processes, as allowed by the Wiener–Khinchin theorem [25, 26]. The latter states that the expectation of the squared modulus of the Fourier transform of a stationary process is the Fourier transform of its time correlation function. We can thus apply this result to the case where the stochastic process is a conserved flux, \(\hat{J}(t)\) [the Cartesian indices have been omitted for clarity], and generalize Eq. (9) to the finite-frequency regime as follows:

$$\begin{aligned} \begin{aligned} S_{\mathcal {T}}(\omega )&= \frac{1}{{\mathcal {T}}} \left\langle \left| \int _0^{\mathcal {T}}\,{\widehat{\!\!{J}}}(t) \mathrm {e}^{i\omega t}\mathrm{d}t \right| ^2 \right\rangle \\&= 2\mathfrak {Re} \int _0^{\mathcal {T}}\left\langle \,{\widehat{\!\!{J}}}(t) \,{\widehat{\!\!{J}}}(0) \right\rangle \mathrm {e}^{i\omega t}\mathrm{d}t + \mathcal {O}({\mathcal {T}}^{-1}). \end{aligned} \nonumber \\ \end{aligned}$$
(29)

More generally, when several fluxes interact with each other, one can define the cross-spectrum of the conserved fluxes as the Fourier transform of the cross time-correlation functions:

$$\begin{aligned} \begin{aligned} S_{lm}(\omega )&= \int _{-\infty }^\infty \langle \,{\widehat{\!\!{J}}}_l(t) \,{\widehat{\!\!{J}}}_m(0) \rangle \,\mathrm {e}^{i\omega t} \mathrm{d}t \\&= \frac{1}{{\mathcal {T}}} \mathfrak {Re} \left\langle \int _0^{\mathcal {T}}\,{\widehat{\!\!{J}}}_l(t) \mathrm {e}^{-i\omega t}\mathrm{d}t \int _0^{\mathcal {T}}\,{\widehat{\!\!{J}}}_m(t) \mathrm {e}^{i\omega t}dt \right\rangle \\&\quad + \mathcal {O}({\mathcal {T}}^{-1}). \end{aligned} \end{aligned}$$
(30)

Onsager’s coefficients can be thus expressed as

$$\begin{aligned} \Lambda _{lm} =\frac{\Omega }{2 k_B} S_{lm}(\omega =0). \end{aligned}$$
(31)

As we shall see, one can leverage on the Wiener–Khinchin theorem and the gauge invariance principles to obtain good estimates (i.e. within  10% accuracy) of the transport coefficients with relatively short EMD simulations (i.e. 10–100 ps). In practice, this result is based on a particular spectral method named cepstral analysis of time series, which we describe below.

5.2 Periodograms and power spectra

Let us focus on one specific flux. In MD simulation, we shall have it as a (discrete time) sample, here denoted with the calligraphic font, of the flux process:

$$\begin{aligned} \mathscr {J}_n \equiv \mathscr {J}(n\epsilon ), \quad n=1,\ldots ,N-1. \end{aligned}$$
(32)

Here \(\epsilon \) is the sampling period, in general a multiple of the timestep of the simulation \(\Delta _t\), so that \(N\epsilon \) is the total length of the simulation. The discrete Fourier transform of the flux time series is defined by

$$\begin{aligned} \tilde{\mathscr {J}}_{k}=\sum _{n=0}^{N-1} \mathrm {e}^{ 2\pi i\frac{kn}{N}} \mathscr {J}_n, \end{aligned}$$
(33)

for \(0 \le k \le N-1\). The sample power spectrum \(\mathscr {S}_k\) (i.e. the periodogram) is defined as

$$\begin{aligned} \mathscr {S}_{k}=\frac{\epsilon }{N} \left| \tilde{\mathscr {J}}_{k} \right| ^2, \end{aligned}$$
(34)

and, for large N, it is an unbiased estimator of the power spectrum of the process, as defined in Eq. (29), evaluated at

$$\begin{aligned} \omega _k= {\left\{ \begin{array}{ll} 2\pi \frac{k}{N\epsilon } &{} \text {for } k\le \frac{N}{2} \\ -\omega _{k-\frac{N}{2}} &{} \text {for } k > \frac{N}{2}, \end{array}\right. } \end{aligned}$$
(35)

namely

$$\begin{aligned} \langle \mathscr {S}_k \rangle = S(\omega _k). \end{aligned}$$
(36)

Since \(\mathscr {J}_{n} \in \mathbb {R}\), we have

$$\begin{aligned} \tilde{\mathscr {J}}_k =\tilde{\mathscr {J}}^*_{N-k}, \mathscr {S}_k=\mathscr {S}_{N-k}, \end{aligned}$$
(37)

which, in the continuous limit, amounts to saying that the power spectrum is an even function of frequency: \(\mathscr {S}(\omega ) = \mathscr {S}(-\omega )\). Thanks to this last point, periodograms are usually reported for \(0\le k\le \frac{N}{2}\) and their Fourier transforms are evaluated as discrete cosine transforms. The space autocorrelations of conserved currents \(\jmath ({\varvec{r}},\Gamma (t))\) are usually short-ranged. Therefore, in the thermodynamic limit, the corresponding fluxes \(J(\Gamma (t)) = \Omega ^{-1} \int _\Omega \jmath ({\varvec{r}},\Gamma (t))\, \mathrm{d} {\varvec{r}}\) can be considered sums of (almost) independent identically distributed (iid) stochastic variables: according to the central-limit theorem their equilibrium distribution is Gaussian. Generalizing this argument allows us to conclude that any conserved-flux process is Gaussian as well. The flux time series is in fact a multivariate stochastic variable that, in the thermodynamic limit, results from the sum of (almost) independent variables, thus tending to a multivariate normal deviate. In particular,

  • for \(k=0\) or \(k=\frac{N}{2}\), \(\tilde{ \mathscr {J}}_k \sim \mathcal {N}\left( 0, \frac{N}{\epsilon }S(\omega _k) \right) \in \mathbb {R}\);

  • for \(k\notin \left\{ 0,\frac{N}{2}\right\} \), \(\mathfrak {Re}\tilde{\mathscr {J}}_k\) and \(\mathfrak {Im}\tilde{\mathscr {J}}_k\) are independent and both \(\sim \mathcal {N}\left( 0, \frac{N}{2 \epsilon }S(\omega _k) \right) \).

Here \(\mathcal {N} (\mu ,\sigma ^2)\) indicates a normal deviate with mean \(\mu \) and variance \(\sigma ^2\). In conclusion, in the large-N (i.e. long-time) limit the periodogram of the time series reads:

$$\begin{aligned} \mathscr {S}_{k} = S \left( \omega _k\right) {\mathscr {\xi }}_{k}, \end{aligned}$$
(38)

the \(\mathscr {\xi } \sim \frac{1}{2\ell } \chi ^2_{2\ell }\) being independent random variables, where \(\chi ^2_{2\ell }\) is the chi-square distribution with \(2\ell \) degrees of freedom and \(\ell \) is the number of independent samples of the current (for instance, \(\ell =3\) when the 3 equivalent Cartesian components of the flux are considered). In particular, \(\bigl \langle \mathscr {\xi }\bigr \rangle =1\) and \(\mathrm {var} \bigl ( \mathscr {\xi }\bigr ) = 1/\ell \).

Equation (38) shows that \(\mathscr {S}_{k=0}\) is an unbiased estimator of the zero-frequency value of the power spectrum, i.e. that \(\langle {\mathscr {S}_0} \rangle = S(0)\), and, through Eq. (31), of the Onsager coefficient we are after. However, this estimator is not consistent, i.e. its variance does not vanish in the large-N limit: The multiplicative nature of the statistical noise makes it difficult to disentangle it from the signal. A way to solve the problem is to apply the logarithm to Eq. (38) to turn the multiplicative noise into an additive one by defining the log-periodogram, \(\mathscr {L}_{k}\), as

$$\begin{aligned} \begin{aligned} \mathscr {L}_{k}&= \log \mathscr {S}_{k} = \log \left( S(\omega _k) \right) + \log ( {\mathscr {\xi }}_k). \\ \end{aligned} \end{aligned}$$
(39)

The quantities \(\log ( {\mathscr {\xi }}_k)\) are iid stochastic variables whose statistics is well known: their mean and variance are simply expressed in terms of the digamma and trigamma functions, \(\psi (\ell )\) and \(\psi ^\prime (\ell )\), respectively [27]. Furthermore, whenever the number of significant (inverse) Fourier components of \(\log (S(\omega ))\) is much smaller than the length of the time series, applying a low-pass filter to Eq. (39) would result in a reduction of the power of the noise, without affecting the signal. To exploit this idea, we define the cepstrum of the time series as the inverse Fourier transform of its sample log-spectrum [28]:

$$\begin{aligned} \mathscr {C}_{n} = \frac{1}{N}\sum _{k=0}^{N-1} \mathscr {L}_{k} \mathrm {e}^{-2\pi i\frac{kn}{N}}. \end{aligned}$$
(40)

A generalized central-limit theorem for Fourier transforms of stationary time series ensures that, in the large-N limit, these coefficients are a set of independent (almost) identically distributed zero-mean normal deviates [29, 30]. It also follows that:

$$\begin{aligned} \begin{aligned} C_n \equiv \langle \mathscr {C}_{n} \rangle&= \frac{1}{N}\sum _{k=0}^{N-1} \log \bigl (S(\omega _k) \bigr ) \mathrm {e}^{-2\pi i\frac{kn}{N}}. \end{aligned} \end{aligned}$$
(41)

Figure 3 confirms the typical behaviour of cepstral coefficients calculated from the low-frequency region of the periodogram. We see that only a few coefficients are in fact substantially different from zero, within statistical uncertainty.

Fig. 3
figure 3

First 40 cepstral coefficients for the log-spectrum of the charge flux in an ab initio MD simulation of molten KCl. They are significantly non vanishing only up to \(n \approx 15\). The vertical dashed line indicates the number of coefficients to retain according to the AIC

Full size image

Let us indicate by \(P^*\) the smallest integer such that \(C_n \approx 0\) for \(P^* \le n \le N-P^*\). By limiting the Fourier transform of the sample cepstrum, Eq. (40), to \(P^*\) coefficients, we obtain an efficient estimator of the zero-frequency component of the log-spectrum, \(\mathscr {L}_0^*\), whose expectation and variance are

$$\begin{aligned} L_0^*\equiv \langle \mathscr {L}_{0}^{*} \rangle&= \log (S_{0}) + \psi (\ell ) - \log (\ell ), \end{aligned}$$
(42)
$$\begin{aligned} \text {var}(\mathscr {L}_{0}^{*})&=\psi ^\prime (\ell )\frac{4P^{*}-2}{N}. \end{aligned}$$
(43)

We thus see that the logarithm of the Onsager coefficient we are after can be estimated from the cepstral coefficients of the flux time series through Eqs. (42-43), and that the resulting estimator is always a normal deviate whose variance depends on the specific system only through the number of these coefficients, \(P^*\). Notice that the absolute error on the logarithm of the conductivity directly and nicely yields the relative error on the conductivity itself. The efficacy of this approach obviously depends on our ability to estimate the number of coefficients necessary to keep the bias introduced by the truncation to a value smaller than the statistical error, while maintaining the magnitude of the latter at a prescribed acceptable level. In Ref. [8], it has been proposed to estimate \(P^*\) using the Akaike’s information criterion [31], even if other more advanced model selection approaches may be more effective [32]. A plot of \(L_0^*\) vs \(P^*\) is shown in Fig. 4. We immediately realize that \(P^*\) returned by the AIC (vertical dashed line) is indeed capable to find the correct value for the log-spectrum—within statistical error—of lowest variance. Furthermore, thanks to the convective invariance principle described in Sect. 4, the cepstral analysis can be extended to multicomponent systems [7]. This multivariate cepstral method has been recently applied to calculate the ab initio thermal conductivity of water at planetary conditions from trajectories as short as a few tens of picoseconds [33]. It has been also shown that the multivariate cepstral analysis is able to substantially reduce the statistical error affecting the estimate of thermal conductivity even for one-component systems: it can in fact decorrelate the finite-frequency power spectrum of non-diffusive fluxes (like the mass flux, or the adiabatic electronic flux in ab initio simulations) from the heat flux power spectrum, which will have its total power considerably reduced, and whose low-frequency portion will be easier to analyse [7]. The statistical tools for time-series analysis presented in this Section have been implemented in the open-source code SporTran, which is freely downloadable from the GitHub repository https://github.com/lorisercole/sportran [34].

Fig. 4
figure 4

The \(\omega =0\) component of the log-spectrum of the charge flux vs \(P^*\) for an ab initio MD simulation of molten KCl. The expectation value (dots, Eq. (42)) and its uncertainty (shaded area, square root of Eq. (43)) are displayed. The vertical dashed line indicates the optimal \(P^*\) predicted with the AIC

Full size image

6 Ab initio heat transport in insulators

Until only a few years ago, it was thought that ab initio MD simulations were in general unsuitable to a GK theory of thermal transport, since the continuous electronic density makes any decomposition of the energy flux into local, atomic contributions fully arbitrary [35]. This consideration clashes with a reductionist picture whereby a fundamental description of the microscopic interactions should be in principle suitable to describe heat transport in the very general GK theory, as well. Once again, this apparent inconsistency stands on the misconception that the definition of microscopic fluxes must be unique, and it is thus solved thanks to the gauge-invariance principle. Apart from correcting such a misconception, the gauge-invariance principle proves to be also a rigorous mathematical tool to derive a well-defined expression (out of the infinitely many possibilities!) for the energy flux directly from DFT, with no ad hoc approximation tailored on the considered physical system [36]. The starting point is the standard DFT definition of the total energy in terms of the Kohn–Sham (KS) eigenvalues \(\varepsilon _v\), eigenfunctions \(\phi _v({\varvec{r}})\), and density \(n({\varvec{r}}) = \sum _v |\phi _v({\varvec{r}})|^2\) [37]:

$$\begin{aligned} E_{{\scriptscriptstyle DFT}}= & {} \frac{1}{2}\sum _{n}M_{n}V_{n}^{2} + \frac{\mathtt {e}^2}{2}\sum _{n,m\ne n}\frac{ Z_{n}Z_{m}}{|{\varvec{R}}_{n}-{\varvec{R}}_{m}|} \nonumber \\&+ \sum _{v}\varepsilon _{v}-\frac{\mathtt {e}^2}{2}\int \frac{n({\varvec{r}})n({\varvec{r}}')}{|{\varvec{r}}-{\varvec{r}}'|}d{\varvec{r}}\mathrm{d}{\varvec{r}}' \nonumber \\&+\int \left( \epsilon _{{\scriptscriptstyle XC}}[n]({\varvec{r}})-\mu _{{\scriptscriptstyle XC}}[n]({\varvec{r}})\right) n({\varvec{r}})\mathrm{d}{\varvec{r}}, \end{aligned}$$
(44)

where \(\mathtt {e}\) is the electron charge, \(\epsilon _{\scriptscriptstyle XC}[n]({\varvec{r}})\) is a local exchange-correlation (XC) energy per particle defined by the relation \( \int \epsilon _{\scriptscriptstyle XC}[n]({\varvec{r}})n({\varvec{r}}) \mathrm{d}{\varvec{r}}=E_{\scriptscriptstyle XC}[n]\), the latter being the total XC energy of the system, and \( \mu _{\scriptscriptstyle XC}({\varvec{r}}) = \frac{\delta E_{\scriptscriptstyle XC}}{\delta n({\varvec{r}})}\) is the XC potential [in this Section and in Sect. 7 we drop the hat on top of classical observables to avoid confusion with quantum operators]. From this we can write a DFT energy density as [38]

$$\begin{aligned} \begin{aligned} E_{{\scriptscriptstyle DFT}}&= \int e_{{\scriptscriptstyle DFT}}({\varvec{r}})\mathrm{d}{\varvec{r}},\\ e_{{\scriptscriptstyle DFT}}({\varvec{r}})&= e_{el}({\varvec{r}})+e_{{\scriptscriptstyle Z}}({\varvec{r}}), \end{aligned} \end{aligned}$$
(45)

where

$$\begin{aligned} e_{el}({\varvec{r}})&=\mathfrak {Re} \sum _{v}\phi _{v}^{*}({\varvec{r}})\bigl (H_{{\scriptscriptstyle KS}}\phi _{n}({\varvec{r}})\bigr ) \nonumber \\&\quad - \frac{1}{2}n({\varvec{r}})v_{{\scriptscriptstyle H}}({\varvec{r}}) +\left( \epsilon _{\scriptscriptstyle XC}({\varvec{r}}) - \mu _{\scriptscriptstyle XC}({\varvec{r}}) \right) n({\varvec{r}}), \end{aligned}$$
(46)
$$\begin{aligned} e_{{\scriptscriptstyle Z}}({\varvec{r}})&= \sum _{n}\delta ({\varvec{r}}-{\varvec{R}}_{n}) \left( \frac{1}{2}M_{n}V_{n}^{2}+w_{n}\right) , \end{aligned}$$
(47)
$$\begin{aligned} w_{n}&=\frac{\mathtt {e}^2}{2}\sum _{m\ne n}\frac{ Z_{n}Z_{m}}{|{\varvec{R}}_{n}-{\varvec{R}}_{m}|}, \end{aligned}$$
(48)

\(H_{\scriptscriptstyle KS}\) is the instantaneous self-consistent Kohn–Sham Hamiltonian, and \(v_{\scriptscriptstyle H}= \mathtt {e}^2 \int \mathrm{d}{\varvec{r}}' \frac{n({\varvec{r}}')}{|{\varvec{r}}-{\varvec{r}}'|}\) is the Hartree potential. An explicit expression for the DFT energy flux,

$$\begin{aligned} {\varvec{J}}^{\scriptscriptstyle E}_{\scriptscriptstyle DFT}= \frac{1}{\Omega } \int {\varvec{r}} \dot{e}_{{\scriptscriptstyle DFT}}({\varvec{r}}) \mathrm{d}{\varvec{r}}, \end{aligned}$$
(49)

obtained by computing the first moment of the time derivative of the energy density in Eqs. (45–48), results in a number of terms, some of which are either infinite or ill-defined in PBC, since the position operator is not periodic. Leveraging on gauge invariance and on a careful breakup and refactoring of the various harmful terms, as explained in Refs. [5, 39], we can recast Eq. (49) in a form suitable to PBC, whose final expression is

$$\begin{aligned} {\varvec{J}}^{\scriptscriptstyle E}_{{\scriptscriptstyle DFT}} ={\varvec{J}}^{{\scriptscriptstyle H}} + {\varvec{J}}^{{\scriptscriptstyle Z}} + {\varvec{J}}^{0} + {\varvec{J}}^{{\scriptscriptstyle KS}} + {\varvec{J}}^{{\scriptscriptstyle XC}}, \end{aligned}$$
(50)

where

$$\begin{aligned}&{\varvec{J}}^{{\scriptscriptstyle H}} = \frac{1}{4\pi {\mathrm \Omega }\mathtt {e}^2}\int \nabla v_{{\scriptscriptstyle H}}({\varvec{r}}) \dot{v}_{{\scriptscriptstyle H}}({\varvec{r}}) \mathrm{d}{\varvec{r}}, \end{aligned}$$
(51)
$$\begin{aligned}&{\varvec{J}}^{{\scriptscriptstyle Z}} = \frac{1}{{\mathrm \Omega }} \sum _{n} \left[ {\varvec{V}}_{n}\left( \frac{1}{2}M_{n}V_{n}^{2} + w_{n}\right) \right. \nonumber \\&\qquad \qquad \left. + \sum _{m\ne n}({\varvec{R}}_{n} - {\varvec{R}}_{m}) \left( {\varvec{V}}_{m} \cdot \frac{\partial w_{n}}{\partial {\varvec{R}}_{m}} \right) \right] , \end{aligned}$$
(52)
$$\begin{aligned}&{\varvec{J}}^{0} = \frac{1}{{\mathrm \Omega }} \sum _{n} \sum _{v}\left\langle \phi _{v} \left| (\hat{{\varvec{r}}}-{\varvec{R}}_{n})\left( {\varvec{V}}_{n}\cdot \frac{\partial \hat{v}_{0}}{\partial {\varvec{R}}_{n}}\right) \right| \phi _{v}\right\rangle , \end{aligned}$$
(53)
$$\begin{aligned}&{\varvec{J}}^{{\scriptscriptstyle KS}} = \frac{1}{{\mathrm \Omega }} \mathfrak {Re} \sum _{v} \langle {\varvec{\bar{\phi }}}_v^c | \hat{H}_{{\scriptscriptstyle KS}}+\varepsilon _v | \dot{\phi }_v^c \rangle , \end{aligned}$$
(54)
$$\begin{aligned}&J_{\alpha =x,y,z}^{{\scriptscriptstyle XC}} = {\left\{ \begin{array}{ll} 0 &{} \mathrm {(LDA)} \\ -\frac{1}{{\mathrm \Omega }} \int n({\varvec{r}}) \dot{n}({\varvec{r}}) \frac{\partial \epsilon ^{{\scriptscriptstyle GGA}} ({\varvec{r}})}{\partial (\partial _\alpha n)} \mathrm{d}{\varvec{r}} &{} \mathrm {(GGA)}. \end{array}\right. } \end{aligned}$$
(55)

Here \(\hat{v}_0\) is the bare, possibly non-local, (pseudo-)potential acting on the electrons and

$$\begin{aligned} \begin{aligned} |{\varvec{\bar{\phi }}}_v^c\rangle&= \hat{P}_c \,\hat{{\varvec{r}}} \,|\phi _v \rangle , \\ |\dot{\phi }_v^c\rangle&= \dot{\hat{P}}_v \,|\phi _v\rangle \end{aligned} \end{aligned}$$
(56)

are the projections over the empty-state manifold of the action of the position operator over the vth occupied orbital, and of its adiabatic time derivative [40,41,42], respectively, \(\hat{P}_v\) and \(\hat{P}_c = 1 - \hat{P}_v\) being the projector operators over the occupied- and empty-states manifolds, respectively. Both these functions are well defined in PBC and can be computed, explicitly or implicitly, using standard density-functional perturbation theory [42]. These rather complicated formulae have been implemented in the open-source code QEHeat of the Quantum ESPRESSO project which has been publicly released [43].

7 Ab initio charge transport in ionic conductors

7.1 Ab initio charge transport in ionic conductors

While in metals the current is carried by delocalised conduction electrons, in electronic insulators the electrons are bound to follow adiabatically the ionic motion and no charge transport can occur if the ion positions are frozen. Nonetheless, when ions are allowed to move, like in ionic liquids, charge can be displaced. Daily life examples range from simple salted water, to liquid electrolytes employed in Li-ion batteries, or to the molten salts (NaCl, KCl, etc.) used as heat exchangers in power plants.Footnote 1 Due to their large electronic bandgap, ionic liquids are in general transparent to visible light and possess a negligible fraction of “free”, conduction electrons. When the quantum nature of the electrons is considered, a question arises about a proper definition for the charge flux, \({\varvec{J}}(t)\), to employ in the GK expression for the electrical conductivity, \(\sigma \):

$$\begin{aligned} \sigma&= \frac{\Omega }{3k_B T} \int _0^\infty \langle \varvec{{J}}(t) \cdot \varvec{{J}}(0) \rangle \, dt. \end{aligned}$$
(57)

Any static partition of the instantaneous total charge to express \({\varvec{J}}\) into atomic contributions is in fact totally arbitrary.

The solution comes from the modern theory of polarisation (MTP), which provides a definition for the polarisation \({\varvec{P}}\) valid for extended systems [45,46,47]. Just like the electronic Hamiltonian and ground state, \({\varvec{P}}\) depends on time through the nuclear coordinates only: we can thus apply the chain rule and write \(\varvec{J}(t)\) as a sum of atomic contributions

$$\begin{aligned} \begin{aligned} \varvec{{J}}(t)&= \dot{{\varvec{P}}}(t) = \frac{\mathtt {e}}{\Omega } \sum _{i=1}^{N_{at}} \mathrm {Z}^*_{i}(t) \cdot \varvec{V}_i(t), \end{aligned} \end{aligned}$$
(58)

where \({\varvec{V}}_i(t)\) is the atomic velocity of the ith atom and \(\mathrm {Z}^*_i(t)\) is a time-dependent tensor whose entries, \(Z^*_{i\alpha \beta } \equiv \tfrac{1}{\mathtt {e}}\tfrac{\partial P_\alpha }{\partial R_{i\beta }}\), are the derivative of the cell polarisation along the Cartesian direction \(\alpha \) with respect to the atomic displacement of atom i along \(\beta \). The dynamical charges \(\mathrm {Z}_i\) are called “Born effective-charge tensors”, and can be computed in perturbation theory [43], or by finite force differences when introducing a finite electric field \({\varvec{E}}\) in accordance to the MTP [48]. Born charge tensors are in general strongly dependent on the atomic positions, and they display large fluctuations along a AIMD trajectory, even for strongly ionic systems [49]. For these reasons, the charge flux defined in Eq. (58) is in general fluctuating and does not vanish identically. Therefore, there is no apparent reason for the GK integral in Eq. (57) to be exactly zero. How come, then, that the electrical conductivity of, say, pure water is zero? And why, instead, under the same microscopic formalism, does a salted water solution have a non-vanishing \(\sigma \)? Another strange “coincidence” reported in the literature [50] is that by replacing in Eq. (58) the time-dependent Born charge tensors with the predefined constant integer charges—the oxidation numbers of the atoms—to define a new flux

$$\begin{aligned} {\varvec{J}}'(t) = \frac{\mathtt {e}}{\Omega } \sum _{i=1}^{N_{at}} q_{S(i)} {\varvec{V}}_i(t) \end{aligned}$$
(59)

leads to the same GK integral, Eq. (57), that is obtained via the alternative definition \({\varvec{J}}\) in Eq. (58).Footnote 2

It seems, therefore, that two main aspects must be understood to answer the fundamental questions raised above: first, the presence of integer numbers suggests some underlying quantisation in the charge transport process; second, the insensitivity of \(\sigma \) to the choice of one of the two different definition \({\varvec{J}} \ne {\varvec{J}}'\) hints at a manifestation of the gauge invariance principle of transport coefficients, described in Sect. 3. We analyse in detail these two aspects in what follows.

7.2 Quantisation of charge transport

In 1983, David Thouless showed that for a quantum system in PBC, with non-degenerate ground state evolving adiabatically under a cyclic, slowly varying Hamiltonian \(\hat{H}(0) = \hat{H}({\tau })\), the total displaced dipole, \(\Delta {\varvec{\mu }}\) is quantised, i.e. the time integral of the current over one cycle \({\tau }\) is equal to a triplet of integers, multiplied by the side L of the cell, which we shall assume cubic for simplicity [51]:

$$\begin{aligned} \Delta {\varvec{\mu }}= \Omega \int _0^{\tau }{\varvec{J}}(t) \mathrm{d}t = \mathtt {e} L {\varvec{Q}} \end{aligned}$$
(60)

with \({\varvec{Q}} = (Q_x,Q_y,Q_z) \in \mathbb {Z}^3\). In AIMD simulations, atomic trajectories are identified by paths in the nuclear configuration space (NCS). The electronic Hamiltonian, \(H(\mathbf {R}(t))\), depends parametrically upon time through the nuclear positions, at time t, \(\mathbf {R}(t) = \left\{ \varvec{R}_i(t) \right\} \), with \(i = 1, \dots , N_{at}\). A path \(\mathcal {C}\) in NCS whose end points are one the the periodic image of the other can be thus considered as a cycle for the electronic Hamiltonian. Therefore, we can thus employ Thouless’ theorem and show, under very general hypothesis, that the electric dipole displaced along \(\mathcal {C}\) is itself partitioned into atomic contributions:

$$\begin{aligned} \Delta {\varvec{\mu }}_\mathcal {C} = \mathtt {e} L \sum _{i=1}^{N_{at}} q_S(i) {\varvec{n}}_i, \end{aligned}$$
(61)

where \({\varvec{n}}_i = (n_{x},n_{y},n_{z})_i\) is the triplet containing the number of cells spanned, along \(\mathcal {C}\), by atom i in x, y and z directions, and \(q_S(i)\) are integer constants which only depend on the species S(i) of atom i, and are shown to coincide with the oxidation numbers suggested by chemical intuition [49]. This provides a quantum-mechanical definition of oxidation numbers where their integerness arises naturally (and not just as an approximation of some real charge), and which identifies them as intrinsically dynamical quantities [49, 52].

7.3 Gauge invariance of electrical conductivity

We are now ready to combine the theory of quantisation of charge transport and the gauge-invariance principle of transport coefficients to answer the questions raised at beginning of this Section. Let us consider a physical path in the NCS from an initial configuration I to the configuration F, as a result, e.g. of an AIMD simulation. Then we elongate the path fictitiously up to the point \(I'\) which is the replica (periodic image) of the initial point I sharing with the point F the same cell of the nuclear configuration space, as depicted by the dashed line in Fig. 5.

Fig. 5
figure 5

Paths in a two-dimensional NCS representing atomic trajectories. \(I'\) is a periodic image of the starting point I belonging to the same periodic cell of the physical endpoint F. The path \(II'\) is the concatenation of the physical trajectory IF with the path \(I'F\) entirely belonging to one cell

Full size image

Due to the additivity of integrals, the electric dipole displaced along the physical path IF reads

$$\begin{aligned} \Delta {\varvec{\mu }}_{IF}&\equiv \int _{IF} \mathrm{d}{\varvec{\mu }} = \Delta {\varvec{\mu }}_{II'} + \Delta {\varvec{\mu }}_{I'F} \;. \end{aligned}$$
(62)

Since the open path \(I'F\) entirely belongs to one cell, \(\Delta {\varvec{\mu }}_{I'F}\) is a bounded quantity. Therefore, thanks to gauge-invariance, to evaluate \(\sigma \) we only need to consider

$$\begin{aligned} \Delta {\varvec{\mu }}_{II'} = \int _{II'} d{\varvec{\mu }} = \mathtt {e} L \sum _{i=1}^N q_{S(i)} {\varvec{n}}_{i} \end{aligned}$$
(63)

since

$$\begin{aligned} \sigma \propto \lim _{{\mathcal {T}}\rightarrow \infty } \frac{\langle |\Delta {\varvec{\mu }}_{IF}({\mathcal {T}})|^2 \rangle }{2{\mathcal {T}}} = \lim _{\tau \rightarrow \infty } \frac{\langle \left| \Delta {\varvec{\mu }}_{II'}({\mathcal {T}})\right| ^2 \rangle }{2{\mathcal {T}}}. \end{aligned}$$
(64)

By the same token, the electric dipole displaced from I to F by to the flux \({\varvec{J}}^\prime (t)\) defined in Eq. (59) is

$$\begin{aligned} \Delta \varvec{\mu ^\prime }_{IF} = \Delta \varvec{\mu }_{II'} +\mathtt {e} \sum _{i=1}^N q_{S(i)} \underbrace{\int _{I'}^F d\varvec{R}_i}_{\text {bounded}} . \end{aligned}$$
(65)

Therefore, we can conclude that

$$\begin{aligned} {\begin{matrix} \sigma ' \propto \lim _{t \rightarrow \infty } \frac{\langle \left| \Delta {\varvec{\mu ^\prime }}_{IF}({\mathcal {T}})\right| ^2 \rangle }{2{\mathcal {T}}}&= \lim _{{\mathcal {T}}\rightarrow \infty } \frac{\langle \left| \Delta {\varvec{\mu }}_{II'}({\mathcal {T}})\right| ^2 \rangle }{2{\mathcal {T}}} \end{matrix}} \end{aligned}$$
(66)

which, by comparison with Eq. (64), proves the equivalence of the electrical conductivities obtained via Eqs. (58) and (59). This is shown in Fig. 6 for an ab initio MD simulation of molten KCl.

Fig. 6
figure 6

Mean square dipole displaced along a trajectory of molten KCl. The slope of the asymptotic behaviour is proportional to the electrical conductivity. Calculations employing the two different definitions Eqs. (58) and (59) share the same slope, and thus the same electrical conductivity. The difference \({\varvec{J}}' - {\varvec{J}}\) is instead non-diffusive (zero slope). The fits of the curves at asymptotic times are also displayed (dashed lines)

Full size image

As explained in detail in Ref. [49], this result is grounded on the hypothesis, which we dubbed strong adiabaticity, that any closed paths in the NCS can be shrunk to a point without closing the electronic gap: this implies—see Eq. (61)—that charge transport can occur only through a net displacement of the ions, as it is typical of stoichiometric ionic conductors. The breach of strong adiabaticity may instead dictate a non-trivial charge-transport regime where charge may be adiabatically transported even in the absence of a net ionic displacement. As shown in Ref. [53] this non-trivial behaviour is intertwined with the presence, in non-stoichiometric electrolytes, of dissolved yet localised electrons, whose displacement is to a large extent uncorrelated to that of the ions. By the same token discussed in Sect. 3.2, we remark that all the conclusions drawn in this Section are independent of the (macroscopic) size L of the system. In fact, even though all the derivation of Eq. (66) is done at finite size, the use of charge fluxes which are well defined in PBC ensures that the GK integrals of their correlation functions are well defined and non-vanishing identically, for any L.

8 Conclusions

To conclude, we believe that we managed to show how the invariance principles of transport coefficients can be employed, within the general Green–Kubo theory of linear response, to demystify some common misconceptions based on the groundless assumption that the microscopic conserved fluxes must be uniquely defined. We have further deployed the power of invariance principles in the construction of numerical tools for the statistical analysis of the time series of the fluxes, produced via equilibrium molecular dynamics simulations. These tools provide accurate values for the transport coefficients from the relatively short trajectories which are accessible to ab initio MD simulations. We also showed how to design an ab initio heat flux suitable to extended systems in PBC, from which the ab initio thermal conductivity can be computed without any further system- or phase-specific assumptions. Finally, we have applied the gauge invariance principle to the ab initio charge transport in ionic systems. We showed that each ion can be associated with a well-defined integer and time-independent charge, and the exact ab initio electrical conductivity can be obtained by replacing, with these integer charges, the time-dependent Born tensors, which enter the definition of the charge flux and whose calculation is a computationally expensive and quite abstruse task. In this way, we recover the classical (Faraday’s) idea of atomic contributions to charge transport [54] and provided a theoretically sound definition to the concept of oxidation states in ionic liquid insulators. We are confident that the results here exposed will be a valid aid for both theorists and practitioners aiming at deeper insights and more efficient implementations about the calculation of transport coefficients from molecular dynamics simulations.