Influence of hydrodynamic interactions on stratification in drying mixtures

Coatings or films formed by drying^1,2 are relevant to many technologies, including latex paints,^3,4 inkjet printing,⁵ and polymer nanocomposites.⁶ Such films are often comprised of multiple components including nanoparticles or colloids,⁷ polymers,⁸ and surfactants.^4,9 It is desirable to control the distribution of components within the dried film, e.g., to ensure the dispersion of nanoparticles^6,10 or to impart antimicrobial properties to a surface.¹¹ Recent experiments showed that quickly dried mixtures of small and large spherical colloids can vertically stratify into layers.^12–14 Unexpectedly, the smaller colloids accumulated at the solvent-air interface and pushed the larger colloids down into the film. Computer simulations qualitatively agreed with the experiments^12–15 and formed the same small-on-top stratification in ternary and polydisperse colloid mixtures¹⁶ and in polymer–polymer and colloid–polymer mixtures.¹⁷ The stratification was more pronounced for larger size ratios and faster drying speeds, although it was found to persist even at moderate drying rates.^13,15,17

Stratification has been modeled theoretically as a diffusion process^7,18 characterized by the film Péclet number for each component, Pe_i = $Hv$/D_i, where H is the initial film height, $v$ is the speed of the solvent–air interface, and D_i is the diffusion coefficient of component i. When Pe_i ≪ 1, diffusion dominates over advection from drying, and a nearly uniform distribution of component i is expected. When Pe_i ≫ 1, component i accumulates toward the drying interface. Trueman et al. developed a model for drying colloid mixtures¹⁹ that qualitatively predicted a maximum in the number of large colloids on top of the film when Pe₁ < 1 and Pe₂ > 1 (here 1 is the smaller colloid and 2 is the larger colloid).^20,21 However, their model does not predict the small-on-top stratification that occurs when Pe₁ > 1 and Pe₂ > 1.^19,20 Fortini et al. hypothesized that such stratification resulted from an osmotic pressure imbalance,^12,16 while Howard et al.^15,17 and Zhou et al.²² proposed models using the chemical potential as the driving force. Computer simulations of polymer mixtures and the Howard et al. model agreed quantitatively,¹⁷ and qualitative agreement was found between the Zhou et al. model and experiments for colloid mixtures in certain parameter regimes.^14,21 The Zhou et al. model is equivalent to the low-density limit of the Howard et al. model.¹⁵ 

For computational and theoretical convenience, many previous studies have employed a free-draining approximation for hydrodynamic interactions to model stratification. This approach incorporates the Stokes drag on each particle but neglects other interactions. The computer simulations^12–17 implicitly modeled the solvent using Langevin or Brownian dynamics methods that treated the solvent as a quiescent viscous background, while the theoretical models^{12,15–17,22} assumed that the particle mobilities were given by the Einstein relation. However, Sear and Warren pointed out in a recent analysis²³ that neglecting hydrodynamic interactions overpredicts the phoretic velocity of a single large colloid in an ideal polymer solution, which must be proportional to the extent of stratification in this simple model. They argued that this discrepancy is caused by the omission of solvent backflow,²³ leading to an incorrect hydrodynamic mobility.²⁴ Their argument is qualitatively supported by differences between the experiments and implicit-solvent simulations of Fortini et al.¹² To our knowledge, though, there has been no direct test of how hydrodynamic interactions influence stratification at finite concentrations or size ratios between components, which are more challenging to analyze theoretically.^23,24

In this article, we systematically demonstrate the influence of hydrodynamic interactions on stratification using a drying polymer mixture as a model system. Computer simulations are ideal tools for such a study because, unlike in experiments, it is possible to artificially remove hydrodynamic interactions between particles in the simulations. Molecular dynamics simulations of drying were performed using both explicit and implicit solvent models. We show that although the implicit-solvent model stratifies for the drying conditions investigated, the explicit-solvent model does not due to the presence of hydrodynamic interactions. Our simulations show that it is critical to incorporate hydrodynamic interactions into models and simulations in order to reliably predict stratification in drying mixtures.

