On the microscopic origin of Soret coefficient minima in liquid mixtures

Temperature gradients induce mass separation in mixtures in a process called thermodiffusion and quantified by the Soret coefficient. The existence of minima in the Soret coefficient of aqueous solutions was controversial until fairly recently, where a combination of experiments and simulations provided evidence for the existence of this physical phenomenon. However, the physical origin of the minima and more importantly its generality, e.g. in non-aqueous liquid mixtures, is still an outstanding question. We report the existence of a minimum in liquid mixtures of non-polar liquids modelled as Lennard-Jones mixtures, demonstrating the generality of this phenomenon. The Soret coefficient minimum originates from a coincident minimum in the thermodynamic factor, and hence denotes a maximimzation of non-ideality mixing conditions. We explain the microscopic origin of this effect in terms of the atomic coordination structure of the mixtures.

Introduction: Thermal gradients induce the transport of colloids in suspension (thermophoresis) and concentration gradients in liquid mixtures and solutions (thermodiffusion) [1][2][3][4][5][6][7]. Thermodiffusion was discovered in the 19 th century: first observed by Ludwig [8], and later systematically investigated by Soret [9]. The Soret coefficient, S T , measures the mass separation of mixtures and suspensions in a thermal field. It is becoming a central property to characterise the non-equilibrium response of soft matter to thermal fields, for example, to quantify thermal transport in micro-scale devices that analyse the interactions of proteins and small molecules in biological liquids. [10] Experimental and computational studies have advanced significantly in recent years, but several outstanding questions remain. One such question is the microscopic origin of the forces driving the phenomenology observed in thermodiffusion measurements. Aqueous solutions feature a particularly rich phenomenology. For example, experiments and simulations report changes in the thermophilicity of alkali halide salt solutions with temperature. The solutions are thermophobic (S T > 0) at high temperatures and thermophilic (S T < 0) at low temperatures [11][12][13]. Interestingly, Gaeta et al. [12] reported minima in the Soret coefficients of NaCl (aq) and KCl (aq) as a function of concentration. These experiments were performed with thermogravitational columns, and the minima could not be reproduced using state-ofthe-art thermal diffusion forced Rayleigh scattering techniques, which circumvent convection effects [13]. Hence the S T minimum has remained controversial for many years. However, this situation has changed with recent experiments and computer simulations of LiCl (aq) , which support the existence of a minimum in the Soret coefficient [14,15]. Very recently, minima at high concentrations (∼ 2 M) were observed for thiocyanate (NaSCN (aq) * o.gittus18@imperial.ac.uk † f.bresme@imperial.ac.uk and KSCN (aq) ) and acetate (CH 3 COOK (aq) ) salt solutions [16,17], giving further impetus to the investigation of the physical origin of the Soret coefficient minima. We note that understanding the physical variables controlling the minimum is potentially relevant to explain the thermoelectric response of aqueous solutions as well (see ref. [15]).
In addition to aqueous electrolyte solutions, minima in the Soret coefficient with composition were observed in mixtures of polar fluids: ethanol/water [5,18,19], dimethyl-sulfoxide/water [20] and acetone/water [20,21]. In all these systems (electrolyte solutions and polar fluid mixtures), one of the components is water. This observation might suggest that the Soret coefficient minima are interlinked with water as a solvent and, therefore, its specific thermal transport properties. Indeed, molecular simulations of atomistic (non-polar) Lennard-Jones (LJ) binary mixtures at supercritical conditions do not offer evidence for the existence of minima in the Soret coefficient with composition [22]. However, some experiments of non-polar or weakly polar liquid mixtures reported maxima/minima in the Soret coefficient (e.g. cyclohexane/cis-decaline) and in some cases, accounting for an extrapolation to infinite dilution, a weak extrema can be inferred (e.g. toluene/1,3-dichlorobenzene) [23]. That work made no attempt to explain the microscopic origin of the extrema in S T , but crucially highlights the importance of the thermodynamic factor, Γ, a key quantity determining the heat of transport [15,24].
To investigate the existence of Soret coefficient minima in non-polar mixtures, and to probe the microscopic origin of such minima, we have performed computer simulations of the simplest liquid binary mixture, modelled with the Lennard-Jones model, which accounts for dispersion interactions. Advancing the discussion below, we show for the first time that a minimum in the Soret coefficient at a specific composition, and constant temperature and pressure, can be observed in simple non-polar liquid mixtures, hence showing that the minimum in the Soret coefficient is a completely general physical phenomenon.
We demonstrate that the minimum in the systems investigated here emerges from a balance of the mutual diffusion coefficient and the thermal diffusion coefficient, which vary in the opposite direction upon changing the composition of the mixtures. Furthermore, we correlate the minimum to the thermodynamic factor and changes in the solvation structure of the particles. Finally, we show that existing theories of thermodiffusion either do not reproduce the minimum and/or predict inaccurate S T values. Methods: We investigate the non-equilibrium behaviour of liquid-liquid binary Lennard-Jones mixtures at constant temperature T , pressure P , and different mole fractions x 1 = N 1 /(N 1 + N 2 ), where N i is the number of particles of type i. Inter-particle interactions were modelled using the LJTS potential, which is the LJ potential V LJ ij (r) = 4 ij [(σ ij /r) 12 − (σ ij /r) 6 ] truncated and shifted at a cutoff radius of r c = 2.5σ, V LJTS ij (r) = (V LJ ij (r) − V LJ ij (r c ))θ(r c − r) with θ being the Heaviside step function. All the particles have the same diameter and mass, but the interactions between of particles of type "1" or "2" are different. The parameters = 22 and σ = σ 11 = σ 22 together with the mass of each particle m 1 = m 2 = m define the usual LJ units. The energy scale was defined in terms of the high-boiling component, with 11 / 22 = 0.6. The Lorentz-Berthelot combining rules were used for the inter-species interaction parameters: 12 = √ 11 22 = √ 0.6 and σ 12 = (σ 11 +σ 22 )/2 = σ. Experiments of liquid-liquid mixtures have reported extrema in S T , therefore we target subcritical conditions for the LJ mixtures. Liquid-liquid mixtures were modelled at temperature T = 0.62 k −1 B (k B is the Boltzmann constant) and pressure P = 0.46 σ −3 , below the critical point [25] for both species. The entire 0 ≤ x 1 ≤ 1 composition domain is therefore expected to be subcritical. We note that along this isobar-isotherm, the thermodynamically stable phase for x 1 = 0 is a solid and we estimate the liquidus to be at x 1 ≈ 0.2 (see sec. 1.1 in the SI for further discussion). Simulations below this mole fraction correspond to a metastable liquid-liquid mixture. Nevertheless, as we will show below, the observed S T minimum is safely within the liquid-liquid mixture portion of the phase diagram.
We performed a variety of equilibrium molecular simulations (EMS), molecular dynamics (MD) and Monte Carlo methods, as well as non-equilibrium molecular dynamics (NEMD) simulations to calculate S T and related quantities. From NEMD, S T was evaluated at the stationary state characterised by zero net mass flux (J 1 = 0) as where x i and w i are the mole and mass fractions of species i (related by w i = x i m i /(x 1 m 1 + x 2 m 2 ) where m i is the mass of species i = 1, 2). EMS methods calcu-late, from [26] where D 12 is the mutual diffusion coefficient and D T is the thermal diffusion coefficient. Onsager's phenomenological coefficients L αβ were calculated from MD simulations in the NVE ensemble, using the Green-Kubo (GK) integral formulas and taking into account the enthalpy terms for the primed coefficient L 1q (see the SI). The chemical potential µ 1 , and subsequently the specific chemical potential µ s,1 = µ 1 /m 1 , was calculated from MD simulations using a free energy perturbation (FEP) method at constant T and P . (∂µ 1 /∂x 1 ) P,T was then calculated from the numerical derivative of µ 1 .
(∂µ 1 /∂x 1 ) P,T was also obtained from Kirkwood-Buff solution theory via two different methods of evaluating the Kirkwood-Buff integrals (KBIs): (1) from particle number fluctuations in grand canonical Monte Carlo (GCMC) simulations and (2) from the extrapolation to infinite system size of finite-volume KBIs [27], which were in turn calculated using radial distribution functions (RDFs) from MD simulations in the NVT ensemble. Thus, three EMS methods were used to calculate (∂µ 1 /∂x 1 ) P,T and subsequently the thermodynamic factor: FEP, KBI(RDF) and KBI(GCMC). Combining these with the GK calculations for the phenomenological coefficients give three corresponding "equilibrium" routes to S T , GK+FEP, GK+KBI(RDF) and GK+KBI(GCMC), in addition to NEMD simulations. In addition, other thermophysical properties such as the equation of state, diffusion coefficients and viscosities were calculated from equilibrium-MD simulations in the NVE, NVT and NPT ensembles. Details about these simulations and calculations are given in the SI. All simulations were performed using the software package LAMMPS [28] (v. 3 March 2020).
Results & discussion: Fig. 1(a) contains the main result of this letter: S T features a minimum as a function of composition at x min(S T ) 1 ∼ 0.5. Fitting cubic functions to the data give x min(S T ) 1 = 0.5 ± 0.1 for all four methods: NEMD, GK+FEP, GK+KBI(RDF) and GK+KBI(GCMC). The thermodiffusion response at the minimum is significantly enhanced with respect to diluted mixtures, by ∼30-40% relative to x 1 = 0.1, 0.9, and by ∼60-80% when compared to the extrapolated value of S T at x 1 = 1. For all compositions, S T < 0 which indicates that species 1 (the low-boiling component) is thermophilic and preferentially collects in the hot region. The S T values calculated from NEMD and all three EMS methods are in excellent agreement: all values agree within their statistical uncertainties. Accepting that the Green-Kubo integrals are sufficiently well converged (as shown in sec. 2.4 of the SI) the agreement indicates that the NEMD results are within the linear regime. This is expected since ∇T T −1 σ 0.007 1. Further verification of linear response is provided in sec. 2.2.2 of the SI for two mole fractions (x 1 = 0.5, 0.9). S T is determined by D 12 and D T (see eq. 2), which monotonically increase and decrease with x 1 , respectively ( fig. 1(b)). Thus, the S T minimum arises from a balance of D T and D 12 , as opposed to being carried through only by one of the transport coefficients.
The Soret coefficient S T can be written in terms of the phenomenological coefficients and the thermodynamic factor Γ, as [26] where m 1 = m 2 ⇔ x 1 = w 1 . The analysis of the differ-ent contributions to the RHS of eq. 5 offers microscopic insight into the mechanisms determining the minimum in S T . We find that the ratio L 1q /L 11 is essentially constant for all compositions ( fig. 1(d)). This suggests that the minimum arises from the w 1 (∂µ 1 /∂w 1 ) P,T term. As shown in fig. 1(c), Γ features a distinctive minimum at x min(Γ) 1 ∼ 0.5. All three EMS methods (FEP, KBI(RDF) and KBI(GCMC)) predict x min(Γ) 1 values in good agreement with each other. The minimum in S T is connected to the minimum in the thermodynamic factor, and therefore the minimum signals the composition at which the mixture features the largest non-ideality, max |1 − Γ|. We note that the S T and Γ minima coincide nearly exactly, due to the weak composition dependence of L 1q /L 11 .
The results presented above indicate that the non-ideal contribution dominates at conditions near the minimum. This can be visualised by splitting Γ into its ideal (id) and excess (ex) parts (see sec. 1.3 in the SI), showing as expected that Γ ex , is responsible for the minimum in Γ. To gain insight into the microscopic origins of the minimum in Γ and S T , we turn to Kirkwood-Buff theory, which connects Γ to the structural properties of the binary mixture, where ρ N is the total number density of the mixture. The Kirkwood-Buff integral (KBI) G ij is defined as the spatial integral over g µV T ij (r) − 1, and quantifies the excess (or deficiency) of species j around i. g µV T ij (r) is the pair correlation function in the grand canonical ensemble. The KBIs can be expressed in terms of the excess coordination numbers n ex ij (r ) = n ij (r ) − n id ij (r ) where n ij is the total coordination number, obtained from an integral over the pair correlation function, and n id ij = (4π/3)ρ N,j r 3 c is the ideal correlation number. ρ N,j is the number density of species j. Hence, the KBIs are given by ρ N,j G ij = lim r →∞ n ex ij (r ) = n ex,∞ we take the first-order approximation to the thermodynamic factor, Γ (1) = 1 − f . As shown in fig. 2(c), Γ (1) results in underestimations of 1-35% across the range of compositions, with larger errors for more non-ideal mixtures, but nevertheless provides insight into the relative importance of the f ij terms. f 12 features a maximum at x 1 ∼ 0.5 ( fig. 2(b)), indicating that the cross-species contribution is responsible for the minimum in Γ (1) and Γ. |(f 11 + f 22 )/f 12 | = 0-0.41 making the cross-species contribution much more significant, and ∼3-4 times larger in the region of the Γ (1) and Γ minima. Consequently, the phenomenology of Γ (1) and Γ are primarily determined by f 12 . Thus, the composition dependence of n ex,∞ 12 ( fig. 2(a)), which represents a net depletion of species 2 around 1 relative to the ideal state, and increases monotonically with x 1 , is the primary microscopic origin of the minimum in Γ and therefore S T .
We have shown above that the Soret coefficient of a simple binary mixture of non-polar liquids features a minimum as a function of composition, and linked this minimum to non-ideal effects, which result in a minimum in the thermodynamic factor and distinctive changes in the structural properties of the mixture, as reflected in the KBIs. Now we examine the accuracy of existing theoretical approaches to predict the S T minimum reported above. We note that existing theoretical models do not accurately predict S T in general. [29][30][31][32][33][34][35] In some cases, even the sign of S T is not predicted correctly. [30][31][32][33][34][35] Especially in earlier works, the discrepancies can, at least in part, be attributed to inaccuracies in experimentally determined properties. For example the Haase [32,36] and Kempers [32] models are very sensitive to partial molar properties, and therefore the equation of state used. [32,33] In computer simulations all the required quantities can be accurately calculated, and as shown here using an exact model, all the theories examined herein feature noticeable deviations from the S T values obtained by direct simulation. Previous simulations have shown that the theories are accurate only in a very limited number of cases, even for simple LJ mixtures and hard-sphere mixtures. [22,34,35] We have tested the standard theoretical models against our simulation data. We calculate the Soret coefficient according to the models of Haase [32,36]  ∼400-500% and ∼300-400% respectively. The Artola-Rousseau-Galliéro/Prigogine model underestimates |S T | by ∼100-110%, predicting values ∼ 10 −1 k B −1 . Furthermore, the model predicts the wrong sign: S ARG/P T > 0 or straddles 0 when accounting for the associated uncertainties.
S H T , S K T and S SF T all possess a minimum because they contain x 1 (∂µ 1 /∂x 1 ) P,T = k B T Γ in the denominator (see the SI). The x 1 (∂µ 1 /∂x 1 ) P,T term in S SF T originates directly from the phenomenological equations for thermodiffusion from linear non-equilibrium thermodynamics (LNET), which the model uses as a starting point for its derivation. For S K T , the x 1 (∂µ 1 /∂x 1 ) P,T term arises naturally from the statistical thermodynamics approach employed by Kempers. Originally an educated guess, S H T can be derived more rigorously within the framework of Kempers [32]. In contrast, S ARG/P T does not predict a minimum: it features a weak concentration dependence, and generally decreases with increasing x 1 . S ARG/P T contains k B T (as opposed to k B T Γ) in the denominator, which is only valid for ideal mixtures. Indeed the Artola-Rousseau-Galliéro/Prigogine model does not explicitly consider a concentration gradient along the reaction coordinate; doing so would result in a similar "non-ideality" term in the denominator. [39] Closing remarks: We close this letter with a discus-sion of our results in the context of recent literature. While the S T minimum examined in this work originates from a coincident minimum in Γ, for LiCl solutions, Γ increases monotonically in the concentration range of the S T minimum. [15] Furthermore, simulations have highlighted the importance of ion solvation structure on the existence of the minimum in LiCl (aq) [40]. In contrast, local structural changes in the LJ mixture are minor (see sec. 1.5 in the SI). We note that many simple salts are thought to have a D 12 minimum at low concentrations (∼ 10 −1 mol dm −3 ) often attributed to ion-pair formation and solute-solvent association. [16,[41][42][43] Recent experiments observed S T minima in NaSCN (aq) , KSCN (aq) , K 2 CO 3(aq) and CH 3 COOK (aq) that are carried through D T . [16,17] S T minima have been observed in mixtures of polar organic solvents with water. Ethanol/water and acetone/water mixtures have S T minima at x 1 ∼ 0.6 and x 1 ∼ 0.5, respectively, roughly coincident with minima in D 12 and Γ for both mixtures. [5, 18-21, 44, 45] Hence, the physical origin of the minima in aqueous solutions and mixtures might be quite different to the one reported here for non-polar liquids, since the LJ mixtures do not feature extrema in either D 12 or D T , which instead change monotonically with composition (see fig. 1(b)).
LJ mixtures are more representative of mixtures of non-polar organic solvents, in which intermolecular interactions are governed by van der Waals forces. At 25°C, cyclohexane/cis-decalin possesses a S T maximum at x 1 ≈ 0.2 but Γ ≈ 1.0 for the entire 0 ≤ x 1 ≤ 1 range (there is a very shallow Γ minimum at x 1 ≈ 0.8), indicative of a different origin compared to the LJ mixture. Also in contrast with the LJ mixture, cyclohexane/benzene features a monotonic increase of S T with x 1 at 25°C, but does have a strong minimum in Γ at x 1 ≈ 0.5. [23] It is evident that the physical origins of S T minima/maxima must be considered on a case-by-case basis for different mixtures, even those belonging to the same class of mixture (aqueous solution, non-polar organic solvent mixtures, etc.). We provide here a proof of principle and demonstrate that S T minima can exist in even the simplest of mixtures, as exemplified by a binary LJ mixture in which the components differ by only the interaction parameter . The concentration dependence of the S T is, at least for simple mixtures, typically attributed to the cross-interactions between unlike particles [22] -a notion sustained in recent reviews [7,46,47]. For dense supercritical LJ mixtures with cross-interactions given by 12 = k 12 √ 11 22 it was found that S T (x 1 ) ≈ bx 1 +c, with the slope b controlled by k 12 . [22] Greater |1 − k 12 | values resulted in greater |b| values. In this work, we identify a mixture with k 12 = 1 that nevertheless features a strong composition dependence and more complex phenomenology (the S T minimum). Clearly, further work is required to explain the composition dependence of S T in liquid mixtures and different thermodynamic conditions.

Supporting Information
On the microscopic origin of Soret   For an ideal solution the phase diagram can be estimated using the equations 2,3 ln x s,1 where R is the gas constant; T m,i and ∆H fus,i are the melting temperature and enthalpy of fusion, respectively, of species i; and x l,1 and x s,1 are the mole fractions of the liquidus and solidus for end-member 1. Eqs. 1&2 further assume that there is no difference in isobaric heat capacity between the supercooled liquid and pure solid of species 1 at T . 3 The predicted phase diagram is shown in fig. 1 The mixture is nearideal at x 1 ≈ 0.21 with thermodynamic factor Γ ≈ 0.8, and the assumption of ideal mixing is therefore expected to be a good approximation. The S T and Γ minima at x 1 ∼ 0.5 (see the main text) are safely within the liquid-liquid mixture region of the phase diagram. We split Γ into its excess and ideal gas contributions, Γ ex and Γ id respectively ( fig. 2).

Ideal and excess contributions to the thermodynamic factor
, the total number density ρ N and its derivative, while Γ ex = x 1 (∂µ ex 1 /∂x 1 ) P,T /(k B T ) depends on the excess chemical potential µ ex 1 and therefore the inter-particle interactions in the system. Γ id monotonically decreases with increasing x 1 , while Γ ex possesses a minimum at x min(Γ ex ) 1 = 0.4 ± 0.1 as determined by fitting a cubic function. Thus, the minimum arises through Γ ex , although its position is shifted by ∼ 0.1 (to x 1 ∼ 0.5) due to the Γ id contribution. 4

Theoretical models
We calculate the Soret coefficient according to the models of Haase 4,5 (S H T ), Kempers 5 (S K T ) Shukla and Firoozabadi 6 (S SF T ), and Artola, Rousseau and Galliéro 7 (S ARG T ): where h i , u i and v i are the partial molar enthalpy, internal energy and volume of species i = 1, 2. R is the gas constant. The Kempers model was derived by considering a non-isothermal two-bulb system, with the main assumption that the stationary state is the macroscopic state with the maximum number of microstates. Within this model, S K T corresponds to the centreof-volume frame of reference, while Haase's earlier educated guess S H T can be derived in the centre-of-mass frame. 5 h 0 i and S 0 T correspond to an ideal gas state at the same temperature (calculated from kinetic theory), and capture the kinetic contribution to S T . The Shukla- In these two kinetic models, thermal diffusion is described as a coupled diffusion-activated process, for which the elementary process can be summarized as a swap between two particles of different species along a temperature gradient. The activation fig. 3), where D 0 i is a constant and the self-diffusion coefficients D i have been corrected for finite-size effects (see sec. 2.6 below). The activation energy transported in molecular motion is better identified with self-diffusion than viscous flow; the use of the latter was historically motivated by the scarcity of self-diffusion data, even for pure components, and justified by the similar activation energies expected from Eyring's rate theory applied to liquids. 10 By using self-diffusion data, an approximation to ∆U c ). While permissible for some mixtures, this approach is unsuitable for others and has been criticised. 11,13,14 In this work, we use the method originally

