Temperature gradients induce mass separation in mixtures in a process called thermal diffusion and are quantified by the Soret coefficient ST. Thermal diffusion in fluid mixtures has been interpreted recently in terms of the so-called (pseudo-)isotopic Soret effect but only considering the mass and moment of inertia differences of the molecules. We demonstrate that the first moment of the molecular mass distribution, the mass dipole, contributes significantly to the isotopic Soret effect. To probe this physical effect, we investigate fluid mixtures consisting of rigid linear molecules that differ only by the first moment of their mass distributions. We demonstrate that such mixtures have non-zero Soret coefficients in contrast with ST = 0 predicted by current formulations. For the isotopic mixtures investigated in this work, the dependence of ST on the mass dipole arises mainly through the thermal diffusion coefficient DT. In turn, DT is correlated with the dependence of the molecular librational modes on the mass dipole. We examine the interplay of the mass dipole and the moment of inertia in defining the isotopic Soret effect and propose empirical equations that include the mass dipole contribution.

Thermal diffusion is the process by which a temperature gradient induces a concentration gradient in a mixture. It was first observed1 in the 19th century by Ludwig and later systematically studied2 by Soret. It was not until the 1910s, however, that thermal diffusion was “theoretically discovered,” and one can speculate that it was overlooked by the original pioneers of kinetic theory because the Soret coefficient vanishes for Maxwellian molecules (point particles that repel each other with a force Fr−5, where r is the inter-particle distance).3 These equations of Enskog and Chapman already implied that components of a mixture could be separated based solely on their mass,3 and in 1919, Chapman was the first to suggest4 the use of thermal diffusion for separating isotopes. Isotope fractionation by thermal diffusion was first experimentally achieved in 1938 by Clusius and Dickel5 and subsequently received considerable theoretical treatment by Furry, Jones, and Onsager, among others.3,6

The pseudo-isotopic Soret effect refers to the contribution to the Soret coefficient ST from the mass and mass distribution of the particles. Using binary mixtures of isotopologues of benzene and cyclohexane, Debuschewitz and Köhler7 demonstrated that the Soret coefficient could be split into two additive contributions,
ST=STchem+STiso.
(1)
The chemical contribution STchem depends on the intermolecular and intramolecular interactions (i.e., the potential energy surface), while the isotopic contribution STiso depends only on the mass and mass distribution of the molecules. Simulations have verified that Eq. (1) holds reasonably well for binary Lennard-Jones mixtures in both the liquid and dense supercritical states even if slight couplings have been observed.8–12 The application of Eq. (1) to the interpretation of ST data for molecular mixtures is now standard procedure.13–16 
To date, all models (empirical and theoretically derived) represent STiso as a function of total masses of the particles Mi and their moments of inertia Ii. For binary mixtures,
STiso=STiso(M1,M2,I1,I2),
(2)
and for mixtures consisting of species with the same mass,
M1=M2=MSTiso=STiso(I1,I2;M),
(3)
where three assumptions are implicit in Eqs. (2) and (3): (1) the effect of the internal mass distribution can be modeled using only the moment of inertia; (2) rigid body moments of inertia are applicable to flexible molecules; and (3) rigid body moments of inertia Ii (3 × 3 matrices) can be reduced to scalars Ii. Assumptions (2) and (3) are justified for rigid linear molecules, such as diatomic molecules. Small polyatomic molecules have typically been characterized by a single principal moment of inertia Ii,k or Ii=3/kIi,k1.7,17–23 A recent study on halobenzene/n-alkane mixtures used STiso(M1,M2) because the moment of inertia is not well defined for flexible alkanes.18 To the best of our knowledge, assumption (1) has not been tested or discussed in the literature.
How best to characterize the internal mass distribution of a molecule? The moments μn of the mass distribution are a natural extension to using the total mass and moment of inertia. For a rigid body with mass density ρ(r),
μn=Vρ(r)rnd3r,
(4)
where rn = rr ⊗ ⋯ ⊗ r (n times) and ⊗ denotes the Kronecker product of two matrices. The zeroth moment μ0 = M is the total mass of the body (scalar). The first moment μ1 = D is the mass dipole (3 × 1 vector), which gives the center of mass when normalized by the total mass, D/M = rcom. The second moment μ2 corresponds to the moment of inertia I (3 × 3 matrix). Higher order moments, μn>2, do not affect rigid body dynamics, and thus, the mass distribution of a rigid body can be completely represented by μn=0,1,2, i.e., by M, rcom, and I. While M and I appear frequently in models of thermal diffusion, the impact of the mass dipole on the Soret coefficient of molecular mixtures has not been considered. Indeed, the mass dipole was introduced only very recently in the context of thermal diffusion when studying the thermal orientation (TO) and thermophoresis of anisotropic colloids.24–26 We show here that the mass dipole is part of a more general description of the mass distribution, which adds additional contributions to the formulation summarized by Eq. (2).

In this work, we investigate the effect of the mass dipole on the Soret coefficient of isotopic mixtures of rod-like molecules. The model we employ allows for the systematic change in the moments of the mass distribution, making it possible to isolate the M, D, and I contributions to the Soret coefficient. We provide a proof of principle for the mass dipole contribution in isotopic mixtures. However, the main conclusions of this work are applicable to non-isotopic mixtures as well, conveniently through Eq. (1).

We investigate the thermal diffusion of mixtures of rod-like molecules in which the components differ only by their internal mass distribution. Following our previous studies,24,25 the molecules were modeled using the shish-kebab model, a rigid chain of tangent spherical monomers of effective diameter σ, of length N = 7. The mass mi of monomer i is given by a point mass at its interaction site. All intermolecular interactions were modeled using the Weeks-Chandler-Andersen (WCA)53 potential,
VWCA(r)=4εσr12σr6+εif r21/6σ,0if r>21/6σ,
(5)
where r is the distance between two monomers and ɛ represents the strength of the purely repulsive pairwise interaction. These parameters, ɛ and σ, together with the unit of mass m are used to define the usual Lennard-Jones units, which are adopted in this work.
The total mass of each molecule was set to M = 56m. Noting that the total mass has been held constant, we report the reduced mass dipole d = D/M, which has a more intuitive physical meaning than D: it is the displacement of the center of mass from the geometric center, d=dûd=rcomrg. The geometric center rg was chosen as the center of symmetry of the molecule (i.e., the inversion center given the Dh point group). Since we consider freely rotating rigid bodies, their rotation is characterized by the moments of inertia I about the axes through the center of mass. The moments of mass for this model (Fig. 1) are given by (for the molecule oriented along the x axis),
M=i=17mi=56m,d=1Mi=17mi(xixg),I=i=17mi(xixcom)2,
(6)
where xi is the distance of monomer i along the molecular axis. There are no analytical solutions to Eq. (6) for a given set {N, M, d, I} without applying additional constraints. We use numerical optimization methods to generate the mass distributions that give the desired {N, M, d, I}. Each mixture is characterized by the mass dipoles di and moments of inertia Ii of species i = 1, 2. By convention, we assign species 1 according to I1 > I2, and d1 > d2 if I1 = I2.
FIG. 1.

The shish-kebab model of length N = 7.

FIG. 1.

The shish-kebab model of length N = 7.

Close modal

We performed both equilibrium and non-equilibrium molecular dynamics (NEMD) simulations of equimolar (mole fractions x1 = x2 = 0.5) binary mixtures targeting temperature T = 3.0 εkB1 and pressure P = 1.121 ɛσ−3, which correspond to a monomer number density of ρN,mon = 0.398 σ−3 (molecule number density ρN,m = ρN,mon/7 = 0.0568 σ−3) and volume fraction ϕ = 0.208. The molecules were modeled as rigid bodies; their translational and rotational degrees of freedom were integrated using the method of quaternions.27 A timestep of δt = 0.002τ was used for all simulations. All simulations were performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)28 (v. 3 March 2020).