We performed nonequilibrium molecular dynamics simulations of drying polymer mixtures using two bead–spring models: one with an explicitly resolved solvent (Sec. II A) and a corresponding one with an implicit solvent (Sec. II B). While the explicit model included full hydrodynamic interactions, the implicit model only incorporated the Stokes drag on each polymer, thus omitting hydrodynamic interactions between polymers. We have chosen to investigate the influence of hydrodynamic interactions on stratification using a polymer mixture rather than the colloid mixture originally studied by Fortini et al.¹² because the polymer mixtures are more computationally amenable to the analysis in Sec. III. However, the results and conclusions presented here extend to other types of drying mixtures. Details of both polymer models are presented in the following.

Polymers were represented as linear chains of monomers with diameter σ and mass m immersed in an explicit solvent of beads of equal size and mass. All monomers and solvent particles interacted with the Lennard-Jones potential,

where ε is the interaction energy and r is the distance between particles. The monomer–monomer interactions were made purely repulsive by truncating the potential at its minimum and shifting it to zero.²⁵ The solvent–solvent and monomer–solvent interactions included the attractive part of U_nb and were truncated at 3.0 σ with the energy and forces smoothed to zero starting from 2.5 σ. Bonds between monomers were represented by the finitely extensible nonlinear elastic potential,²⁶ 

with the standard parameters r₀ = 1.5 σ and κ = 30 ε/σ². This choice prevents unphysical bond crossing because the equilibrium bond length is approximately 0.97 σ at a temperature of T = 1.0 ε/k_B, where k_B is Boltzmann’s constant. This model corresponds to good solvent conditions for the polymers.

The drying process was simulated using the method developed by Cheng and co-workers.^27–29 The polymer solution was supported by a smooth structureless substrate modeled by a Lennard-Jones 9–3 potential at the bottom of the simulation box,

where z is the distance between the particle and the substrate and ε_w = 2 ε is the strength of the interaction. Interactions were truncated for z > 3.0 σ. The solvent vapor above the liquid film was confined by an additional potential of the same form as U_w at the top of the simulation box, which was made purely repulsive by truncating it for z > (2/5)^1/6σ. The solvent was evaporated by deleting a small number of randomly chosen solvent particles from the top 20 σ of the vapor. To maintain temperature control, monomers and solvent in a slab of height 20 σ above the substrate were weakly coupled to a Langevin thermostat with T = 1.0 ε/k_B and friction coefficient 0.1 m/τ,^30–32 where $τ=mσ2/ε$ is the unit of time.

All simulations were performed using the HOOMD-blue simulation package (version 2.2.2) on multiple graphics processing units^33,34 with a time step of Δt = 0.005 τ. The solvent coexistence densities were 0.664 σ⁻³ and 0.044 σ⁻³ for the liquid and vapor phases, respectively. The simulation box was periodic in the x and y directions with a length of L = 50 σ. The height of the box was L_z = 540 σ, where H = 500 σ was the initial film height and the remaining 40 σ was filled with solvent vapor. We investigated a mixture of short polymers of M₁ = 10 monomers and long polymers of M₂ = 80 monomers at an initial monomer density of roughly 0.06 σ⁻³ each. The simulations consisted of 730 000 solvent particles, 900 long chains, and 7100 short chains, resulting in a total of 873 000 particles. By deleting one solvent particle every 1.25 τ, we obtained an evaporative flux of 3.2 × 10⁻⁴/σ²τ. The drying rate for this model solvent is estimated to be considerably faster than for solvents like water that are often used in experiments, but is close to the slowest rate that was computationally accessible to us. However, the simulated films were also much thinner than in experiments, resulting in realistic Péclet numbers.

The implicit-solvent model was constructed by matching the structure of the polymer chains in the explicit solvent at infinite dilution. Interactions between the monomers and the solvent swelled the polymers beyond the size of a chain with only the monomer–monomer interactions, so additional effective interactions were required. We measured the distance distribution between monomers separated by two bonds along a chain surrounded by solvent and a chain without solvent. We fit the negative logarithm of the ratio of the explicit and implicit distribution with a spline potential,³⁵ 

The fit resulted in parameters n = 3.8, ε_s = 0.68 ε, r_s = 1.32 σ, and r_c = 2.0 σ. This potential, modeling the soft repulsion between monomers due to the solvent, was added to the bare monomer–monomer interactions. We validated this effective model by measuring the polymer structure over a range of concentrations and compositions up to a total monomer number density of 0.375 σ⁻³, finding overall good agreement between the implicit and explicit models (see the supplementary material).

