Control of Stratification in Drying Particle Suspensions via Temperature Gradients

A potential strategy for controlling stratification in a drying suspension of bidisperse particles is studied using molecular dynamics simulations. When the suspension is maintained at a constant temperature during fast drying, it can exhibit"small-on-top"stratification with an accumulation (depletion) of smaller (larger) particles in the top region of the drying film, consistent with the prediction of current theories based on diffusiophoresis. However, when only the region near the substrate is thermalized at a constant temperature, a negative temperature gradient develops in the suspension because of evaporative cooling at the liquid-vapor interface. Since the associated thermophoresis is stronger for larger nanoparticles, a higher fraction of larger nanoparticles migrate to the top of the drying film at fast evaporation rates. As a result, stratification is converted to"large-on-top". Very strong"small-on-top"stratification can be produced with a positive thermal gradient in the drying suspension. Here we explore a way to produce a positive thermal gradient by thermalizing the vapor at a temperature higher than that of the solvent. Possible experimental approaches to realize various thermal gradients in a suspension undergoing solvent evaporation, and thus to produce different stratification states in the drying film, are suggested.


I. INTRODUCTION
The drying of colloidal suspensions has been studied for several decades. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] Recently, drying-induced stratification phenomena in polydisperse colloidal mixtures have attracted great attention, [6,7,[18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] as they point to a quick, facile, one-pot method of depositing layered multifunctional coating films on a surface. In a particle suspension undergoing drying, the vertical distribution of particles is controlled by the Péclet number, Pe = Hv e /D, where H is the thickness of the suspension film, v e is the receding speed of the liquid-vapor interface during evaporation, and D is the diffusion coefficient of the particles. [9,34] The Péclet number characterizes the competition between diffusion and evaporation-induced particle migration. When Pe 1, the particles build up near the interface and their final distribution in the dry film may develop gradients, while for Pe 1, the particles diffuse fast enough to mitigate evaporative effects and are expected to be uniformly distributed in the deposited film. [9] In the case of a suspension of a bidisperse mixture of particles made from the same material but having different diameters, d l and d s , the final distribution of particles is determined by two Péclet numbers, Pe l and Pe s , for the large and small particles, respectively. If the Stokes-Einstein relationship holds, then Pe l /Pe s = d l /d s > 1. When Pe l > 1 > Pe s , Trueman et al. found the so-called "large-on-top" stratification, [12,13] where the larger (smaller) particles are enriched (depleted) near the receding interface. Recently, Fortini et al. discovered the counterintuitive "small-on-top" stratification in the regime of Pe l > Pe s 1, i.e., when drying is extremely rapid. [18,21] Since then, a number of experimental, [19,24,27,29,31,33] theoretical, [20,25,26] and simulation [22,23,28,30,32] studies have been reported on the stratification phenomena in drying suspensions of polydisperse particles and their physical mechanisms. The idea of diffusiophoresis being responsible for "small-on-top" stratification is widely supported. [18,20,22,25,26,30] In this picture, when Pe s 1, the smaller particles congregate near the receding interface during evaporation and their distribution develops a gradient that decays into the drying film. Further, when the volume fraction of the smaller particles, φ s , is above certain threshold that depends on Pe s , this gradient generates a diffusiophoretic force that is strong enough to push the larger particles out of the interfacial region. Consequently, the larger particles are depleted near the interface, resulting in "small-on-top" stratification.
The key ingredient of the diffusiophoretic model is that the cross-interaction between the large and small particles has asymmetric effects on the phoretic drift of particles and drives the larger ones away from the interfacial region faster than the smaller ones. [20,25] Therefore, the size asymmetry, quantified as α = d l /d s , is a crucial parameter that controls the outcome of stratification, with larger α favoring "small-on-top" stratification. Martín-Fabiani et al. studied a system with the smaller particles coated with hydrophilic shells and explored the effect of changing the pH of the initial dispersion. [19] In a dispersion with low pH, α is large enough to lead to "small-ontop" stratification. When the pH is raised, α is reduced as the hydrophilic shells swell substantially, and stratification is suppressed.
The approach of Martín-Fabiani et al. can be used for systems where the particle size can be tuned with external stimuli. [19] However, other possible approaches of controlling stratification for systems with fixed particle sizes have rarely been explored. In a previous work, [30] we used molecular dynamics (MD) modeling to study drying suspensions of a binary mixture of nanoparticles and found that for fast evaporation rates, the solvent can develop a negative temperature gradient toward the interface because of evaporative cooling effect. This temperature gradient induces thermophoresis, in which the larger particles are pushed more strongly into the interfacial region where the temperature is lower and the solvent density is higher. The competition between thermophoresis generated by evaporative cooling and diffusiophoresis can thus suppress "small-on-top" stratification at ultrafast drying rates or even turn the stratification into "large-on-top". [30] This discovery further indicates that thermophoresis, with a controlled thermal gradient other than the naturally occurring evaporative cooling, may be used to control stratification. In this paper, we employ MD modeling to test this idea in detail and demonstrate that stratification in a drying suspension can be controlled on demand with a temperature gradient imposed on the system, i.e., via controlled thermophoresis.

II. METHODS
We performed MD simulations on a suspension of a bidisperse mixture of nanoparticles. [30] The solvent is modeled explicitly as beads of mass m and interacting with each other via a Lennard-Jones (LJ) potential, U LJ (r) = 4 (σ/r) 12 − (σ/r) 6 − (σ/r c ) 12 + (σ/r c ) 6 , where r is the center-to-center distance between beads, is an energy scale, σ is a length scale, and the potential is truncated at r c = 3σ. The nanoparticles are modeled as spheres with a uniform distribution of LJ beads at a mass density 1.0m/σ. [35,36] The large nanoparticles (LNPs) have diameter d l = 20σ and mass m l = 4188.8m, and the small nanoparticles (SNPs) have diameter d s = 5σ and mass m s = 65.4m. The size ratio is α = 4. The nanoparticle-nanoparticle interactions are given by an integrated form of the LJ potential for two spheres with a Hamaker constant, A nn , characterizing the interaction strength. [35,36] In this study, A nn = 39.48 . To ensure that nanoparticles are well dispersed in the initial suspension, the nanoparticle-nanoparticle interactions are rendered purely repulsive by truncating them at 20.574σ, 13.085σ, and 5.595σ for the LNP-LNP, LNP-SNP, and SNP-SNP pairs, respectively. The nanoparticle-solvent interactions are described by a similar integrated form of the LJ potential with a Hamaker constant A ns = 100 and a cutoff length d/2 + 4σ, where d is the nanoparticle diameter. [37] The nanoparticle-solvent interactions thus have attractive tails, which make the effective diameter of a nanoparticle larger than its nominal diameter. [30]. The size ratio is defined here based on the nominal diameters of LNPs and SNPs. If their effective diameters are used, then the size ratio is about 3.4.
The entire system consists of ∼ 7 × 10 6 LJ beads, 200 LNPs, and 6400 SNPs. The system is placed in a rectangular simulation cell of dimensions L x × L y × L z , where L x = L y = 201σ, and L z = 477σ. The liquid-vapor interface is in the x-y plane, in which periodic boundary conditions are imposed. In the initial suspension, the thickness of the liquid film is about 304σ. The volume fractions of LNPs and SNPs in the initial dispersion are φ l = 0.068 and φ s = 0.034, respectively. Along the z-axis, all particles are confined in the simulation cell by two walls at z = 0 and z = L z . The particle-wall interaction is given by a LJ-like 9-3 potential, 3 , where the interaction strength W = 2.0 , h is the distance between the particle center and the wall, and h c is the cutoff length of the potential. For the solvent beads, D W = 1σ and h c = 3σ (0.8583σ) at the lower (upper) wall. With these parameters, the liquid solvent completely wets the lower wall while the upper wall is purely repulsive. For the nanoparticles, both walls are repulsive with D W = d/2 and h c = 0.8583D W , where d is the nanoparticle diameter.
To model evaporation of the solvent, a rectangular box of dimensions L x × L y × 20σ at the top of the simulation cell is designated as a deletion zone and a certain number (ζ) of vapor beads of the solvent in this zone are removed every τ , where τ = σ(m/ ) 1/2 is the reduced LJ unit of time. In this paper, two evaporation rates ζ = 30 and ζ = 5 are adopted. At these rates, the liquid-vapor interface retreats during evaporation at almost a constant speed, v e . The value of v e is determined for each evaporating suspension by directly computing the location of the interface as a function of time. The diffusion coefficients of nanoparticles are calculated with direct, independent simulations and the results are D l = 3.61 × 10 −3 σ 2 /τ for LNPs and D s = 2.11 × 10 −2 σ 2 /τ for SNPs at the initial volume fractions of nanoparticles prior to evaporation (see Supporting Information). The ratio D s /D l = 5.8 is higher than α = 4, the value expected from the Stokes-Einstein relation, which may be due to the finite concentrations of nanoparticles. [38] With values of D l , D s , v e , and H determined, the Péclet numbers for LNPs and SNPs, Pe l and Pe s , are computed for each evaporating system.
The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [39] was employed for all simulations reported here. A velocity-Verlet algorithm with a time step δt = 0.01τ was used to integrate the equation of motion. We have performed tests to confirm that the results reported here remain unchanged if a smaller time step is used. In the thermalized zone(s) specified for each system, a Langevin thermostat with a small damping rate Γ = 0.01τ −1 was used for the solvent beads. We have confirmed that this weak damping is strong enough to ensure a constant temperature in each thermalized liquid zone.
All results are presented below in the LJ units. Here we provide a rough mapping of these units to real ones in Table I by mapping the LJ solvent adopted in this paper to a liquid with a critical point similar to water, a typical solvent used in drying experiments. [7]. The details of this FIG. 1. Schematics of three types of thermalizations during solvent evaporation: (a) Only a thin layer of the liquid solvent adjacent to the bottom wall is thermalized at T l ; (b) All liquid and vapor are thermalized at T l ; (c) A thin layer of the liquid solvent adjacent to the bottom wall is thermalized at T l while the vapor zone at some distance away from the equilibrium liquid-vapor interface is thermalized at Tv. We set T l = 1.0 /kB and Tv can be higher or lower than T l to create a thermal gradient.
mapping are provided in the Supporting Information.

III. RESULTS AND DISCUSSION
Our goal is to demonstrate that a temperature gradient and the associated thermophoretic effect can be used to control stratification in a drying suspension of a polydisperse mixture of nanoparticles. We have previously shown that particles of different sizes have different thermophoretic responses to a thermal gradient. [30] In our previous work, only a thin layer of the liquid solvent adjacent to the bottom wall is thermalized at T l during evaporation, as shown in Fig. 1(a). Because of evaporative cooling at the liquid-vapor interface, a negative temperature gradient develops and its magnitude is larger for faster evaporation rates. The negative thermal gradient induces a positive gradient of the solvent density toward the interface, which generates a driving force to transport nanoparticles into the interfacial region. [40,41] The thermophoretic driving force is stronger for larger particles. The Soret coefficient, S T , can be used to characterize the strength of thermophoretic motion with respect to diffusive motion of particles. We have performed independent simulations of thermophoresis at A ns = 100 and found that for the LNPs, S T ∼ 0.1 K −1 , while for the SNPs, their thermophoretic response is extremely weak and S T is almost 0 (see Supporting Information). As a result, for very fast evaporation relatively more LNPs than SNPs are driven toward the interface in a drying bidisperse suspension. [30] The thermophoresis caused by evaporative cooling competes with the diffusiophoresis that leads to "small-on-top" stratification at fast drying rates, which is why only weak "small-on-top" stratification was observed in our previous simulations. [30] In certain cases the "small-on-top" stratification expected by the existing theory [25] was even converted to "largeon-top" in the presence of strong evaporative cooling. [30] Based on the physical picture depicted above, it is natural to investigate the effects of a controlled thermal gradient on stratification in a drying suspension. In this paper, we explore this idea by comparing three types of thermalization schemes as sketched in Fig. 1. The Scheme A is the same as in our previous work in which only a 10σ thick layer of the liquid solvent at the bottom of the suspension is thermalized at T l [ Fig. 1(a)]. [30] Evaporative cooling leads to a negative temperature gradient in the suspension toward the interface. In Scheme B, all solvent beads in the simulation cell are thermalized at T l [ Fig. 1(b)] and thus there are no thermal gradients during evaporation. In Scheme C, in addition to a liquid layer of thickness 10σ thermalized at T l near the bottom wall, the vapor beads with z-coordinates between L z − 150σ and L z are coupled to a thermostat with a target temperature T v [ Fig. 1(c)]. In this way, a thermal gradient is generated in the suspension with its direction and magnitude controlled by the difference between T v and T l , the thickness of the film, and the strength of evaporative cooling (i.e., the evaporation rate). For all systems studied in this paper, For Scheme A, the systems are labeled as T l 1.0 ζ y where the subscript y denotes the value of ζ. For Scheme B, T 1.0 ζ y is used to emphasize that the entire system is maintained at 1.0 /k B during evaporation. For Scheme C, the systems are labeled as T l 1.0 T v x ζ y , where x indicates the value of T v . All systems studied are listed in Table II.
We also studied systems with ζ = 5 and T v < T l , which show negative thermal gradients in the suspension and thermophoresis similar to those in T l 1.0 ζ 30 and T l 1.0 ζ 5 where evaporative cooling occurs. However, we observed condensation of droplets in the vapor phase if T v is made lower than the temperature at the liquid-  vapor interface in Scheme A with the same ζ. Despite this unwanted effect, cooling the vapor at a temperature lower than that of the suspension could be one experimental approach to apply a negative thermal gradient for systems that evaporate slowly or for which the effect of evaporative cooling is not as strong as that for the model LJ liquid employed in our simulations. The last six systems in Table II with T v varying from 0.75 /k B to 1.1 /k B are included in the Supporting Information.
In the main text we focus on the first five systems in Table II. Snapshots of the first five nanoparticle suspensions in Table II  Pe s > 1, [20] the thermophoresis associated with the negative temperature gradient works against diffusiophoresis and transports more LNPs into the interfacial region. As a result, the two systems exhibit "large-on-top" stratification.
When all solvent beads in the simulation cell are thermalized during evaporation, the temperature in the entire system is constant and no thermal gradients are produced. Thermophoresis is thus suppressed and only diffusiophoresis remains. The expected outcome is "smallon-top" stratification for Pe l Pe s > 1. The results from T 1.0 ζ 30 and T 1.0 ζ 5 confirm this prediction, as shown in Figs. 2(c) and (d). For example, comparing the last snapshot for T l 1.0 ζ 5 (the second row of Fig. 2) and that for T 1.0 ζ 5 (the fourth row of Fig. 2), the transition from "large-on-top" to "small-on-top" is visible after the thermal gradients and the associated thermophoresis are inhibited, especially from the distribution of LNPs in the drying films. This transition is verified quantitatively by an order parameter of stratification, which is discussed below (see Fig. 4).
The last row of Fig. 2 shows the snapshots for T l 1.0 T v 1.2 ζ 5 . In this system, the vapor about 23σ above the initial liquid-vapor interface prior to evaporation are thermalized at Consequently, there is a positive temperature gradient in the liquid solvent along the film's normal direction toward the interface. The solvent density develops a negative gradient and the accompanied thermophoresis drives LNPs toward the substrate. As a result, thermophoretic and diffusiophoretic effects are in synergy and strong "small-on-top" stratification is generated, which is apparent in Fig. 2(e) where the LNPs are enriched near the substrate during drying.
To understand quantitatively the stratification phenomena in drying particle suspensions, we plot the temperature and density profiles in Fig. 3. The local temperature T (z) at height z is computed from the average kinetic energy of the solvent beads in the spatial bin [z − 2.5σ, z + 2.5σ]. [42] The temperature profiles in the top row of Fig The local density of solvent or nanoparticles is com- represents the number of particles in the spatial bin [z − 0.5σ, z + 0.5σ] and m i is the particle mass. A nanoparticles straddling several bins is partitioned based on its partial volume in each bin. When computing the solvent density, the volume occupied by the nanoparticles is subtracted. The second row of Fig. 3 shows the solvent density as a function of height and the profiles exhibit gradients in accordance with the thermal gradients. Particularly, a positive (negative) thermal gradient generates a negative (positive) density gradient for the solvent and the stronger the thermal gradient, the stronger the density gradient. This correlation results from the fact that local thermal equilibrium is always maintained even at the fastest evaporation rates adopted in our simulations. [42] The density profile of the solvent affects the receding speed, v e , of the liquid-vapor interface. The data in Table II show that at the same ζ, the value of v e is slightly lower for T l 1.0 ζ y under Scheme A than for T 1.0 ζ y under Scheme B. For T l 1.0 ζ y , the evaporative cooling causes a positive gradient of the solvent density. The average solvent density is thus higher for T l 1.0 ζ y than for T 1.0 ζ y , as shown in Figs. 3(b), (f), (j), and (n). As a result, the liquid-vapor interface recedes more slowly in T l 1.0 ζ y than in T 1.0 ζ y at the same ζ. The density profiles for LNPs and SNPs are shown in the bottom two rows of Fig. 3, respectively. These profiles demonstrate the phoretic response of the nanopar-ticles to the thermal gradients (equivalently, the density gradients of the solvent induced by the thermal gradients) as well as the effects of the evaporation rate. For all simulations discussed here, the evaporation rates are high enough such that Pe l Pe s > 1. The corresponding fast receding liquid-vapor interface tends to trap both LNPs and SNPs just below the interface. If no other factors are at play, this effect combined with a large enough φ s is expected to yield "small-on-top" stratification via the diffusiophoresis mechanism as suggested by Sear and collaborators [18,25] and Zhou et al.. [20] This scenario is indeed the case for T 1.0 ζ 30 and T 1.0 ζ 5 , as shown in the third and fourth columns of Fig. 3 where there are no thermal gradients. The diffusiophoresis model also implies that the degree of "small-on-top" stratification is enhanced when the evaporation rate is increased. [18,20] However, as shown later, T 1.0 ζ 5 actually exhibits stronger "small-on-top" stratification than T 1.0 ζ 30 , even though the evaporation rate is increased six fold in the latter system. This discrepancy may be partially due to the small thickness of the suspension film studied in our simulations, which is limited by the available computational resources. The effect of film thickness on stratification is explored in a separate study. [43] Another reason may be that when the evaporation rate is increased, the drying time is shortened and there is less time for LNPs to diffuse out of the interfacial region via diffusiophoresis. As a result, "small-on-top" stratification is weakened when the evaporation rate is enhanced beyond certain threshold. This trend indicates that "small-on-top" stratification is most enhanced at some Pe l and is diminished if Pe l is increased further, which is consistent with two recent reports. [28,43] When only a thin layer of solvent beads at the bottom wall is thermalized, the temperature in the vicinity of the liquid-vapor interface decreases because of evaporative cooling effect. The resulting enhancement of the solvent density at the interface leads to thermophoretic drift of nanoparticles with the effect more significant for larger particles. This physical picture explains the observations for T l 1.0 ζ 30 and T l 1.0 ζ 5 . In these two systems, the SNPs are found to accumulate at the surface of the evaporating suspension as Pe l Pe s > 1 [Figs. 3(d) and (h)]. However, a significant accumulation of LNPs is found just below the enriched surface layer of SNPs, as shown in Figs. 3(c) and (g). The net outcome is actually "large-on-top" stratification, which will be confirmed later with an order parameter quantifying stratification (see Fig. 4). Furthermore, the degree of "large-on-top" stratification is stronger for T l 1.0 ζ 5 than for T l 1.0 ζ 30 , indicating a delicate competition between diffusiophoresis and thermophoresis. The lower evaporation rate in T l 1.0 ζ 5 suppresses both processes but it appears that diffusiophoresis favoring "small-on-top" is mitigated slightly more, creating stronger "large-on-top" for T l 1.0 ζ 5 . In our previous work, [30] we obtained a state diagram of stratification with systems all thermalized with Scheme A (i.e., a thin layer of liquid solvent contact-  ing the substrate is thermalized at T l = 1.0 /k B ) and only observed weak "small-on-top" stratification at values of Pe s and φ s far exceeding the critical values predicted by the diffusiophoretic model of Zhou et al.. [20] The presence of thermophoresis at fast evaporation rates may help understand the discrepancy between the simulations and the theory. [30] Indeed, when thermophoresis is suppressed, systems that are driven into the "large-ontop" regime by thermophoresis can be turned into (usually weak) "small-on-top". Examples are the transition from T l 1.0 ζ 30 to T 1.0 ζ 30 and that from T l 1.0 ζ 5 to T 1.0 ζ 5 . To achieve strong "small-on-top" stratification, a natural idea is to enable thermophoresis that works in conjunction with diffusiophoresis. This cooperation requires a thermal gradient during evaporation that is opposite to the one induced by evaporative cooling. To realize this, we thermalize the vapor zone from L z − 150σ to L z at a temperature T v > T l . The data in the fifth column of Fig. 3  It is expected that for systems thermalized with Scheme C and T v < T l , a negative thermal gradient develops in the liquid solvent, similar to the evaporative cooling case in Scheme A. Consequently, systems under Scheme C with T v < T l could display "large-ontop" stratification as long as the thermal gradient is large enough. These cases are in fact observed and discussed in detail in the Supporting Information, where some complications are noted related to droplet condensation in a vapor that is thermalized at low temperatures. Even for T v T l , the thermal gradient in the drying suspension can still be negative if evaporative cooling is strong enough. This is the case for T l To quantify stratification, we define an order parameter using the full density profiles of nanoparticles. [30] The mean heights of LNPs and SNPs are computed as The order parameter of stratification is then computed as (2 z l − 2 z s )/H(t), i.e., the mean separation between LNPs and SNPs nor- malized by H(t)/2, where H(t) is the instantaneous thickness of the suspension. In the equilibrium suspension prior to evaporation, both z l and z s are very close In Fig. 4 the order parameter of stratification is plotted against the extent of drying, quantified as (H(0) − H(t))/H(0), for the first five systems listed in Table II. It is clear that T 1.0 ζ 30 and T 1.0 ζ 5 exhibit "small-on-top" stratification since diffusiophoresis dominates while thermal gradients and thermophoresis are absent. The extent of stratification is slightly stronger for T 1.0 ζ 5 , though it dries more slowly. "Large-on-top" is observed for T l 1.0 ζ 30 and T l 1.0 ζ 5 and is again stronger for T l 1.0 ζ 5 that has a smaller evaporation rate. Although thermophoresis is much weaker for T l 1.0 ζ 5 because of the reduced evaporation rate, diffusiophoresis favoring "small-on-top" is suppressed even more when evaporation is slowed down and the delicate interplay of the two phoretic processes leads to stronger "large-on-top" stratification for T l 1.0 ζ 5 . A dramatic "small-on-top" state is clearly demonstrated in Fig. 4 for T l 1.0 T v 1.2 ζ 5 . Note that in the equilibrium suspension, φ l = 2φ s . If in the final dry film all SNPs were on top of all LNPs (i.e., complete stratification) but each group is uniformly distributed in its own region, then z l = H(t)/3 and z s = 5H(t)/6, yielding (2 z l − 2 z s )/H(t) = −1. As shown in Fig. 4, the order parameter of stratification reaches a minimal value around −0.5 for T l 1.0 T v 1.2 ζ 5 , indicating that the vertical distribution of the binary mixture of nanoparticles is substantially segregated in the drying film with SNPs on top of LNPs. This outcome is visually apparent in Fig. 2(e) as well.
The stratification order parameter used here and in Ref. [30] is based on the average position of nanoparti-cles, which is the first moment of their density profile in the entire drying film. This order parameter describes the systems studied here well and the identification of a stratified state is consistent with the classification based on the overall trend of the density profile of nanoparticles in the bulk of the drying film. Namely, "smallon-top" stratification generally corresponds to a negative gradient of the density profile of LNPs and a positive or nearly zero gradient of the density profile of SNPs from the bottom of the film to the top, while "large-on-top" is the other way around. However, this order parameter may not be applicable to oscillating density profiles or systems with only local stratification. In these more complicated situations, some other characteristics of the nanoparticle distribution including the higher moments of the density profile may be necessary to classify stratification. For all of our simulations, there is always a layer enriched with SNPs at the top of the drying film because Pe s > 1. However, it is misleading to call all these systems "small-on-top". Instead, information on the nanoparticle distribution in the film below this SNPrich surface layer should be taken into account as well. The order parameter used here fulfills this goal and yields a more consistent classification scheme for the outcome of stratification.
Evaporative cooling is a natural effect in a fast drying liquid. If a particle suspension is placed on a substrate that is kept at a constant temperature and the suspension undergoes very fast solvent evaporation, then a temperature lower than that of the substrate is expected at the evaporating interface, resulting in a negative thermal gradient in the suspension. T l 1.0 ζ 30 and T l 1.0 ζ 5 studied here are set up to mimic such situations. However, it is challenging to maintain a constant temperature or induce a positive thermal gradient along the normal direction toward the interface in a drying suspension, especially when the evaporation rate is high. One possible approach is to dissolve a gas (e.g., N 2 , Ar, He, or CO 2 ) into the solvent (e.g., water). Beaglehole showed that heating a water film with a dissolved gas from above or below produces very different temperature distributions within the liquid. [44,45] When heated from below, a fairly uniform temperature is found throughout the liquid. However, when the liquid is heated from above, a temperature gradient develops in it with the temperature higher at the liquid-vapor interface. Then it may be possible to study the effect of solvent evaporation on the particle distribution in a drying film under isothermal conditions and positive thermal gradients, similar to Scheme B and C.
In most experiments, films are much thicker than those studied here with MD and evaporation rates are much lower by a factor of 10 4 to 10 5 for drying at room temperature, about 45% of the critical temperature of water. In these systems, evaporative cooling is negligible and heat transfer is fast enough to make temperature uniform throughout a drying film. [46,47] To mimic this situation, Scheme B is used to maintain an isothermal drying film by coupling all solvent beads including vapor to a weak Langevin thermostat with a small damping rate, Γ = 0.01τ −1 . To address whether hydrodynamic interactions are screened in Langevin dynamics, [48] we run an additional simulation for T 1.0 ζ 30 with the Langevin thermostat replaced by a pairwise thermostat based on dissipative particle dynamics (DPD) with a weak friction coefficient γ = 0.05m/τ . [49] With the DPD thermostat, local momentum conservation is preserved throughout the simulation box and hydrodynamic interactions are expected to be fully captured. The results with the DPD thermostat are very close to those discussed here with the Langevin thermostat. These results are included in the Supporting Information. Under Scheme A and C, local momentum conservation is fulfilled away from the thermalized zones. All tests indicate that the Langevin thermostat adopted here is weak enough such that the viscosity of the LJ liquid is only weakly altered and the screening effect on hydrodynamic interactions is negligible.
In all simulations discussed thus far, the temperature of the thermalized liquid zone is always 1.0 /k B . The highest temperature used for the thermalized vapor zone is 1.2 /k B , which is close to the critical temperature, T c , of the LJ solvent with r c = 3.0σ. Furthermore, all simulations start with systems in which the nanoparticles are uniformly dispersed and then the evaporation process and the thermal gradient are imposed simultaneously. With this approach the interplay between diffusiophoresis and thermophoresis is investigated. To ensure that the physical mechanism of controlling stratification via thermophoresis is applicable to systems with both liquid and vapor temperatures way below T c , we run an additional simulation for RT l 0.9 T v 1.0 ζ 5 , i.e., with the bottom layer of the solvent adjacent to the lower wall thermalized at 0.9 /k B while the vapor zone above the liquid-vapor interface thermalized at 1.0 /k B . For RT l 0.9 T v 1.0 ζ 5 , the system is first relaxed under the imposed thermal gradient, which causes the LNPs to drift toward the lower wall via thermophoresis. The SNPs are still uniformly dispersed in the film as they are almost irresponsive to the thermal gradient. Then the solvent is evaporated from the relaxed system. The results for RT l 0.9 T v 1.0 ζ 5 are discussed in detail in the Supporting Information and fully consistent with the idea that thermophoresis from a positive thermal gradient from the bulk of a film to its drying front works in synergy with diffusiophoresis to enhance "small-on-top", while a negative thermal gradient works against diffusiophoresis to suppress "small-on-top" and promote "large-on-top".

IV. CONCLUSIONS
In this paper we focus on how stratification can be controlled in a drying suspension of a bidisperse mixture of nanoparticles via MD simulations with an explicit solvent model. We demonstrate that a thermal gradient and the induced thermophoresis can be used to alter stratification from "large-on-top" all the way to strong "small-on-top". This strategy is based on the observation that particles of different sizes in a suspension have different responses to a thermal gradient. In particular, larger particles experience a larger driving force that transports them into cooler regions where the solvent density is higher. For A ns = 100 adopted here, the smaller nanoparticles show little or even no response to a thermal gradient. When a suspension undergoes fast drying and only a thin layer of the solvent adjacent to the substrate is thermalized at T l , mimicking an experimental situation where the substrate supporting the suspension is maintained at a constant temperature during solvent evaporation, a negative temperature gradient develops in the suspension because of the evaporative cooling effect that makes the temperature at the evaporating interface to drop below T l . A larger fraction of the larger nanoparticles are driven into the interfacial region via the thermophoresis induced by this thermal gradient. As a result, the fast drying suspensions display "large-on-top" stratification instead of "small-on-top" expected by the diffusiophoresis model in which the suspension is assumed to be isothermal during evaporation. [18,20,25] Interestingly, when the entire suspension is maintained at T l during drying by thermalizing all solvent beads in the simulation cell, they do exhibit "small-on-top" stratification at fast evaporation rates, consistent with the prediction of the diffusiophoresis model. [18,20,25] However, the degree of stratification is found to be weak, probably due to the fact that φ s is small and the liquid film is thin for the simulations reported here. When a positive thermal gradient is induced in the suspension by thermalizing the vapor at a temperature sufficiently higher than T l , all larger nanoparticles are propelled toward the substrate. In this case, the synergy between thermophoresis and diffusiophoresis is underlying the observation of strong "small-on-top" stratification. Our results thus reveal a potentially useful strategy of controlling stratification via a regulated thermal gradient in a drying suspension of polydisperse particles.
The film thickness in our simulations prior to evaporation is about 300σ ∼ 105 nm. For a temperature difference 0.1 /k B across the film, the magnitude of the thermal gradient is about 0.5 K/nm if we take /k B ∼ 550 K as in Table I. This thermal gradient is several orders of magnitude larger than a typical experimental value. However, the suspensions simulated here are at temperatures not far from the critical temperature of the solvent, which allows the evaporation process of the solvent to be fast enough that can be modeled with MD. As a result, the evaporation rates in the MD simulations are also much higher than those in experiments. Nevertheless, as already discussed in Ref. [30], such high evaporation rates are needed to drive a sub-micron thin film suspension of nanoparticles into the "small-on-stop" regime (i.e., Pe s > 1), demonstrated in silico with our simulations. It is an interesting challenge if such a scenario can be realized experimentally, for example, by bringing the suspension close to the critical point of its solvent, as it may be relevant to the fabrication of multilayered thin film coatings. Because of high evaporation rates and the resulting strong evaporative cooling, large thermal gradients are needed in order to control stratification in drying thin films. For films with micrometer to millimeter thickness as in many experiments, [7] evaporation rates and thermal gradients smaller by a factor of 10 4 to 10 5 than those employed in MD simulations, i.e., those within the typical experimental range, will be sufficient to drive the systems into the same physical regime where thermophoresis is comparable to diffusiophoresis. In this sense, the results obtained here from MD simulations with thin films are scalable to much thicker films studied in experiments.

S1. Diffusion Coefficients of Nanoparticles:
The diffusion coefficients of the large nanoparticles (LNPs) and small nanoparticles (SNPs) were determined by an independent simulation. A suspension of LNPs and SNPs with the volume fractions (φ l = 0.060 and φ s = 0.033) close to those in the initial suspension discussed in the main text was prepared. The suspension has a cubic shape with edge length 202.6σ. Periodic boundary conditions were used in all three directions. The mean square displacements of both LNPs and SNPs as a function of time are shown in Fig. S1. The data show a clear transition from ballistic regime at short times to diffusive regime at long times. The diffusion coefficients are D l = 3.61 × 10 −3 σ 2 /τ for LNPs and D s = 2.11 × 10 −2 σ 2 /τ for SNPs. A simulation for a cubic simulation cell with edge length 101.3σ gives similar results. S2. Thermophoresis of Nanoparticles: To understand the thermophoresis of nanoparticles, we performed additional simulations for suspensions of only SNPs or only LNPs at volume fractions close to those in the suspension of the mixture of SNPs and LNPs. Each suspension was first equilibrated at T = 1.0 /k B . Then a thermal gradient was introduced into the system by thermalizing a top region of the liquid solvent and all vapor at T H = 1.0 /k B while thermalizing a layer of the solvent adjacent to the bottom wall at T L , as shown in Fig. S2. Two values of T L , 0.9 /k B and 0.7 /k B , were used to generate a thermal gradient with different magnitudes in the direction perpendicular to the liquid-vapor interface. The average position of nanoparticles in each system is plotted in Fig. S3. Since T L is lower than the initial temperature of the equilibrium suspension, the liquid contracts and the liquid-vapor interface recedes when the thermal gradient is imposed. Our data show that for the nanoparticle-solvent interaction with A ns = 100 , the SNPs first move toward the substrate because of the contraction of the liquid solvent. After this transient phase, the SNPs do not show any response to the imposed thermal gradient and their average position remains almost constant with time. However, the LNPs show a clear thermophoretic response to the thermal gradient and drift toward the cooler region where the liquid density is higher. These independent simulations demonstrate that for the parameters used in this paper, the LNPs exhibit strong thermophoresis while the SNPs exhibit no thermophoresis. Thermophoresis can be characterized by the Soret coefficient S T ≡ D T /D, where D T is the particle's thermal diffusion coefficient and D is its diffusion coefficient. We can estimate D T from v T = −D T ∇T , where v T is the drift velocity of the particles under the given thermal gradient ∇T . For a temperature difference of 0.3 /k B imposed on a liquid film with thickness about 100σ, ∇T 3 × 10 −3 /(k B σ). The LNPs exhibit positive thermophoresis and migrate about 30σ by average toward the cooler region during 5 × 10 4 τ (the dark brown curve in Fig. S3). The drift velocity v T is about −6 × 10 −4 σ/τ . The thermal diffusion coefficient of the LNPs is thus D T = −v T /∇T 0.2σ 2 k B /(τ ). Since for the LNPs, the diffusivity D = 3.61 × 10 −3 σ 2 /τ , the Soret coefficient of the LNPs is S T = D T /D 55k B / . If we take /k B 550 K as for water, then S T 0.1 K −1 for the LNPs. For a temperature difference 0.1 /k B across the film (the orange curve in Fig. S3), a similar analysis yields S T 110k B / 0.2 K −1 . Therefore, S T is at the order of 0.1 K −1 for the LNPs in our simulations. Since the SNPs do not exhibit much response to an imposed thermal gradient, we can effectively take their Soret coefficient to be 0.

S3. Additional Simulations and Results: Simultaneously Imposed Thermal Gradient and Evaporation
We ran 6 additional simulations in which the systems were thermalized with Scheme C with T l = 1.0 /k B and T v = 1.1 /k B , 1.05 /k B , 1.0 /k B , 0.9 /k B , 0.85 /k B , and 0.75 /k B at an evaporation rate set by ζ = 5, as listed in Table 1 Fig. S6(c).
For T l 1.0 T v 1.05 ζ 5 , the evaporative cooling in the solvent is slightly stronger than the heating from the vapor thermalized at T v = 1.05 /k B . The temperature profile has a very weak negative gradient in the evaporating suspension, as shown in Fig. S5(e). Accordingly, the density profile of the solvent has a very weak positive gradient, as shown in Fig. S5(f). Both LNPs and SNPs accumulate below the receding interface as Pe l > Pe s > 1 [see Figs. S5(g) and S5(h) and Figs. S6(a) and S6(b)]. The distribution of nanoparticles in this system does not stratify in the early stage of drying, though there is a transition to weak "large-on-top" at late times, as shown by the stratification order parameter in Fig. S6(c).

SI-2
Although for T l 1.0 T v 1.0 ζ 5 the vapor zone is coupled to a thermostat with a target temperature 1.0 /k B , which is equal to the temperature in the thermalized liquid layer at the bottom of the suspension, the evaporative cooling leads to a negative temperature gradient in the suspension, as shown in Fig. S5(i). Consequently, the solvent density increases from its bulk value in the bottom thermalized liquid layer to a higher value near the liquid-vapor interface, as shown in Fig. S5(j). The associated thermophoresis drives the LNPs toward the interfacial region [see Figs. S5(k) and S6(a)], though the SNPs are still accumulated right below the liquid-vapor interface as Pe s > 1 [see Figs. S5(l) and S6(b)]. The net outcome for T l 1.0 T v 1.0 ζ 5 is "large-on-top", as shown in Fig. S6(c). T l 1.0 T v 1.2 ζ 5 has been discussed in detail in the main text. Fig. S6 shows the results of the average position of nanoparticles normal to the interface and the order parameter of stratification for T l This comparison clearly shows the transition from strong to weak "small-on-top" stratification, to almost no stratification, and finally to "large-on-top" when T v − T l is reduced from 0.2 /k B to 0.
T l 1.0 ζ 5 has a weak negative thermal gradient in the evaporating suspension because of evaporative cooling. Via the associated thermophoretic process, the LNPs are pushed toward the liquid-vapor interface where the solvent density is higher. Thermophoresis overpowers diffusiophoresis where the LNPs are pushed away from the interface by the concentration gradient of the SNPs. As a result, T l 1.0 ζ 5 exhibits "large-on-top" stratification. Motivated by this observation, we ran T l 1.0 T v 0.75 ζ 5 that has a large |T v − T l |. Our original goal was to induce a large negative thermal gradient in the evaporating suspension at a relatively small v s such that the accompanying thermophoresis could generate strong "large-on-top" stratification in the drying suspension. However, as the initial vapor phase is at equilibrium with a liquid at T l = 1.0 /k B , cooling the vapor down to T v = 0.75 /k B leads to droplet condensation in the vapor phase, as shown in Fig. S4(d). Eventually, a liquid film of the solvent beads is formed at the top of the simulation cell. Because of the condensation, the vapor density near the liquid-vapor interface decreases very quickly and as a result, the solvent evaporates at a rate much higher than the target rate set by ζ = 5. The actual evaporation rate in T l The average position of nanoparticles during drying is plotted in Fig. S7 as a function of the extent of drying for the 5 systems discussed in the main text (see Table 2), as well as T l 1.0 T v 0.75 ζ 5 . It is clear that the LNPs have a large thermophoretic response while the SNPs only show changes in their average position in the late stage of drying, when the thermal gradient is varied. The behavior of SNPs is predominately affected by the receding speed of the liquidvapor interface. The variation of their average position in the late stage of drying under different thermal gradients is due to the change in the distribution of LNPs in the drying film. For example, when the LNPs are concentrated in a region, the SNPs are driven out of the same region because of crowding.
In Fig. S8, the order parameter of stratification is plotted as a function of the extent of drying for T l 1.0 ζ 30 , T l 1.0 ζ 5 , and T l 1.0 T v 0.75 ζ 5 . T l 1.0 T v 0.75 ζ 5 exhibits "large-on-top" stratification with an amplitude very close to that in T l 1.0 ζ 30 for reasons discussed previously. The reason that T l 1.0 ζ 5 shows even stronger "large-on-top" stratification than T l 1.0 ζ 30 is discussed in the main text.
In Fig. S9, we plot the average position of nanoparticles and the average separation between LNPs and SNPs for the 3 systems with T v < T l , i.e., T l 1.0 T v 0.75 ζ 5 , T l 1.0 T v 0.85 ζ 5 , and T l 1.0 T v 0.9 ζ 5 . While the SNPs show very little response to a thermal gradient, the distribution of LNPs is sensitive to the strength of the thermal gradient. When T v approaches T l from below, the negative thermal gradient in the evaporating suspension becomes smaller in magnitude and the driving force for LNPs to migrate to the interfacial region decreases. However, the receding speed of the liquid-vapor interface also decreases as T v is increased toward T l , and the diffusiophoretic driving force that pushes LNPs way from the interfacial region thus decreases too. The complicated interplay of these two factors makes "large-on-top" stratification stronger when T v is increased toward T l as in T l

S4. Additional Simulations and Results: Evaporation Imposed after Relaxation under a Given Thermal Gradient
In all simulations that have been discussed so far in the main text and this Supporting Information, the evaporation process and the thermal gradient are imposed simultaneously. All these systems start from a state in which the nanoparticles are uniformly dispersed in the suspension prior to evaporation and thermophoresis. In this way, the interplay between evaporation-induced nanoparticle migration and thermophoretic motion is investigated. In this section, we report an additional simulation in which the thermal gradient is first imposed on the nanoparticle suspension. The system is relaxed under the given thermal gradient. Then the solvent is evaporated out of the relaxed system.
The system is designated as RT l 0.9 T v 1.0 ζ 5 , which is thermalized under Scheme C with T l = 0.9 /k B and T v = 1.0 /k B , both well below the critical temperature of the Lennard-Jones solvent. The letter R in the system label indicates SI-3 that the system is first relaxed under the imposed thermal gradient. After relaxation, the solvent is evaporated at a rate ζ = 5. In Fig. S10, the snapshots of RT l 0.9 T v 1.0 ζ 5 are shown. The snapshot at t = 0 in Fig. S10 is for the relaxed system under the imposed thermal gradient. This gradient is apparent in the temperature profile (the blue curve) shown in Fig. S11(a). The solvent density develops a negative gradient and decreases toward the liquid-vapor interface, which can be easily seen from the density profile (the blue curve) in Fig. S11(b). Because of thermophoresis discussed in detail in Sec. S2, the LNPs by average migrate toward the lower part of the suspension where temperature is lower. This thermophoretic migration is difficult to see from the snapshot at t = 0 in Fig. S10 but is confirmed by the density profile of LNPs plotted in Fig. S11(c), where the blue curve is for the system at t = 0, right before the evaporation process is initiated. The overall negative gradient of the blue curve in Fig. S11(c) indicates that, as the outcome of thermophoresis, the LNPs are more concentrated in the lower region of the suspension in the relaxed state at t = 0. This result is also reflected in the average position of LNPs at t = 0 plotted in Fig. S12(a). As the SNPs do not exhibit any thermophoretic responses, their distribution in the relaxed suspension is not affected by the imposed temperature gradient, as shown in Figs. S10, S11(d), and S12(b). In terms of the nanoparticle distribution, the relaxed suspension is already in a "small-on-top" state, as indicated by the negative value of the stratification order parameter at t = 0 shown in Fig. S12(c) .
During the evaporation process of RT l 0.9 T v 1.0 ζ 5 , a weak negative temperature gradient emerges in the drying suspension as shown in Fig. S11(a) because of evaporative cooling, even though T v > T l . Correspondingly, the solvent density develops a weak positive gradient and increases slightly near the liquid-vapor interface, as shown in Fig. S11(b). The outcome of thermophoresis during evaporation is therefore to drive the LNPs toward the liquid-vapor interface, which is apparent from the density profiles of LNPs shown in Fig. S11(c). At t = 0 just prior to evaporation, the density profile of LNPs [the blue curve in Fig. S11(c)] exhibits a negative gradient as discussed previously. At t = 2.4 × 10 5 τ when 1.2 million solvent beads are evaporated, the density profile of LNPs [the purple curve in Fig. S11(c)] has developed a positive gradient, indicating that by average the LNPs migrate toward the interface during evaporation because of evaporative cooling and the associated thermophoresis. Although the SNPs do not respond to a thermal gradient, they accumulate just below the liquid-vapor interface as Pe s > 1. This accumulation can be seen from the density profiles of SNPs plotted in Fig. S11(d).
For RT l 0.9 T v 1.0 ζ 5 , because the imposed thermal gradient is outweighed by evaporative cooling, the temperature profile in the drying suspension has a negative gradient from the bottom of the suspension to the liquid-vapor interface. The resulting thermophoresis makes the LNPs to drift toward the interface and their accumulation is even more significant than the enrichment of SNPs near the receding liquid-vapor interface. As a consequence, the "smallon-top" stratification in the relaxed system is weakened during drying, as shown by the plot of the stratification order parameter vs time in Fig. S12(c).

S5. Additional Simulations and Results: DPD Thermostat
In all simulations discussed previously, a Langevin thermostat that is weakly coupled to the solvent beads has been used to control temperature. In Scheme A and Scheme C, a layer of solvent adjacent to the bottom wall is thermalized in a drying suspension. Additionally, the vapor zone is also thermalized in Scheme C. In these schemes, the most part of the system including the liquid-vapor interface has no thermostats and the dynamics is controlled by the interactions between the particles (i.e., the Hamiltonian). Local momentum conservation is preserved away from the thermalized zones. The hydrodynamic interactions between nanoparticles are captured with this approach.
In Scheme B, all solvent beads in the system are coupled to a Langevin thermostat in order to achieve an isothermal system during evaporation. To address whether the observed phenomena including stratification may be affected by the Langevin thermostat employed in the simulation, though the strength of the thermostat is weak as we use a small damping rate, 0.01τ −1 , we performed an additional simulation in Scheme B with the temperature of the system controlled by a thermostat based on dissipative particle dynamics (DPD) [for detail, see https://lammps.sandia. gov/doc/pair_dpd.html]. In the DPD thermostat, two particles experience an equal but opposite viscous drag force that depends on the difference between the velocities of the two particles. [49] For the entire system, all these drag forces cancel each other and the net drag is 0. The DPD thermostat is thus momentum-conserving and appropriate for capturing the hydrodynamic interactions in a suspension. We simulated the drying process of T 1.0 ζ 30 , where the temperature is kept at 1.0 /k B using a DPD thermostat with a small drag coefficient γ = 0.05m/τ . The results, including the snapshots of the system, the temperature and density profiles, and the order parameter of stratification, are shown in Figs. S13, S14, and S15, respectively. In summary, all results obtained with the DPD thermostat are very close to those obtained with the Langevin thermostat. Therefore, it is safe to use a weak Langevin thermostat to control temperature in molecular dynamics simulations of particle suspensions.

S6. Mapping between Lennard-Jones and Real Units
All results in this paper are reported in the reduced Lennard-Jones (LJ) units based on (energy unit), m (mass unit), and σ (length unit). Here we provide details on how the rough mapping between the LJ and real units given in Table 1 of the main text is derived. The critical temperature of the LJ fluid adopted here is estimated to be about 1.18 /k B . [42,50] By mapping this temperature to the critical temperature of water, 647 K, we obtain /k B 550 SI-4 K and 7.6 × 10 −21 J. The critical pressure of the LJ fluid used here as the solvent is estimated to be around 0.12 /σ 3 . [51,52] By mapping this pressure to the critical pressure of water, 22 MPa, we obtain σ 0.35 × 10 −9 m. The critical density of the LJ fluid used in this paper is around 0.31m/σ 3 . [42,[50][51][52] By mapping this density to the critical density of water, 322 kg/m 3 , we obtain m 4.5 × 10 −26 kg. After the mapping of , m, and σ to real values is determined, the mapping of the other LJ units to the real units can be derived, which is summarized in Table 1 of the main text. However, the rough mapping only serves as a crude guidance to expedite the comparison of the results reported in this paper with experimental measurements. The mapping is not expected to be accurate as the LJ fluid model adopted here is not anticipated to be an accurate model of water. For example, in this mapping the LJ liquid density is about 0.64m/σ 3 672 kg/m 3 at temperature 1.0 /k B 550 K. [15] However, the water density at 550 K is about 750 kg/m 3 . The dynamic viscosity of the LJ liquid used here is about 1.0m/(τ σ) 1.5 × 10 −4 Pa s at 1.0 /k B 550 K. [37,53] However, the viscosity of water at 550 K is about 30% lower at about 1.0 × 10 −4 Pa s.