1. Non-equilibrium molecular dynamics simulations

Boundary-driven NEMD simulations of the mixtures were performed at an average density of ρN,mon = 0.4σ−3. An elongated (tetragonal) simulation cell of dimensions (Lx, Ly, Lz) = (50, 50, 100)σ was used under 3D periodic boundary conditions. Two thermostatting regions, hot and cold, were located in the center and edges of the simulation, respectively (see Fig. 2). The thermostatting regions had a width Δz = 7σ and extended over the entire (x, y) plane such that the temperature gradients were generated along the z-direction. For the thermostatting procedure, a Langevin thermostat with a time constant (damping parameter) of 1τ was applied to the hot and cold regions at temperatures Th and Tc, respectively. The thermostat temperatures were set to Tc = 2.5 εkB1 and Th = 3.5 εkB1, unless specified otherwise. The Langevin thermostats do not conserve linear momentum, so the system’s center-of-mass velocity was subtracted from each particle at every timestep in order to ensure linear momentum conservation. In the stationary state, this setup results in two equal but opposite temperature gradients and therefore in equal and opposite heat fluxes such that the system is completely periodic. In order to avoid potential artifacts associated with the dynamics of the Langevin-thermostated particles,29,30 only data pertaining to molecules with a center of geometry at a minimum distance of (N/2)σ = 3.5σ from either thermostatting region were used for analysis. For each system, ten statistically independent replicas were generated, each consisting of an initial te ≥ 5 × 104 τ to establish the stationary state followed by production runs of equal length between 10–30 × 104 τ. Owing to the small signal/noise ratio, it is not possible to obtain well-converged mole fraction profiles in the transient regime preceding the stationary state. The time required to reach the stationary state is of the order of the characteristic time ts = l2/D, where l is the distance between the two thermostats and D is the self-diffusion coefficient. For our NEMD simulations, te/ts ∼2–5, using the smaller of the two self-diffusion coefficients (corresponding to species i = 1, 2) for each mixture.

FIG. 2.

Representative (a) temperature T and monomer number density ρN,mon profiles and (b) mole fraction, x1 and x2, profiles for the NEMD simulations. The blue (cold) and red (hot) indicate the location of the thermostatting regions in the simulation cell. The profiles correspond to the (I1 = I2, d1, d2) = (802, 2.7σ, 0) mixture.

FIG. 2.

Representative (a) temperature T and monomer number density ρN,mon profiles and (b) mole fraction, x1 and x2, profiles for the NEMD simulations. The blue (cold) and red (hot) indicate the location of the thermostatting regions in the simulation cell. The profiles correspond to the (I1 = I2, d1, d2) = (802, 2.7σ, 0) mixture.

Close modal

2. Equilibrium simulations

Equilibrium molecular dynamics (EMD) simulations of the mixtures were performed in the NVT ensemble. A cubic simulation cell of length L ≈ 56.05σ containing 104 molecules was used. The temperature was controlled by the Nosé–Hoover chain thermostat with three chains and a time constant of 10τ. Sampling consisted of 10–30 replicas for each system; each replica was equilibrated for 103 τ followed by 2 × 104 τ of production.

Simulations of the pure fluid (same particle model but with the mass distribution given by mi = m for all monomers) in the NPT ensemble were performed to determine the phase coexistence diagram. These simulations used cubic cells with 3 × 104 molecules. Temperature (pressure) was controlled using a Nosé-Hoover chain thermostat (barostat) with three chains and a time constant of 1τ (10τ). To ensure that hysteresis effects31 were not significant, two replicas were performed for each (P, T) state point; only systems where the S2 order parameter (see Sec. III A) of the two replicas agreed to within their statistical uncertainties were included. Along each isobar, the two replicas were generated using upward (downward) branches starting from a low-density isotropic (high-density nematic) system with gradually decreasing (increasing) temperature. Each replica consisted of 0.4–12 × 104 τ of equilibration starting from a nearby state point, followed by a 2 × 104 τ production run.

The Soret coefficient ST was evaluated from the NEMD simulations at the zero-mass-flux (J1 = 0) stationary state as
ST=1w1w2w1TJ1=0=1x1x2x1TJ1=0,
(7)
where xi and wi are the mole and mass fractions of component i. The local gradients ∇T and ∇x1 were determined by fitting a straight line to their profiles within a range of ±10σ around the selected state point at T=3.0εkB1. Verification of linear response is given in the supplementary material.
The mutual diffusion coefficient D12 was calculated from
D12=L11ρ(1w1)Tμs,1w1P,T,
(8)
where μs,1 is the specific chemical potential of component 1. For ideal mixtures,32 such as the isotopic mixtures considered here with the same intermolecular interactions between different components,
μs,1w1P,T=kBTw1[M1w1(M1M2)].
(9)
The phenomenological coefficients L11 were calculated using the Green–Kubo integral formula
L11=V3kBlimt0tJ1(t)J1(0)dt,
(10)
where V is the volume of the simulation cell and the factor of 3 in the denominator averages the contributions from each spatial dimension. The mass flux is given by
J1=1Vi=1N1j=17mj,ivj,i,
(11)
where v is the velocity and the sums run over all seven monomers of each of the N1 molecules of species 1. Equation (10) was evaluated from the equilibrium-NVT simulations using a correlation time of t′ = 500τ for the upper limit of the integral, which is sufficient for well-converged integrals (see the supplementary material).

Unlike in previous studies,33–35 the thermal diffusion coefficient DT could not be evaluated from equilibrium simulations using the Green–Kubo approach. Owing to the rigid body constraint forces, reliable per-particle virial stress tensor-based instantaneous heat fluxes could not be obtained even when using the recently derived centroid formulation36,37 (and using only the x component38,39). The formulation of instantaneous heat flux in rigid body dynamics has long been established.40 Further work is required to implement this capability in LAMMPS. Instead, DT was calculated from DT = STD12. This approach relies on the fact that finite-size effects associated with ST from NEMD and D12 from equilibrium-NVT simulations are small and additionally cancel. It has recently been shown that the hydrodynamic correction of Yeh and Hummer41 for self-diffusion coefficients can also be applied to D12.42 This correction amounts to an increase in D12, D1, and D2 of <4% (see the supplementary material), which is smaller than the uncertainties associated with D12, ST, and therefore DT. Regarding ST from NEMD, extrapolating the simulation cell lengths Lx = Ly → ∞ decreases ST by 4% (see the supplementary material). Thus, finite-size effects in ST and D12 are expected to at least partially cancel when calculating DT. This is because finite-size effects in STDT/D12 include those from D12; based on our finite-size analysis of D12 and ST, we expect finite-size effects of 1% in DT, which is smaller than its associated uncertainties. Thus, we use and report the uncorrected D12, D1, and D2 values, which give more internally consistent values for ST, D12, and DT. We note that using the infinite-size D12, D1, and D2 does not change any of the qualitative trends reported in this work (see the supplementary material).

Self-diffusion coefficients were calculated using the Einstein relation. Rotational diffusion coefficients and shear viscosities were calculated from their Green–Kubo equations using an upper integration limit of 200τ and 500τ, respectively. (See the supplementary material). As with D12, these transport properties were calculated from equilibrium-NVT simulations that use a global Nosé–Hoover chain thermostat. Global thermostats have been shown to give transport properties that are statistically indistinguishable from those of the NVE ensemble.30 Especially given the weak coupling strength (large time constant), the thermostat is not expected to significantly affect the correlation functions and transport properties calculated in this work.