The polymer long-time dynamics in the explicit solvent were matched in the implicit solvent using Langevin dynamics simulations.^30–32 This technique incorporates the effects of Brownian motion and Stokes drag from the solvent, but neglects hydrodynamic interactions between monomers. The monomer friction coefficients were adjusted for each polymer to give the same polymer center-of-mass diffusion coefficient at infinite dilution as in the explicit solvent. This approach gives the correct long-time dynamics, but distorts the internal relaxation modes of the polymers.¹⁷ We measured diffusion coefficients in the explicit solvent of D₁ = 0.0342 σ²/τ for the M₁ = 10 polymers and D₂ = 0.0071 σ²/τ for the M₂ = 80 polymers from the polymer center-of-mass mean squared displacement in a cubic box with edge length L = 50 σ, giving 2.92 m/τ and 1.76 m/τ for the friction coefficients, respectively (see the supplementary material for details).

The liquid–vapor interface was modeled by the repulsive part of a harmonic potential.^15,17 (The complete form of the potential can be found in Ref. 15.) In order to closely match the explicit-solvent simulations, the spring constant and position of the interface were adjusted to obtain similar initial density profiles to the explicit-solvent model. We used spring constants of 0.1 ε/σ² and 0.18 ε/σ² for the monomers of the short and long chains, respectively. The position of the interface was measured throughout the explicit-solvent simulations and used directly in the implicit-solvent simulations. The minimum of the potential was offset for the long polymers by −1 σ.

We performed 25 drying simulations for both models and measured the monomer density profiles in the film. Figure 1 shows the average profiles at three times: before evaporation, at a film height of roughly 3H₀/4, and at a film height close to H₀/2. We used the average interface speed and the diffusion coefficients at infinite dilution to estimate the film Péclet numbers of Pe₁ ≈ 7 and Pe₂ ≈ 33. (Throughout this discussion, 1 denotes the shorter polymer, while 2 is the longer polymer.) Initially, the polymers were nearly uniformly distributed in the film. When drying began, both the long and short polymers accumulated at the moving interface, as expected from their Péclet numbers. However, there was a noticeable difference in the distribution of chains within the film between the two models. Consistent with the previous simulations and the Howard et al. theoretical model,¹⁷ the implicit-solvent model stratified with the shorter polymers on top of the longer polymers [Fig. 1(b)]. In stark contrast, the explicit-solvent model showed no small-on-top stratification [Fig. 1(a)] and in fact more long polymers than short polymers accumulated immediately below the liquid–vapor interface.

The qualitatively different behavior of the two models suggests a significant difference in the relative migration speeds of the polymer chains, which was previously shown to give rise to stratification for the implicit-solvent model.^12,15,17,22 At isothermal conditions, the relative migration speed Δu predicted by the Howard et al. model¹⁵ is

where u₁ is the migration speed of the short polymers, μ_i is the chemical potential of component i computed from the equilibrium free-energy functional, and $Mi$ is the effective mobility for component i. We previously estimated $Mi$ from the equilibrium diffusion coefficient D_i using the Einstein relation, $Mi=Di/kBT$⁠, within the free-draining approximation of hydrodynamics.^15,17 The Howard et al. model was previously shown to give quantitative agreement with implicit-solvent simulations of a polymer mixture.¹⁷ Stratification of the short polymers requires that Δu > 0, which occurs when the ratio of the chemical potential gradients is larger than the ratio of mobilities. Our simulations and Eq. (5) then suggest several possible explanations for the different drying behavior:

• temperature gradients from evaporative cooling affect μ_i or the migration velocities in the explicit-solvent model, but such gradients are absent from the implicit-solvent model,

• differences in the chemical potentials of the two models lead to different diffusive driving forces,

• differences in the equilibrium diffusion coefficients change $Mi$ for the two models, and/or

• hydrodynamic interactions mediated by the explicit solvent alter the migration velocities.

We probed each of these effects in turn. As will be shown, we found that only hydrodynamic interactions (D) had a significant enough effect to explain the lack of stratification in the explicit-solvent model.