Thermal conductivity
It is well established that the Soret effect reduces the thermal conductivity of the mixture. 17 From LNET, the thermal conductivity λ is given in terms of the phenomenological coefficients L αβ and L αβ as where the enthalpic terms in the primed coefficients cancel exactly to give the same thermal conductivity as for the unprimed coefficients. λ can therefore be split into two contributions: (1) the thermal conductivity in the absence of coupling effects, λ 0 = L qq /T 2 or λ 0 = L qq /T 2 , and (2) the mass-heat coupling term, δλ = −L 1q L q1 /(L 11 T 2 ) or δλ = −L 1q L q1 /(L 11 T 2 ), such that λ = λ 0 + δλ = λ 0 + δλ . The Soret effect should always decrease λ, and thus We show in fig. 6 the thermal conductivity of the mixture, and the contribution due to heat-mass coupling. The thermal conductivities calculated via NEMD and equilibrium-MD

Simulation details
The Soret coefficient S T and thermal conductivity λ were calculated from boundary-driven non-equilibrium molecular dynamics simulations (NEMD) in the stationary state. An elongated (tetragonal) simulation cell of dimensions (L x , L y , L z ) = (20, 20, 30)σ was used, with 3D periodic boundary conditions. Two thermostatting regions, hot and cold, were located in the centre and edges of the simulation, respectively (see fig. 7). The thermostatting regions had a width ∆z = 3σ and extended over the entire (x, y) plane, such that the temperature gradients were generated along the z-direction. For the thermostatting procedure, a simple velocity rescaling algorithm was used to maintain the hot and cold thermostatting regions at temperatures T h and T c , respectively. The velocities of all particles in each region were rescaled, every timestep, by a factor α = K t /K c where K t and K c are the target and current kinetic energies of the region. This velocity rescaling procedure does not conserve linear momentum, so the system's centre-of-mass velocity was subtracted from each particle at every time step in order to ensure linear momentum conservation. A timestep of δt = 0.002 τ was used.  For each system, the average density, mole fraction and thermostat temperatures were chosen to give P zz within ±0.001 σ −3 of the target pressure P = 0.46 σ −3 , where P zz is the pressure tensor component parallel to the heat flux vector. Relatively small (by simulation standards) temperature differences ∆T = T h − T c ≈ 0.04 k −1 B and resulting gradients ∇T ≈ 2.9×10 −3 k −1 B σ −1 were used to accurately target the thermodynamic state. For each system, 10 statistically independent replicas were generated, each consisting of an initial 2 × 10 4 τ to establish the stationary state, followed by 2 × 10 6 τ for data collection.
In the stationary state, this set-up results in two equal but opposite temperature gradients, and therefore in equal and opposite heat fluxes, such that the system is completely periodic.
The heat flux across the system, J q = (0, 0, ±J q ), can be obtained from the continuity where A = L x ×L y is the cross-sectional area of the simulation cell, δt is the timestep, and ∆U is the internal energy exchanged at each timestep. The employed simple rescaling thermostat only changes the kinetic energy such that ∆U = ∆K. The factor of 2 in the denominator accounts for the two heat fluxes (equal magnitude and opposite direction) generated in this setup. The thermal conductivity λ was then calculated using Fourier's Law: where ∇T is the local temperature gradient. The Soret coefficient S T was calculated using eq. 1 in the main text, again employing local values of the gradients ∇T and ∇x 1 . Local densities ρ, mol fractions x 1 and temperatures T were determined from a 1σ bin close to the centre of each NVE compartment. The position of the bin was chosen such that T = 0.62 k −1 B . The local values of ∇T and ∇w 1 were determined by fitting straight lines to the temperature and mole fraction profiles within a range of ±2.5σ around the selected state point.