All statistical uncertainties reported in this work refer to the 95% confidence interval of the mean, unless stated otherwise.

The simulated fluid exhibits an isotropic–nematic phase transition, and we are interested in the behavior of the mixtures in the isotropic phase. Choosing the aspect ratio of the molecules is a compromise: the maximum mass dipole achievable is determined by the length N of the shish-kebab molecule, but the nematic phase is shifted to lower volume fractions for greater aspect ratios. It is therefore necessary to map out the phase coexistence diagram to determine which thermodynamic conditions correspond to the isotropic phase.

The uniaxial orientational order parameter S2 is zero for an isotropic fluid and unity for a perfectly aligned system. S2 is the ensemble average of the largest positive eigenvalue P2 of the Q tensor,
Qαβ=1Nmi=1Nm32ûi,αûi,β12δαβ,α,β=x,y,z,
(12)
S2=P2(nû)=P2(cosθn)=32cosθn12,
(13)
where the director n is the eigenvector associated with P2 and û=ûd is the unit vector along the long molecular axis. Numerical data for the estimated isotropic–nematic coexistence temperatures and upper/lower bounds on the coexistence densities are given in the supplementary material.

The phase diagram of the pure fluid is shown in Fig. 3. In classical systems, ensemble averages over thermodynamic observables are independent of the mass and mass distribution of the particles, and the phase diagram is therefore applicable to the isotopic mixtures under equilibrium conditions. All the NEMD simulations correspond to thermodynamic conditions, T ≥ 2.0 εkB1 and ϕ < 0.24, that are safely within the isotropic region of the phase diagram. Each set of NEMD profiles corresponds to the equation of state (EOS) at the simulated pressure, as defined by the component (Pzz) parallel to the heat flux vector. The inset of Fig. 3 shows that NEMD simulations reproduce the equilibrium EOS at P = 1.121ɛσ−3 and therefore that local equilibrium is fulfilled.

FIG. 3.

Phase diagram and equation of state (EOS). T, ρN,mon, and ϕ are the temperature, monomer number density, and volume fraction, respectively. The data are color coded according to the uniaxial (orientational) order parameter S2. The dashed isotropic–nematic coexistence line is to guide the eye. The inset shows the EOS at pressure P = 1.121ɛσ−3, as predicted by equilibrium-NPT and NEMD simulations. Symbols: I = isotropic and N = nematic.

FIG. 3.

Phase diagram and equation of state (EOS). T, ρN,mon, and ϕ are the temperature, monomer number density, and volume fraction, respectively. The data are color coded according to the uniaxial (orientational) order parameter S2. The dashed isotropic–nematic coexistence line is to guide the eye. The inset shows the EOS at pressure P = 1.121ɛσ−3, as predicted by equilibrium-NPT and NEMD simulations. Symbols: I = isotropic and N = nematic.

Close modal

In this section, we consider mixtures with I = I1 = I2 but d1d2 for which all current models of thermal diffusion predict ST = 0 (see, e.g., the reviews Refs. 1316 and 43 and the references contained therein). Figures 4(a-i) and 4(b-i) contains the main result of this work: mixtures with components that differ only by their mass dipole have non-zero Soret coefficients. Indeed, the mole fraction profiles for one such mixture are shown in Fig. 2, providing clear evidence for species separation in the thermal field. In all cases, ST > 0, which indicates that the species with the greater mass dipole is thermophobic and preferentially collects in the cold region. Assigning an interaction strength typical of simple molecular fluids, ɛ/kB = 102 K, gives ST = 1–15 × 10−3 K−1 at T = 300 K (T = 3.0 εkB1) for the mass dipole contribution [Figs. 4(a-i) and 4(b-i)]. This is of comparable magnitude to the entire (pseudo-)isotopic Soret effect in many molecular mixtures.7,17,18

FIG. 4.

Transport coefficients for the isotopic mixtures with (a) I = I1 = I2 and d2 = 0 and (b) I = I1 = I2 = 802 and d1d2. The (i) Soret coefficient ST, (ii) thermal diffusion coefficient DT, and (iii) mutual diffusion coefficient D12 as a function of d1d2. di and Ii are the mass dipole and moment of inertia of species i = 1, 2. The solid lines show equations that model the data: see the main text for details.

FIG. 4.

Transport coefficients for the isotopic mixtures with (a) I = I1 = I2 and d2 = 0 and (b) I = I1 = I2 = 802 and d1d2. The (i) Soret coefficient ST, (ii) thermal diffusion coefficient DT, and (iii) mutual diffusion coefficient D12 as a function of d1d2. di and Ii are the mass dipole and moment of inertia of species i = 1, 2. The solid lines show equations that model the data: see the main text for details.

Close modal

We probe the phenomenology of STDT/D12 by considering the thermal diffusion coefficient DT and mutual diffusion coefficient D12. As shown in Fig. 4(a-iii), D12 decreases with increasing I for mixtures with I = I1 = I2 and d2 = 0 but is essentially constant with d1. Furthermore, the constant value of D12 determined by fitting to the D12(d1I = 802, d2 = 0) data accurately predicts the D12 values for mixtures with d2 ≠ 0 as well [Fig. 4(b-iii)], indicating that D12 does not significantly depend on either mass dipole, d1 or d2, at least relative to the statistical uncertainties associated with the data.

With regard to DT, all the data for the I = I1 = I2 and d2 = 0 mixtures can be accurately described by a power law dependence DT=a(d1d2)k [Fig. 4(a-ii)], indicating that DT is only weakly dependent on the moment of inertia. The exponent k=kd2=0=1.6±0.1 was determined by fitting the equation ln DT = k ln(d1d2) + ln a to the d2 = 0 data (see the supplementary material). Considering mixtures with d2 ≠ 0, fits to the equation DT=a(d1d2)kd2=0 accurately reproduces the simulation data [Fig. 4(b-ii)]. Thus, the effect of changing d2 while keeping d1d2 constant can be captured by the prefactor a, which increases with increasing d2, as shown explicitly in the inset of Fig. 8(b) in the supplementary material.

Finally, we return to ST, which we predict with the equation
ST[d1,d2,I=I1=I2]=DT[d1,d2]D12[I]=a(d1d2)kcI[I],
(14)
where k=kd2=0 describes the whole dataset and a and cI were determined by the fits to DT and D12, respectively. We show in Figs. 4(a-i) and 4(b-i) that Eq. (14) accurately predicts ST and captures its weak dependence on I, which can be discerned for large values of d1d2, as shown in the inset of Fig. 4(a-i).

We additionally calculate D12 using the Darken approximation DMSDMSDarken=x2D1+x1D2, where DMS is the Maxwell–Stefan diffusion coefficient related to the Fickian D12 via D12 = ΓDMS, where Γ is the thermodynamic factor (Γ = 1 for the ideal mixtures considered here). The Darken approximation D12Darken predicts values in excellent agreement with the (non-approximate) D12 values [Fig. 5(a)], indicating that displacement cross correlations between particles do not play an important role in determining D12. Thus, the behavior of D12 can be rationalized in terms of the self-diffusion coefficients D1 and D2. As shown in Fig. 5(b), D12Darken features a weak but discernible dependence on d1 (and therefore d1d2), which was masked by the greater uncertainties associated with D12 calculated from the Green–Kubo integral formula. D1 and D2 also depend weakly on the mass dipole d1 [Figs. 5(c) and 5(d)]. Unsurprisingly, the impact of d1 on D1 is greater than on D2 with ΔD1D2 ∼ 1–4, where ΔDi = Di(d1) − Di(d1 = 0.5σ). Thus, the weak dependence of D12 on the mass dipoles d1 and d2 can be primarily attributed to the effect of di on Di. (See the supplementary material for additional data on the diffusion coefficients and related quantities).