Solvent evaporation can lead to cooling at fast drying rates.²⁷ For our explicit-solvent model, the temperature is expected to be fixed at T = 1.0 ε/k_B by the Langevin thermostat close to the substrate and, at pseudo-steady state, to decrease linearly in distance from the substrate with a slope proportional to the evaporative flux. The temperature profile T(z) measured in the simulations, shown in Fig. 2(a), is consistent with this expectation. Here, the local temperature T(z) is defined by

where the sum is taken over the N_k particles in a slab k centered around z with thickness 2Δz = 1 σ and $vj$ is the speed of particle j. On the other hand, there were no temperature gradients in the implicit-solvent model because the solvent was treated as an isothermal viscous background.

Temperature gradients can give rise to mass flux (i.e., the Soret effect),³⁶ and the mobility is temperature dependent through the viscosity.³⁷ It has recently been shown that temperature gradients can influence the stratification of nanoparticles.³⁸ In order to test how the presence of temperature gradients influenced stratification in the explicit-solvent model, we removed the Langevin thermostat and instead coupled all particles to a dissipative particle dynamics (DPD) thermostat^32,39,40 to eliminate temperature gradients and performed a simulation at the same drying rate. The DPD thermostat applies pairwise random and dissipative forces to all particles and so conserves momentum and preserves hydrodynamic interactions. We set the DPD friction coefficient to be 0.5 m/τ, which has been shown to have a negligible effect on the viscosity of a fluid of nearly hard spheres.⁴⁰ Figure 2(a) shows that although there was originally roughly 12% evaporative cooling in the explicit-solvent model, the DPD thermostat maintained a constant temperature throughout the film.

We then compared the distribution of polymers in the film with and without evaporative cooling at the same drying rate. The temperature gradient in the film led to a corresponding gradient in the local total density (solvent and polymer), which was higher in colder regions [Fig. 2(b)]. This in turn decreased the total film height compared to the isothermal case for an equal amount of the solvent evaporated. To account for this difference, we compared the monomer density profiles at equal film height. The monomer density profiles for the M₁ = 10 and M₂ = 80 polymers are shown in Figs. 2(c) and 2(d) at film height 360 σ, when the temperature gradient was most pronounced. These monomer density profiles were indistinguishable with and without evaporative cooling. Accordingly, we excluded temperature effects as a possible reason for the different stratification behavior in the implicit- and explicit-solvent models.

In the absence of pressure and temperature gradients, diffusive flux is driven by chemical potential gradients within the regime of linear response.⁴¹ A mismatch in the density or composition dependence of the chemical potentials could result in different stratification behavior for the explicit- and implicit-solvent models. We accordingly measured the relevant chemical potentials for both models using test insertion methods. We first generated trajectories of polymer mixtures in cubic periodic simulation boxes over a wide range of mixture compositions (see the supplementary material), saving 100 independent configurations for each composition.

For the explicit-solvent mixtures, the chemical potential of the solvent μ₀ was calculated by Widom’s insertion method,⁴² 

where β = 1/k_BT, $μ0id$ is the chemical potential of an ideal gas of particles at the same temperature and density as the solvent in the mixture, and ΔU is the change in the potential energy on insertion of a test particle. The ensemble average is taken over configurations and insertions. We performed 5 × 10⁵ insertions per configuration, which we found to be sufficient to obtain a converged value for μ₀.

The chemical potential of the polymers was estimated by the chain increment method⁴³ for both the explicit- and implicit-solvent models. The incremental chemical potential μ⁺(M) to grow a chain of length M + 1 from a chain of length M is

where the ensemble average is now taken over configurations and test insertions onto the ends of chains of length M. The incremental chemical potential μ⁺(M) converges to a constant value μ⁺(∞) for sufficiently large M.⁴³ The greatest deviations from μ⁺(∞) occur for a single monomer, μ⁺(0), which can instead be measured by Widom insertion.⁴² The chemical potential of the chain was then approximated as⁴⁴ 

where $μiid$ is the chemical potential of an ideal gas of chains of length M_i at the same chain density and temperature as in the mixture.

