Thermal diffusion in binary mixtures: transient behavior and transport coefficients from equilibrium and non-equilibrium molecular dynamics

Equilibrium and non-equilibrium molecular dynamics simulations are combined to compute the full set of coefficients that appear in the phenomenological equations describing thermal transport in a binary mixture subject to a constant thermal gradient. The Dynamical Non-Equilibrium Molecular Dynamics approach (D-NEMD) is employed to obtain the microscopic time evolution of the density and temperature fields, together with that of the mass and energy fluxes. D-NEMD enables to study


Introduction
Thermal transport in simple dense fluids is an interesting phenomenon whose description in the framework of classical statistical mechanics has been elegantly reviewed, in particular, in a vintage chapter by K. E. Gubbins 1 . In this paper, we shall focus on the coupled mass and energy diffusion in a binary mixture subject to a constant thermal gradient and consider the evaluation of the thermal transport coefficients via computer simulation. It is well known 2-4 that a system subjected to a constant thermal gradient will be driven to a non-equilibrium steady state such that energy (heat) flows at a constant rate through it, while the mass flux stops and a constant concentration gradient is established. This thermal diffusion process, known as the Ludwig-Soret effect, was first observed in 1856 in sulfur-solfate solutions and has since been identified experimentally in a broad range of systems, ranging from gases to polymer mixtures [5][6][7] . The experimental characterization of the Ludwig-Soret effect is, however, challenging 8 . In mutual diffusion, in fact, thermally driven flows are considerably smaller than, for example, density driven flows and therefore difficult to measure accurately. Molecular Dynamics (MD) based simulations have then been used 4,[9][10][11][12][13][14][15][16][17][18][19][20][21] , together with purely theoretical work, to investigate thermal transport and compute, in particular, the Soret coefficient. Important technical issues, such as the proper definition of the microscopic estimator for the heat flux, have been clarified by these calculations 12 , but they are limited to the steady state. As discussed more in detail in the 2 following, the approach that we employ, the D-NEMD method 22 , makes it possible to study also the onset of the transport phenomena and the evolution of hydrodynamic fields towards the steady state, if the latter exists. We shall exploit this feature to study the Ludwig-Soret effect focusing at first on the transient response of the system as described by the relaxation of the time dependent number and temperature density fields and the mass and energy currents. Following this relaxation in time, the transient transport mechanisms and the time scales for the different fluxes will be investigated in some detail. Calculations in the steady state will then be combined with the evaluation, from the time autocorrelation function of the mass current at equilibrium, of the mutual diffusion coefficient to compute the Onsager thermal transport coefficients 23,24 that characterize the Ludwig-Soret effect.
The paper is organized as follows. We begin by summarizing the key phenomenological equations, indicating the relevant transport coefficients and how they will be calculated from simulations. This is followed by a description of the D-NEMD approach and the definition of the estimators for the microscopic fields and fluxes studied. A smoothing scheme adopted to reduce the noise on the estimated currents is also discussed. We conclude the Methods and Model section providing details on the adopted model (an equimolar mixture of Lennard-Jones Argon and Krypton) and on the set-up chosen to simulate the constant thermal gradient. Results are then discussed, followed by some conclusions.