FIG. 5.

(a) The mutual diffusion coefficient D12 vs the Darken approximation D12Darken. (b) D12Darken, (c) ΔD1, and (d) ΔD2 as a function of the mass dipole, d1, of species 1. ΔDi = DiDi(d1 = 0.5σ) is the change in self-diffusion coefficient Di of species i = 1, 2 relative to the corresponding mixture with d1 = 0.5σ.

FIG. 5.

(a) The mutual diffusion coefficient D12 vs the Darken approximation D12Darken. (b) D12Darken, (c) ΔD1, and (d) ΔD2 as a function of the mass dipole, d1, of species 1. ΔDi = DiDi(d1 = 0.5σ) is the change in self-diffusion coefficient Di of species i = 1, 2 relative to the corresponding mixture with d1 = 0.5σ.

Close modal

The mass dipole also affects the rotational diffusion of the rod-like molecules. We show in Fig. 6(a) that the rotational diffusion coefficient Dr,1 decreases with increasing mass dipole d1. Dr,i is the rotational diffusion coefficient of species i = 1, 2. The rotational diffusion coefficient decreases with increasing moment of inertia, as expected. The translational and rotational diffusion coefficients are coupled with 15 > (Di/Dr,i)/(σ2 rad−2) > 13 for all the mixtures considered in this section. As shown in Fig. 6(b), the ratio D1/Dr,1 increases with the mass dipole d1 but is statistically independent of I within the statistical uncertainties of our data. We note that comparing to previous simulations44 of systems of the same WCA shish-kebab model but at the different temperature T = 1 εkB1, the (N = 7, ρN,m = 0.0568 σ−3) state point studied in this work is expected to be in the semidilute regime, where the translational and rotational diffusion are known to be coupled.45 

FIG. 6.

Rotational and translational diffusion coefficients and their correlation with the thermal diffusion coefficient DT. (a) Dr,1 as a function of the mass dipole d1 of species 1. (b) The ratio D1/Dr,1 as a function of d1. (c) DT vs −(Dr,1Dr,2). (d) DT vs −(D1D2). Dr,i and Di are the rotational and translational diffusion coefficients, respectively, of species i = 1, 2.

FIG. 6.

Rotational and translational diffusion coefficients and their correlation with the thermal diffusion coefficient DT. (a) Dr,1 as a function of the mass dipole d1 of species 1. (b) The ratio D1/Dr,1 as a function of d1. (c) DT vs −(Dr,1Dr,2). (d) DT vs −(D1D2). Dr,i and Di are the rotational and translational diffusion coefficients, respectively, of species i = 1, 2.

Close modal

We conclude from Fig. 4 that the Soret coefficient depends on the mass dipoles mainly through DT. We show in Figs. 6(c) and 6(d) that DT is positively correlated with the differences in rotational diffusion coefficients, −(Dr,1Dr,2), and translational diffusion coefficients, −(D1D2). We expect the rotational correlation reflects librational modes that influence the thermal transport and therefore DT.

Hence, to probe the microscopic mechanism associated with librational modes and its dependence on the molecular internal degrees of freedom, we calculate the power spectra of various velocity autocorrelation functions (VACFs). We consider the “all-atom” VACF, Catom, as well as the rotational and translational center-of-mass VACFs of the rigid molecules, Crot and Ctrans, respectively,
Catom(t)=vi(t)vi(0)vi2(0),
(15)
Crot(t)=ωj(t)ωj(0)ωj2(0),
(16)
Ctrans(t)=vCOM,j(t)vCOM,j(0)vCOM,j2(0),
(17)
where vi is the velocity of monomer i, ωj is the angular velocity of rigid molecule j, vCOM,j is the center-of-mass velocity of rigid molecule j, and t is the elapsed time from an arbitrary starting time. The power spectrum is then given by the Fourier transform of the corresponding VACF,
IX(ν)=limtttCX(t)ei2πνtdt,
(18)
where IX is the intensity and ν is the frequency.

The power spectra of species 1 in the I = I1 = I2 = 802 and d2 = 0 mixtures are shown in Fig. 7. For large mass dipoles, Iatom features a peak at 0.5τ1 indicative of a librational mode. Increasing the mass dipole increases the frequency of the librational mode. This librational mode is also reflected in Irot. For low mass dipoles, a single peak is observed at ∼0.2–0.3τ−1, which contains contributions from both the libration and a lower frequency mode. Upon increasing d, the librational mode is blue-shifted and two distinct peaks are observed.

FIG. 7.

Power spectra of species 1 of the (I1 = I2, d1, d2) = (802, d1, 0) mixtures. The intensity IX vs frequency ν for the (a) “all-atom” (Iatom) (b) rotational (Irot) and (c) translational center-of-mass (Itrans) power spectra. In (a) and (b), the vertical solid lines show the librational mode frequencies νlib. In (b), the dashed-dotted (-.) lines denote the fitted spectra Ifit and the inset shows the individual fitted functions for the (I1 = I2, d1, d2) = (802, 2.7σ, 0) mixture.

FIG. 7.

Power spectra of species 1 of the (I1 = I2, d1, d2) = (802, d1, 0) mixtures. The intensity IX vs frequency ν for the (a) “all-atom” (Iatom) (b) rotational (Irot) and (c) translational center-of-mass (Itrans) power spectra. In (a) and (b), the vertical solid lines show the librational mode frequencies νlib. In (b), the dashed-dotted (-.) lines denote the fitted spectra Ifit and the inset shows the individual fitted functions for the (I1 = I2, d1, d2) = (802, 2.7σ, 0) mixture.

Close modal
In order to disentangle contributions from the different modes, we fit the rotational power spectra Irot using the following equations:
Ifit=IL+IDHO,
(19)
IL=aL1+ννLcL2,
(20)
IDHO=Dγ(ν02ν2)2+γ2ν2,
(21)
where the Lorentzian-like function IL captures the low frequency peak at ν ∼ 10−2τ−1. In spectroscopy, librations are typically modeled as damped harmonic oscillators (DHO). In Eq. (21) for the DHO, D is the oscillator intensity constant, γ is the damping constant, and ν0 is its frequency (band position). Using the I = I1 = I2 = 802 and d2 = 0 mixtures as an example, we show in Fig. 7(b) that fitting Eq. (19) accurately reproduces Irot. Extracting the librational frequencies νlib = ν0 from the fits, it is evident that νlib increases with d1 (Fig. 7).

In Fig. 8, we show that DT is highly correlated with the difference in librational frequencies of species 1 (νlib,1) and species 2 (νlib,2). The effect of the mass dipole on νlib can be large; we observe νlib,1νlib,2 values of up to 100% of νlib,2. We fit the data to the equation DT = b(νlib,1νlib,2), demonstrating a direct proportionality relationship with DT. Given their statistical uncertainties, all data points agree with the fitting except for the (I1 = I2, d1, d2) = (402, 2.8σ, 0) mixture, which corresponds to the greatest νlib,1νlib,2 value of all the mixtures considered in this section. This suggests that DT may saturate at sufficiently large νlib,1νlib,2 values, beginning from 0.5τ1. We conclude that the librational frequency is sensitive to the mass dipole and that the mass dipole contribution to DT is highly correlated with the modification of the librational mode and therefore the short-time dynamics of the rod-like molecules.

FIG. 8.

Thermal diffusion coefficient DT as a function of νlib,1νlib,2, where νlib,i is the librational frequency of species i. The solid line shows the fit to DT = b(νlib,1νlib,2). Error bars correspond to maximum/minimum values based on the sensitivity of the fitting procedure for Eqs. (19)(21).