In order to improve convergence of the ensemble averages, test particles were inserted onto the chain ends at positions consistent with the bond length distribution.^44,45 The test particle position was drawn with random orientation relative to the chain end at distance r distributed according to r² exp(−βϕ(r)) in the range 0.7 σ ≤ r ≤ 1.3 σ, where ϕ(r) is the interaction potential between bonded monomers given by Eqs. (1) and (2). An additional analytical contribution to the chemical potential was then included to account for the weighted sampling.⁴⁴ To measure μ⁺(10) and μ⁺(80), we performed 100 trial insertions per chain end in each configuration and, finding that μ⁺(10) ≈ μ⁺(80) within statistical uncertainty, took μ⁺(∞) ≈ μ⁺(10). The monomer excess chemical potential, μ⁺(0), was measured using 5 × 10⁵ random insertions per configuration, as for the solvent.

In our implicit-solvent theoretical model [Eq. (5)],^15,17 the diffusive driving force was given by the gradient of μ_i. However, this expression must be modified to account for the presence of the solvent in the explicit model. If the mixture is approximately incompressible, the solvent and polymer concentrations are not independent. Assuming that the volume occupied by one polymer chain of length M_i is roughly equal to M_i solvent particles, the chemical potential μ_i is replaced in Eq. (5) by the exchange chemical potential Δμ_i = μ_i − M_iμ₀.^46,47 For the implicit model, the solvent is effectively incorporated into the free energy of the polymers, giving Δμ_i = μ_i and recovering Eq. (5), as expected. Batchelor showed that the diffusion process resulting from gradients of Δμ_i is equivalent to sedimentation of the solute (polymers) in a force-free solvent.^47,48

Although Δμ_i was determined across a wide range of mixture compositions (Fig. 3), we found that it was essentially only a function of the total monomer density ρ_m. As expected, Δμ_i increased for larger ρ_m. This is in qualitative agreement with our analysis of stratification in mixtures of hard-sphere chains, where the excess chemical potential per monomer depended primarily on the total monomer packing fraction and more weakly on the total chain number density.¹⁷ Most importantly, Δμ_i is in excellent agreement between the two models for both the short and long polymers. We then expect both the explicit- and implicit-solvent models to give equivalent driving forces for diffusion and stratification.

Having found good agreement in the exchange chemical potential between the two models, we considered dynamic effects in Eq. (5). The polymer mobility $Mi$ relates the diffusive driving force on the polymers to the migration velocity. In our theoretical model, we assumed that $Mi$ could be estimated from the equilibrium diffusion coefficient D_i. In our simulations, D_i for the implicit-solvent model was matched to the explicit-solvent model in the limit of infinite dilution, but the models may deviate at higher concentrations. Differences in the concentration-dependence of D_i, and hence $Mi$⁠, may accordingly influence the relative migration speeds of the components during drying.

We measured D_i for both the explicit- and implicit-solvent models across a wide range of concentrations and compositions from the mean-squared displacement of the polymer centers of mass. The values of all measured diffusion coefficients are available in the supplementary material. At low total monomer densities, the agreement between the two models was good, as expected from the model fitting, with larger discrepancies at higher monomer densities. Over the range of compositions relevant to the drying simulations, the agreement between the two models was quite good, with a deviation of at most 20% for the short polymers and 30% for the long polymers relative to the explicit-solvent model.

According to Eq. (5), a larger value of $M2/M1$⁠, and so approximately D₂/D₁, should result in a larger difference in migration velocities for the same chemical potential gradient and lead to stronger stratification. For all concentrations and compositions investigated, D₂/D₁ was larger for the explicit-solvent model (Fig. 4). This suggests that stronger stratification might be expected for the explicit model on the basis of the diffusion coefficients, which is the opposite of the observed behavior. Accordingly, we concluded that differences in the equilibrium diffusion coefficients were not likely the reason for the lack of stratification in the explicit-solvent model.

Given that temperature gradients, chemical potentials, and equilibrium diffusion coefficients do not explain the difference in stratification in the implicit- and explicit-solvent models, it remains now to test the influence of hydrodynamic interactions. Hydrodynamic interactions are inherently present in the explicit-solvent model, but are absent from the implicit-solvent Langevin dynamics simulations. The recent analysis by Sear and Warren demonstrated how solvent backflow reduces the phoretic motion of a large colloid in an ideal polymer solution.²³ Brady showed from a rigorous microscopic perspective how this fundamentally results from a neglect of hydrodynamic interactions, which modify the mobility.²⁴ 