Linear response
We demonstrate in fig. 8 that the magnitude of temperature gradients used in this work are within the linear regime. For these additional simulations, production runs of length 1.

Finite-size effects
We show in fig. 9 that the lateral simulation cell length L ⊥ = L x = L y has a significant effect on S T . While all but one (L ⊥ /σ = 10, 20 for We do not observe appreciable finite-size effects for the thermal conductivity. This is consistent with the mechanism of thermal transport in molecular liquids (and liquid mixtures), which is dominated by collisions between nearest neighbors, setting a characteristic length scale for heat transport at ∼ 1σ.

Partial molar properties
Partial molar properties z i were calculated using the equations z 1 = Z +(1−x 1 )(∂Z/∂x 1 ) P T N2 and z 2 = Z − x 1 (∂Z/∂x 1 ) P T N2 , where Z is the corresponding extensive property. Z(x 1 ) were calculated from MD simulations in the N P T ensemble. At each selected x 1 , two additional simulations were performed at x 1 ± 0.01, and (∂Z/∂x 1 ) P T N2 was then calculated by fitting a straight line through these three points.

Onsager's phenomenological coefficients
The phenomenological coefficients L αβ were calculated using the Green-Kubo integral formula where V is the volume of the simulation cell, and the factor of 3 in the denominator averages the contributions from each spatial dimension. In order to calculate L qq , L 1q , L q1 and L 11 , the expressions for the heat flux J q and mass flux J 1 in terms of microscopic quantities are required. These are: where the sum runs over all N 1 particles of species 1, and in the case of two-body interactions the Irving-Kirkwood formula for heat flux is where U i is the per-particle internal energy which can be split into potential and kinetic energy K i = (1/2)m i v 2 i contributions; S i is the per-particle stress tensor; F ij is the force exerted on particle i by particle j; and r ij = r j − r i where r i is the position vector of particle i.
The unprimed coefficients were calculated from simulations in the NVE ensemble. A correlation time of t c = 10τ was used for the upper limit of eq. 10. Selecting this integration limit is a compromise between sampling efficiency and minimizing the resulting truncation dotted lines show J q (t) · J 1 (0) and L q1 error. We show in fig. 11 that 10 τ is sufficiently long to achieve a well-converged integral, while the exhaustive extent of our sampling is reflected in the associated uncertainty.
The primed coefficients L αβ were then obtained from L αβ using the formulas where h s,i is the specific molar enthalpy of component i = 1, 2. These formulas can be derived by comparing expressions for entropy production in LNET, or trivially from the Green-Kubo 16 integral formula (eq. 10) noting that J q = J q − (h s,1 − h s,2 )J 1 . 19,23 The difference between J q and J q corresponds to the heat transported due to diffusion.
We show in fig. 12 the values of coefficients L αβ and L αβ . Consistent with Onsager's reciprocal relations, L 1q = L q1 and L 1q = L q1 to within their associated uncertainties. We therefore average over both reciprocal coefficients to give a single value each for L 1q = L q1

Chemical potentials
The chemical potential of species i was calculated in terms of its ideal (id) and excess (ex) contributions, µ i = µ id i + µ ex i . The ideal gas contribution µ id i is given by the formula where ρ N,i = x i ρ N and Λ i are the number density and thermal de Broglie wavelength of species i, respectively. The densities were obtained from MD simulations in the N P T ensemble. By convention we set Λ 2 = h/ √ 2πm 2 k B T = 1 ⇒ ln(Λ 3 2 ) = 0 (note that Λ 1 = Λ 2 since m 1 = m 2 ) for T = 0.62 k −1 B , and in doing so fix the energy scale by assigning a value to Planck's constant h in Lennard-Jones units. This amounts to a constant shift in µ i and does not affect (∂µ i /∂x 1 ) P,T .
The excess chemical potential µ ex i was calculated from MD simulations at constant T and P using a free energy perturbation (FEP) method. A single particle was inserted into the simulation cell via (n − 1) small "steps" along a reversible alchemical thermodynamic path, analogous to slowly "growing" the particle. The Gibbs free energy change for this particle insertion, ∆ N +1 N G ex , is given by where V is the volume of the simulation cell, ∆ ζi+1 ζi V(ζ, r) = V ζi+1 (ζ, r) − V ζi (ζ, r) and ζ is a coupling parameter that connects the reference (N ) and perturbed (N + 1) systems according to V(ζ, r) = ζV N +1 (r) + (1 − ζ)V N (r), with ζ taking values from 0 to 1. In order to avoid singularities when ζ → 0, the FEP simulations were performed using a soft-core version 24 of the Lennard-Jones potential (LJSC) truncated and shifted at a cutoff radius of r c = 2.5σ such that V LJSCTS ij (r c ))θ(r c − r) with θ being the Heaviside step function. Values of n = 1 and scaling constant α = 0.5 were used. For ζ = 1, eq. 19 reduces to the standard Lennard-Jones potential, and in the limit ζ → 0 no work is required to change from α = 0 to 0 < α < ∞ (i.e. the two potentials have equivalent initial states). Therefore, the soft-core version gives the same free energy difference as the standard LJ (or LJTS) potential. 24 Each FEP simulation used a cubic simulation cell containing 5000 particles, and additionally the single particle being inserted. The systems were equilibrated for 20τ in the NVT ensemble, then for 2 × 10 3 τ in the NPT ensemble. The particle insertion was then performed in steps of 0.01ζ over a 10 4 τ production run. A timestep of 0.002τ was used. Temperature (pressure) was controlled using the Nosé-Hoover chain thermostat (barostat), with 3 chains, and a time constant of 1τ (5τ ). Sampling consisted of 1700-2000 replicas for µ ex 1 , and 500-600 replicas for µ ex 2 . The chemical potentials and the derivative (∂µ 1 /∂x 1 ) P,T = (∂µ id 1 /∂x 1 ) P,T +(∂µ ex 1 /∂x 1 ) P,T are shown in fig. 13. (∂µ id 1 /∂x 1 ) P,T and (∂µ ex 1 /∂x 1 ) P,T were evaluated at x 1 by fitting a straight line through the three points within x 1 ± 0.01 and x 1 ± 0.05, respectively.

