Modeling Solution Drying by Moving a Liquid-Vapor Interface: Method and Applications

A method of simulating the drying process of a soft matter solution with an implicit solvent model by moving the liquid-vapor interface is applied to various solution films and droplets. For a solution of a polymer and nanoparticles, we observe “polymer-on-top” stratification, similar to that found previously with an explicit solvent model. Furthermore, “polymer-on-top” is found even when the nanoparticle size is smaller than the radius of gyration of the polymer chains. For a suspension droplet of a bidisperse mixture of nanoparticles, we show that core-shell clusters of nanoparticles can be obtained via the “small-on-outside” stratification mechanism at fast evaporation rates. “Large-on-outside” stratification and uniform particle distribution are also observed when the evaporation rate is reduced. Polymeric particles with various morphologies, including Janus spheres, core-shell particles, and patchy particles, are produced from drying droplets of polymer solutions by combining fast evaporation with a controlled interaction between the polymers and the liquid-vapor interface. Our results validate the applicability of the moving interface method to a wide range of drying systems. The limitations of the method are pointed out and cautions are provided to potential practitioners on cases where the method might fail.

To model the drying process of a soft matter solution, a key challenge is the treatment of the solvent. To capture factors that may be important in an evaporation process, such as hydrodynamic interactions between solutes [30,34], evaporation-induced flow in the solution (e.g., capillary flow in an evaporating droplet showing the coffeering effect) [39], and instabilities during drying including Rayleigh-Bénard and Bénard-Marangoni instabilities [8,[40][41][42], it is ideal to include the solvent explicitly in a computational model [15,24,25,[30][31][32]34,[36][37][38]. However, such models are computationally extremely expensive as the solvent particles significantly outnumber the solutes at realistic volume fractions. Usually, systems containing several million particles or more have to be considered even for a few hundred nanoparticles [24,43]. Because of the large system size, only fast evaporation rates can be modeled using this approach [15,24,25,[30][31][32]34,[36][37][38]. Considering the limitations of explicit solvent models, it is natural to explore alternative approaches that are computationally more efficient and able to quickly yield results that are at least qualitatively reasonable. One such effort is to model the solvent as an uniform and viscous medium in which the dispersed solutes are moving around. This approach leads to various implicit solvent models of soft matter solutions [30,34,38,44].
Recently, implicit solvent models have been applied to study the evaporation process of particle suspensions and polymer solutions and reveal many interesting phenomena induced by drying [26][27][28][29][30][33][34][35]38]. Fortini et al. used Langevin dynamics simulations based on an implicit solvent model to study the drying of a bidisperse colloidal suspension film and demonstrated the counterintuitive "small-on-top" stratification in which the smaller particles are predominately distributed on top of the larger particles after drying [26,27]. Tatsumi et al. used a similar model to investigate the role of evaporation rates on stratification [35]. Howard et al. and Statt et al. employed this approach to simulate the drying of colloidal suspensions [30], colloidal mixtures [28], polymer-colloid mixtures [29], polymer-polymer mixtures [29,34], and polydisperse polymer mixtures [33], and observed stratifying phenomena as well. Recently, we demonstrated that comparable stratification could be observed for colloidal suspensions in both explicit and implicit solvent models that are carefully matched [38].
In the implicit solvent approach to modeling drying, all the solutes are confined by a potential barrier, which represents the confinement effect of the liquid-vapor interface between the solvent and its vapor. One simple form of the confining potential is a harmonic potential. In a previous work, Tang and Cheng analyzed the capillary force experienced by a small spherical particle at a liquid-vapor interface and showed that this harmonic approximation is reasonable for many situations [45]. A rigorous physical foundation was thus established for the implicit solvent approach. In this paper, we review the general method of using an implicit solvent model to study the drying process of a soft matter solution. A careful implementation of the model has removed certain undesirable effects occurring at the liquid-vapor interface in several previous studies [26,28,29]. We further apply the method to various systems including solution films and droplets. Our results indicate that this approach can be applied to solutions with a good solvent where the solutes are initially well dispersed. Through the current study, the advantages and possible deficiencies of the implicit solvent approach are revealed and summarized.
The paper is organized as follows. In Section 2, the model and simulation method are introduced. The applications of the method to various systems are presented in Section 3. We briefly summarize the method and findings in Section 4.

Model and Simulation Methodology
In an implicit solvent model, the solvent is treated as a uniform, viscous, and isothermal background [38]. The motion of a particle in this background is typically described by a Langevin equation that includes Stokes' drag as a damping term. The damping rate can be chosen according to the effective viscosity of the solvent, the diffusion coefficient of the particle, and the particle size. To model the liquid-vapor interface, a potential barrier is used to confine all the particles in the liquid phase. As shown previously [45], a harmonic potential can be employed as the barrier. The evaporation process of the solvent, during which the liquid-vapor interface recedes, can be mimicked by moving the location of the confining potential's minimum. Below we discuss these two main ingredients of the implicit solvent approach of modeling solution drying.

Langevin Dynamics
In an implicit solvent, the equation of motion of particles is given by the following Langevin Equation [46,47], where m is the mass of the i-th particle, r i is its position vector, t is time, f ij is the force on the particle from its interaction with the j-th particle, F D i is Stokes' drag, and F R i is a random force. Stokes' drag can be expressed as where ξ is the friction (drag) coefficient of the particle. To maintain the system at a constant temperature T, the random force needs to follow the constraint set by the fluctuationdissipation theorem, where · stands for an ensemble average, k B is the Boltzmann constant, δ ij is the Kronecker delta, and δ(t − t ) is the Dirac delta function.
In the dilute limit, the diffusion coefficient D of the particle is related to the friction coefficient ξ through which is known as the Einstein relation. For a small particle with radius R in a flow with a low Reynolds number, Stokes' law states that where η is the viscosity of the fluid. This yields the Stokes-Einstein relation, The friction coefficient ξ is usually written in terms of a damping time Γ as ξ ≡ m/Γ. As a result, Stokes' drag becomes If Stokes' law holds, then Γ = m/(6πηR). If we further assume the particle has a uniform mass density ρ, then m = 4 3 πR 3 ρ and the damping time Γ becomes The implication of this relationship is that for a bidisperse mixture of particles of size ratio α ≡ R l /R s . where R l is the radius of the larger particles and R s is the radius of the smaller particles, the damping time of the larger particles should be α 2 times of that of the smaller ones in order for the Stokes-Einstein relation to hold for both.

Moving Interface Method
When the solvent evaporates from a solution, the liquid-vapor interface recedes. To mimic this process, the location of the minimum of the potential barrier that is used to confine all the particles (in general, solutes) in the liquid solvent is moved toward the solution phase at speed v e . For evaporation at a constant rate, the position of the liquid-vapor interface at time t is given by where H(0) is the initial thickness of a suspension film or the initial radius of a spherical droplet. Similarly, H(t) is the thickness of a drying film or the radius of a drying droplet at time t. For a particle whose center is within distance R from the liquid-vapor interface, where R is a radius parameter, the particle experiences a force normal to the interface. Otherwise, the particle does not interact with the interface. Therefore, the force exerted by the interface on the particle is where k s is a spring constant, ∆z n is the distance from the center of the particle to the instantaneous location of the liquid-vapor interface (which is flat for a film but curved for a droplet), and θ is the contact angle of the solvent on the surface of the particle. Mathematically, ∆z n = z n − H(t), where z n is the particle's coordinate along the z-axis for a flat interface with the bottom of the film at z = 0 or along the radial direction for a spherical interface with the center of the droplet at z = 0. The minimum of the confining potential is thus located at H(t) − R cos θ. In Equation (9), a positive (negative) value indicates that the force is toward the liquid solvent (vapor phase). Some ambiguities exist in the literature regarding the physical interpretation of k s . Pieranski proposed that k s = 2πγ with γ being the interfacial tension of the liquid-vapor interface [48]. However, this expression completely neglects capillary effects. Recently Tang and Cheng [45] analyzed the capillary force exerted on a spherical particle of radius R at a liquid-vapor interface when the particle is displaced out of its equilibrium location and showed that the linear approximation in Equation (9) generally works well but the spring constant should be understood as [45,49] k s = 2πγ ln(2L/R) (10) where L is the lateral span of the liquid-vapor interface. In the simulations reported here, we set k s = 3.0 /σ 2 unless otherwise noted. We also typically set θ = 0, in which case the potential barrier representing the liquid-vapor interface becomes the right half of a harmonic potential. Such a potential ensures that all the solutes are confined in the solution phase. The contact angle θ can also be adjusted to tune the interaction between the solute and liquid-vapor interface. For example, a nonzero θ can be used for systems with attractive solute-interface interactions [50].

Drying of Solution Films of Polymer-Nanoparticle Mixtures
Evaporation of a mixed solution of polymer chains and nanoparticles has been studied via molecular dynamics (MD) simulations with both explicit [15] and implicit [28] solvent models. Here we apply the moving interface method to study a mixture of 3200 linear polymer chains, each consisting of 100 connected beads of mass m and diameter σ, and a varying number of nanoparticles. All the beads interact through a Lennard-Jones (LJ) 12-6 potential, where r is the separation of two beads, sets the strength of interaction, and r c is the cutoff of the potential. To model a good solvent, in this study the LJ interactions are truncated at r c = 2 1/6 σ, rendering the potentials purely repulsive. Alternatively, one can increase the cutoff to include attractive interactions. However in this case the temperature of the solution must be above its theta temperature. Otherwise, phase separation will occur in the initial dilute solution.
Adjacent beads on a chain are connected by a spring described by a finitely extensible nonlinear elastic (FENE) potential, where R o = 1.5σ and k = 30 /σ 2 [51]. The LJ term in the FENE potential is truncated at 2 1/6 σ. With these parameters, the model describes polymer chains in a good solvent with a root-mean-square radius of gyration R g 7.2σ for 100-bead chains. A nanoparticle is modeled as a uniform sphere of LJ beads at a mass density of 1.0 m/σ 3 [52,53]. The interaction between a nanoparticle and a monomer bead on a polymer chain is given by an integrated LJ potential [52], where R is the nanoparticle radius, r is the center-to-center distance between the nanoparticle and monomer, and A np is a Hamaker constant setting the interaction strength. Here nanoparticles with two different radii, R m = 10σ or R s = 2.5σ, are studied. To facilitate the discussion below, we call the former medium nanoparticles (MNPs) and the latter small nanoparticles (SNPs). In both cases, we set A np = 100 /σ 2 . The nanoparticle-polymer interaction is purely repulsive with the potential in Equation (13) truncated at 10.858σ for MNPs and 3.34σ for SNPs. The nanoparticle-nanoparticle interaction is given by an integrated LJ potential for two spheres [52], with where R 1 and R 2 are the radii of the nanoparticles, r is the distance between their centers, and A nn is a Hamaker constant. Here we use A nn = 39.48 [52]. The nanoparticlenanoparticle interaction is also purely repulsive with the potential truncated at 20.574σ for R 1 = R 2 = 10σ or 5.595σ for R 1 = R 2 = 2.5σ. In this manner, the nanoparticles and polymer chains are guaranteed to be well dispersed in the implicit solvent prior to evaporation. The mixed solution of the nanoparticles, consisting of either 171 MNPs or 10,944 SNPs, and polymer chains is placed in a rectangular simulation cell with dimensions of L x × L y × L z , where L x = L y = 200σ and L z = 800σ. Periodic boundary conditions are used in the xy plane. Along the z direction, the solution is confined from below by a wall located at z = 0. The interaction between a particle (either a nanoparticle or a polymer bead) and the wall is given by a LJ 9-3 potential, with w = 2.0 , a being the particle radius, h as the shortest distance between the particle center and the wall. The particle-wall potential is truncated at h c = 0.858a to make the wall purely repulsive. Prior to evaporation, the solution is confined from above by the liquid-vapor interface, which is modeled as a potential barrier as in Equation (9) with H(0) = 800σ. The contact angle θ = 0 is used for both polymer beads and nanoparticles. All the simulations reported here are performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [54]. A velocity-Verlet algorithm with a time step of dt = 0.005τ is used to integrate the equation of motion. All the particles in the system are thermalized at a temperature of T = 1.0 /k B with a Langevin thermostat. The damping time Γ is set to 1τ for the polymer beads, 6.25τ for SNPs, and 100τ for MNPs, respectively. The drying process is modeled by moving the liquid-vapor interface downward to H f = 200σ at a constant speed v e . Simulations are run for v e = 0.024σ/τ and 0.006σ/τ.
After equilibration, each initial solution is a uniform mixture of the polymer chains and nanoparticles. The initial volume fraction is φ n,0 = 0.0224 for both SNPs and MNPs. The initial number density of the polymer beads is ρ p,0 = 0.01σ −3 . The size of a polymer chain can be quantified by its root-mean-square radius of gyration, R g . The 100-bead chains have swollen conformations in the implicit solvent adopted here, similar to those in an explicit LJ solvent [55], and R g 7.2σ. Therefore, R m > R g > R s . The diffusion coefficients of the polymer chains and nanoparticles in each mixed solution are calculated with independent MD simulations. For the polymer chains, the diffusion coefficients, D p 0.00915σ 2 /τ in the mixtures with MNPs and D p 0.00949σ 2 /τ in the mixtures with SNPs. The small difference in D p is likely caused by the different surface area of the nanoparticles in each mixture. At the same volume fraction, the total surface area of SNPs is 4 times that of MNPs. The diffusion coefficient of MNPs is D m 0.0138σ 2 /τ while of SNPs is D s 0.0971σ 2 /τ. Note that even though R m > R g , the computed results show D m > D p , indicating that the diffusion coefficients do not follow the Stokes-Einstein relation.This is not surprising as in the implicit solvent adopted here, the viscous damping is applied to each polymer bead. The polymer dynamics thus follows the Rouse model instead of the Zimm model [56] and a polymer chain freely drained by the implicit solvent cannot be treated as a solid object with a size equal to its R g . With the damping time of Γ p = 1τ for each bead, the diffusion coefficient of a 100-bead chain is estimated to be k B T/(100ξ p ) = Γ p k B T/(100m) = 0.01σ 2 /τ, which is very close to the computed value of D p . For the nanoparticles, D s /D m 7, which is slightly larger than the value ( 5.8) for the same sized nanoparticles in an explicit solvent [38]. Both are larger than the size ratio R m /R s = 4. As discussed in Section 2.1, the diffusion coefficient of a particle following the Langevin dynamics is D = Γ n k B T/m n in the dilute limit, where Γ n is the damping time of the particle and m n is its mass. The expected diffusion coefficient in the dilute limit for SNPs is 0.0955σ 2 /τ, which is close to the calculated value in the polymer-SNP mixture. For MNPs, the diffusion coefficient in the dilute limit is expected to be 0.0239σ 2 τ, which is 1.7 times the computed value in the polymer-MNP mixture.
Using the above results, the Péclet number of each species can be computed as Pe = H(0)v e /D, where H(0) is the initial film thickness, v e is the receding speed of the liquid-vapor interface, and D is the diffusion coefficient. The Péclet number thus characterizes the competition between evaporation-induced migration and particle diffusion. The values of Pe for the polymer-nanoparticle mixtures studied here are summarized in Table 1, where Pe n is the Péclet number for the nanoparticles while Pe p is for the polymer chains.
The snapshots of the mixed solutions undergoing drying are shown in Figure 1 for various times during evaporation and the corresponding density profiles of the nanoparticles and polymer chains are shown in Figure 2. Prior to evaporation, all the chains and nanoparticles are uniformly distributed in the solution. After drying, the polymer is enriched near the descending interface for all four systems studied, which feature two nanoparticle sizes and two evaporation rates. Because of the repulsion with the polymer chains, the nanoparticles are expelled from this skin layer of polymer and become accumulated just below the skin layer. As a result, a stratified state is created, similar to the results previously obtained with an explicit solvent model [15]. In the implicit solvent simulation of Howard et al., similar stratification is also observed with nanoparticles concentrated below a concentrated layer of polymer chains [29]. However, a layer of nanoparticles is also found on top of the polymer skin layer in those simulations (see Figures 6 and 7 of Ref. [29]). This is an effect of the potential barrier used in Ref. [29] to represent the liquid-vapor interface, where the contact angle was set to 90 • for both nanoparticles and polymer chains. The resulting potential is thus attractive for both species and the attraction can be quite strong for nanoparticles, leading to their adsorption at the interface. Since here we use θ = 0 for all the solutes, this effect is removed, as evident in Figure 1 (see the results for the polymer-MNP mixtures).
Although R m > R g , the results in Figure 1 for the polymer-MNP mixtures cannot be classified as "small-on-top" stratification [26][27][28]36,57] as the polymer chains actually diffuse more slowly than the nanoparticles. In each mixed solution, the polymer chains have a larger Péclet number and should be effectively treated as the "larger" species. However, they are always found on top as a skin layer after drying. When the evaporation rate is reduced, the skin layer becomes thicker, signaling enhanced stratification. Such "polymer-on-top" stratification is also found in mixtures with SNPs that have a much smaller radius (see the bottom row of Figure 1) and in other reports [29]. A question naturally arises: Does "polymer-on-top" stratification always occur in all polymer-particle mixtures that undergo a drying process? Or equivalently, can "particle-on-top" stratification be realized in drying polymer-particle mixtures? Although these questions are still open at this point, a careful examination of the results shown in Figure 1 has offered certain clues. In particular, for the polymer-SNP mixtures some SNPs are observed on-top-of the polymer skin layer during drying. Since R g /R s 3 and D p /D s 0.1, the results indicate that the polymer-SNP mixtures may have a tendency of exhibiting "small-on-top" stratification (in this case, "particle-on-top") but are not there yet as the size ratio of 3 is still relatively small. Furthermore, when the evaporation rate is reduced by a factor of 4 from v e = 0.024σ/τ to 0.006σ/τ, the number of SNPs above the polymer skin layer slightly increases, though the skin layer thickens in this case. These results indicate that to realize "particle-on-top" stratification, a even larger ratio between the chain size and the nanoparticle radius may be needed, which can be achieved by using longer chains. A slower evaporation rate may further favor "particle-on-top" stratification. Work along this line will be reported in the future.   The distribution of the nanoparticles and polymer chains shown in Figure 1 can be quantified by their density profiles along the direction of drying, as shown in Figure 2. The local number density of polymer beads is defined as ρ p (z) = n p (z)/ L x L y ∆z , where n p (z) is the number of polymer beads with their z-coordinates between z − ∆z/2 and z + ∆z/2. The local mass density of nanoparticles is defined as ρ n (z) = n n (z)m n / L x L y ∆z , where n n (z) is the partitioned contribution to the slice parallel to the xy plane and from z − ∆z/2 to z + ∆z/2 from all the nanoparticles straddling that slice, weighted by their intersection volume with the slice, and m n is the nanoparticle mass. For the results shown in Figure 2, ∆z = 2σ. To facilitate comparison, the local number density of polymer beads is normalized by their initial average number density (0.01σ −3 ) and the local mass density of nanoparticles is normalized by their initial average mass density (0.0224mσ −3 ). The normalized density profiles (ρ) in Figure 2 reveal several features hidden in the snapshots. First of all, as the evaporation rate is reduced, the peak height of the density profile decreases for both nanoparticles and chains but the range where the density profile exhibits gradients is widened. This indicates that the perturbation caused by the receding interface propagates farther as it take longer time to achieve a certain stage of drying at slower evaporation rates. Secondly, at the same evaporation rate and stage of drying, the peak height as well as the gradient range and magnitude of the polymer density are similar in the two mixtures containing nanoparticles with different diameters but the larger nanoparticles develop a higher peak in their density profile. Furthermore, the gradient of the density profile occurs in a smaller spatial range and as a result, the corresponding density profile has a steeper gradient for the larger nanoparticles. The density profiles in Figure 2 and the snapshots in Figure 1 therefore both point to a more dramatic "polymer-on-top" stratification when the nanoparticles are larger than the polymer chains. Finally as shown in the bottom row of Figure 2, there is an excess of SNPs near the liquid-vapor interface in the equilibrium solution as they are smaller than the polymer chains and can get closer to the interface [36], which leads to a weak peak in the density profile of SNPs above the polymer skin layer during evaporation. The presence of such SNPs near the receding interface is also visible in the snapshots shown in the bottom row of Figure 1.

Drying of Suspension Droplets of Bidisperse Mixtures of Nanoparticles
Stratification has mostly been discussed in the context of drying films [26][27][28][29][34][35][36]58]. However, stratification can also occur in drying droplets under appropriate conditions [59][60][61][62], which may lead to fast procedures of making core-shell structures. In industry, spray drying processes are frequently practiced, where droplet drying is a critical step [63]. In one such process, a particle suspension is injected from a nozzle or an injector into a flowing gas. The liquid jet of the suspension then breaks into droplets, which further dry in the hot gas into clusters of particles (i.e., solutes). The drying of a single droplet was recently studied using the Leidenfrost effect: a droplet is levitated on a hot substrate by its own vapor and then let dry [64]. Here we use the moving interface method to study the drying of a suspension droplet of a mixture of nanoparticles of two different radii, motivated by the possibility of creating a core-shell cluster with one type of nanoparticles in the outside shell while the other type in the inner core. In this context, a bidisperse nanoparticle mixture stratifies radially into either "small-on-outside" or "large-on-outside".
The droplet contains 200 large nanoparticle (LNPs) of radius R l = 20σ and 102,400 SNPs of radius R s = 2.5σ, which are initially confined by a spherical potential barrier inside a sphere with radius H(0) = 2000σ. Their initial volume fractions are φ l = φ s = 0.0002. Although the liquid-vapor interface of a droplet is curved, we still adopt Equation (9) for the particle-interface interaction [65]. We also set the contact angle θ = 0 for both LNPs and SNPs. The evaporation process is mimicked by decreasing the radius of the droplet at a constant speed, v e . The instantaneous radius at time t since the start of the evaporation is thus R d = H(0) − v e t. The nanoparticle-nanoparticle interactions are given by Equation (14) with A nn = 39.48 . All these interactions are purely repulsive with the corresponding potentials truncated at 40.571σ, 23.086σ, and 5.595σ for the LNP-LNP, LNP-SNP, and SNP-SNP pairs, respectively. A Langevin thermostat is used to maintain the temperature of the system at T = 1.0 /k B . To conserve the Stokes-Einstein relation, the damping time is set to 1600τ for LNPs and 25τ for SNPs according to Equation (8). Four independent drying simulations are performed with v e ranging from 1.8 × 10 −2 σ/τ to 1.44 × 10 −4 σ/τ.
The snapshots of the droplets at various stages of drying are shown in Figure 3 and the density profiles of nanoparticles are shown in Figure 4. The density of nanoparticles is defined as ρ(r) = n i m i /[V(r + ∆r) − V(r)], where n i is the number of i-type nanoparticles in a spherical shell from r to r + ∆r, m i is the mass of one i-type nanoparticle, and V(r) = 4 3 πr 3 is the volume of a sphere of radius r. For a nanoparticle straddling multiple shells, its mass is partitioned to each shell in proportion to its intersection volume with that shell. During evaporation, the droplet radius R d decreases and the average density of nanoparticles increases as drying progresses. To facilitate comparison, the local density ρ(r) is scaled as ρ(r) = ρ(r)/β 3 , where β = H(0)/R d . The scaled density profiles are plotted in Figure 4.
)"( %"!' %*$! %#)# !""" #%+ +#$ %*%! %#$( !""" !"" #$" %%"" %$$" After equilibration and prior to evaporation, the droplet is a uniform mixture of LNPs and SNPs (see the first column of Figure 3) with an average total nanoparticle density of 4 × 10 −4 m/σ 3 . The diffusion coefficients are found to be D l = 0.0414σ 2 /τ for LNPs and D s = 0.378σ 2 /τ 2 for SNPs with independent MD simulations. The ratio of D s /D l is about 9.1, which is just slightly larger than the size ratio R l /R s = 8. Since the initial volume fraction of the nanoparticles in the droplet is very low, the computed result of the diffusion coefficient is very close to its value in the dilute limit, which is 0.382σ 2 /τ for SNPs and 0.0477σ 2 /τ for LNPs, respectively. With D l and D s , the Péclet number is then computed for both LNPs and SNPs at a given evaporation rate (v e ) and the results are included in Table 2.  During evaporation, both SNPs and LNPs are first enriched near the liquid-vapor interface, though the degree of enrichment is lower at slower drying rates. At high drying rates (e.g., v e = 1.8 × 10 −2 σ/τ and 1.8 × 10 −3 σ/τ), because of the similar physical mechanism leading to "small-on-top" stratification in drying films of bidisperse particles [26,36,57], SNPs form a concentrated shell near the interface while LNPs are pushed out of this region and into the interior of the droplet. In the final state, a "small-on-outside" cluster is clearly observed (see the third and fourth row of Figure 3 and the bottom row of Figure 4). The simple model discussed here thus points to the possibility of creating core-shell clusters of particles by drying suspension droplets rapidly. Real spray drying processes are of course more complicated with many factors that are not captured by our simple model, such as air invasion, cracking, and buckling [64,66]. Despite these limitations, our results indicate that increasing drying rates may lead to new strategies of controlling the structure of the resulting clusters or creating new structures.
When the drying rate is reduced by lowering v e to 6 × 10 −4 σ/τ (with Pe l = 29 and Pe s = 3.2), "large-on-outside" stratification is observed, as shown in the second rows of Figures 3 and 4. During drying, both SNPs and LNPs are accumulated at the receding liquid-vapor interface but the enrichment of LNPs is more significant. Eventually, LNPs are more enriched at the surface of droplet. When v e is further reduced to 1.44 × 10 −4 σ/τ (with Pe l = 7 and Pe l = 0.76), a uniform distribution of LNPs and SNPs is found in the final dried droplet (see the first rows of Figures 3 and 4). In this case, evaporation is slow and SNPs are almost uniformly distributed in the droplet in the entire process of drying. LNPs are first accumulated near the receding interface in the early stage of drying but eventually become uniformly distributed too. The simulation results thus demonstrate that the distribution of nanoparticles transitions from uniform to "large-on-outside" to "small-on-outside" as the evaporation rate is increased, which is consistent with the prediction of the diffusive model of stratification in drying films proposed by Zhou et al. [57].