We performed a simple simulation to test the effect of hydrodynamic interactions on polymer migration. We constructed a simulation box with L = 40 σ and L_z = 80 σ. Half of the box (z > 0 σ) was filled with pure solvent, while the other half (z < 0 σ) was filled with a solution of long and short polymers. A semipermeable membrane at z = 0 σ separated the two compartments, and walls were placed at z = ±40 σ. Interactions with the walls were given by Eq. (3) truncated at distances greater than 3 σ with ε_w = 1 ε, while the monomer interactions with the membrane were modeled by Eq. (3) truncated at (2/5)^1/6 σ with ε_w = 1 ε. The solvent was allowed to exchange across the membrane, coming to equal chemical potential, but the polymers were confined to z < 0 σ. For the explicit-solvent model, there were 78 720 solvent particles in the box and the initial total monomer density of the polymers in the compartment was ρ_m = 0.2 σ⁻³ (12 800 monomers) with half of the monomers belonging to each type of polymer (ρ₁ = ρ₂ = 0.1 σ⁻³). The total density of the solution in z < 0 σ roughly corresponded to the measured equilibrium density (see the supplementary material). The polymers in the explicit solvent were used as starting configurations for the implicit-solvent model.

We subsequently removed the membrane and allowed the polymers to diffuse. At infinite dilution, the polymers would be expected to diffuse at roughly the same rate in both models because their center-of-mass diffusion coefficients were matched at these concentrations. At finite concentrations, however, hydrodynamic interactions between the polymers, qualitatively associated with solvent backflow, may alter the diffusive flux out of equilibrium.²⁴ Since the initial gradient in Δμ_i across the membrane should be comparable for both models based on Fig. 3, the diffusive driving force is the same in the explicit and implicit solvent, and any differences in the measured flux should be due to hydrodynamic interactions out of equilibrium. Figure 5 shows the fraction of monomers which moved into the opposite side of the box (z > 0 σ) at a given time, averaged over 25 independent simulations. At the steady state, half of the monomers, summed over both short and long polymers, should have z > 0 σ.

It is clear that the dynamics of the implicit-solvent model are faster than the explicit-solvent model, in agreement with the considerations in Refs. 23 and 24. This is also consistent with the measured total monomer density profiles in the drying films, computed from Fig. 1, that show more monomers near the liquid–vapor interface for the explicit-solvent model. Moreover, the long (M₂ = 80) polymers and short (M₁ = 10) polymers responded differently to the same density gradient between the two models. In the implicit-solvent model, the long polymers migrated faster than the short polymers, consistent with Eq. (5) and stratification in the drying simulations. On the other hand, the long polymers migrated much slower than the short polymers in the explicit-solvent model, explaining the lack of stratification in the drying simulations.

We demonstrated the influence of hydrodynamic interactions on stratification in a drying polymer mixture using explicit-solvent and implicit-solvent computer simulations. The implicit-solvent model predicted stratification in agreement with the previous simulations and theoretical considerations.¹⁷ However, no such stratification was found for the explicit-solvent model under the drying conditions considered. Despite good mapping of the equilibrium bulk properties (chemical potential, diffusion coefficient) between the explicit- and implicit-solvent models, hydrodynamic interactions out of equilibrium were shown to alter the polymer diffusion in a way that is consistent with the lack of stratification. Our analysis directly tests and confirms that implicit-solvent simulations and theoretical models lacking hydrodynamic interactions are not capable of quantitatively predicting stratification, in agreement with the analysis of Sear and Warren²³ for the special case of a large colloid in an ideal polymer solution. In future, hydrodynamic interactions must be incorporated into any simulations aiming to study stratification, e.g., through explicit-solvent molecular dynamics or with an appropriate mesoscale simulation method.

While there are currently no experiments available for stratification in polymer mixtures, colloid mixtures have been shown to stratify in experiments.^12–14 The extent of colloid stratification in the experiments¹² and recent explicit-solvent computer simulations³⁸ appears weaker than predicted by models lacking hydrodynamic interactions, consistent with our analysis. Larger chemical potential gradients on the larger component would be required to increase the stratification, which could be induced by, for example, additional cross-interactions between the components.

See supplementary material for the single-chain structure for both models and equilibrium bulk properties of mixtures.

We gratefully acknowledge the use of computational resources supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University. Financial support for this work was provided by the Princeton Center for Complex Materials, a U.S. National Science Foundation Materials Research Science and Engineering Center (Award No. DMR-1420541).