Self-diffusion coefficients and shear viscosity
The self-diffusion coefficients D i of species i = 1, 2 were calculated from the average mean square displacement (MSD) using the Einstein relation where r i is the position vector of a particle of species i, t is the elapsed time from arbitrary starting time t 0 , and d = 3 is the number of spatial dimensions. D i was therefore calculated by fitting to the equation |r i (t) − r i (0)| 2 = 6tD i , excluding the first 10 τ of data to ensure only the diffusive regime was sampled (as opposed to the ballistic regime). In order to account for finite-size effects, the "infinite-size" diffusion coefficient D i,0 was calculated by extrapolation to L −1 = 0, where L is the length of the cubic simulation cell. This finite-size analysis is shown in fig. 14, and follows from the equation derived by Yeh and Hummer 25 using a simple hydrodynamic model of a particle surrounded by a solvent of viscosity η YH in a periodically replicated simulation cell, where D i,PBC are the finite-size diffusion coefficients calculated from our MD simulations, and ξ is a dimensionless constant equal to 2.837297 for a cubic simulation cell with 3D periodic boundary conditions. 25 The same expression was obtained earlier 26 by Dünweg and Kremer using a closely related derivation. For each state point, the i = 1 and i = 2 data sets were simultaneously fit to eq. 21, giving a single viscosity η YH for the mixture.
The shear viscosity η was additionally calculated using the Green-Kubo integral formula where α = β such that P αβ are the off-diagonal elements of the pressure tensor. Results were averaged over (α, β) = (x, y), (x, z) and (y, z). A correlation time of t c = 10τ was used for the upper limit of the integral, which is sufficient to obtain well-converged integrals (see fig. 15).
In order to accurately target a specific thermodynamic states when calculating D i,PBC , the NVT ensemble was sampled using a temperature-control algorithm that alters dynamics compared to the Newtonian dynamics of the NVE ensemble. It is therefore necessary to check whether the employed global Nosé-Hoover (three chains) thermostat affects the computed self-diffusion coefficients. We show in fig. 14