FIG. 8.

Thermal diffusion coefficient DT as a function of νlib,1νlib,2, where νlib,i is the librational frequency of species i. The solid line shows the fit to DT = b(νlib,1νlib,2). Error bars correspond to maximum/minimum values based on the sensitivity of the fitting procedure for Eqs. (19)(21).

Close modal

In this section, we examine how the mass dipoles affect the Soret coefficient of mixtures with non-zero moment of inertia contributions. We consider mixtures with components that differ only in their moments of inertia (I1I2, d1 = d2) with I1, I2 = 10–2702 and ratios I1/I2 = 1–27 and then mixtures with I1I2 andd1d2.

First we consider the simplest case of mass-symmetric molecules: I1I2 and d1 = d2 = 0. The Soret coefficients for these mixtures are shown in Figs. 9(a) and 9(b). In all cases ST > 0, indicating that the component with the greater moment of inertia (species 1) is thermophobic. For molecular mixtures, the pseudo-isotopic Soret effect is usually modeled by the empirical equation
ST=STchem+STiso=STchem+aMM1M2M1+M2+bII1I2I1+I2,
(22)
where STiso was originally taken by analogy to the description of gaseous mixtures3,19,46 and adjusted for liquids by the prefactors aM and bI.13,20 For the mixtures considered here, Eq. (22) reduces to
ST=bII1I2I1+I2,
(23)
which is also retrieved from the theoretical model47 of Villain-Guillot and Würger (derived by considering velocity fluctuations in a hard-bead model with elastic collisions) when setting M1 = M2. As shown in Fig. 9(a), Eq. (23) is in relatively good agreement with our simulation results. However, it is unable to predict all data points within their associated uncertainties. Previous simulations of binary isotopic mixtures of Lennard-Jones dumbbells with I1/I2 = 1–40 have shown very good agreement with the moment of inertia contribution given by Eq. (23).8,12 These simulations did not account for the mass dipole, but this contribution is likely to be small with d1d2 ≤ 1σ for the dumbbells. More generally, Eq. (22) has been able to accurately model the experimental Soret coefficients of relatively low molecular mass non-polar mixtures7 and isotope mixtures.7,19–23 However, M1 + M2 and I1 + I2 in these mixtures were almost constant; the terms in Eq. (22) could also be written in terms of absolute differences M1M2 and I1I2. Indeed, Wittko and Köhler found that the isotopic substitution C6H12 to C6D12 results in a constant shift in ST, irrespective of the nature of the other component.17 This trend is not reproduced by Eq. (22). Thus, we also fit our simulation data to
ST=bI(I1I2),
(24)
which similarly reproduces the ST values relatively well, as shown in Fig. 9(b).
FIG. 9.

The Soret coefficient ST as a function of (a) the relative (I1I2)/(I1 + I2) and (b) absolute (I1I2) differences in moment of inertia Ii of species i = 1, 2 for mixtures with mass dipoles d1 = d2 = 0; (c) (I1I2) for mixtures with d1 = d2; and (d) (d1d2) for mixtures with (I1, d1, I2, d2). In (d), the solid lines show fits to Eq. (25), and the label units for Ii and di are 2 and σ, respectively.

FIG. 9.

The Soret coefficient ST as a function of (a) the relative (I1I2)/(I1 + I2) and (b) absolute (I1I2) differences in moment of inertia Ii of species i = 1, 2 for mixtures with mass dipoles d1 = d2 = 0; (c) (I1I2) for mixtures with d1 = d2; and (d) (d1d2) for mixtures with (I1, d1, I2, d2). In (d), the solid lines show fits to Eq. (25), and the label units for Ii and di are 2 and σ, respectively.

Close modal

We emphasize that Eqs. (23) and (24) are empirical equations, and how accurately they model the moment of inertia contribution must be considered on a case-by-case basis for each mixture. Additionally, bI is typically fit simultaneously with multiple other coefficients,7,17 making it difficult to isolate the moment of inertia term. Isotopic substitution changes both the mass and inertia of a molecule. Thus, even for isotopic mixtures, the moment of inertia must be considered together with the total mass (and the mass dipole). Only recently, in simulations of isotopic mixtures of Lennard-Jones dumbbells with M1 = M2, has the moment of inertia contribution been studied independently from total mass contribution.12 The shish-kebab model used in this work allows us to isolate the moment of inertia contribution from both the total mass and mass dipole contributions.

Next, we consider mixtures with I1I2 and d1 = d2 ≠ 0. As shown in Fig. 9(c), increasing d1 = d2 while holding I1 and I2 = 102 constant decreases ST, and this decrease is more pronounced when I1 (i.e., I1I2) and ST are larger. In terms of its relative magnitude, ST decreases by 15% when increasing d1 = d2 from 0 to 2σ for the entire 100 ≤ (I1I2)/(2) ≤ 260 range (15% is comparable to the uncertainties for I1I2 = 1002). This demonstrates that the moment of inertia contribution is coupled with the mass dipole even when d1 = d2.

Finally, we consider mixtures with I1I2 and d1d2. As shown in Fig. 9(d), ST can be tuned by the mass dipoles d1 and d2 when holding I1 and I2 constant. In all cases, ST increases with d1d2, and the mass dipole contribution can both enhance or compete with the moment of inertia contribution. All current models of thermal diffusion predict that when the components differ by only their internal mass distribution, the component with the greater moment of inertia is thermophobic (ST > 0). Crucially, we demonstrate that by increasing −(d1d2), the competing mass dipole contribution causes a reversal in sign of ST and species 1 (I1 > I2) becomes thermophilic (ST < 0), or can lead to the inhibition of thermal diffusion (ST = 0).

Following from the power law dependence of the mass dipole contribution [Eq. (14)] described in Sec. III B, we model the I1I2 and d1d2 mixtures using the equation
ST=a(d1d2)|d1d2|k1+cI,
(25)
where k=kd2=0=1.6±0.1 as determined in Sec. III B. The prefactor a accounts for the strength of the mass dipole contribution including the coupling with the moments of inertia, while the constant cI(I1, I2) represents only the inertia contribution. As shown in Fig. 9(d), Eq. (25) is in excellent agreement with our simulation data; all data points agree to within their statistical uncertainties. a depends on both I1 and I2, indicating that the mass dipole contribution is coupled with the moments of inertia.
In this section, we propose an empirical equation for the isotopic Soret effect that includes the mass dipole contribution. Building on the success of equations of the form STiso=STM(M1,M2)+STI(I1,I2) for modeling experimental7,17,19–23,48 and simulation8 data, we assume that contributions from the different moments of the mass distribution are only weakly coupled and seek an equation of the form
STiso(M1,M2,d1,d2,I1,I2)=STM(M1,M2)+STd(d1,d2)+STI(I1,I2).
(26)
All mixtures considered in this work correspond to M1 = M2 and thus STM=0. However, we note that STM(M1M2)/(M1+M2), STM(M1M2), and STM(M11M21) have been used to model a range of molecular mixtures.7,17,18,48 For STd, we propose
STd(d1,d2)=(a1ds+a2)(d1d2)|d1d2|k1,
(27)
where ds = min{d1, d2} is the smaller mass dipole, a1 and a2 are constants, and k=kd2=0=1.6±0.1 as determined in Sec. III B. The (a1ds + a2) term stems from the observation that a in Eq. (14) is approximately linear with ds (d2 = ds for the I1 = I2 mixtures considered in Sec. III B), as evident from the inset of Fig. 8(b) in the supplementary material. For STI, we test both Eq. (23) (model A) and Eq. (24) (model B). Thus, for the isotopic mixtures studied in this work, the empirical equations are model A,
STA=STd(d1,d2)+bII1I2I1+I2,
(28)
and model B,
STB=STd(d1,d2)+bI(I1I2).
(29)
For the mass dipole contribution STd, the parameters a1 and a2 were fit using only the I1 = I2 mixtures. Analogously, for the moment of inertia contribution STI, the parameters bI were fit using only the d1 = d2 mixtures. The (I1I2, d1d2) mixtures were not used in the fitting procedure, and STX=A,B values for these mixtures are therefore true predictions.