Macroscopic model and relevant coefficients
To set the stage and introduce some notation, let us recall that the key fields for describing thermal diffusion in a binary mixture are J e , the energy current density, and J 1 the diffusion current density of species 1 relative to the center of mass frame of reference. For a binary mixture J 2 = −J 1 , where the subscripts 1, 2 are the indexes of the species. In the following 1 will refer to Kr and 2 to Ar. Adopting a hydrodynamic point of view, the divergence 3 of the energy current and species diffusion current drives the time evolution of the energy and species mass density fields, respectively (the other continuity equations relating the time derivatives of the total mass density to the divergence of the total momentum, and the time derivative of the momentum density to the divergence of the stress tensor are not relevant in this work). In 1971, Trimble and Deutch 25 , derived the linearized hydrodynamic equations for a two-component mixture starting from a microscopic description based on Kubo's linear response theory. They showed, in particular, that the currents of interest here can be expressed as where κ is the coefficient of thermal conductivity (as in the Fourier law), D is the mutual diffusion coefficient (as in Fick's law) and D T the coefficient of thermal diffusion (as in Soret's effect). The equations above refer to the case, considered in the following, of an isotropic system. In the more general case, e.g. liquid crystals, the thermal conductivity is in fact a tensor. Note that here and in the following we adopt the notation of Trimble and Deutch to identify the thermal transport coefficients and the relevant fields (see, for example, the discussion in the Appendix of Trimble and Deutch 25 for a review of possible alternative definitions). In the same paper 25 , expressions for the transport coefficients in terms of time correlation functions were obtained.
The connection between the coefficients computed via time correlation functions at equilibrium and quantities directly accessible to experiments, however, is problematic for mixtures because it requires auxiliary calculations of thermodynamic quantities such as partial enthalpies and derivatives of the chemical potential 26,27 . The exception is the mutual diffusion coefficient D for which a straightforward Green-Kubo expression is established 25 .
Moreover, useful relationships exist for the thermal transport coefficients and energy current 4 densities in the, non-equilibrium, stationary states 25 . When combined, as discussed below, these relationships make it possible to compute the thermal conductivity and diffusion coefficients more straightforwardly. In this work, we propose to take advantage of this observation to define a suitable computational approach. To see how, we recall that for a system subject to a temperature gradient the stationary state is such that where the coefficient λ T is defined implicitly in the last equation. The thermal and mutual diffusion coefficients are related, in the definition of the Soret coefficient of the mixture 6 , by S T = D T /D. This coefficient, accessible also in experiments, provides a measure of the relative strength of thermally induced and concentration driven diffusion. In the following, we shall consider an equimolar mixture in which the thermal gradient is directed along the z axis. In this case, the Soret coefficient can be measured directly, in stationary conditions, where χ 1 is the mole fraction of species 1.
Based on the observations above, all the thermal transport coefficients can be computed from atomistic simulations by combining equilibrium and non-equilibrium molecular dynamics simulations. Our approach can be summarized as follows. Equilibrium MD will be used to obtain the mutual diffusion coefficient D by computing the time correlation function of the mass current of species 1 in the system. Non-equilibrium MD will instead be used to simulate the binary mixture under the constant thermal gradient and access quantities in the steady state. D-NEMD will be employed to characterize the transient and to verify the establishment of the steady state by following the evolution of the density and temperature fields, and of the associated currents, from the appearance of the thermal gradient to sta-5 tionarity. In the stationary condition, a direct measurement of the energy flux at a known value of the thermal gradient will provide us with an estimate of the thermal coefficient , while the ratio of the slope of the density (related to the mole fraction) and temperature fields will measure the Soret coefficient. From these calculations we can

Simulation methods and microscopic expressions of observables
As mentioned in the Introduction, the non-equilibrium molecular dynamics simulation will be performed using the D-NEMD approach 22 . This method implements numerically the Onsager principle on the regression of fluctuations 23,24 and has been applied to a variety of problems [28][29][30][31][32][33][34][35] , including a preliminary investigation of the Ludwig-Soret effect 20 . In this approach, the binary mixture is described at a microscopic level as a system of N α particles of each species (α = 1, 2 is the species index) in the phase space Γ ≡ {r i 1 , p i 1 , r i 2 , p i 2 , i 1 = 1, . . . , N 1 , i 2 = 1, . . . , N 2 }, where r iα , p iα are the positions and momenta of particle i of species α. The time evolution of the system is governed, at equilibrium, by the (timeindependent) Hamiltonian H 0 (Γ), and time evolution operator S 0 (t) such that Γ t = S 0 (t)Γ, where we have indicated with Γ t the phase space point at time t, and for convenience of notation, we assume that the initial time of the evolution is t 0 = 0, using the notation According to the prescriptions of statistical mechanics, a generic macroscopic field, indicated as O(r, t), can be obtained as the ensemble average of the corresponding microscopic observable over phase space. If we indicate this microscopic observable asÔ(r, Γ) ≡ The last equation, known as the Kubo-Onsager relationship, provides the formal basis for D-NEMD 36,37 . It states that the time evolution of a macroscopic field can be obtained as the average over the ensemble at the initial time of the time evolved microscopic observable.
The time evolution of the microscopic observable can be obtained via MD for quite general dynamical systems. If the probability density f (Γ, 0) can be sampled by simulation (e.g. when it corresponds either to equilibrium or to stationary conditions) then, the relationship above generates the D-NEMD algorithm to calculate the fields.
This algorithm can be summarized in three steps (see also figure 1): (1) sample a set of initial conditions from f (Γ, 0) (e.g. at equilibrium using Monte Carlo or molecular dynamics driven by H 0 );  (2) evolve these initial configurations under the non-equilibrium dynamics of the system and compute the microscopic observable along each trajectory; (3) compute the macroscopic field as the average of the microscopic observable over the trajectories.
As mentioned in the previous subsection, we shall study the Ludwig-Soret effect in an equimolar Ar-Kr mixture, with constant thermal gradient directed along the z axis of the system.
The macroscopic fields relevant for this phenomenon are the number density field of the species and the temperature field. The corresponding microscopic observables are defined as (in all equations below α = 1, 2)n andT In the equation above,p iα (t) is the momentum in the center of mass reference, i.e.p iα (t) = where v cm is the center of mass velocity, V the volume of the system, and k B is the Boltzmann constant. In all calculations discussed in the following, v cm = 0 and so we drop the tilde on the momenta to simplify the notation. The brackets in eq. (8) 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 an average with respect to the initial ensemble at t 0 = 0 while, consistent with eq. (6), the time evolution is that of the system subject to the thermal gradient. A detailed description of the equilibrium sampling and non-equilibrium dynamics are given in sections Simulation set-up and Results and Discussion. To the density and temperature fields are associated the mass current for each speciesĴ and the energy current 13,38 where φ(r iα , r j β ) is the (pair) potential between particle i of species α and particle j of species β, and the prime indicates that self interactions must not be considered. F iαj β = −∇ iα φ(r iα , r j β ) is the force on particle i of species α due to particle j of species β.
As mentioned above, see eq. (3), given the symmetry used in our simulation, the ratio of the slopes of n 1 (r, t) (which can be related to the mole fraction) and T (r, t) in the steady state will allow us to measure the Soret coefficient, while the ratio of the constant energy flux over the imposed constant thermal gradient, see second line of eq. (2), will provide an estimate of λ T . Finally, the mutual diffusion coefficient (to be computed from equilibrium molecular dynamics simulations) is defined as 25 9 The estimate of the thermal transport coefficients outlined above requires a precise measurement of the currents. The averages of these quantities, however, present a significant level of noise that hides the behavior of the system both in the transient and in the steady state, and prevents converging the D-NEMD calculation with reasonable computational effort 20 . The problem affects all the currents, but, as shown in the following, it is particularly critical for the mass currentĴ 1 of species 1. It is, however, possible to mitigate this difficulty by observing that the typical time scales over which the currents show significant changes are considerably longer than those corresponding to the natural microscopic evolution of the system. One can then take advantage of this fact to smooth the, slow varying, noisy signal where τ w is the time width of the averaging window (nb. in the time intervals τ ∈ [t 0 , t 0 + τ w ] and τ ∈ [t f − τ w , t f ], the limits of integration and the normalization are adjusted accordingly to the available window, centered around τ ). The stability of this smoothing procedure can be verified by varying the width of time window over which the average is taken along each trajectory. Note that this approach still preserves the full average over the whole set of non-equilibrium initial conditions and is therefore very different from attempts to obtain non-equilibrium properties from long time averaging along a single trajectory.
Then particle velocities are rescaled in such a way that, within each slab, the (local) kinetic energy matches the target temperature while the (local) momentum does not change. More in detail, velocity rescaling is performed as follows. At each timestep, the center of mass velocity of the particles in the slab is computed as where the index M can refer either to a slab in the cold region, i.e. M = 1, 2 or 3 or to a slab in the hot region, i.e M = 28, 29 or 30 . The (local) instantaneous kinetic energy, excluding the center of mass motion, in the M -th slab is then calculated as and used to obtain the scaling factor where N M (t) = α iα∈M 1 is the total number of particles in the M -th slab at the time t and T M the target temperature. The particles are finally assigned new momenta defined as As mentioned above, this procedure ensures that, layer by layer, the thermostat does not influence the momentum density and, as a consequence, the total momentum of the whole 13

Conclusions
In this paper, we have combined equilibrium and non-equilibrium simulations to obtain the full set of phenomenological coefficients appearing in the Onsager description of thermal transport. The transient behavior of the system was investigated in detail, to the best of our knowledge for the first time, using the D-NEMD approach to non-equilibrium simulations.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 It presents a picture with a clear time scale separation between a fast regime, t ≈ 50 in which the mixture behaves like a single component fluid, reaching a steady state for the temperature and overall density profiles, and a slow regime in which the behavior of the species differentiates and the Ludwig-Soret effect builds up. The time evolution of the Ar (lighter species in our mixture), in particular, is non trivial, following the profile of the average density for short times, building a definite oscillation and then inverting its overall slope to relax to its linear stationary profile on the longer time scale. In the stationary regime the momentum current of the fluid vanishes while a steady heat flux from the hot to the cold thermostat is established, as illustrated by the asymptotic behavior of the energy current.