Kirkwood-Buff integrals from grand canonical Monte Carlo
KBIs were calculated from grand canonical Monte Carlo simulations. For each simulation, the temperature of the ideal gas reservoir was set to T = 0.62 k −1 B , and the input chemical potentials µ 1 and µ 2 were determined from the free-energy perturbation simulations described in sec. 2.5. In all cases, a cubic simulation cell with length L = 20σ was used. Each MC step, 100 trial displacement and 100 trial exchanges (insertions or deletions with equal probability) were attempted. For the displacement moves, a maximum translation distance of 1.0σ was allowed. Trial moves were accepted/rejected using the standard Metropolis criterion. Each replica was first equilibrated in an NVT MD simulation at the target x 1 and ρ N for 200 τ , followed by 10 5 MC steps, and finally a production run of 10 7 MC steps. Sampling consisted of 400-500 statistically independent replicas for each state point.
The Kirkwood-Buff integrals G ij were calculated from the particle number fluctuations according to where N i is the number of particles of species i and V is the volume of the simulation cell.
The average pressures are within ±0.008 σ −3 of P = 0.46 σ −3 , corresponding to differences in number density on the order of δρ N ∼ 10 −4 σ −3 relative to those calculated from NPT MD simulations (see fig. 1(a)). Differences in G ij and subsequently Γ due to the slightly different thermodynamic states being sampled are expected to be insignificant compared to their associated statistical uncertainties.