We show in Fig. 10 that both models are in relatively good agreement with our simulation data. For the d1 = d2 mixtures, STX=STI have maximum absolute errors (MAEs) of 7.2 × 10−2 kBɛ−1 and 7.6 × 10−2 kBɛ−1 for models A and B, respectively. These are 4 times greater than the MAE of 1.8 × 10−2 kBɛ−1 for the I1 = I2 mixtures for which STX=STd. Thus, Eqs. (28) and (29) more accurately model the mass dipole contribution compared to the moment of inertia contribution. For the (I1I2, d1d2) mixtures, models A and B have MAEs of 4.0 × 10−2 kBɛ−1 and 3.0 × 10−2 kBɛ−1, respectively. We note that while a full range of mass dipoles (di ≤ 3 for this particle model) has been explored, the simulated (I1I2, d1d2) mixtures correspond to a more limited range of moments of inertia [see Fig. 9(d)]. Using the same range of (I1, I2) as the d1 = d2 mixtures would presumably result in higher MAEs for the (I1I2, d1d2) mixtures. For the mixtures studied in this work, the accuracy of the empirical models is limited by the (in)accuracy of the STI term rather than the assumption that the moments of the mass distribution are uncoupled.

FIG. 10.

Empirical models (a) A [Eq. (28)] and (b) B [Eq. (29)] for the isotopic Soret effect. (i) ST vs STX for model X = A, B, where STX and ST are the predicted and reference (simulated) Soret coefficients, respectively. (ii) The error STXST as a function of the mass dipole contribution STd,X and moment of inertia contribution STI,X of model X = A, B. Data are color coded according to the error STXST.

FIG. 10.

Empirical models (a) A [Eq. (28)] and (b) B [Eq. (29)] for the isotopic Soret effect. (i) ST vs STX for model X = A, B, where STX and ST are the predicted and reference (simulated) Soret coefficients, respectively. (ii) The error STXST as a function of the mass dipole contribution STd,X and moment of inertia contribution STI,X of model X = A, B. Data are color coded according to the error STXST.

Close modal

Next, we consider the validity of the weak coupling approximation. As shown in Fig. 4(a-iii), D12 increases (ST decreases) by 20% when decreasing I = I1 = I2 from 1002 to 402, which corresponds to the largest change in the mass dipole contribution due to the coupling with the moments of inertia. The a1 and a2 coefficients in STd [Eq. (27)] include an average over the effect of I on D12; we therefore expect a maximum error of 10% in STd due to the coupling with I1 and I2. With regard to the moment of inertia contribution (as discussed in Sec. III C), ST decreases by 15% when increasing d1 = d2 from 0 to 2σ for the (I2 = 102, d1 = d2) mixtures. Extrapolating this result to the entire range of mass dipoles di ≤ 3σ and moments of inertia, Ii ≤ 2602, we expect coupling effects in STI of ≲30% for the mixtures considered in this work.

The molecules with di ≠ 0 exhibit thermal orientation49 (TO): they adopt an average orientation with respect to the heat flux vector, quantified by cosθi=ûd,iûJq for species i = 1, 2, where ûd,i and ûJq are the unit vectors in the directions of the mass dipole and heat flux, respectively. Representative cosθz,i=ûd,iẑ profiles from the NEMD simulations are shown in Fig. 11(a). In all cases, ⟨cos θi⟩ ≥ 0, indicating that on average the mass dipole points toward the cold source (i.e., the heavy side points toward the cold source). This is consistent with simulations of mass-asymmetric diatomic molecules for which the heavier atom preferentially orients toward the cold source49 and follows the trend in the isotopic Soret effect of binary mixtures where the heavier component migrates toward the cold source (assuming the components are otherwise identical). As shown in Fig. 11(b-i), ⟨cos θ1⟩/|∇T| increases and then begins to saturate with d1. Holding d1 constant, ⟨cos θ1⟩/|∇T| increases with decreasing I = I1 = I2. ⟨cos θ1⟩/|∇T| also depends on d2 as demonstrated for the I = 802 mixtures; increasing d2 decreases ⟨cos θ1⟩/|∇T|. However, the effects of I and d2 on ⟨cos θ1⟩ are an order of magnitude smaller than its dependence on d1.

FIG. 11.

Thermal orientation of the rod-like molecules. (a) The average orientation ⟨cos θz,i⟩ of species i = 1, 2 as a function of position z in the NEMD simulations. The blue (cold) and red (hot) indicate the location of the thermostatting regions in the simulation cell. The profiles correspond to the same system shown in Fig. 2. (b) Average orientation ⟨cos θ1⟩ of species 1 for different mixtures: (i) ⟨cos θ1⟩/|∇T| vs the mass dipole d1 and (ii) the Soret coefficient ST vs ⟨cos θ1⟩ for different temperature gradients ∇T.

FIG. 11.

Thermal orientation of the rod-like molecules. (a) The average orientation ⟨cos θz,i⟩ of species i = 1, 2 as a function of position z in the NEMD simulations. The blue (cold) and red (hot) indicate the location of the thermostatting regions in the simulation cell. The profiles correspond to the same system shown in Fig. 2. (b) Average orientation ⟨cos θ1⟩ of species 1 for different mixtures: (i) ⟨cos θ1⟩/|∇T| vs the mass dipole d1 and (ii) the Soret coefficient ST vs ⟨cos θ1⟩ for different temperature gradients ∇T.

Close modal

Does thermal orientation have an effect on the Soret coefficient? Based on the general formalism of linear non-equilibrium thermodynamics32 (LNET), as a coupling effect, it does. For rigid rod-like colloids in the dilute regime, the TO effect decreases ST = DT/D primarily by increasing the diffusion coefficient D of the colloid.25 TO has also been shown to impact the Soret coefficient of Janus particles.26 However, we show in Fig. 11(b) that for the molecular mixtures studied here the impact of TO is smaller than the uncertainties associated with our ST data. As exemplified by three different mixtures, ST does not show a significant dependence on ⟨cos θi⟩. (For each mixture, greater/smaller |⟨cos θi⟩| values were achieved by increasing/decreasing ∇T in additional NEMD simulations; see the supplementary material for details). The |⟨cos θ⟩/∇T| values observed in this work and in other49 molecular fluids are ∼2–3 orders of magnitude smaller than those observed in the colloidal regime.24–26 This suggests that TO effects large enough to significantly impact ST may be more difficult to achieve in molecular mixtures. However, it is worth pointing out that thermal polarization, which emerges from the thermal orientation of polar molecules, can give rise to measurable electrostatic fields in molecular fluids.50 As a coupling effect, it leads to a reduction in the thermal conductivity of the polar fluid, as described by LNET.50,51 Thermal polarization can play an important role in determining the thermoelectric response of aqueous solutions.52 The impact of thermal polarization on the thermal diffusion of polar molecular mixtures remains an open question.

We have uncovered the mass dipole contribution to the isotopic Soret effect, showing that mixtures of components that differ only by the first moment of their mass distributions can have non-zero Soret coefficients. To the best of our knowledge, all current models of thermal diffusion describe the isotopic contribution STiso in terms of only the total mass (the zeroth moment) and moments of inertia (the second moment) of the components. Our results demonstrate that other moments of the mass distribution must be included for a complete description of the (pseudo-)isotopic Soret effect in fluid mixtures.