Drying of Solution Droplets of a Polymer Blend
The moving interface method can also be used to study the drying process of polymer solution droplets. In this section we focus on the solution of a polymer blend, which is a symmetric mixture of short polymer chains of type A and B, with 2048 chains of each type. Each chain is linear and consists of 12 beads of mass m that are bonded by the FENE potential in Equation (12) [51]. In addition to the bonded interaction, all the nonbonded pairs of beads interact via a LJ potential with a cutoff distance of 2 1/6 σ, i.e., a purely repulsive potential. The systems are kept at T = 1.0 /k B via a Langevin thermostat with a damping time of Γ = 10τ. As a result, at low concentrations the polymer chains adopt swollen conformations as in a good solvent. To model an incompatible polymer blend, the strength of the self interaction is set to AA = BB = while the cross repulsion is stronger with AB = 2.0 . The critical value of AB for a symmetric blend to phase separate depends on its density. For a melt of a symmetric mixture of linear N-bead chains at a density of 0.85 m/σ 3 , Grest et al. found that it phase separates into an ordered phase at AB (1 + 3.4/N) [67]. Therefore, the mixtures studied here are expected to phase separate as the packing density of the systems approaches the melt density.
In the initial state prior to drying, all the chains are uniformly dispersed in a sphere with radius 100σ, as shown in Figure 5a. All the polymer beads are confined in this sphere with a spherical potential barrier as described in Section 3.2 with H(0) = 100σ. The potential given in Equation (9) with k s = 3.0 /σ 2 and R = 2σ is used for the interaction between a polymer bead and the liquid-vapor interface. The initial density of each droplet is about 0.012 m/σ 3 . During drying, the radius of the droplet is reduced to H f = 24σ at a rate of v e = 1.52 × 10 −2 σ/τ. At this final radius, the packing density of the polymer beads is about 0.85 m/σ 3 , which is the melt density of the blend under an external pressure of about 5 /σ 3 [67]. The entire drying process thus lasts 5000τ. Figure 5. (a) A solution droplet of a polymer blend prior to drying. Depending on the interaction with the receding liquid-vapor interface, polymeric particle with different morphologies are obtained after drying: (b) A Janus particle is produced with θ = 0 for both polymer A (red) and B (blue) after drying for 5 × 10 3 τ followed by relaxation for 2 × 10 5 τ; (c) A core-shell particle is produced with θ = 0 for polymer A (red) while θ = π/2 for polymer B (blue) after drying for 5 × 10 3 τ followed by relaxation for 5 × 10 4 τ. For clarity, the droplets after drying and relaxation are sliced through the center to show the interior. Figure 5 shows that the morphology of the resulting polymeric particle depends on the drying conditions and the polymer-interface interaction controlled by the contact angle θ in Equation (9). In Figure 5b, θ = 0 for both polymers. The two polymers phase separate after fast drying as the cross repulsion (A-B) is sufficiently stronger than the self repulsion (A-A and B-B) in the two components [67]. Domains of each component are clearly visible in the polymeric particle. After relaxation, a Janus particle is produced with each component occupying half a sphere.
In Figure 5c, θ = 0 for polymer A while θ = π/2 for polymer B, indicating that the liquid-vapor interface is repulsive for polymer A but attractive for polymer B. The latter thus has a tendency to be adsorbed at the liquid-vapor interface. As a result, the B-type chains are enriched at the droplet surface during drying while all the A-type chains are pushed into the interior of the sphere. In the final state, with a radius of 24σ, a core-shell distribution can be clearly identified. The core-shell structure remains completely stable during the subsequent relaxation period. The results thus indicate that the solute-interface interaction is an important factor affecting the structures produced by drying [50] and can be used to produce polymeric particles with different morphologies.