Kirkwood-Buff integrals from molecular dynamics (NVT ) simulations
For an infinitely large and open three-dimensional system, the KBI for mixture components i and j is defined as where r is the radial distance between particles i and j. The infinite-size KBIs were calculated from a finite-size analysis of finite-volume KBIs defined for finite and open subvolumes embedded in a reservoir. 28 For spherical symmetry w(r) = 4πr 2 (1 − 3x/2 + x 3 /2) valid for x < 1, where x = r/(2R) and R is the radius of the subvolume. g µV T ij is the pair correlation function in the thermodynamic limit. It has been shown that G V ij scales linearly with R −1 and that the infinite-size G ij can be obtained by extrapolating the linear regime to R −1 = 0. [28][29][30] We calculate the RDFs from simulations of N = 10 5 to N = 5 × 10 6 particles in the NVT ensemble. Because fluctuations transform between ensembles, g µV T ij cannot be formally replaced with g N V T ij . However, for a sufficiently large system, the correlation lengths of particle density fluctuations are small compared to the linear dimension of the simulation cell, and local correlations are expected to be well reproduced. In addition to a large system size, The data shown corresponds to mole fraction x 1 = 0.5 and a system size of N = 5 × 10 6 particles. Symbols: circles (· · • · ·) and downward triangles (· · · ·) denote the use of pair correlation function g ij (uncorrected) and the corrected g GV ij , respectively, although these cannot be easily distinguished on the scale of the plots.
we apply the tail correction 30 of Ganguly and van der Vegt to ensure the correct asymptotic limit (lim r→∞ g ij = 1) of the RDF.
where N j is the total number of particles j in the simulation cell of volume V . However, we note that in our simulations, this correction does not have an appreciable effect on the linear regime and subsequently the extrapolated G ij (see fig. 16). Other corrections 31,32 to the RDFs for calculating KBIs have been proposed; the correction of Ganguly and van der Vegt was found to be the most accurate for a WCA system, relative to a larger reference system (L = 80σ). 29 We show in fig. 17 a finite-size analysis of the infinite-size KBIs. The final values for G ij in the main text were taken from the largest system size, and the associated uncertainties were estimated from the convergence with system size.