For the isotopic mixtures of rigid linear molecules examined in this work, the dependence of ST on the mass dipole arises mainly through the thermal diffusion coefficient DT. In turn, DT is correlated with both long-time dynamics—differences in rotational and translational diffusion coefficients, and short-time dynamics—the modification of a librational mode. Regarding the latter, greater mass dipoles give rise to higher frequency librations (of frequency νlib) with changes of up to 100%, and we find that DTνlib,1νlib,2. Through the self-diffusion coefficients Di=1,2, the mutual diffusion coefficient D12 features a very weak dependence of the mass dipoles, which can be attributed primarily to the effect of di on Di. However, D12 does depend significantly on the moment of inertia, which together with DT(d1, d2) affects how ST varies with the mass dipoles. Indeed, in mixtures with d1d2 and I1I2, the mass dipole contribution depends on both moments of inertia.

The mass dipole contribution can both enhance or compete with the moment of inertia contribution, giving rise to new phenomenology. For example, it is possible to design isotopic mixtures with M1 = M2 where the component with a greater moment of inertia is thermophilic. Additionally, the moment of inertia contribution depends on the mass dipoles, even when d1 = d2 and the mass dipole contribution vanishes. Our results show that changes in the internal mass distribution can be used to tune the Soret coefficient, including its sign and the thermophilicity of the components.

Building on the success of previous empirical models, we show that the mass dipole contribution can be incorporated as an additional additive contribution to the (pseudo-)isotopic Soret effect. This approach is valid when the contributions are only weakly coupled. For the mixtures considered in this work, the weak coupling assumption leads to estimated errors of ≲10% and ≲30% for the mass dipole and moment of inertia contributions, respectively. Indeed, the accuracy of the models is limited by the (in)accuracy of the moment of inertia term rather than the weak coupling assumption. Overall, the proposed empirical equations reproduce our simulation data well.

Our work highlights the importance of internal degrees of freedom, such as the mass dipole, in determining the thermodiffusion response of binary mixtures. Further work is required to assess the magnitude of the mass dipole contribution in other molecular mixtures and determine how it affects the phenomenology of these mixtures under the influence of thermal fields.

Supplementary material for this article. (1) Verification of linear response for the Soret coefficients, ST, calculated from the NEMD simulations. (2) Finite-size analysis of the Soret coefficients. (3) The phenomenological coefficient L11, shear viscosity η, and rotational diffusion coefficient Dr from their Green–Kubo relations. (4) Self-diffusion coefficients from the Einstein relation, related quantities, and finite-size analysis. (5) Estimated isotropic–nematic coexistence conditions. (6) Fitting empirical equations for the thermal diffusion coefficient DT.

We thank the Leverhulme Trust for Grant No. RPG-2018-384. We gratefully acknowledge a Ph.D. studentship (Project Reference No. 2135626) for O.R.G. sponsored by ICL’s Chemistry Doctoral Scholarship Award, funded by the EPSRC Doctoral Training Partnership Account (Grant No. EP/N509486/1). We also acknowledge Imperial College Research Computing Service, http://doi.org/10.14469/hpc/2232.

The authors have no conflicts to disclose.

Oliver R. Gittus: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (equal). Fernando Bresme: Conceptualization (supporting); Methodology (supporting); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
C.
Ludwig
, “
Diffusion zwischen ungleich erwärmten orten gleich zusammengesetzter lösungen
,” in
Sitzungsberichte der Bayerischen Akademie der Wissenschaften
(
Mathematisch-Naturwissenschaftliche Klasse
,
1856
), Vol.
539
.
2.
C.
Soret
, “
Sur l’état d’équilibre que prend, du point de vue de sa concentration, une dissolution saline primitivement homogène, dont duex parties sort portèes à des temperatures différentes
,”
C. R. Arch. Sci. Phys. Nat., Genève
2
,
48
61
(
1879
).
3.
R. C.
Jones
and
W. H.
Furry
, “
The separation of isotopes by thermal diffusion
,”
Rev. Mod. Phys.
18
,
151
224
(
1946
).
4.
S.
Chapman
, “
XIII. The possibility of separating isotopes
,”
London, Edinburgh, Dublin Philos. Mag. J. Sci.
38
,
182
186
(
1919
).
5.
K.
Clusius
and
G.
Dickel
, “
Neues verfahren zur gasentmischung und isotopentrennung
,”
Naturwissenschaften
26
,
546
(
1938
).
6.
W. H.
Furry
,
R. C.
Jones
, and
L.
Onsager
, “
On the theory of isotope separation by thermal diffusion
,”
Phys. Rev.
55
,
1083
1095
(
1939
).
7.
C.
Debuschewitz
and
W.
Köhler
, “
Molecular origin of thermal diffusion in benzene + cyclohexane mixtures
,”
Phys. Rev. Lett.
87
,
055901
(
2001
).
8.
G.
Galliéro
,
B.
Duguay
,
J.-P.
Caltagirone
, and
F.
Montel
, “
Thermal diffusion sensitivity to the molecular parameters of a binary equimolar mixture, a non-equilibrium molecular dynamics approach
,”
Fluid Phase Equilib.
208
,
171
188
(
2003
).
9.
P.-A.
Artola
,
B.
Rousseau
, and
G.
Galliero
, “
A new model for thermal diffusion: Kinetic approach
,”
J. Am. Chem. Soc.
130
,
10963
10969
(
2008
).
10.
D.
Reith
and
F.
Müller-Plathe
, “
On the nature of thermal diffusion in binary Lennard-Jones liquids
,”
J. Chem. Phys.
112
,
2436
2443
(
2000
).
11.
G.
Galliéro
, “
Thermal diffusion in Lennard–Jones fluids in the frame of the law of the corresponding states
,”
Fluid Phase Equilib.
224
,
13
22
(
2004
).
12.
H.
Hoang
and
G.
Galliero
, “
Predicting thermodiffusion in simple binary fluid mixtures
,”
Eur. Phys. J. E
45
,
42
(
2022
).
13.
D.
Niether
and
S.
Wiegand
, “
Thermophoresis of biological and biocompatible compounds in aqueous solution
,”
J. Phys.: Condens. Matter
31
,
503003
(
2019
).
14.
W.
Köhler
and
K. I.
Morozov
, “
The Soret effect in liquid mixtures—A review
,”
J. Non-Equilib. Thermodyn.
41
,
151
197
(
2016
).
15.
P.-A.
Artola
and
B.
Rousseau
, “
Thermal diffusion in simple liquid mixtures: What have we learnt from molecular dynamics simulations?
,”
Mol. Phys.
111
,
3394
3403
(
2013
).
16.
S.
Wiegand
, “
Thermal diffusion in liquid mixtures and polymer solutions
,”
J. Phys.: Condens. Matter
16
,
R357
R379
(
2004
).
17.
G.
Wittko
and
W.
Köhler
, “
Universal isotope effect in thermal diffusion of mixtures containing cyclohexane and cyclohexane-d12
,”
J. Chem. Phys.
123
,
014506
(
2005
).
18.
B.
Pur
,
W.
Köhler
, and
K. I.
Morozov
, “
The Soret effect of halobenzenes in n-alkanes: The pseudo-isotope effect and thermophobicities
,”
J. Chem. Phys.
152
,
054501
(
2020
).
19.
J.
Schirdewahn
,
A.
Klemm
, and
L.
Waldmann
, “
Thermodiffusion in D2-HT und anderen wasserstoffgemischen
,”
Z. Naturforsch. A
16
,
133
144
(
1961
).
20.
W. M.
Rutherford
, “
Effect of mass distribution on the isotopic thermal diffusion of substituted benzenes
,”
J. Chem. Phys.
81
,
6136
6139
(
1984
).
21.
W. M.
Rutherford
, “
Effect of mass distribution on the isotopic thermal diffusion of benzene
,”
J. Chem. Phys.
86
,
5217
(
1987
).
22.
W. M.
Rutherford
, “
Isotopic thermal diffusion of carbon disulfide in the liquid phase
,”
J. Chem. Phys.
86
,
397
399
(
1987
).
23.
W. M.
Rutherford
, “
Effect of carbon and hydrogen isotopic substitutions on the thermal diffusion of benzene
,”
J. Chem. Phys.
90
,
602
603
(
1989
).
24.
J.
Olarte-Plata
,
J. M.
Rubi
, and
F.
Bresme
, “
Thermophoretic torque in colloidal particles with mass asymmetry
,”
Phys. Rev. E
97
,
052607
(
2018
).
25.
O. R.
Gittus
,
J. D.
Olarte-Plata
, and
F.
Bresme
, “
Thermal orientation and thermophoresis of anisotropic colloids: The role of the internal composition
,”
Eur. Phys. J. E
42
,
90
(
2019
).
26.
F.
Bresme
,
J. D.
Olarte-Plata
,
A.
Chapman
,
P.
Albella
, and
C.
Green
, “
Thermophoresis and thermal orientation of Janus nanoparticles in thermal fields
,”
Eur. Phys. J. E
45
,
59
(
2022
).
27.
T. F.
Miller
,
M.
Eleftheriou
,
P.
Pattnaik
,
A.
Ndirango
,
D.
Newns
, and
G. J.
Martyna
, “
Symplectic quaternion scheme for biophysical molecular dynamics
,”
J. Chem. Phys.
116
,
8649
8659
(
2002
).
28.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in’t Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
, and
S. J.
Plimpton
, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
29.
C. D.
Daub
,
P.-O.
Åstrand
, and
F.
Bresme
, “
Thermo-molecular orientation effects in fluids of dipolar dumbbells
,”
Phys. Chem. Chem. Phys.
16
,
22097
22106
(
2014
).
30.
J. E.
Basconi
and
M. R.
Shirts
, “
Effects of temperature control algorithms on transport properties and kinetics in molecular dynamics simulations
,”
J. Chem. Theory Comput.
9
,
2887
2899
(
2013
).
31.
A.
Kowaguchi
,
P. E.
Brumby
, and
K.
Yasuoka
, “
Phase transitions and hysteresis for a simple model liquid crystal by replica-exchange Monte Carlo simulations
,”
Molecules
26
,
1421
(
2021
).
32.
S. R.
de Groot
and
P.
Mazur
,
Non-Equilibrium Thermodynamics
(
Dover
,
1984
).
33.
O. R.
Gittus
and
F.
Bresme
, “
On the microscopic origin of Soret coefficient minima in liquid mixtures
,”
Phys. Chem. Chem. Phys.
25
,
1606
1611
(
2023
).
34.
A.
Perronace
,
G.
Ciccotti
,
F.
Leroy
,
A. H.
Fuchs
, and
B.
Rousseau
, “
Soret coefficient for liquid argon-krypton mixtures via equilibrium and nonequilibrium molecular dynamics: A comparison with experiments
,”
Phys. Rev. E
66
,
031201
(
2002
).
35.
N. A. T.
Miller
,
P. J.
Daivis
,
I. K.
Snook
, and
B. D.
Todd
, “
Computation of thermodynamic and transport properties to predict thermophoretic effects in an argon-krypton mixture
,”
J. Chem. Phys.
139
,
144504
(
2013
).
36.
D.
Surblys
,
H.
Matsubara
,
G.
Kikugawa
, and
T.
Ohara
, “
Application of atomic stress to compute heat flux via molecular dynamics for systems with many-body interactions
,”
Phys. Rev. E
99
,
051301
(
2019
).
37.
D.
Surblys
,
H.
Matsubara
,
G.
Kikugawa
, and
T.
Ohara
, “
Methodology and meaning of computing heat flux via atomic stress in systems with constraint dynamics
,”
J. Appl. Phys.
130
,
215104
(
2021
).
38.