Drying of Solution Droplets of Diblock Copolymers
Block copolymers can also be employed to produce structured polymeric nanoparticles [68]. For a bulk system of a block copolymer with incompatible blocks, there are well-known ordered structures such as lamellar, cylindrical, and spherical phases depending on the fraction of the two components on the chain [69]. In this section, we use a setup almost identical to the one in Section 3.3 to study the drying of a solution droplet of symmetric diblock copolymer chains. Here, the number ratio between monomers in the A blocks and those in the B blocks is 1:1 and we vary the chain length N from 12 to 96. The total number of monomers is fixed at 49,152 and the number of chains is thus 49,152/N. The bonded interaction is given by the FENE potential in Equation (12). All the monomers, if not directly bonded, interact with each other via a LJ 12-6 potential as in Equation (11). The nonbonded interaction is purely repulsive. The interaction strength in the LJ potential is for the A-A and B-B pairs. The strength of A-B interactions is given by AB , which is varied from 2 to 8 . A Langevin thermostat with a damping time of Γ = 10τ is used to control the temperature at 1.0 /k B . In the initial state of each solution (e.g., see Figure 6a for N = 24), all the polymer chains are confined in a sphere with radius 100σ by a spherical potential barrier described by Equation (9) with k s = 3.0 /σ 2 and R = 1.0σ. During a drying period that lasts 5000τ, the radius of the droplet is reduced to H f = 24σ, yielding a packing density of 0.85 m/σ 3 for the polymer beads in the final droplet.
Similar to the case of a polymer blend, the results in Figure 6 for N = 24 show that the structure of the resulting polymeric particle of diblock chains depends on the blockinterface interaction and the strength of the block-block repulsion. The polymer-interface interaction is determined by the contact angle θ in Equation (9). In Figure 6b,c, θ = 0 for both blocks and the interface therefore appears neutral for all the monomers. The value of AB at which the order-disorder transition (ODT) occurs depends on the total density of the system. Grest et al. [67] found that for a symmetric diblock with N = 20 at a density of 0.85 m/σ 3 , ODT occurs at AB 3.9 . In the current systems with N = 24, the critical value of AB is expected to be slightly smaller. As shown in Figure 6b, with AB = 2.0 , only small domains of each type of blocks are observed after drying. No significant growth of the domain size is observed after relaxation. When the cross repulsion between the two blocks is increased to AB = 8.0 , bringing the system deeply into the ordered-phase region in the phase diagram [69], the blocks start to aggregate during drying, as shown in Figure 6c. After relaxation, stripes of different blocks are clearly visible because of the incompatibility of the blocks, resembling the lamellar phase in bulk diblock copolymers [69]. (blue) blocks is produced with θ = 0 for both blocks and AB = 2.0 ; (c) A particle with stripes of A and B blocks is produced with θ = 0 for both blocks and AB = 8.0 ; (d) A particle with a mixture core enclosed by a bilayer shell (i.e., a layer of B blocks on the outside and a layer of A blocks on the inside) is produced with θ = 0 for block A (red), θ = π/2 for block B (blue), and AB = 8.0 . In each case, the droplet is dried to a radius of 24σ during 5 × 10 3 τ, followed by relaxation for 2 × 10 5 τ.
For clarity, the droplets in (b-d) are sliced through the center to show the interior.
In Figure 6d, nonneutral block-interface interactions are adopted with θ = 0 for block A while θ = π/2 for block B. As a result, monomers in the B blocks have a tendency to adsorb at the liquid-vapor interface while the A blocks are repelled by the interface. After drying, the resulting particle is enveloped by a layer of B blocks adsorbed at the liquid-vapor interface, with a layer of A blocks inside that is bonded to the B blocks at the surface. The core region of the final particle is filled with a mixture of diblock chains that tend to phase separate. In a larger drop filled with longer diblock chains, an onion-like particle with layers of different blocks alternating radially is expected to form. Studies on this and other interesting morphologies will be reported in the future. Figure 7 shows the effect of chain length (N) on the particle morphology, while the length ratio of the two blocks is kept at 1:1. For all these systems, AB = 2.0 and θ = 0 for both blocks, indicating that the interface is neutral and repulsive for both components. After fast drying to H f = 24σ, the resulting polymeric particle has a disordered structure at N = 12 but phase separation occurs for longer chains, resulting in patchy polymeric articles. After a relaxation process in which the radius of the confining spherical potential is fixed, rough stripe-like domains are formed in the polymeric particles with each domain dominated by one type of blocks. Each domain can be regarded as a patch. The number of patches is thus reduced in the relaxation process. Furthermore, the number of patches decreases as N is increased, as shown in Figure 7. This can be understood by noticing that the domain size is larger when the block length is longer. As N increases, the critical value of AB at which ODT occurs for the diblock copolymers decreases. For example, at a melt density of 0.85/σ 3 , ODT occurs at AB = 1.85 for N = 40 and 1.28 for N = 100 [67]. Since we fix AB at 2.0 , the system enters and moves deeper into the ordered region of the phase diagram as the block size increases. This explains the trend that the stripes and patches grow larger and their boundaries sharpen as N increases. For N = 96, the final polymeric nanoparticle after relaxation has a surface pattern that features just a few large stripes and patches. Our results indicate that the chain length is a useful parameter to tune for controlling the surface pattern, including the number of patches, of a polymeric particle produced via drying a solution droplet of diblock copolymers. Experimentally, the drying of solution droplets discussed above can be realized if a solution is broken into droplets that are suspended and dried in air [68]. Another technique that can be used to produce polymeric particles or nanoparticle clusters from solution droplets is flash nanoprecipitation (FNP) invented by Johnson and Prud'homme [70,71], in which a solute or a mixture of solutes (e.g., drug molecules, nanoparticles, and polymers) is first dissolved in a solvent and then the solution is rapidly mixed with a non-solvent for the solute(s). In this process, the solution is broken into small droplets, which are dispersed in the non-solvent. As the solvent and non-solvent are miscible, the solvent is quickly extracted from the droplets, which shrink rapidly, and the solutes are compressed into particles or clusters by the surrounding non-solvent. This process is quite similar to the drying of a solution droplet. During drying, the solvent leaves the droplet via an evaporation process. In FNP, the solvent leaves the droplets via diffusion into the nonsolvent. In both cases, the surface of the droplet recedes, which is either a liquid-vapor interface or an interface between the droplet and the non-solvent, and the droplet shrinks. Therefore, it is not surprising that polymeric particles with morphologies similar to some of those in Figures 5 and 6 were produced using FNP [72]. The moving interface method discussed in this paper can thus be applied to droplets undergoing FNP and to address questions such as how the mixing rate of the solvent and non-solvent affects the structure of the resulting particles.