Melting points of the pure components
The melting temperatures T m of the pure components i = 1, 2 were determined using the direct coexistence method in the N P z T ensemble. This involves preparing a solid-liquid interface at a given (T, P ) state point and observing whether the system completely melts, freezes, or remains in coexistence. In MD simulations, cooling a homogeneous liquid below its melting point usually results in a metastable supercooled liquid; freezing is not observed unless prohibitively long simulations are performed. Likewise, heating a solid slightly above its melting point results in a superheated solid. The presence of the interface lowers the 24 kinetic barrier for melting/freezing.
First, the densities of the solid and liquid phases were determined from NPT simulations.
Initial configurations were prepared by joining two half-boxes, one of the FCC-solid the other of the liquid, each containing N = 5324 particles and equilibrated in the NVT ensemble. The combined N = 10648 system was then equilibrated for 20τ , ensuring that the solid phase did not melt by restraining each solid particle to its equilibrium position r 0 with a harmonic potential, V(r) = k(r − r 0 ) 2 /2 of force constant k = 10 σ −2 . A 4 × 10 4 τ production run in the N P z T ensemble was then performed, applying a barostat only to the direction parallel to the surface normal (the z-direction). A Nosé-Hoover chain barostat was used, with 3 chains, and a time constant of 4 τ . Likewise, a Nosé-Hoover chain thermostat was used, also with 3 chains, and a time constant of 1 τ . Sampling consisted of an initial 10 statistically independent replicas, followed by an additional 40 replicas if the initial set did not all either melt or freeze. In all cases, a homogeneous phase was observed by the end of the production run.
The simulations were performed in the N P z T ensemble, and as such the simulation cell dimensions (L x and L y ) were chosen to be consistent with the density of the solid at the (T, P ) state point. Otherwise, the solid phase would possess internal stress, 33,34 corresponding to a higher free energy, and resulting in the overestimation of the coexistence pressure and/or underestimation of the melting temperature.
Each replica was determined to have either frozen or melted using the Q 6 (l = 6) bondorientational order parameter. The Steinhardt order parameters Q l were introduced to characterize local orientational order in atomic structures, 35 and are given by Y lm (θ(r ij ), φ(r ij )) (29) for each particle i, where Y lm (θ, φ) are the spherical harmonics, θ and φ are the polar angles of "bond" vector r ij between i and neighbour j, and N n is the number of nearest neighbours to particle i. Q l is therefore a rotationally invariant non-negative amplitude. Q l adopt welldefined values for high-symmetry structures; for a perfect FCC crystal (F m3m) and N n = 12, Q 6 = 0.575. 36 We show in fig. 18 that the Q 6 order parameter can be used to distinguish between the FCC solid and liquid phases. The replica was determined to have frozen if Probability of freezing P f as a function of temperature T for species 1 (left) and 2 (right).
Solid lines show the fitted sigmoid-like functions (eq. 30).
Q 6 > 0.496, or melted if Q 6 < 0.375. The average Q 6 was taken over the last 10 τ of the trajectory.
In the thermodynamic limit, a solid melts at T > T m , and a liquid freezes for T < T m .
However, a finite-size system may stochastically melt or freeze with probabilities P m and P f = 1 − P m respectively. For each component, T m was determined by fitting P f (T ) to the sigmoid-like function ( fig. 18): where d controls the sharpness of the probability profile, with δ 10−90 = 2.178d being the width of the interval where P f goes from 0.1 to 0.9. The melting (coexistence) temperature is defined by P m (T m ) = P f (T m ) = 0.5. The melting temperatures are shown in table 1 alongside other coexistence properties (coexistence densities and enthalpies were determined by interpolating data from the NPT simulations).