Due to an implementation detail in LAMMPS the y and z components of the heat flux for rigid bodies are highly unphysical and should not be used.

39.
LAMMPS Development Team
, LAMMPS documentation, https://docs.lammps.org/compute_heat_flux.html (accessed 1 July 2022).
40.
D. J.
Evans
and
S.
Murad
, “
Thermal conductivity in molecular fluids
,”
Mol. Phys.
68
,
1219
1223
(
1989
).
41.
I.-C.
Yeh
and
G.
Hummer
, “
System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions
,”
J. Phys. Chem. B
108
,
15873
15879
(
2004
).
42.
S. H.
Jamali
,
L.
Wolff
,
T. M.
Becker
,
A.
Bardow
,
T. J. H.
Vlugt
, and
O. A.
Moultos
, “
Finite-size effects of binary mutual diffusion coefficients from molecular dynamics
,”
J. Chem. Theory Comput.
14
,
2667
2677
(
2018
).
43.
M.
Eslamian
and
M. Z.
Saghir
, “
A critical review of thermodiffusion models: Role and significance of the heat of transport and the activation energy of viscous flow
,”
J. Non-Equilib. Thermodyn.
34
,
97
131
(
2009
).
44.
D. M.
Heyes
, “
Translational and rotational diffusion of rod shaped molecules by molecular dynamics simulations
,”
J. Chem. Phys.
150
,
184503
(
2019
).
45.
M.
Doi
and
S. F.
Edwards
,
The Theory of Polymer Dynamics
(
Clarendon Press
,
1986
).
46.
J. M.
Kincaid
,
E. G. D.
Cohen
, and
M.
López de Haro
, “
The Enskog theory for multicomponent mixtures. IV. Thermal diffusion
,”
J. Chem. Phys.
86
,
963
975
(
1987
).
47.
S.
Villain-Guillot
and
A.
Würger
, “
Thermal diffusion in a binary liquid due to rectified molecular fluctuations
,”
Phys. Rev. E
83
,
030501
(
2011
).
48.
S.
Hartmann
,
W.
Köhler
, and
K. I.
Morozov
, “
The isotope Soret effect in molecular liquids: A quantum effect at room temperatures
,”
Soft Matter
8
,
1355
1360
(
2012
).
49.
F.
Römer
,
F.
Bresme
,
J.
Muscatello
,
D.
Bedeaux
, and
J. M.
Rubí
, “
Thermomolecular orientation of nonpolar fluids
,”
Phys. Rev. Lett.
108
,
105901
(
2012
).
50.
O. R.
Gittus
,
P.
Albella
, and
F.
Bresme
, “
Polarization of acetonitrile under thermal fields via non-equilibrium molecular dynamics simulations
,”
J. Chem. Phys.
153
,
204503
(
2020
).
51.
F.
Bresme
,
A.
Lervik
,
D.
Bedeaux
, and
S.
Kjelstrup
, “
Water polarization under thermal gradients
,”
Phys. Rev. Lett.
101
,
020602
(
2008
).
52.
S.
Di Lecce
and
F.
Bresme
, “
Thermal polarization of water influences the thermoelectric response of aqueous solutions
,”
J. Phys. Chem. B
122
,
1662
1668
(
2018
).
53.
J. D.
Weeks
,
D.
Chandler
, and
H. C.
Andersen
,
J. Chem. Phys.
54
,
5237
5247
(
1971
).
Published open access through an agreement with JISC Collections

Supplementary Material