Temperature gradients induce mass separation in mixtures in a process called thermal diffusion and are quantified by the Soret coefficient *S*_{T}. 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 *S*_{T} = 0 predicted by current formulations. For the isotopic mixtures investigated in this work, the dependence of *S*_{T} on the mass dipole arises mainly through the thermal diffusion coefficient *D*_{T}. In turn, *D*_{T} 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.

## I. INTRODUCTION

Thermal diffusion is the process by which a temperature gradient induces a concentration gradient in a mixture. It was first observed^{1} in the 19th century by Ludwig and later systematically studied^{2} 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 *F* ∝ *r*^{−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 suggest^{4} the use of thermal diffusion for separating isotopes. Isotope fractionation by thermal diffusion was first experimentally achieved in 1938 by Clusius and Dickel^{5} and subsequently received considerable theoretical treatment by Furry, Jones, and Onsager, among others.^{3,6}

*S*

_{T}from the mass and mass distribution of the particles. Using binary mixtures of isotopologues of benzene and cyclohexane, Debuschewitz and Köhler

^{7}demonstrated that the Soret coefficient could be split into two additive contributions,

^{8–12}The application of Eq. (1) to the interpretation of

*S*

_{T}data for molecular mixtures is now standard procedure.

^{13–16}

*M*

_{i}and their moments of inertia

*I*

_{i}. For binary mixtures,

*I*_{i}(3 × 3 matrices) can be reduced to scalars

*I*

_{i}. 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

*I*

_{i,k}or $Ii=3/\u2211kIi,k\u22121$.

^{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.

*μ*_{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*

*r*^{⊗n}=

**⊗**

*r***⊗ ⋯ ⊗**

*r***(**

*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}=

**is the mass dipole (3 × 1 vector), which gives the center of mass when normalized by the total mass,**

*D***/**

*D**M*=

*r*_{com}. The second moment

*μ*_{2}corresponds to the moment of inertia

**(3 × 3 matrix). Higher order moments,**

*I*

*μ*_{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*,

*r*

_{com}, and

**. While**

*I**M*and

**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.**

*I*^{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

**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).**

*I*## II. METHODS

### A. Particle model

^{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

*m*

_{i}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,

*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.

*M*= 56

*m*. 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

**: it is the displacement of the center of mass from the geometric center, $d=d\u22c5u\u0302d=rcom\u2212rg$. The geometric center**

*D*

*r*_{g}was chosen as the center of symmetry of the molecule (i.e., the inversion center given the

*D*

_{∞h}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),

*x*

_{i}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

*d*

_{i}and moments of inertia

*I*

_{i}of species

*i*= 1, 2. By convention, we assign species 1 according to

*I*

_{1}>

*I*

_{2}, and

*d*

_{1}>

*d*

_{2}if

*I*

_{1}=

*I*

_{2}.

### B. Simulation details

We performed both equilibrium and non-equilibrium molecular dynamics (NEMD) simulations of equimolar (mole fractions *x*_{1} = *x*_{2} = 0.5) binary mixtures targeting temperature *T* = 3.0 $\epsilon kB\u22121$ 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 (*L*_{x}, *L*_{y}, *L*_{z}) = (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 *T*_{h} and *T*_{c}, respectively. The thermostat temperatures were set to *T*_{c} = 2.5 $\epsilon kB\u22121$ and *T*_{h} = 3.5 $\epsilon kB\u22121$, 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 *t*_{e} ≥ 5 × 10^{4} *τ* to establish the stationary state followed by production runs of equal length between 10–30 × 10^{4} *τ*. 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 *t*_{s} = *l*^{2}/*D*, where *l* is the distance between the two thermostats and *D* is the self-diffusion coefficient. For our NEMD simulations, *t*_{e}/*t*_{s} ∼2–5, using the smaller of the two self-diffusion coefficients (corresponding to species *i* = 1, 2) for each mixture.

#### 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 10^{4} 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 10^{3} *τ* followed by 2 × 10^{4} *τ* of production.

Simulations of the pure fluid (same particle model but with the mass distribution given by *m*_{i} = *m* for all monomers) in the *NPT* ensemble were performed to determine the phase coexistence diagram. These simulations used cubic cells with 3 × 10^{4} 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 effects^{31} were not significant, two replicas were performed for each (*P*, *T*) state point; only systems where the *S*_{2} 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 × 10^{4} *τ* of equilibration starting from a nearby state point, followed by a 2 × 10^{4} *τ* production run.

### C. Transport properties from simulations

*S*

_{T}was evaluated from the NEMD simulations at the zero-mass-flux (

*J*_{1}= 0) stationary state as

*x*

_{i}and

*w*

_{i}are the mole and mass fractions of component

*i*. The local gradients ∇

*T*and ∇

*x*

_{1}were determined by fitting a straight line to their profiles within a range of ±10

*σ*around the selected state point at $T=3.0\epsilon kB\u22121$. Verification of linear response is given in the supplementary material.

*D*

_{12}was calculated from

*μ*

_{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,

*L*

_{11}were calculated using the Green–Kubo integral formula

*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

**is the velocity and the sums run over all seven monomers of each of the**

*v**N*

_{1}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 *D*_{T} 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 formulation^{36,37} (and using only the *x* component^{38,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, *D*_{T} was calculated from *D*_{T} = *S*_{T}*D*_{12}. This approach relies on the fact that finite-size effects associated with *S*_{T} from NEMD and *D*_{12} from equilibrium-*NVT* simulations are small and additionally cancel. It has recently been shown that the hydrodynamic correction of Yeh and Hummer^{41} for self-diffusion coefficients can also be applied to *D*_{12}.^{42} This correction amounts to an increase in *D*_{12}, *D*_{1}, and *D*_{2} of $<4%$ (see the supplementary material), which is smaller than the uncertainties associated with *D*_{12}, *S*_{T}, and therefore *D*_{T}. Regarding *S*_{T} from NEMD, extrapolating the simulation cell lengths *L*_{x} = *L*_{y} → ∞ decreases *S*_{T} by $\u223c4%$ (see the supplementary material). Thus, finite-size effects in *S*_{T} and *D*_{12} are expected to at least partially cancel when calculating *D*_{T}. This is because finite-size effects in *S*_{T} ≡ *D*_{T}/*D*_{12} include those from *D*_{12}; based on our finite-size analysis of *D*_{12} and *S*_{T}, we expect finite-size effects of $\u223c1%$ in *D*_{T}, which is smaller than its associated uncertainties. Thus, we use and report the uncorrected *D*_{12}, *D*_{1}, and *D*_{2} values, which give more internally consistent values for *S*_{T}, *D*_{12}, and *D*_{T}. We note that using the infinite-size *D*_{12}, *D*_{1}, and *D*_{2} 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 *D*_{12}, 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.

## III. RESULTS AND DISCUSSION

### A. Equation of state and phase diagram

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.

*S*

_{2}is zero for an isotropic fluid and unity for a perfectly aligned system.

*S*

_{2}is the ensemble average of the largest positive eigenvalue

*P*

_{2}of the

**tensor,**

*Q***is the eigenvector associated with**

*n**P*

_{2}and $u\u0302=u\u0302d$ 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 $\epsilon kB\u22121$ 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 (*P*_{zz}) 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.

### B. Mass dipole contribution and its microscopic origin

In this section, we consider mixtures with *I* = *I*_{1} = *I*_{2} but *d*_{1} ≠ *d*_{2} for which all current models of thermal diffusion predict *S*_{T} = 0 (see, e.g., the reviews Refs. 13–16 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, *S*_{T} > 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, *ɛ*/*k*_{B} = 10^{2} K, gives *S*_{T} = 1–15 × 10^{−3} K^{−1} at *T* = 300 K (*T* = 3.0 $\epsilon kB\u22121$) 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}

We probe the phenomenology of *S*_{T} ≡ *D*_{T}/*D*_{12} by considering the thermal diffusion coefficient *D*_{T} and mutual diffusion coefficient *D*_{12}. As shown in Fig. 4(a-iii), *D*_{12} decreases with increasing *I* for mixtures with *I* = *I*_{1} = *I*_{2} and *d*_{2} = 0 but is essentially constant with *d*_{1}. Furthermore, the constant value of *D*_{12} determined by fitting to the *D*_{12}(*d*_{1}; *I* = 80*mσ*^{2}, *d*_{2} = 0) data accurately predicts the *D*_{12} values for mixtures with *d*_{2} ≠ 0 as well [Fig. 4(b-iii)], indicating that *D*_{12} does not significantly depend on either mass dipole, *d*_{1} or *d*_{2}, at least relative to the statistical uncertainties associated with the data.

With regard to *D*_{T}, all the data for the *I* = *I*_{1} = *I*_{2} and *d*_{2} = 0 mixtures can be accurately described by a power law dependence $DT=a(d1\u2212d2)k$ [Fig. 4(a-ii)], indicating that *D*_{T} is only weakly dependent on the moment of inertia. The exponent $k=kd2=0=1.6\xb10.1$ was determined by fitting the equation ln *D*_{T} = *k* ln(*d*_{1} − *d*_{2}) + ln *a* to the *d*_{2} = 0 data (see the supplementary material). Considering mixtures with *d*_{2} ≠ 0, fits to the equation $DT=a(d1\u2212d2)kd2=0$ accurately reproduces the simulation data [Fig. 4(b-ii)]. Thus, the effect of changing *d*_{2} while keeping *d*_{1} − *d*_{2} constant can be captured by the prefactor *a*, which increases with increasing *d*_{2}, as shown explicitly in the inset of Fig. 8(b) in the supplementary material.

*S*

_{T}, which we predict with the equation

*a*and

*c*

_{I}were determined by the fits to

*D*

_{T}and

*D*

_{12}, respectively. We show in Figs. 4(a-i) and 4(b-i) that Eq. (14) accurately predicts

*S*

_{T}and captures its weak dependence on

*I*, which can be discerned for large values of

*d*

_{1}−

*d*

_{2}, as shown in the inset of Fig. 4(a-i).

We additionally calculate *D*_{12} using the Darken approximation $DMS\u2248DMSDarken=x2D1+x1D2$, where *D*_{MS} is the Maxwell–Stefan diffusion coefficient related to the Fickian *D*_{12} via *D*_{12} = Γ*D*_{MS}, 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) *D*_{12} values [Fig. 5(a)], indicating that displacement cross correlations between particles do not play an important role in determining *D*_{12}. Thus, the behavior of *D*_{12} can be rationalized in terms of the self-diffusion coefficients *D*_{1} and *D*_{2}. As shown in Fig. 5(b), $D12Darken$ features a weak but discernible dependence on *d*_{1} (and therefore *d*_{1} − *d*_{2}), which was masked by the greater uncertainties associated with *D*_{12} calculated from the Green–Kubo integral formula. *D*_{1} and *D*_{2} also depend weakly on the mass dipole *d*_{1} [Figs. 5(c) and 5(d)]. Unsurprisingly, the impact of *d*_{1} on *D*_{1} is greater than on *D*_{2} with Δ*D*_{1}/Δ*D*_{2} ∼ 1–4, where Δ*D*_{i} = *D*_{i}(*d*_{1}) − *D*_{i}(*d*_{1} = 0.5*σ*). Thus, the weak dependence of *D*_{12} on the mass dipoles *d*_{1} and *d*_{2} can be primarily attributed to the effect of *d*_{i} on *D*_{i}. (See the supplementary material for additional data on the diffusion coefficients and related quantities).

The mass dipole also affects the rotational diffusion of the rod-like molecules. We show in Fig. 6(a) that the rotational diffusion coefficient *D*_{r,1} decreases with increasing mass dipole *d*_{1}. *D*_{r,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 > (*D*_{i}/*D*_{r,i})/(*σ*^{2} rad^{−2}) > 13 for all the mixtures considered in this section. As shown in Fig. 6(b), the ratio *D*_{1}/*D*_{r,1} increases with the mass dipole *d*_{1} but is statistically independent of *I* within the statistical uncertainties of our data. We note that comparing to previous simulations^{44} of systems of the same WCA shish-kebab model but at the different temperature *T* = 1 $\epsilon kB\u22121$, 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}

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

*C*

_{atom}, as well as the rotational and translational center-of-mass VACFs of the rigid molecules,

*C*

_{rot}and

*C*

_{trans}, respectively,

*v*_{i}is the velocity of monomer

*i*,

*ω*_{j}is the angular velocity of rigid molecule

*j*,

*v*

_{COM,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,

*I*

_{X}is the intensity and

*ν*is the frequency.

The power spectra of species 1 in the *I* = *I*_{1} = *I*_{2} = 80*mσ*^{2} and *d*_{2} = 0 mixtures are shown in Fig. 7. For large mass dipoles, *I*_{atom} features a peak at $\u223c0.5\tau \u22121$ indicative of a librational mode. Increasing the mass dipole increases the frequency of the librational mode. This librational mode is also reflected in *I*_{rot}. 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.

*I*

_{rot}using the following equations:

*I*

_{L}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*=

*I*

_{1}=

*I*

_{2}= 80

*mσ*

^{2}and

*d*

_{2}= 0 mixtures as an example, we show in Fig. 7(b) that fitting Eq. (19) accurately reproduces

*I*

_{rot}. Extracting the librational frequencies

*ν*

_{lib}=

*ν*

_{0}from the fits, it is evident that

*ν*

_{lib}increases with

*d*

_{1}(Fig. 7).

In Fig. 8, we show that *D*_{T} 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 $\u223c100%$ of *ν*_{lib,2}. We fit the data to the equation *D*_{T} = *b*(*ν*_{lib,1} − *ν*_{lib,2}), demonstrating a direct proportionality relationship with *D*_{T}. Given their statistical uncertainties, all data points agree with the fitting except for the (*I*_{1} = *I*_{2}, *d*_{1}, *d*_{2}) = (40*mσ*^{2}, 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 *D*_{T} may saturate at sufficiently large *ν*_{lib,1} − *ν*_{lib,2} values, beginning from $\u223c0.5\tau \u22121$. We conclude that the librational frequency is sensitive to the mass dipole and that the mass dipole contribution to *D*_{T} is highly correlated with the modification of the librational mode and therefore the short-time dynamics of the rod-like molecules.

### C. Coupling with the moment of inertia contribution

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 (*I*_{1} ≠ *I*_{2}, *d*_{1} = *d*_{2}) with *I*_{1}, *I*_{2} = 10–270*mσ*^{2} and ratios *I*_{1}/*I*_{2} = 1–27 and then mixtures with *I*_{1} ≠ *I*_{2} and*d*_{1} ≠ *d*_{2}.

*I*

_{1}≠

*I*

_{2}and

*d*

_{1}=

*d*

_{2}= 0. The Soret coefficients for these mixtures are shown in Figs. 9(a) and 9(b). In all cases

*S*

_{T}> 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

^{3,19,46}and adjusted for liquids by the prefactors

*a*

_{M}and

*b*

_{I}.

^{13,20}For the mixtures considered here, Eq. (22) reduces to

^{47}of Villain-Guillot and Würger (derived by considering velocity fluctuations in a hard-bead model with elastic collisions) when setting

*M*

_{1}=

*M*

_{2}. 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

*I*

_{1}/

*I*

_{2}= 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

*d*

_{1}−

*d*

_{2}≤ 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 mixtures

^{7}and isotope mixtures.

^{7,19–23}However,

*M*

_{1}+

*M*

_{2}and

*I*

_{1}+

*I*

_{2}in these mixtures were almost constant; the terms in Eq. (22) could also be written in terms of absolute differences

*M*

_{1}−

*M*

_{2}and

*I*

_{1}−

*I*

_{2}. Indeed, Wittko and Köhler found that the isotopic substitution C

_{6}H

_{12}to C

_{6}D

_{12}results in a constant shift in

*S*

_{T}, 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

*S*

_{T}values relatively well, as shown in Fig. 9(b).

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, *b*_{I} 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 *M*_{1} = *M*_{2}, 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 *I*_{1} ≠ *I*_{2} and *d*_{1} = *d*_{2} ≠ 0. As shown in Fig. 9(c), increasing *d*_{1} = *d*_{2} while holding *I*_{1} and *I*_{2} = 10*mσ*^{2} constant decreases *S*_{T}, and this decrease is more pronounced when *I*_{1} (i.e., *I*_{1} − *I*_{2}) and *S*_{T} are larger. In terms of its relative magnitude, *S*_{T} decreases by $\u223c15%$ when increasing *d*_{1} = *d*_{2} from 0 to 2*σ* for the entire 100 ≤ (*I*_{1} − *I*_{2})/(*mσ*^{2}) ≤ 260 range (15% is comparable to the uncertainties for *I*_{1} − *I*_{2} = 100*mσ*^{2}). This demonstrates that the moment of inertia contribution is coupled with the mass dipole even when *d*_{1} = *d*_{2}.

Finally, we consider mixtures with *I*_{1} ≠ *I*_{2} and *d*_{1} ≠ *d*_{2}. As shown in Fig. 9(d), *S*_{T} can be tuned by the mass dipoles *d*_{1} and *d*_{2} when holding *I*_{1} and *I*_{2} constant. In all cases, *S*_{T} increases with *d*_{1} − *d*_{2}, 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 (*S*_{T} > 0). Crucially, we demonstrate that by increasing −(*d*_{1} − *d*_{2}), the competing mass dipole contribution causes a reversal in sign of *S*_{T} and species 1 (*I*_{1} > *I*_{2}) becomes thermophilic (*S*_{T} < 0), or can lead to the inhibition of thermal diffusion (*S*_{T} = 0).

*I*

_{1}≠

*I*

_{2}and

*d*

_{1}≠

*d*

_{2}mixtures using the equation

*a*accounts for the strength of the mass dipole contribution including the coupling with the moments of inertia, while the constant

*c*

_{I}(

*I*

_{1},

*I*

_{2}) 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

*I*

_{1}and

*I*

_{2}, indicating that the mass dipole contribution is coupled with the moments of inertia.

### D. An empirical equation for the isotopic Soret effect

^{7,17,19–23,48}and simulation

^{8}data, we assume that contributions from the different moments of the mass distribution are only weakly coupled and seek an equation of the form

*M*

_{1}=

*M*

_{2}and thus $STM=0$. However, we note that $STM\u221d(M1\u2212M2)/(M1+M2)$, $STM\u221d(M1\u2212M2)$, and $STM\u221d\u2212(M1\u22121\u2212M2\u22121)$ have been used to model a range of molecular mixtures.

^{7,17,18,48}For $STd$, we propose

*d*

_{s}= min{

*d*

_{1},

*d*

_{2}} is the smaller mass dipole,

*a*

_{1}and

*a*

_{2}are constants, and $k=kd2=0=1.6\xb10.1$ as determined in Sec. III B. The (

*a*

_{1}

*d*

_{s}+

*a*

_{2}) term stems from the observation that

*a*in Eq. (14) is approximately linear with

*d*

_{s}(

*d*

_{2}=

*d*

_{s}for the

*I*

_{1}=

*I*

_{2}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,

*a*

_{1}and

*a*

_{2}were fit using only the

*I*

_{1}=

*I*

_{2}mixtures. Analogously, for the moment of inertia contribution $STI$, the parameters

*b*

_{I}were fit using only the

*d*

_{1}=

*d*

_{2}mixtures. The (

*I*

_{1}≠

*I*

_{2},

*d*

_{1}≠

*d*

_{2}) 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 *d*_{1} = *d*_{2} mixtures, $STX=STI$ have *maximum* absolute errors (MAEs) of 7.2 × 10^{−2} *k*_{B}*ɛ*^{−1} and 7.6 × 10^{−2} *k*_{B}*ɛ*^{−1} for models A and B, respectively. These are $\u223c4$ times greater than the MAE of 1.8 × 10^{−2} *k*_{B}*ɛ*^{−1} for the *I*_{1} = *I*_{2} 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 (*I*_{1} ≠ *I*_{2}, *d*_{1} ≠ *d*_{2}) mixtures, models A and B have MAEs of 4.0 × 10^{−2} *k*_{B}*ɛ*^{−1} and 3.0 × 10^{−2} *k*_{B}*ɛ*^{−1}, respectively. We note that while a full range of mass dipoles (*d*_{i} ≤ 3 for this particle model) has been explored, the simulated (*I*_{1} ≠ *I*_{2}, *d*_{1} ≠ *d*_{2}) mixtures correspond to a more limited range of moments of inertia [see Fig. 9(d)]. Using the same range of (*I*_{1}, *I*_{2}) as the *d*_{1} = *d*_{2} mixtures would presumably result in higher MAEs for the (*I*_{1} ≠ *I*_{2}, *d*_{1} ≠ *d*_{2}) 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.

Next, we consider the validity of the weak coupling approximation. As shown in Fig. 4(a-iii), *D*_{12} increases (*S*_{T} decreases) by $\u223c20%$ when decreasing *I* = *I*_{1} = *I*_{2} from 100*mσ*^{2} to 40*mσ*^{2}, which corresponds to the largest change in the mass dipole contribution due to the coupling with the moments of inertia. The *a*_{1} and *a*_{2} coefficients in $STd$ [Eq. (27)] include an average over the effect of *I* on *D*_{12}; we therefore expect a maximum error of $\u223c10%$ in $STd$ due to the coupling with *I*_{1} and *I*_{2}. With regard to the moment of inertia contribution (as discussed in Sec. III C), *S*_{T} decreases by $\u223c15%$ when increasing *d*_{1} = *d*_{2} from 0 to 2*σ* for the (*I*_{2} = 10*mσ*^{2}, *d*_{1} = *d*_{2}) mixtures. Extrapolating this result to the entire range of mass dipoles *d*_{i} ≤ 3*σ* and moments of inertia, *I*_{i} ≤ 260*mσ*^{2}, we expect coupling effects in $STI$ of ≲30% for the mixtures considered in this work.

### E. Thermal orientation

The molecules with *d*_{i} ≠ 0 exhibit thermal orientation^{49} (TO): they adopt an average orientation with respect to the heat flux vector, quantified by $\u27e8cos\theta i\u27e9=u\u0302d,i\u22c5u\u0302Jq$ for species *i* = 1, 2, where $u\u0302d,i$ and $u\u0302Jq$ are the unit vectors in the directions of the mass dipole and heat flux, respectively. Representative $\u27e8cos\theta z,i\u27e9=u\u0302d,i\u22c5z\u0302$ 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 source^{49} 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 *d*_{1}. Holding *d*_{1} constant, ⟨cos *θ*_{1}⟩/|∇*T*| increases with decreasing *I* = *I*_{1} = *I*_{2}. ⟨cos *θ*_{1}⟩/|∇*T*| also depends on *d*_{2} as demonstrated for the *I* = 80*mσ*^{2} mixtures; increasing *d*_{2} decreases ⟨cos *θ*_{1}⟩/|∇*T*|. However, the effects of *I* and *d*_{2} on ⟨cos *θ*_{1}⟩ are an order of magnitude smaller than its dependence on *d*_{1}.

Does thermal orientation have an effect on the Soret coefficient? Based on the general formalism of linear non-equilibrium thermodynamics^{32} (LNET), as a coupling effect, it does. For rigid rod-like colloids in the dilute regime, the TO effect decreases *S*_{T} = *D*_{T}/*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 *S*_{T} data. As exemplified by three different mixtures, *S*_{T} 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 other^{49} 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 *S*_{T} 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.

## IV. CONCLUSIONS

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 *S*_{T} on the mass dipole arises mainly through the thermal diffusion coefficient *D*_{T}. In turn, *D*_{T} 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 *D*_{T} ∝ *ν*_{lib,1} − *ν*_{lib,2}. Through the self-diffusion coefficients *D*_{i=1,2}, the mutual diffusion coefficient *D*_{12} features a very weak dependence of the mass dipoles, which can be attributed primarily to the effect of *d*_{i} on *D*_{i}. However, *D*_{12} does depend significantly on the moment of inertia, which together with *D*_{T}(*d*_{1}, *d*_{2}) affects how *S*_{T} varies with the mass dipoles. Indeed, in mixtures with *d*_{1} ≠ *d*_{2} and *I*_{1} ≠ *I*_{2}, 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 *M*_{1} = *M*_{2} where the component with a greater moment of inertia is thermophilic. Additionally, the moment of inertia contribution depends on the mass dipoles, even when *d*_{1} = *d*_{2} 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

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Sitzungsberichte der Bayerischen Akademie der Wissenschaften*

_{12}

*n*-alkanes: The pseudo-isotope effect and thermophobicities

_{2}-HT und anderen wasserstoffgemischen

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.