Summary and Discussion
In this paper, we have reviewed the method of modeling particle suspensions, polymer solutions, and their mixtures using an implicit solvent with the liquid-vapor interface mimicked by a potential barrier confining all the solutes in the solvent. Their drying process can be studied with the moving interface method, in which the location of the liquid-vapor interface, i.e., the location of the confining potential barrier's minimum, is moved in a prescribed manner. The evaporation rate can be tuned by varying the speed at which the interface is moved. Various evaporation patterns, including drying films and droplets, can be realized with an appropriate choice of the way in which the interface (i.e., the equipotential surface of the confining potential) is moved. For example, the interface is flat and is translated along its perpendicular direction when a film is dried while it is spherical and is shrunk radially for a drying droplet.
With the moving interface method, we have studied the drying behavior of a mixed solution of polymer chains and nanoparticles, a suspension droplet of a bidisperse mixture of nanoparticles, a solution droplet of a polymer blend, and a solution droplet of diblock copolymer chains. A rich set of structures are formed after drying, including stratified films, core-shell clusters (i.e., radially stratified clusters) of nanoparticles, Janus polymeric particles, core-shell polymeric particles, and patchy polymeric particles. These structures are consistent with those observed previously in explicit solvent simulations and experiments, indicating that the moving interface method with an implicit solvent model can be used in certain situations to yield realistic results. Since the solvent is not treated explicitly, such method has the advantage of significantly reducing the number of particles needed to describe a physical system and thus allows the modeling of much larger systems over longer times and at slow evaporation rates that may be directly comparable with those used in experiments.
Caution needs to be taken about where the moving interface method can be applied. In this method, the solvent is treated as a uniform, viscous, and isothermal medium. Therefore, the solvent remains as a background during drying and does not exhibit any flow. The moving interface method is thus only applicable to situations when the solvent flow is not a crucial factor. For drying films, spherical droplets, and cylindrical droplets, this condition can be satisfied and the systems can be modeled with the moving interface method. Examples include a liquid film generated by dip coating [73], respiratory droplets suspended in air [17], solution droplets created by a spray drying process [74], and a polymer fiber during electrospinning [75]. However, for systems in which the solvent develops flow patterns during drying, the moving interface method cannot be directly employed. One example is the drying of a sessile droplet on a substrate, in which a capillary flow emerges during evaporation, transporting solutes to the edge of the droplet and leading to the famous coffee-ring effect if the peripheral of the droplet is pinned [39,76]. Since the capillary flow is not captured by an implicit solvent model, it is impossible to produce the coffee-ring deposits with the simplest moving interface method. However, in these situations the moving interface method can be combined with other techniques such as lattice-Boltzmann method that is able to describe a flow field to study the drying process [77]. Furthermore, there are still many scenarios in which the flow field of the solvent is not a dominant factor and the moving interface method discussed here can be useful because of its computational efficiency.
In the current implementation of the moving interface method, the liquid-vapor interface is moved to simulate the drying process but its shape remains unchanged. Therefore, the method cannot be directly applied to systems where the liquid-vapor interface develops instabilities during evaporation. For example, when a spin-coated polymeric film is dried, it may develop Rayleigh-Bénard-Marangoni convective instabilities that cause the evolution of surface morphology of the film [42]. To capture these effects, the moving interface method needs to be extended to include a model for the liquid-vapor interface that allows the interface to deform and exhibit various instabilities during the drying process.
Another potential issue with the moving interface method is regarding the dynamics of polymer chains in an implicit solvent. When the Langevin dynamics is applied to each monomer, a polymer chain essentially follows the Rouse dynamics with the implicit solvent effectively draining through the chain. It is understood that in a dilute polymer solution, the chain dynamics are better described by the Zimm model, where hydrodynamic interactions make a chain and the solvent in its pervaded volume to behave as a solid object moving through the surrounding solvent [56]. The previous simulations by Statt et al. have revealed that an implicit solvent model of polymer solutions ignoring hydrodynamic interactions could yield outcomes of evaporation even qualitatively inconsistent with those from an explicit solvent model that is matched to the implicit one as far as the equilibrium solution properties are concerned [34]. To address this issue, new implicit solvent models need to be developed to yield polymer dynamics matching those in an explicit solvent. Such models can then be used with the moving interface method to study the drying process of polymer-containing soft matter solutions and produce results that agree with those from experiments and explicit-solvent simulations.